Study of the Dimensionality Reduction Approach for the Efficient Simulation and Optimization Design of Terahertz Lens With Electrically-Large Dimension

Dielectric lenses are widely used in terahertz imaging and communication systems to focus or collimate the Gaussian beam by adjusting the phase distribution of wave front. However, due to the high frequency and short wavelength of the terahertz band, the focusing lenses used in the application systems are usually electrically large, which bring the extremely difficulty to carry out efficient electromagnetic (EM) simulation and optimization design. In this paper, a dimensionality reduction concept was introduced to achieve efficient design and optimization of electrically-large terahertz lenses, which have symmetric structures and are commonly used in quasi-optical systems, such as circular lenses and cylindrical lenses. To precisely solve the EM problem with reduced dimension, a two-dimensional moment method (2D-MOM) for homogeneous dielectric targets was studied and successfully developed by solving the surface coupled integral equation discretized with appropriate basis and test functions. Then, the dimensionality reduction approach with the combination of the ray-tracing method (RTM) and the 2D-MOM was developed for the shaped design of terahertz lens with high efficiency. A 0.3THz lens with diameter 10cm was designed with the proposed approach as an example, with its pattern measured by a terahertz field scanning platform. It’s found that, with the dimensional reduction, the unknowns can be reduced more than 1500 times and the memory required can be reduced more than 2.5 million times, as compared to the traditional 3D-MOM simulation. And the simulation results agree well with the experiments, which both demonstrate the greatly improved performance of the lens designed by the proposed approach, as compared to the standard lens.


