Coherent Fourier Optics Model for the Synthesis of Large Format Lens Based Focal Plane Arrays

—Future sub-millimetre imagers are being developed with large focal plane arrays (FPAs) of lenses to increase the field of view (FoV) and the imaging speed. A full-wave electromagnetic analysis of such arrays is numerically cumbersome and time-consuming. This paper presents a spectral technique based on Fourier Optics combined with Geometrical Optics for analysing, in reception, lens based FPAs with wide FoVs. The technique provides a numerically efficient methodology to derive the Plane Wave Spectrum (PWS) of a secondary Quasi Optical component. This PWS is used to calculate the power received by an antenna or absorber placed at the focal region of a lens. The method is applied to maximize the scanning performance of imagers with monolithically integrated lens feeds without employing an optimization algorithm. The derived PWS can be directly used to define the lens and feed properties. The synthesized FPA achieved scan losses much lower than the ones predicted by standard formulas for horn based FPAs. In particular, a FPA with scan loss below 𝟏 𝐝𝐁 while scanning up to ±𝟏𝟕. 𝟓° ( ~ ± 𝟒𝟒 beam-widths) is presented with directivity of 𝟓𝟐𝐝𝐁𝐢 complying with the needs for future sub-millimetre imagers. The technique is validated via a Physical Optics code with excellent agreement.

In this paper, we propose, the characterization of wide field of view imagers via the derivation of their plane wave spectrum (PWS) in reception. The approach simplifies the design of the lens focal plane array since both the lens shape and the feed radiation properties can be derived directly from the PWS without the need of using an optimization algorithm. The optimal radiation pattern of an antenna feed can then be directly derived by applying a conjugate field match condition [12]. In the case of absorbers, their optimal angular response can be derived by linking the PWS to an equivalent Floquet mode circuit as in [13].
In [14], a numerical evaluation of the incident PWS in a reflector system was described. A much simpler approach using Fourier Optics (FO) was proposed in [13] and [15]. Over a limited applicability domain, the later approach leads to analytical expressions for the PWS for specific geometries for broadside or slightly squinted incident angles. In this work, we extend the FO approach for quasi-optical systems with multiple components and wide-angle applications by combining it with a numerical GO based technique in reception. The analyses in [13] and [15] were aiming to focal plane arrays of bare absorbers. Therefore, the derived PWS has not taken into account the quadratic dependence of the focal field phase. Here, to properly include the coupling between two QO components in the PWS field representation, especially for off-focuses cases, the quadratic phase is efficiently introduced by applying a local phase linearization around the observation point in the focal plane.
The developed technique is then applied to the synthesis of a wide field of view imager complying with the needs for future sub-millimetre imagers for security applications [8], [16]− [18]. For these applications antenna gains of about of 50 to 60dBi are required [8] with about 100 × 100 beams.
Various solutions have been proposed in the past to improve the scanning performance of quasi-optical systems either using Gaussian horn feeds combined with shaped reflector or lens antennas [19]− [22] (with most of the cases over sizing the radiating aperture) and/or determining an optimum focal surface [23], [24], where the array elements are placed [17]; or by using array clusters of feeds to achieve a conjugate field match condition with the focal plane field [25]− [27]. This work considers a relatively simple FPA architecture based on a lens array. All the lens feeds are placed over a flat surface, enabling monolithic integration at high frequencies. The surface shape of the lenses is linked directly to the phase of the incident PWS, while the radiation of the lens feeds is matched to the amplitude of the PWS via a Gaussian model approximation. For simplicity, the main reflector aperture is modelled as a symmetric non-oversized parabola. The obtained performances, validated via a conventional PO analysis, show significantly lower scan loss than it would be obtained by placing Gaussian horns in the optimal focal surface of such reflector as in [17]. The proposed technique could be easily extended to more practical reflector implementation (e.g. a Dragonian dual reflector [28]) by linking the PWS derivation to a GO field propagation in the reflector system and adjust accordingly the lens surfaces, as well as in combination with oversized shaped surfaces.
The paper is structured as follows. Section II describes the proposed FO/GO methodology to derive a PWS field representation in a multi-cascade quasi-optical system, while Section III extends this technique to wide-angle optics. In Section IV, the methodology is applied to a Fly's eye lens array, and Section V presents an application case. Concluding remarks are given in Section VI.

