Compared with traditional spherical and aspherical optical surfaces, freeform optical surface offers more degrees of design freedom, and it can be used in the design of imaging systems with more advanced system specifications, better imaging performance, more compact structure and novel functions. During freeform imaging system design and optimization, high imaging performance is an important design target. In addition, the freeform surfaces should be easier to be tested and fabricated. Interferometric surface testing is one of the most accurate methods for freeform surface and it is now increasingly used. To decrease the testing difficulty, the sag difference between the freeform surface and the base sphere or base conic should be as small as possible. For rotationally symmetric systems, the sag difference can be controlled easily and efficiently by using circular-domain orthogonal polynomial surfaces such as Zernike polynomials surface and Q2D polynomials surface. However, for nonrotationally symmetric freeform systems, as rectangular field-of-view is often used, the freeform surfaces often have rectangular aperture, the ability of circular-domain orthogonal polynomial surfaces is limited. Therefore, it is necessary to establish a design method to control the testing difficulty of freeform surface with rectangular aperture.A design method of imaging systems using square-domain orthogonal polynomials freeform surface is proposed. Two kinds of square-domain orthogonal polynomials of Chebyshev polynomials and Legendre polynomials are analyzed and used. The inner product of the surface sag difference using orthogonal surfaces is related to the weighted square sum of the polynomial coefficients. For Chebyshev polynomials, as its weight function is a complicated function of x and y, it is not straightforward to use this property to control the sag difference. However, the sag difference can be controlled by constraining the sum of sag difference around the margins of the rectangular aperture to be zero. In addition, piston and tilt terms in orthogonal surface description should be zero. This can be controlled by constraining linear combinations of surface coefficients to be zero during optimization. For Legendre polynomials surface type, the constraint on the aperture margin can still be used. In addition, as the weight function of Legendre polynomials is one, the square sum of the polynomial coefficients can be used directly, which can be integrated into the total merit function during optimization. Detailed mathematical equations for establishing the design constraints and merit functions can be found in Eqs. (8), (9), (15)-(18). Several design examples are used to show the feasibility and effect of the proposed design method. For Chebyshev polynomials freeform surface, a freeform off-axis three-mirror system whose primary and tertiary mirrors are integrated into one surface is designed. Compared with the design using traditional XY polynomials surface without sag difference constraints, the sag difference of the freeform surfaces in the system using Chebyshev polynomials surface is effectively controlled (Fig.2-3, Tab.3). For Legendre polynomials freeform surface, three kinds of off-axis three-mirror systems are designed: a system whose primary and tertiary mirrors are integrated into one surface, a system with the traditional zig-zag structure, and a system with a cylindrical package and real exit pupil. The design constraints on the sag difference at aperture margins, and the constraints on square sum of surface coefficients are used. Compared with the design using traditional XY polynomials surface without sag difference constraints, the sag difference of the freeform surfaces in the system using Legendre polynomials surface is effectively controlled (Fig.4, Fig.6, Fig.8, Tab.4-6). For the commonly used rectangular surface aperture in freeform imaging system, a design method of freeform imaging system using square-domain orthogonal polynomials freeform surface is proposed. Chebyshev polynomials freeform surface and Legendre polynomials freeform surfaces are used and discussed. Based on the mathematical properties of the two kinds of polynomials, the mathematical constraints on the sag difference at the margins of the rectangular aperture and the constraints on the square sum of the polynomial coefficients are derived. Several design examples are given to show the feasibility and effect of the proposed design method. The design results show that, using the proposed design method and square-domain orthogonal Chebyshev and Legendre polynomials, the surface sag difference between the freeform surface and the base surface can be reduced effectively, and the testing difficulty can be reduced. The proposed method can be used in the design and development of all kinds of freeform imaging systems, and can be easily implemented in optical design software and other computing platforms and environments.