I. INTRODUCTION
Terahertz (THz) wave, which is generally referred to the spectrum from 0.1 to 10 THz (or more strictly with lower bound limited to be 0.3 THz), has unique properties due to its special position in the electromagnetic (EM) spectrum, which lies in the gap between the electronics and photonics [1]- [3].
The associate editor coordinating the review of this manuscript and approving it for publication was Wei E.I. Sha .
Unlike optical and infrared radiation, THz wave offers the property of being able to 'see through' obscuring materials such as clothing, cardboard, plastics, and wood with relatively little loss [4]. Compared to microwave and lower radio frequency wave, comparatively broad band-width and high resolution are available in THz band [5]. The above unique advantages make THz technologies promising for plenty of applications in imaging [6], [7], communication [8], [9], and other fields [10], [11]. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Since the lower bound of the THz band overlapped with the submillimeter wave and its upper bound overlapped with the infrared, both the technologies borrowed from the electronics and the photonics play important roles in this band. And the quasi-optical transmission and radiation naturally becomes one of the key technologies in this special band [12], [13], which is based on the concept of fairly well collimated beam propagating in free space. In THz quasi-optics, the beam can be focused and aligned using curved reflectors or dielectric lenses [14]- [16]. For a multiple cascaded quasi-optics, the occlusion problem must be carefully considered if the reflectors are implemented, which makes the design quite complicated. In contrast, dielectric lenses can be implemented to shape the THz beam with inline topologies with the occlusion problem easily avoid and are favored in lots of applications [17]. On the other hand, limited by the comparatively low power available from the THz source in nowadays' technology status, it is usually to implement quasi-optics with large aperture [18] to achieve high beam efficiency to satisfy special applications.
In THz band, it is a great challenge for the simulation and optimization of the focusing lens with electrically-large aperture. The design theory of conventional hyperbolic planoconvex lens widely used in THz quasi-optics is based on the thin lens approximation [19], in which the THz waves or rays are assumed to propagate parallelly within the lens. This will result in large deviation especially for the design of large-aperture lens with thickness far more than the working wavelength. Additionally, for the electrically-large lens, the accurate analysis and design with full-wave EM simulation will cost tremendous computational complexity and resources.
As an example, for simulating a hyperbolic plano-convex lens with a aperture diameter of 10 cm at 0.3THz with the method of moments (MOM) as shown in Fig.1, it is necessary to use over 4 million triangular patches with λ/10 dimension (λ is the wavelength) to model the surface of the lens. The amount of the triangular patches is in proportional to the number of unknowns N in MOM simulation (For dielectric problems, the number of unknowns is twice the number of patches), while the memory requirement of MOM is in proportional to N 2 . Hence, one can get that the memory requirement for analyzing such a lens is more than 1000TB, which is nearly inaccessible. Even employing the fast multipole techniques (FMM) in the MOM simulation, in which the memory requirement can be greatly reduced and become in proportional to N 3/2 , the memory requirement is still more than 300GB for such a lens. An alternative way with high efficiency, the high-frequency approximation methods, such as physical optics (PO), geometric optics (GO), are also used to analyze terahertz lens. However, these methods cannot deal with the multiple reflection between the front and back surfaces of THz lens appropriately, and hence lead to large deviation especially as the dielectric constant of the lens material become large. The Body-of-Revolution (BOR) MoM [20]- [22] is also an elegant way of dealing these lenses which are rotational symmetric. However, BOR-MOM cannot be used to solve the problem of cylindrical lenses, which is also commonly used in THz quasi-optics. Now, how to design and optimize the electricallylarge THz lenses and evaluate their EM performance accurately with high efficiency has become an important and challenging problem for the quasi-optical system design in THz band.
In this paper, an efficient simulation and design method based on the concept of dimensionality reduction was proposed to achieve the optimal design and accurate EM simulation for the electrically-large lenses in THz band with high efficiency. Considering THz lenses with symmetric structures commonly used in quasi-optical systems, such as circular lenses and cylindrical lenses, we demonstrated the theoretical rationality to transform a three-dimensional (3D) EM problem for lens simulation into a two-dimensional (2D) one. To precisely solve the EM problem with reduced dimension, the 2D-MOM for homogeneous dielectric objects was developed by solving the surface coupled integral equation discretized with appropriate basis and test functions. The processing method of the singular points was also investigated and successfully developed for the 2D-MOM, which is always essential to ensure the precision of the MOM analysis. Then, an approach with the combination of the ray-tracing method (RTM) and the 2D-MOM was proposed for the optimization and evaluation of THz plano-convex lens with shaped design. And a proof-of-principle quasi-optics including a dielectric circular lens with diameter 10cm was designed by the proposed optimization method and fabricated in 0.3THz band. While another lens was designed based on the standard thin lens model for comparison. For the simulation of such electrically-large lens with 2D-MOM, only 5360 of unknowns are required, which is nearly 1500 times less than the traditional 3D-MOM simulation, and only 480 MB memory is required, which is 2.5 million times less than the traditional 3D-MOM simulation. The patterns of the two lenses designed in different ways were finally measured by a THz field scanning platform. The experimental results agree well with the simulation and both verified the effectiveness of the proposed theory and approach in this paper. This paper is organized as follows. In Section II, the rationality of the dimensionality reduction approach to analyze the focusing lens is discussed theoretically. The theory of 2D-MOM for solving homogeneous dielectric objects is proposed in Section III. In Section IV, the conventional design theory of standard hyperbolic plano-convex lens is briefly introduced, and the optimization design approach for terahertz lens, based on the combination of RTM and 2D-MOM, is also proposed. The 2D-MOM simulation results for the standard hyperbolic plano-convex lens and the optimized lens are also given in Section IV. In Section V, a 0.3THz test platform is built and the simulation results are verified experimentally. Finally, a conclusion is drawn in Section VI.