II. COHERENT FOURIER OPTICS
In this section, a plane wave spectrum (PWS) representation, for the magnitude and phase of the focal field is developed. This PWS is derived for a generic QO component illuminated by a plane wave, using a new coherent FO approach. In [15], the PWS represented only the magnitude of the focal field, since the effort was focused on analysing incoherent detectors. Conversely, including the phase in the PWS is now essential for accurately representing the coupling between multiple QO components, depicted in the scenario shown in Fig. 1, or for evaluating the performance of a QO system with a coherent detection scheme. The phase can be efficiently introduced in the PWS by applying a local linearization as shown in this section.
Let us consider a generic focusing QO component illuminated by a plane wave ⃗⃗ =̂− ⃗⃗ • ⃗ , with wavevector ⃗⃗ . As shown in Fig. 1, an equivalent FO sphere centered at the focus of the component can be used to represent the direct field, ⃗ ( ⃗ ), on the focal plane ( = 0) in terms of a PWS ( [15], [13]): where is the radius of the equivalent FO sphere, ⃗⃗ = sin̂, with being the wave-number of the medium surrounding the focal plane, and ⃗⃗ ( ⃗⃗ ) is the PWS of the direct field. The last quantity can be calculated as follows [13]: where ̂=̂+ √1 − 2 / 2 , and ̂× [ ⃗⃗ ( ⃗⃗ ) ×̂] is the GO field component tangent to the equivalent FO sphere. This GO field is defined over the angular sector subtended by the optical system ( 0 in Fig. 1). This GO field can be calculated analytically [13] when a parabolic reflector or elliptical lens is illuminated by a slightly skewed incident plane wave ( ≤ 11°). For larger illumination angles and for a generic QO component, a numerical GO based approach can be employed [29]. Specifically, the field over its FO sphere can be expressed as follows: where, / represent a scenario involving a transmitting (e.g. lens) or a reflective (e.g. mirror) surface, respectively; ⃗⃗ ( ⃗ ) is the incident field evaluated at the point ⃗ of the QO surface (see Fig. A.1 in the Appendix); ̿ = ⊥̂⊥̂⊥ + ∥̂∥̂∥ and ̿ = ⊥̂⊥̂⊥ + ∥̂∥̂∥ are the transmission and reflection dyads, respectively; ⊥ ( ⊥ ) and ∥ ( ∥ ) are the perpendicular and parallel transmission (reflection) coefficients on the surface, respectively; ̂⊥ /∥ (̂⊥ /∥ ) represents the polarization unit vector of the transmitted (reflected) rays; 1 / and 2 / are the principal radii of curvature of the transmitted/reflected wave fronts; is the length of the GO ray propagating from the QO component to the FO sphere, Fig. 1. The expression of the GO parameters in (3) for the transmission case is provided in the Appendix. As for the detailed derivation of the reflection and refraction cases, the reader is addressed to [30]. The integral in (1) resembles an inverse Fourier transform which relates the spectral field ⃗⃗ to the spatial one, ⃗ ( ⃗ ), except for the presence of the quadratic phase term, − 0 ( ) 2 /(2 ) . As an example for demonstrating the importance of including the quadratic phase term into the PWS representation, let us consider a parabolic reflector with a diameter of = 141.4 0 , and a f-number # = 2. The reflector is assumed illuminated by a polarized plane wave with | 0 | = 1 V/m. The same scenario is going to be analysed throughout this paper. As an example here, an incident  Fig. 2(a). The position of the geometrical flashpoint, ⃗ is also shown. We define the geometrical flash point as the position of the peak of the focalized field over a focal plane assuming that no higher order aberrations are present, i.e. the beam deviation factor (BDF) is 1. Figure 2(b) shows the magnitude of the LHS of (2) along = sin when = 0. By considering the quadratic phase term constant and equal to the one taken at the flashpoint (i.e. − 0 |⃗ ⃗⃗ | 2 /(2 ) ) the spectrum is flat over the reflector spectral domain (solid black line). ). Figures 2(c) and (d), show the magnitude and phase values of the focal field ⃗ , respectively, along the x-axis in a region close to the flash point. The result obtained when assuming constant quadratic phase term (solid black line), as in [15], are compared against a reference solution using a standard PO based code (dotted blue line). It is evident that the magnitude of the focal field is accurately represented, but the phase is not.
To properly represent the phase, we can rewrite (1) as the product of two spatial functions: where ( ⃗ ) = − 0 |⃗ ⃗⃗ | 2 /(2 ) is the quadratic phase term, and Specifically, the CFO spectrum is given by: where * is the convolution operator, and Φ( ⃗⃗ ) is the Fourier transform of the quadratic phase term, which can be expressed analytically as: With reference to the previous example, the grey curve of Fig. 2(b) shows the variation of the magnitude of the coherent FO spectrum. The spectrum is now highly oscillating and shifted with respect to the one of the FO approximation. In Figs. 2(c) and (d), it is shown that both the magnitude and phase of the focal field are accurately calculate using (5). However, due to the oscillations present in the convoluted spectrum (see Fig.  2(b)), the numerical convergence of the integral in (5) is more demanding than the one in (1).
We can simplify the calculation of the coherent FO spectrum by approximating the quadratic phase term with a linear one which accurately represents the field only at the surrounding of a specific position in the focal plane. This position is referred to as the CFO position, ⃗ . This approximation is achieved by introducing a new reference system in the focal plane centred at this position, where ⃗ ′ = ⃗ − ⃗ (Fig. 1), and retaining the first two terms of the Maclaurin series of the quadratic phase argument: where ⃗⃗ = ⃗ . Therefore, the convolution operation in (6) simply results in a shift of the FO spectrum along ⃗⃗ . where the coherent FO spectrum is approximated as follows: When examining Fig. 2(b), we can notice that the approximated coherent FO spectrum (dash red line) has a spectral domain similar to the one calculated by using (6), but without oscillations. In Figs. 2(c) and (d), both magnitude and phase of the spatial field are reported (dash red line), respectively. The agreement of the obtained results with the one of the PO solution (dash blue line) is evident. The grey region shown in these figures corresponds to the applicability region of the approximation (8). This region is defined as a circle, centred on ⃗ , with diameter where a maximum phase error of /8 is allowed in the quadratic phase term: where and # are the diameter and the f-number of the corresponding QO component, respectively. It is worth noting that this applicability region does not depend on the chosen CFO position, ⃗ . Figure 3(a) shows the diameter of this applicability region for a few f-number cases of a parabolic reflector versus its diameter . It can be noted that as the reflector f-number increases, the number of beams that could be analysed using (8) decreases. For comparison, the dashed curves in the figure show the applicability region of the spectrum in (1) when a constant quadratic phase evaluated at ⃗ = ⃗ is considered. In the latter case, applicability region depends on the chosen ⃗ , and the approximation can only be used for a region close to the origin, and small # . The diffractive coupling between a primary QO component and a secondary one, as shown in Fig. 1, can be represented using the PWS in (11). The focal field of this secondary QO component can also be represented using (1) and (2). In this case, the GO field at the FO sphere of the secondary QO component, ⃗⃗ ( ⃗⃗ ), is calculated by propagating each incoming plane wave from the spectrum of the primary QO component to the FO sphere of the secondary component. As a result, ⃗⃗ ( ⃗⃗ ) is the summation of the contribution of each plane wave from the spectrum of the primary component.