II. THEORETICAL RATIONALITY OF DIMENSIONALITY REDUCTION APPROACH
In this section, we discussed the theoretical rationality to transform a 3D EM problem for lens simulation into a 2D one, by analyzing the Gaussian beam solutions of the paraxial wave equations in 3D and 2D cases. Gaussian beam theory can be regarded as the guided wave theory of quasi-optical technology, which is a paraxial approximation derived from the Helmholtz equation.
From EM theory, a simple harmonic wave in free space should satisfy the Helmholtz equation: where k 0 is the wave number in free space, represents any scalar field. In a Cartesian coordinate system, the wave propagation direction is assumed to be positive x direction, then can be written as: Substituting (2) into (1) gives the following expression: According to Goldsmith [12], when the beam is moderately well collimated, which is satisfied by the Gaussian beam in the quasi-optics, the paraxial approximations can be employed as Then (3) can be simplified as The solution for µ is of the form: where λ 0 is the free-space wavelength, w 0 is the radius of the beam waist, w is the beam radius, R is the radius of curvature for the equiphase surface and is the additional phase shift generated during the propagation of the Gaussian beam: For the 2D case, the expression of the paraxial wave equation becomes Solving (11) yields the 2D Gaussian beam expression: where w 2d is the radius of the 2D Gaussian beam, R 2d is the radius of curvature and 2d is the phase shift: By comparing (9) and (14), it is found that the 2D Gaussian beam has exactly same expression of the radius of curvature as that of the 3D Gaussian beam. In a quasi-optics, the function of the lens is to adjust the radius of curvature of the input beam, to realize the beam focusing or collimation as the output beam leaving the lens [23]. Therefore, it is theoretically reasonable to analyze the 3D lens as a 2D problem.
Equation (16) gives the classical transformation formula of the focusing lens [19], where f is the focal length of the standard lens, R 1 and R 2 are the radius of curvature of the output beam at the input and output planes of the lens, respectively.
If the beam is diverging along the propagation direction, the radius of curvature for the beam is positive. While the beam is converging along the propagation direction, the radius is negative. VOLUME 8, 2020 In THz quasi-optics, the standard plano-convex lens with hyperbolic surface was commonly used. However, the design theory of conventional hyperbolic plano-convex lens is based on the thin lens approximation, in which the beam radius at the input and output planes of the lens are assumed to be equivalent. This approximation results in large deviation especially for the design of large-aperture lens with thickness far more than the working wavelength. Hence, the accurate and efficient EM simulation for the performance of the lens is important to evaluate and achieve an optimized design. Based on the concept of dimensionality reduction as discussed in this section, a 2D-MOM was studied and developed in the following section for the simulation of THz lens.

III. MOMENT METHOD FOR SOLVING TWO DIMEN-SIONAL DIELECTRIC OBJECTS
The method of moments (MOM) in electromagnetism was proposed by Harrington [24] in 1968, which is an accurate numerical method and has been developed into one of the main algorithms for computational electromagnetics. MOM is suitable for solving EM problems in open regions such as the antenna radiation problems and EM scattering problems. In addition, MOM is also very effective in solving conductor or homogeneous dielectric problems [25]- [27]. In fact, the 2D-MOM for the analysis of the dielectric objects has been discussed in some previous documents. Based on the surface equivalence principle, Peterson has established the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) for the homogeneous dielectric cylinders in [28]. However, just as Peterson himself has discussed in detail in [28], the numerical stability of the 2D-MOM formulation based on EFIE or MFIE is expected to be improved. Therefore, in this paper, we solved 2D homogeneous dielectric EM problem based on the PMCHW [29]- [31] integral equation, which has better numerical stability.

A. INTEGRAL EQUATIONS FOR 2D DIELECTRIC EM PROBLEMS
As shown in Fig.2, let l denotes the boundary curve of the 2D dielectric object with the outward normal vector n. The electric field E in and the magnetic field H in , defined to be the field due to the impressed source (TE plane wave), are incident on and induce the equivalent current J and the equivalent magnetic current M on l.
The relative permittivity and permeability of the background space and the dielectric body are ε 1 , µ 1 , ε 2 , and µ 2 , respectively. For the moment method to analyze the homogeneous dielectric object, we established the PMCHW integral equations for dielectric lens: (18) where Z m is the wave impedance, L m and K m are two integral operators. And the subscript, m, is equal to 1 or 2, representing the background space and the dielectric space, respectively. In the 2D problem, the integral operators are: where k m is the wave number of the electromagnetic wave in the corresponding region. In the 2D case, the Green's function in free space is the zero-order Hankel function of the second kind [28]: where ρ and ρ are the position vectors of the field point and the source point with respect to the global coordinate origin, respectively. And the gradient of the 2D Green's function is: where H 1 (x) is the one-order Hankel function of the second kind, and ρ SF is the unit direction vector from the source point to the field point:

B. DISCRETIZATION OF 2D INTEGRAL EQUATIONS
In this section we discussed how to choose the appropriate basis function to discretize the integral equations of the dielectric object in 2D-MOM. When dealing with 3D problems, we usually choose the RWG triangular basis functions [32] to discretize the unknowns. For the 2D case, we choose rectangular window function as the basis function, as shown in Fig. 3(a). The boundary curve of the dielectric region is divided into N segments. And the equivalent sources J and M can be approximated by the superposition of sectionalized basis functions: where J j and M j are the values of current and magnetic current, f j and g j are the basis functions defined on the segment j: where µ (x) is a step function, j is the length of segment j, l ∈ 0, j , z is the unit direction vector of positive z direction in the Cartesian coordinate system, and l j is the unit direction vector of the segment j: The next step in MOM is to take a testing procedure. As shown in Fig.3(b), we choose the sectionalized triangle functions ν i and w i as the testing functions: where i is the length of segment i, l ∈ (0, i ) and l i is the unit direction vector of segment i: Equation (19) and (20) are tested respectively with ν i and w i , yielding which can be written in matrix form as: Once the elements of [A] and [b] are determined, we may solve the resulting system of (34) for the unknown column vector [α]:

C. SOLVING THE INTEGRAL KERNELS
In solving (34), the solution of the coefficient matrix [A] is the most complicated. By substituting (19), (20) and (21) into (32) and (33), we can get six different types of integral kernels: In the six integral kernels, since the line divergence of the current basis function is zero, the value of the second integral kernel is also zero, and the other five integral kernels need to be solved. It should be noted that when i = j, there will VOLUME 8, 2020 be singularity problems in the integral kernels. We take the first and the sixth integral kernels as examples to introduce the solution method.

1) FIRST INTEGRAL KERNEL
When i = j, ρ − ρ approaches to zero, the integral of the Hankel function is singular, which corresponds to the diagonal elements in the solution matrix of the first integral kernel.
In order to solve the problem of singular points, we may use the small-argument asymptotic form of the Hankel function on the segment where the source point coincides with the field point: where γ = 1.7810724. As shown in Fig.4, the field point F is integrated over the boundary curve of the dielectric object. Suppose l j is clockwise, n j is the outer normal unit vector, and z points to the paper facing outward. We can get: And d is the vertical distance from the field point to the section line, and it may be a positive value or a negative value: Let's do the singular point integration of the source area first: When d → 0, (41) can be solved directly: The expression of the elements on the diagonal in the solution matrix of the first type of integral kernel is: For the off-diagonal elements in the solution matrix of the first integral kernel, there is no singularity problem. These elements can be directly solved by interpolation method. In order to ensure both the accuracy and speed of the solution, different interpolation precisions need to be selected according to the distance between the field point and the source point. The farther the distance is, the fewer the number of interpolation points is. Similarly, the closer the distance is, the more interpolation points need to be selected. Fig.5 shows the basic interpolation scheme.
The expression of the off-diagonal elements can be written as: where Q is the number of interpolation points.