Coupling of the QO System to Antenna Feeds
Once the PWS of a QO system is derived, the coupling to antenna based feeds can be analysed resorting to a reception formulation [12] where the equivalent Thévenin open circuit voltage of each antenna can be evaluated as follows: ⃗⃗ is the far field radiated to the FO sphere, by the antenna when equivalently fed by a current of 0 ; and is the wave impedance of the medium in which the antenna is embedded.
The power delivered to the load of the receiving antenna can be calculated as being the total power radiated by the antenna when fed with the current 0 . The power delivered to the feed is maximized when its far field is equal to the conjugate of the CFO spectrum. This condition is referred to as the conjugate field match condition.
After calculating the power delivered to the load, one can estimate the aperture efficiency of the entire QO system as = ⁄ , where = 0.5| 0 | 2 / 0 ; 0 is the free space impedance, 0 is the amplitude of the plane wave incident on the main QO component, and is its physical area. By using reciprocity, the electric field, ⃗⃗ , that the same antenna feed would radiate in (  ,  ,  ), at a far distance from the QO system, can be evaluated as follows:

III. WIDE-ANGLE OPTICS
The method reported in Section II can accurately represent the PWS of a QO component within the FO applicability region introduced in [15]. However, this region limits the maximum size of a FPA under analysis. In this section, we extend the CFO method derived previously to cases where the FPA is larger than this applicability region.
For this purpose, we divide a large FPA into sub-regions where at the centre of each sub-region a FO sphere (off-centred) is placed, as shown in Fig. 4. The GO field is then evaluated over the new sphere. The maximum subtended angle of the sphere ( in Fig. 4) is then defined by the region where GO field exists. Once the centre of the reference system is changed to ⃗⃗ , identical steps to the ones described in Section II can be taken to derive the PWS.
The validity region of the FO representation is directly proportional to the radius chosen for the FO sphere [15]. Moreover, the field tangent to an off-focus FO sphere can be approximated by using the GO ray propagation when the surface of the sphere is far from the caustic points of the geometry (where the focal field is maximum). Specifically, the GO representation is accurate when the radius of the off-focus FO sphere is greater than the depth of focus, Δ = 1.77 / , where = /(4 # ) is the Fresnel number. Therefore, to maximize this radius, it is convenient to define as the z-position where the off-focus sphere is tangent to the surface of the QO component (Fig. 4). For this case, the maximum rim angle for the m-th off-focus sphere can be expressed as follows where | ⃗⃗ | indicates the distance of the centre of off-focus FO sphere, in the focal plane, from the QO component origin. As example cases, we considered a parabolic reflector with = 500 0 and # = 4, or one with # = 0.6, illuminated by a plane wave impinging with an angle far from the broadside. The normalized focal field for these two cases is shown in Fig. 5, and compared, with an excellent agreement, with a standard PO based approach.
The FO representation of the focal fields can be derived by performing approximations in the PO radiation integral as described in [15]. Specifically, approximations on magnitude, vector, and phase terms in the integrand. The applicability region for the FO method can then derived by imposing a maximum acceptable value for the error committed in these approximations, for the magnitude and vector cases, and Φ for the phase. By following the same steps as in [15], for the mth off-focus FO, one can define the following validity regions: Figure 3(b) shows the validity region of the off-focus FO for reflectors with different f-numbers assuming = 0.2 and Φ = /8 as in [15]. As it can be easily seen, the broadside FO validity region is larger for greater f-numbers. However, for reflectors with large f-number this region decreases more rapidly as the sphere is farther away from the focus. This is due to the fact that 1/sin grows rapidly when the reflector fnumber is large. It is worth noting that following similar steps, one can derive the FO applicability region in the vertical direction with respect to the focal plane ( in Fig. 4) as described in [30]. This vertical applicability region can be extended further by displacing the center of the equivalent FO sphere in the -direction. This extension leads to the possibility of analyzing non-focal plane arrays such as imaging reflector antennas for satellite communications [31].

IV. FLY'S EYE LENS ARRAY
In this section, it is clarified how the proposed CFO methodology can be applied to FPA based on lens antennas. The geometry of the problem is sketched in Fig. 6.
By extending the applicability region of the FO method, see (16), a large format lens based FPA such as the one in Fig. 6 is divided into several regions. In the middle of each region an off-focus FO sphere is centred. Around the apex position of each lens element, a local phase linearization is performed, see (8), where ⃗ = ⃗ is chosen. As the result, the PWS of the reflector, ⃗⃗ ( ⃗⃗ ), is derived at the surrounding of the lens element. Each plane wave of this spectrum is propagated using a GO approach to a FO sphere defined inside the lens element as shown in the inset of Fig. 6, as:  where ⃗ is the corresponding point on the lens surface, and is the length of the corresponding transmitted GO ray between the lens and FO surface (see inset of Fig. 6).
To derive the PWS of the lens fed by the reflector, the GO fields in (17) are coherently summed: where Ω is the integration domain which is the entire angular region subtended by the off-focus FO sphere of the reflector. When the plane waves impinging on the lens are characterized by small incident angles, i.e. ≤ 11°, the GO field can be approximated (with a 20% maximum error in the field magnitude estimation) as follows [13]:  cos ) where is the eccentricity of the elliptical lens. The condition ≤ 11° and the FO limit given in (16), define the validity region of (19).
To check the validity of the above methodology, let us consider the same reflector geometry described in the previous section but including a focal plane array of elliptical lenses. In defined as # = / (see inset of Fig. 6), and in all the four cases # = 0.6 (i.e., the lens is truncated). In Figs. 7(a) and (b), the lens under analysis is positioned at the focal plane of a reflector at ⃗ = 8 0 # ̂, the plane wave angles of incidence are = 8 0 / and 10 0 / , respectively. The results are compared to those obtained using a standard PO for the left column and multi-surface PO for the rest. The multi-surface PO code is based on the formulation described in [32]. The excellent agreement inside the validity region of the FO is evident. The radius of the central FO applicability region for the discussed parabolic reflector is approximately 10.5 0 # . To demonstrate the necessity of a coherent FO representation for the reflector's focal field, in Fig.  7(a) the lens focal field is also calculated assuming a constant quadratic phase term in the spectrum of the reflector focal field (1). From this figure, it is evident that one commits a large error in analysing the coupling of the lens to the reflector by not accurately describing the quadratic phase term.
In Fig. 7(c)-(d) scenarios which involve off-centred FO spheres in the x-direction are considered. The lens under analysis is positioned in the focal plane of the reflector at ⃗ = 40 0 # ̂ and the reflector is illuminated by a plane wave with an incident angle of = 40 0 / . In Fig. 7(d), the propagation to the lens FO sphere requires the use of the numerical GO, given in (18), since 0 > 11°. The agreement with the multi-surface PO evaluation is very good for all case.