2) SIXTH INTEGRAL KERNEL
As shown in Fig.6, the singular point integration of the source region is also performed first:  where When d → 0, (45) can be solved directly by using (38): For off-diagonal elements, the solution method of the 6th integral kernel is different from other integral kernels, because the 6th type of integral kernel involves the line divergence of the basis function g j : which makes the solution matrix of the 6th integral kernel also have singularity problems when source point and field point are defined on adjacent segments, |i − j| = 1 or |i − j| = N − 1. It cannot be directly solved by interpolation method. These singularities are discussed separately.
When i − j = −1 or i − j = N − 1, d → 0, the expression of 6th integral kernel is: Combined with interpolation for the second term in (50), we can get: 1 (k m ρ(s)−ρ j+1 ) (51) Similarly, when i − j = 1 or i − j = 1 − N , and d → 0, the solution of the 6th type of integral kernel can be obtained: Expression of the rest elements of the solution matrix can be obtained with interpolation method:

D. VERIFICATION OF 2D-MOM
In order to verify the effectiveness of the 2D-MOM developed in this section, we apply the 2D-MOM and the analytical method to calculate the equivalent electric and magnetic current distribution on the surface of a homogeneous 2D dielectric cylinder in 0.3THz respectively, which is illuminated by TE plane wave. The diameter of the dielectric cylinder is ten free-space wavelengths, and the relative permittivity is 2.25.
The results of the 2D-MOM and the analytical method are shown in Fig.7. The calculation results of these two different VOLUME 8, 2020 methods agree well with each other, which demonstrates the effectiveness and accuracy of the proposed 2D-MOM.

IV. RAY TRACING METHOD AND TWO-DIMENSIONAL MOMENT METHOD TO OPTIMIZE TERAHERTZ LENS
In order to simulate the focusing performance of terahertz lens by 2D-MOM, it is necessary to use 2D Gaussian beam as the excitation. In this section, the conventional design theory of standard hyperbolic plano-convex lens was introduced briefly in Part A. Then in part B, we analyzed a cylindrical lens with a section size of 3cm × 3cm by FEKO and 2D-MOM, respectively. In part C, we designed a standard hyperbolic plano-convex circular lens, and 2D-MOM was implemented for the simulation of the standard THz lens. And great deviation was found between the theoretical values and the simulation results. In Part D, an approach with the combination of the RTM and 2D-MOM was proposed for the optimization and evaluation of THz plano-convex lens with shaped design. Fig.8 gives a schematic show of a Gaussian beam being reshaped by a lens, where w in is the input beam waist size, L 1 is the distance from the input beam waist to the lens, D is the lens diameter, t is the lens thickness, w out is the output beam waist size and L 2 is the distance from the lens to the output beam waist. According to the Gaussian beam theory discussed in section II, one can get the radius of the curvatures for the beams at the input and the output planes:

A. DESIGN AND ANALYSIS FOR THE STANDARD HYPERBOLIC PLANO-CONVEX LENS
Under the thin lens approximation, one can get the relationship between the beam waists of the input and output beams as In THz quasi-optics, the standard hyperbolic plano-convex lens with hyperbolic surface was commonly used to realize above beam transformation with the thin lens approximation. The focal length of lens, f , was determined by R 1 and R 2 under the relationship in (16). The 3D equation [33] for the lens surface can be given as where ε r is the relative permittivity of the lens. From (57), the thickness of the lens, t, is found to be:

B. CYLINDRICAL LENS
We first designed a cylindrical lens, which was illuminated by a Gaussian beam with a beam waist radius of 1.91mm. The distance between the beam waist position and the lens was 0.06m, and the output beam was focused at x = 0.1m. The equation of the section curve for the cylindrical lens can be written as where y ∈ [−0.015m, 0.015m].
147140 VOLUME 8, 2020   Fig.9(a) shows the simulation results of the focusing performance for the cylindrical lens using FEKO, Fig.9(b) shows the field distribution of the output beam in the x-y plane at the focusing position.
It is easy to see that the calculation results of 2D-MOM are basically consistent with FEKO, which verifies the accuracy of 2D-MOM, and also shows the effectiveness of 2D-MOM in analyzing the focusing performance of cylindrical lens, which is the problem that cannot be solved by BOR-MOM, highlighting the application value of dimension reduction method. Table 1 gives the parameters of the input beam and the output beam in an actual focusing system, as well as the required size of a standard hyperbolic plano-convex circular lens.

C. CIRCULAR LENS
According to (57), at the same time, the output plane of the lens coincides with the y-z plane in the Cartesian coordinate system, the surface equation of the standard lens can be written as: The equation for the section curve of the standard lens in the x-y plane is: Then, the 2D-MOM developed in section III is implemented for the simulation of the 0.3THz focusing quasi-optics with the standard hyperbolic plano-convex lens. The 2D Gaussian beam given in section II is used as the excitation. The length of the line segment of the lens boundary curve is less than one tenth of the wavelength. The number of the segments, N, equals 2680, so the number of unknowns generated during the 2D-MOM is 5360. The corresponding memory requirement is nearly 480MB. The simulation results of the 2D-MOM for the standard lens are shown in Fig.10. Fig.10(a) shows the 2D distribution of the magnitude of the output electric field after the input beam passes through the standard lens, the red mark point is the theoretical beam focus position. Fig.10(b) gives the magnitude distribution of electric field along the axis of the standard lens, and x = 0 is referred to the output plane of the lens. Fig.10(c) plots the transversal magnitude distribution of electric field along y direction at the distance x = 0.3577m, where the magnitude of the output beam reaches its maximum.
From the results in Fig.10, the waist radius and waist position of the output beam can be evaluated as From Table 1, we can see that the theoretical values of the beam waist radius and waist position for the output beam are: It is found that, the actual focusing performance of the standard lens is significantly different from the theoretical prediction with thin lens approximation. Therefore, it is necessary to optimize the design of the standard lens.