V. WIDE FIELD OF VIEW WITH NON-HOMOGENOUS LENS ANTENNA ARRAYS
It is well known that the scanning capabilities of reflector antennas are limited for large off-broadside angles. Focal plane arrays of homogenous (i.e. identical) horns or lenses have scanning properties proportional to the size of the beam illuminating the focal plane. In [17], formulas to derive the field of view (defined with a 3dB scan loss criterion) were given for opto-mechanical imaging systems. At low frequencies, the use of feed clusters has been proposed to enlarge the field of view [25]- [27]. Here, we investigate, instead, the possibility to enlarge the field of view by properly designing lens based feeds (lens dimension, lens surface and lens feed). The concept is applied to a focal plane array where the elements will be nonhomogenous. The feeds of the lens array are placed over a flat surface to facilitate a monolithic integration at high frequencies.
For lens elements close to the focus of the reflector, the quadratic phase in (1) and the comma phase in the associated reflector CFO spectrum are not significant, and a homogenous lens array can be used with negligible scan penalty.
For mm-and sub-wavelength systems, the use of large fnumbers ( # >1) is common due to their intrinsic larger scanning property [23]. In these cases, the quadratic phase term is the dominant source of error for off-focus lenses and the CFO has a dominating linear phase term. To achieve a conjugate field matching condition, the lens feeds should be laterally displaced along the lens focal plane with respect to the lens focus. For elements even farther away from the centre, the CFO spectrum contains higher order phase terms. These phase terms lead to a widening of the beams impinging on the lens array. To improve the coupling to these distorted fields, one can first enlarge the lens diameters (amplitude matching) and introduce a nonrotationally symmetric lens feed. Secondly, the phase of the distorted CFO spectrum can be matched by reshaping the surface of the lenses. Fig. 8 schematically describes a possible composition of an optimum focal plane array. Here, different regions, filled with different types of lenses, have been identified. As an application case, we consider a scenario compatible with wide-angle QO systems used in the state-of-the-art compact imaging systems [8], [16]− [18] where antenna gains of about of 50 to 60dBi are needed with about 100 × 100 beams.
As the baseline for the design of the FPA, we consider a silicon elliptical lens ( = 11.9) of variable diameter and coated with a standard quarter wavelength matching layer with relative permittivity of = √ = 3.45. The parameters of the considered reflector coupled lens antenna are listed in Table  1. The far field of linearly -polarized lens feeds is modelled via a Gaussian beam as follows: (20) where = sin cos , and = sin sin ; 0 = 1 V/m is a normalization factor; 0 and 0 are chosen in such a way that the antenna far field matches the CFO spectrum at −11dB normalized level. The Gaussian patterned antenna feeds are placed at the lower focus of each elliptical lens. Figure 9 shows the field on the reflector focal plane when 0 (i.e. broadside direction), 15.5, 23.5, 34, and 43.75 beams are scanned. The maximum of the focal field for each considered scanning position is located inside one of the validity region of the central, 1 st , 2 nd and 3 rd off-focus FO sphere located at ⃗⃗ = 18.2 0 # ̂, 32.3 0 # ̂, and 44.4 0 # ̂, respectively. When the reflector is scanning 15.5 beams, the focal field exhibits asymmetric sidelobes, due to the comma phase terms as described in [13], while scanning 23.5, 34 and 43 beams the first two side lobes and the main lobe of the focal field are merged, due to higher order phase errors.
In Fig. 10, the scan loss of this incident focal field is shown (solid grey line). The circle mark represents the number of beams scanned ( = /( / 0 )) through the parabolic reflector before reaching a scan loss of 3dB. The value is obtained by using eq. (3) of [17]. It is worth noting that the incident scan loss curve (solid grey line) calculated here    Table 1. The yellow, green, and orange regions represent the identical (21), displaced (22), and enlarged elements (23) regions, respectively. The value identified by the grey circle symbol shows the number of beams scanned with less than 3dB scan loss [17]. The cross marks indicate the scan loss obtained by using the PO solver of GRASP [33].
In the following subsections, the four FPA regions identified for optimizing the scanning performance of the reflector system are described. In the top row of Figs. 11(a)-(d), the magnitude and phase of the CFO spectrum of the parabolic reflector is shown for several of the cases in Fig. 9.

A. Region 1: Homogenous Lens Array with Identical Feeds
In this region, the diameter of these lenses is chosen as = 2 0 # which roughly corresponds to the width of the main beam of the reflector focal field when looking at the broadside direction. In Fig. 11(a), the CFO spectrum of the lens is compared to the corresponding one calculated from the antenna far field, when the lens element is placed at the reflector focus. It can be noted an excellent matching between the two fields (middle and bottom rows). As a result, the aperture efficiency for the central array element is about 80%. Figure 10 shows the scan loss when an array of homogeneous lenses with identical centred feeds are considered (solid black lines). It is worth noting that for this lens array the scan loss reaches 3dB only after scanning 23 beams. The rapid increase of the loss is due to the phase mismatch between the CFO spectrum and the antenna far field. This phase mismatch is mainly due to the quadratic phase of the reflector focal field. One can calculate the quadratic phase difference is the radius of the reflector FO sphere. Imposing a maximum of /2 phase difference leads to a scan loss of 0.5 dB. Taking this scan loss value as the limit, the maximum number of beams scanned by homogenous lens array (i.e. with identical uniformly spaced feed elements) defines the limit for this region as follows: In Fig. 10, this region is marked with a yellow colour. As expected, at the edge of this region, the identical element array exhibits about 0.5 dB of scan loss. Within the region identified by (21), the architecture of the proposed optimum lens based FPA is also synthesized using identical elements. The scan loss of this array is also shown in Fig. 10 (blue line).

B. Region 2: Homogenous Lens Array with Displaced Feeds
For elements farther away than 1 , see (21), the CFO spectrum exhibits a linear phase as can be seen in Fig. 11(b). One can conjugate match this phase term by displacing the lens feeds laterally in their respective lens focal planes. In this    As shown in Fig. 11(b), for an element 15.5 beams away from the centre of the reflector focal plane both the magnitude and phase of the incident field and the antenna far field are well matched reaching an aperture efficiency of 76%. Fig. 12(a) summarizes the optimum feed position (indicated in Fig. 8) using the procedure described above.
The limit of this region is associated to the higher order phase distortions in the reflector CFO, specifically the comma error. By using the formula derived in [13], for the estimation of the comma phase in the PWS of a parabolic reflector, one can calculate a maximum number of beams scanned by the displaced feeds reaching at most 0.5 dB of scan loss, as follows: In Fig. 10, this region is marked with a green colour. Within this region, the architecture of the proposed optimum lens based FPA is synthesized using the homogenous lens array with displaced feeds. The scan loss of this array is shown in Fig. 10 (blue line). As expected, at the edge of the region identified by (22), this array exhibits about 0.5 dB of scan loss. The performance of the homogenous lens array with displaced feeds is significantly improved with respect to the one with identical elements (black line).Region 3: Non-Homogeneous Lens Array For elements farther away than 2 , see (22), the diameter of the lens elements should increase to compensate the widening of the reflector focal field due to the higher order phase distortions. As shown in Fig. 9, this focal field is asymmetric in this region. We define a larger rim (i.e. diameter) for the lenses in this region by finding the contour of the reflector focal field at a certain level with respect to its maximum, referred here as lens edge taper level. As an example, Fig. 9 shows that a lens element close to edge of the FPA is defined with an edge field taper level of ~− 7 dB. An automatic procedure is established to define the lens rim for every element by initially using a −11 dB edge field taper. However, as mentioned in Sec. III, the FO validity region is also limited in the vertical direction. Therefore, the considered lens heights and consequently their diameters are limited. In the described example scenario, this maximum lens diameter is ~5 0 # . The implemented automatic produce limits the diameters to this number, and consequently, the obtained edge taper levels are reduced at the edge of the array. Fig. 12(b) shows the obtained lens diameters and field edge levels for the considered scenario. The reported edge taper level is for the worst case of the 1D cut over the lens surface, e.g. for scanning in -direction along when = 0. As consequence, the Gaussian beam waists in (20) will be different now in the two main planes.
(a) (b) Fig. 12. The geometrical parameters of the synthesized non-homogeneous lens array. (a) Gaussian feed parameters (black curves), and feed displacement in the lens focal plane (red curve). (b) Diameter of the lens elements (black curve), and edge taper level for each lens for the worst case 1D cut over its surface. The yellow, green, and orange regions represent the identical (21), displaced (22), and enlarged elements (23) regions, respectively. Fig. 11(c) shows, for the lens element located 23.5 beams away from the centre, that the field match between the lens CFO and Gaussian feed is very good, both in magnitude and phase. Figs. 12(a) and (b) summarize the optimum Gaussian feed parameters and lens diameters and for all regions, respectively. By using the formula derived in [17], one can calculate the maximum number of beams scanned in this region with a scan loss below 0.5 dB, as follows: In Fig. 10, this region is marked with orange colour. Within this region, the proposed optimum lens-based array is synthesized using the design steps described in this subsection. As expected, at the edge of this region, the array exhibits about 0.5 dB of scan loss.

C. Region 4: Non-Homogeneous Array with Shaped Lens Surface
For elements farther away than 3 , see (23), the CFO spectrum cannot be matched with a translated non-symmetric Gaussian lens feed. Fig. 11(d) shows a significant difference in phase distribution between the two, leading to about 5dB scan loss for this case. To improve this scan loss, one can reshape the surface of the lens to remove the higher order phase terms on the lens CFO. Specifically, the difference between the phase of the elliptical lens CFO spectrum and the translated nonsymmetric Gaussian lens feed, referred to as the hologram phase, is approximated by a Zernike expansion [34], [35]. The surface of the elliptical lens is then modified using the following expression: where is the modification of the height of the lens (see Zernike approximation of the hologram phase; and = cos is the -component of the wave-vector in the lens material. In the region outside the one identified by (23), the proposed optimum lens-based array is synthesized according to the design steps described in this subsection using enlarged lens elements with modified elliptical surfaces. The scanning performance of this array is shown in Fig. 10 (blue line).
As an example case, the surface of a lens element located at ⃗ = 43.75 0 # ̂ is considered. Firstly, the hologram phase for this example case is calculated. Secondly, this phase is represented by a =30 =30 Zernike expansion, Fig. 13(a). Finally, the required height modification over the elliptical shape is calculated using (24), Fig. 13(b). The required modification of the lens surface is within the specifications given by commercial silicon micro-machining companies [36]. By reshaping the surface of this lens element, the system scan loss is improved from 5 to 1 dB.

D. Validation of the Methodology
In this subsection, the coupling of the described quasi-optical system calculated using the proposed methodology is compared to the one obtained by performing a PO analysis that exploits the reciprocity of the problem and studies it in transmission. In particular, the field radiated outside the lens antenna is obtained by using an in-house developed PO formulation similar to the one described in [37]. Depending on the array element under study, the lens surface is either elliptical or modified elliptical. According to the size of the lens element and its distance from the parabolic reflector, the field is calculated in the lens radiative near field or in the far field region. This field is then provided to the PO solver of GRASP [33] as a tabulated source illuminating the parabolic reflector, to obtain the field radiated by the entire quasi-optical system. In the proposed CFO method, the first-order PO diffraction effects are taken into account; while in GRASP simulation, the diffraction contribution from the edges (using PTD method) are also included. Table 2 compares the aperture efficiency, evaluated with both methods for the four considered examples in Fig. 11. The same excellent agreement can be observed in Fig. 14, where the radiation patterns of the complete quasi-optical system are shown. Moreover, the scan loss obtained by the PO analysis in transmission is shown with cross marks in Fig. 10. Again, the results are very well matched to the ones obtained by the proposed CFO method. It is worth noting, that the CFO derivation provides the lens and feed geometries with a single calculation that lasts about 4 minutes per lens element. In comparison, the PO analysis in transmission takes about 30 minutes in the same computer. Therefore, this second analysis procedure would lead to very long elapsed times to find the optimal lens geometry using iterative simulations. All the simulations were performed by using a single core Intel i7-4790 processor with a clock frequency of 3.6 GHz, Cache and RAM memory of 10MB and 16GB, respectively. The solid lines, and dot marks represent the pattern obtained in transmission by using PO, and reception by using the proposed method, respectively. The former and the inset, illustrating the pattern in the u-v plane, are calculated by using GRASP.

VI. CONCLUSION
Imaging systems at millimetre and sub-millimetre wavelengths are entering a new era with the development of large format arrays of detectors. A fly's eye lens array coupled to absorbers or antennas is a common FPA architecture for such imagers. Typically, such FPAs are coupled to a quasi-optical (QO) system involving reflectors. For large QO systems, a fullwave electromagnetic analysis is not feasible since it is numerically cumbersome and time-consuming.
In this paper, the original Fourier Optics (FO) procedure has been extended to derive the spectrum of the incident field on a reference system centred on antennas located at a large distance from the focus. The procedure, named here "coherent" FO, has been used to express the spectrum of the incident field in realistic cases which include large arrays of lenses within reflectors focal planes. In particular, the methodology can be linked to spectral techniques commonly used for arrays, such as Floquet mode theory, for analysing absorbing mesh grids, and antenna in reception formalism to analyse the performance of antenna feeds in reception. By introducing the off-focus FO approach, the proposed coherent FO representation can be used  to analyse, and design systems composed of tens of thousands of pixels, while the original FO would provide accurate spectra for only a few tens of lenses. The technique can be used to assess the scanning performance of large format lens based FPAs. In particular, by using the developed analysis tool, it was shown that a scan loss lower than the one of the direct field (given by standard formulas in the literature) can be achieved for a wide-angle optics coupled to a lens based FPA. The proposed array is synthesized according to the described design rules, namely field matching between the CFO spectrum and the far field of the lens feed. It is worth noting that in this design process no numerical optimization algorithms were employed.
Here, scan loss of less than 1 dB has been achieved while scanning up to ±17.5° (~± 44 beam-widths) for an example relevant to the state-of-the-art wide-angle imaging systems with reflector f-number of 2 and directivities of 52dBi. Finally, the proposed technique has been validated via a standard Physical Optics based analysis in transmission with excellent agreement.

APPENDIX: GO PROPAGATION THROUGH DIELECTRIC MATERIALS
The EM fields reflected by curved multiple surfaces can be evaluated using a GO formalism as described in [29], [38], [39]. The propagation of GO fields through dielectric surfaces is instead, to our knowledge, not exhaustively treated in the literature. This appendix summarizes the formulas describing the field transmission and propagation, a key aspect for analysing lenses with the proposed CFO formalism.
In particular, the transmitted GO electric field at an observation point, ⃗, inside a dielectric object (Fig. 1A) can be expressed as follows: where ̿ = ⊥̂⊥̂⊥ + ∥̂∥̂∥ is the transmission dyad; ⊥ and ∥ are the perpendicular and parallel transmission coefficients on the surface, respectively; ̂⊥ /∥ , and ̂⊥ /∥ represent the polarization unit vectors of the incident and transmitted rays, respectively; ⃗⃗ ( ⃗ ) =̂− 0̂• ⃗⃗ is the incident plane wave on the dielectric object propagating along ̂ direction; is the distance between the refraction point, ⃗ , and observation point, ⃗; is the propagation constant in the denser medium; 1 and 2 are the principal radii of curvature of the transmitted wave front and can be calculated according to an equation system as follows: where is the relative permittivity of the dielectric object, and ̂ is the normal unit vector at the dielectric interface pointing toward the direction of the impinging wave (Fig. 1A); is the refraction angle; ̂1 , and ̂2 are the principal unit directions of the surface; 1 , and 2 are the principal radii of curvature of the surface. It is worth noting that the expression of the GO transmitted field, (A1), can be derived by asymptotically evaluating the PO surface integral at the interface between the two media. The GO ray contribution corresponds to the stationary phase point of this PO integral. For further details, the reader is addressed to [30], where the generalization to an arbitrary astigmatic incident wave front is also discussed.