D. THE OPTIMAL DESIGN FOR THE THZ LENS
In this section, an optimization design approach for terahertz lens, based on the combination of RTM and 2D-MOM, was studied and developed. We innovatively established a reasonable RTM to simulate the process of Gaussian beam passing through the lens, and the focusing performance of the optimized lens was evaluated by 2D-MOM. Then we need to adjust the parameters to be optimized according to the 2D-MOM simulation results until we get the shape of the focusing lens that we want. Based on this, we can achieve high-efficiency optimized design of terahertz lens with electrically-large dimension.
The basic process of the optimization method is schematically shown in Fig.11.
And the practical design procedure is decomposed and described in detail as following: 1) For the standard lens, the parameters to be optimized should be determined first, which is the starting point  of establishing the ray tracing model. By changing the values of the parameters to be optimized, the crosssection curves of different shapes of the lens can be obtained, which is convenient for optimization; 2) According to geometrical optics approximation, the ray-tracing model is established to simulate the focusing process of the terahertz lens on an incident Gaussian beam. The schematic diagram of the ray tracing model is shown in Fig.12. respectively. The intersection of the rays and the desired focusing plane is H x out i , y out i , where i = 1, 2, 3 · · · Num, Num equals 199. And the distance between the feed and the lens is 0.2m, and the distance between the focal position of the output beam and the lens is 0.3m. It should be noted that in the ray tracing model, we set these rays to be distributed in the beam radius range of the Gaussian beam and the distribution of y ip i is equally spaced: where w is the beam radius of the Gaussian beam at the input plane of the lens. Fig.12(b) is a partial enlarged view of geometric rays passing through the lens in the ray tracing model. According to the shape and the refractive index of the lens, the path of each ray can be traced and solved using the law of refraction.
where phiP in i and phiP out i are the incident angle and exit angle of the rays at the illumination surface of the lens, respectively. Also, phiQ in i and phiQ out i are the incident angle and exit angle of the rays at the dark side of the lens, respectively. And n is the refractive index of the lens.
3) By changing the parameters to optimize the shape of the lens, the weighted sum, D average , of the off-axis distances of all geometric rays at the desired focal plane is minimized: where η is the weight coefficient of these rays. It should be noted that for the RTM, the selection of the weight coefficient of each ray is very important, which is directly related to the optimization results. In a Gaussian beam, the distribution of electric field strength is Gaussian as the off-axis distance changes, so the weight coefficients of these rays are chosen as Gaussian distribution in the optimization process: where w x op i is the beam radius of the Gaussian beam at the output plane of the lens. The smaller D average , we think the better the focusing performance of the lens. 4) After the curve equation of the optimized lens is obtained, the 2D-MOM is used to simulate and analyze the lens to check the focusing performance. Then, the parameters to be optimized can be adjusted according to the simulation results. According to (61), a parameterized quadratic curve was implemented to model the boundary curve of the standard lens to be optimized: Based on RTM, we use GA to optimize these undetermined parameters. a 0 and a 1 are encoded to form chromosomes, and the initial population is generated. The value of a 2 depends on the diameter of the lens. D average is the fitness value of the population in each generation. The smaller D average is, the stronger adaptability of population to environment is. Then, according to the evolutionary rules, the initial population is selected, crossed, and mutated. Finally, the optimal solution is obtained. Fig.13 shows the evolutionary curve and the optimized lens curve.
The equation for the section curve of the optimized lens is: where y ∈ [−0.050435m, 0.050435m].
The simulation results of the optimized lens using 2D-MOM are shown in Fig.14. Fig.14(a) shows the 2D distribution of the magnitude of the output electric field after the input beam passes through the optimized lens, the red mark point is the theoretical beam focus position. Fig.14(b) gives the magnitude distribution of electric field along the axis of the optimized lens, and x = 0 is referred to the output plane of the lens. Fig.14(c) plots the transversal magnitude distribution of electric field along y direction at the distance x = 0.2966m, where the magnitude of the output beam reaches its maximum.
From the results in Fig.14, the waist radius and waist position of the output beam can be evaluated as Comparing (73) and (62), (74) and (63) respectively, we can see that the optimized lens has better focusing effect VOLUME 8, 2020  Table 1, which verified the effectiveness of the proposed approach.
For comparison, another design based on the optimization of the focal length, f , is also performed. In this optimization model, the curved boundary of the lens is chosen as standard hyperbolic, while the focal length f of the lens is set as the parameter to be optimized, to achieve the minimum D average based on the similar ray tracing algorithm introduced in section IV. After optimization, the best focal length, f b , of the standard hyperbolic can be obtained as The value of f can be seen in Table 1. Then we can get the cross-section curve equation of the lens with the optimal focal length: where y ∈ [−0.050435m, 0.050435m]. The focusing performance of the lens with the best focal length is simulated and analyzed by using 2D-MOM, and the results are shown in Fig.15. Fig.15(a) shows the 2D distribution of the magnitude of the output electric field after the input beam passes through the optimized lens with the best focal length, the red mark point is the theoretical beam focus position. Fig.15(b) gives the magnitude distribution of electric field along the axis of the lens, and x = 0 is referred to the output plane of the lens. Fig.15(c) plots the transversal magnitude distribution of electric field along y direction at the distance x = 0.3301m, where the magnitude of the output beam reaches its maximum. From the results in Fig.15, the waist radius and waist position of the output beam can be evaluated as: (77) Comparing the 2D-MOM simulation results of the standard lens and the best focal length lens, it is found that the lens with the best focal length has better focusing effect on the incident Gaussian beam than the standard lens. While compared with the optimized lens with the proposed method, the performance of the lens with the best focal length is still expected to be greatly improved.

V. EXPERIMENT VERIFICATION OF TERAHERTZ LENS
In order to validate the optimization and simulation results based on the concept in Section IV, a 0.3THz test platform was developed to measure the focusing performance of the standard lens and the optimized lens. The test platform mainly consists of the terahertz transceiver, the computercontrolled near-field scanning platform, and the fabricated lens and feed horn with the corresponding beam waist. The terahertz transceiver is composed of microwave vector network analyzer (VNA), frequency multiplier and heterodyne receiver. The schematic diagram of terahertz transceiver is shown in Fig.16.
The VNA provides Ku-band radio frequency (RF) and local oscillator (LO) continuous swept signal source with a bandwidth of 3.33GHz and a fixed frequency difference of 25MHz. The RF signal is upconverted into RF1 of 0.3THz band by frequency multiplier. Then, the RF1 signal is divided into two branches by the coupler, one of which is transmitted to the feed horn to radiate the Gaussian beam directly. The LO signal is also upconverted and divided into two branches. One branch signal is mixed with RF1, subsequently input VNA as reference intermediate frequency (IF) signal. And the other branch signal is mixed with the signal received from the probe, input VNA as the measured IF signal. Then demodulate the reference and the measured IF signals coherently, and transmit the field information of 0.3THz to the computer for data processing. Fig.17 shows a photograph of the test platform. In the experiment, the near-field scanning range covers a 2D area of 60 mm × 60 mm in the y-z plane. And the center of the scanning area is aligned with the center of the lens and the aperture of the feed horn. The terahertz signal is coupled to the feed horn to radiate the Gaussian beam. After focused by the lens, the field pattern of the output Gaussian beam is scanned by the probe on the scanning platform. The field information of the output beam in different planes is determined by changing the distance between the probe and the output plane of the lens. Fig.18 and Fig.19 respectively show the measured results for the quasi-optics with standard hyperbolic lens and optimized lens. Fig.18(a) gives the magnitude distribution of electric field for the output beam along the axis of the standard lens. Fig.18(b), (c) and (d) plot the transversal magnitude distribution of electric field at the distance x = 0.3447m, where the magnitude of the output beam reaches its maximum. From the results in Fig.18, the waist radius and waist position of the output beam from the standard lens can be evaluated as w outsta = 4.15mm (79) Here, the waist radius is computed as the average of the waist radius in x-y and x-z plane. Fig.19(a) gives the magnitude distribution of electric field for the output beam along the axis of the optimized lens. Fig.19(b), (c) and (d) plot the transversal magnitude distribution of electric field at the distance x = 0.2872m, where the magnitude of the output beam reaches its maximum. From the results in Fig.19, the waist radius and waist position of the output beam from the optimized lens can be evaluated as w outopt = 3.28mm (81) It's seen that, the experimental results are basically consistent with the simulation results, which proves the effectiveness of optimization and evaluation methods with the VOLUME 8, 2020  dimensional reduction concept in this paper. With comparison of Fig.14 (b) and Fig.19 (a), a slight difference of the field distribution for the output beam can be observed, which may come from the imperfection of the experimental setup, including the unavoidable reflection of the platform and metal box of the receiver.

VI. CONCLUSION
In this paper, an approach of dimensionality reduction for the efficient simulation and optimization design of terahertz lens with electrically-large dimension was proposed.
The theoretical rationality of the dimensionality reduction for the THz lens with symmetrical structures was firstly demonstrated. A 2D moment method for homogeneous dielectric targets was studied and successfully developed by solving the surface coupled integral equation discretized with appropriate basis and test functions, to precisely solve the EM problem with reduced dimension. Then, the dimensionality reduction approach with the combination of the RTM and the 2D-MOM was developed for the shaped design of THz lens with high efficiency. With the proposed approach, a 0.3THz lens with diameter 10cm was designed as an example, with its pattern measured by a THz field scanning platform. It's found that, with the dimensional reduction, only 5360 of unknowns are required for the simulation of such an electrically-large lens, which is 1500 times less than the traditional 3D-MOM simulation, and only 480MB memory is required, which is 2.5 million times less than the traditional 3D-MOM simulation. And the simulation results agree well with the experiments, which both demonstrate the greatly improved performance of the lens designed by the proposed optimization method, as compared to the standard lens, and verify the effectiveness of the proposed theory and approach.