Application of the semi-analytical Fourier transform to electromagnetic modeling

The Fast Fourier Transform (FFT) algorithm makes up the backbone of fast physical optics modeling. Its numerical effort, approximately linear on the sample number of the function to be transformed, already constitutes a huge improvement on the original Discrete Fourier Transform. However, even this orders-of-magnitude improvement in the number of operations required can fall short in optics, where the tendency is to work with field components that present strong wavefront phases: this translates, as per the Nyquist-Shannon sampling theorem, into a huge sample number. So much so, in fact, that even with the reduced effort of the FFT, the operation becomes impracticable. Finding a workaround that allows us to evade, at least in part, these stringent sampling requirements is then fundamental for the practical feasibility of the Fourier transform in optics. In this work we propose, precisely, a way to tackle the Fourier transform that eschews the sampling of second-order polynomial phase terms, handling them analytically instead: it is for this reason that we refer to this method as the “semi-analytical Fourier transform”. We present here the theory behind this concept and show the algorithm in action at several examples which serve to illustrate the vast potential of this approach. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

(e.g. LPIA) are defined in the spatial domain. As we know, the connection between these two domains (namely, the spatial domain and the spatial-frequency domain) is given by the Fourier transform: it becomes paramount, then, if we are going to combine within the same system methods defined in both domains, to optimize the Fourier-transforming step.
These two basic ideas are added to the field-tracing concept, alongside the tearing and interconnection previously mentioned, to further improve its performance: first, to always solve Maxwell's equations in the domain where it is most efficient to do so, and second, to perform the consequently unavoidable Fourier transforms as efficiently as possible. The above is illustrated in Fig. 1: to realize the smooth combination of different field solvers in different domains, the Fourier transform is needed. Typically, for a complicated modeling task which involves various field solvers and requires multiple times shift between different domains, the performance of the Fourier transform technique significantly influences the efficiency of the simulation. Fig. 1. Illustration of the physical optics modeling of a light path by a field-tracing diagram [5]. P indicates the physical free-space propagation operator in the κ-domain. B andB represent the physical response of an optical element in the ρ-or in the κ-domain respectively. The entire process includes the handling of six electromagnetic field components.
We come then to the question of the different computational implementations of the mathematical operation that is the Fourier transform. The brute-force approach that stems from simply discretizing the Fourier integral is commonly referred to as the Discrete Fourier Transform, or DFT, and presents a complexity of O(N 2 ) for a function with N samples. The computational effort of the operation is already greatly reduced by the Fast Fourier Transform, or FFT, the most famous version of which is the one introduced by Cooley and Tukey [6][7][8], and which requires N log N individual computations.
In optics, however, even the vast computational improvement of the FFT over the DFT often falls short: the Nyquist-Shannon sampling theorem, which stipulates a condition of minimum sampling, can yield an impracticably large sample number even for relatively simple fields, like a spherical wave with moderately high numerical aperture (NA). Only the most paraxial fields can be feasibly Fourier-transformed with the FFT. This caveat for the practicality of the FFT in optics comes as a consequence of the need to resolve the complex amplitude with a 2π-wrapped phase [9,10].
Various authors have proposed methods for the fast calculation of the field distribution in the focal region of a non-paraxial imaging system [11,12]. A series of papers based on the Chirp Z-transform (CZT) is recognized by the optics community [13,14]. All these algorithms are founded on the convenient property whereby the quadratic phase term in the Fourier transform is extracted. It can be used to advantage in terms of numerical effort since, according to convolution theory, a single Fourier-transform operation can be replaced by a pair of Fast Fourier Transform steps which, even together, require a much lower sampling number.
Another useful trick included in this technique is that the Fourier transform of a quadratic phase turns out to also be quadratic. It is demonstrated by L. Mandel and E. Wolf for different diffraction integrals [15]. However, the discussion of this quadratic-phase trick tends to be limited to one-dimensional (x 2 ) or, at most, separable problems (x 2 , y 2 ), without including the cross term xy which constitutes a prominent part of diffraction theory [16,17]. What's more, in some cases the scaling factors of the quadratic phase term are deduced via the paraxial approximation (which mathematically translates as using the Taylor expansion of a spherical phase around its center, instead of the full spherical phase). Consequently, the performance and accuracy of any approaches thus limited may suffer in the face of non-paraxial configurations.
In fact, generalization of the CZT to consider two dimensions and the cross-talk term are been done [18,19]. However, although the end result bears a definite resemblance to the conclusions we present below, the derivation of the CZT is centered on the Fourier-transform kernel of the already-discretized operation, while the method we put forth focuses, instead, on the field term. In 2016, Adrianov proposed another derivation process which starts from the phase feature of the field, but his work is applied to a pulse simulation [20] and, consequently, the algorithm is restricted to one dimension. The author does mention the possibility to extend the approach to two dimensions, but does not elaborate further in this direction. The similarities between the methods here cited and our own proposed technique stem from the use of convolution theory applied to quadratic phase terms.
It should be emphasized that the reference work mainly considers the propagation problems, whereas we concentrate on the FT itself, with the aim of reaching conclusions which can be applicable in other situations. The semi-analytical FT can be employed for all the different situations arising in field tracing in which the FT is needed, e.g. for propagating to grating or for the calculation of the E z component from E x and E y .
In this work, the authors introduce an algorithm, which we have named "semi-analytical Fourier transform", whose aim is to efficiently compute the Fourier transform of a field without reverting to any approximations. We describe the full derivation in Section 2, taking as our starting point the mathematical definition of the Fourier transform and using convolution theory. We then consider the application of the algorithm to numerical simulations. Two fundamental simulation techniques, free-space propagation and calculation of the E z component, for which Fourier analysis is a vital tool, are discussed in relation with the proposed technique. Setting our sights on the more general objective of reducing the sampling effort for fields with general wavefronts, we try out different fitting methods in order to extract a polynomial of the second order from said general wavefronts. Simulation results are presented alongside the mathematical derivations to illustrate the potential of the semi-analytical Fourier transform.

Derivation of the semi-analytical Fourier transform
The electromagnetic field at a certain plane is given by where = 1, . . . , 6 to account for all six (three electric and three magnetic) components. In our notation, r = (x, y, z) and ρ = (x, y) are respectively the position vector and its projection onto the transversal plane. The complex amplitude V (ρ) can be separated into two parts, amplitude, |V (ρ)|, and phase, exp (iγ (ρ)), as shown in Eq. (1). Let us next turn our attention to the phase term exp (iγ (ρ)) in Eq. (1). The relationship between the complex exponential function and the trigonometric functions, namely exp (iγ (ρ)) = cos γ (ρ) + i sin γ (ρ), is well-known. The periodic nature of the trigonometric functions constrains all phase points to the range (−π, π], also referred to as a "2π-wrapped phase". It is important to note that a phase function that is "smooth" (C 1 ) in its unwrapped incarnation becomes jaggedly discontinuous when wrapped. The steeper the slope of the original, unwrapped, phase function, the denser the teeth of its wrapped counterpart become. The numerical computation of the FFT requires a well-sampled field which fulfills the Shannon-Nyquist theorem of at least two sampling points per period of cos γ (ρ) and sin γ (ρ); the periods are inversely proportional to the slope of the unwrapped phase function, i.e., to the gradient of γ (ρ). The largest slope in the phase function restricts the period which must be used in the sampling.
In order to overcome this sampling overkill that plagues the FFT algorithm, we must reduce the maximum gradient of γ (ρ). To this end, we present one workaround in this publication: the semi-analytical Fourier transform. The fundamental notion on which this method rests, as the name itself suggests, is the analytical handling of certain phase terms for which a known analytical solution exists. Linear phases, for instance, benefit from the well-known shift theorem which constitutes a property of the Fourier transform. In this work we shall concentrate on the polynomial of second order.
We rewrite the representation of the field with We extract the polynomial of the second order with real-valued coefficients C, D x and D y . In addition, we combine the amplitude and whatever remains of the phase after extracting the polynomial of the second order, and denote it by U (ρ). Although U (ρ) must still be sampled in the numerical operation, arg [U (ρ)] has a less pronounced local gradient than the full phase function arg [V (ρ)] = γ (ρ). Undoubtedly, by this approach, the sampling effort is reduced. However, that only helps, if the demanded FFT can be performed without the sampling of exp iψ q (ρ) . It means that we must handle this information analytically in the entire process.
Next, we attempt to calculate the Fourier transform of V (ρ) without sampling the quadratic phase term defined in Eq. (5).
Plugging Eq. (2) into the Fourier transform equation, we obtaiñ where κ = k x , k y is the projection of the spatial frequency vector onto the plane transversal to z. From the convolution theorem, we know that where the symbol " * " indicates convolution.
In principle, the term F κ [U (ρ)] must be treated numerically. On the other hand, we find that the Fourier transform of the polynomial of second order F κ exp iψ q (ρ) can be solved with the integral formula With the help of Eq. (8), we deduce the Fourier transform of the polynomial of second order and with constant factors α,C,D x andD y It should be remarked that when C 2 − 4D x D y = 0, the two-dimensional phase function ψ q (ρ) degenerates into a rotated one-dimensional function. For such special cases, we need to do the coordinate transformation and then treat them similarly in one dimension.
Substituting for F κ exp iψ q (ρ) in Eq. (7) and expanding the convolution, we get And by defining the coordinate transform vector β = (β x , β y ) the integral in Eq. (13) can be written as The interpretation of Eq. (15) as the product of a Fourier transform and a polynomial of the second orderṼ has the same form as the field representation in the spatial domain, that is, a numerical spectrum and an analytical quadratic phase. In a numerical implementation of the method here proposed, all the above-presented formulas can be written in the discrete, finite-dimensional version, e.g.
and y n ∈ − ∆Y 2 , ∆Y 2 . According to sampling theory, the sampling number is determined by the extension of V (x m , y n ) andṼ k xm , k yn in both domains, i.e. N x = ∆X/ 2π ∆K x and N y = ∆Y / 2π ∆K y . These numbers refer to the full field V (x m , y n ). Then, considering our chosen field representation, in the case of a strong quadratic phase ψ q (x m , y n ), at the boundary of the effective range there would be a very high local gradient ∇ψ q (x m , y n ), which will result in a large extension in the other domain: this is the reason behind the high sampling effort of the full field. On the other hand, the sampling required for the residual field remains relatively small. It is important to note that, in the semi-analytical Fourier transform, the quadratic phase is always handled analytically. When the semi-analytical Fourier transform is being performed, the first FFT is numerically used on U (x m , y n ), with a much lower sampling number than that required for the full field. The second, inverse, FFT is not only applied to the residual spectrum, but also to an additional phase term. In Eqns. (10)(11)(12) we see that the quadratic phase factors in the spatial-frequency domainD x ,D y and C are inversely proportional to the quadratic factors in the spatial domain. Therefore, when the field possesses a strong quadratic phase in the spatial domain, the corresponding quadratic phase in the spatial-frequency domain will be very weak, which means we can neglect its influence on the sampling of Ũ (κ) exp iφ q (κ) . Therefore, even though an additional FFT is used, the numerical effort of performing two FFT is still smaller than that required for the FFT of the fully sampled, complete field.

The inverse semi-analytical Fourier transform
Analogously to the regular Fourier transform, the semi-analytical FT has its reverse operation, namely the inverse semi-analytical Fourier transform.
The spectrum in the spatial-frequency domain is written as whereφ q (κ) is the polynomial of second order in the κ-domain which has been previously described in Eq. (10). Using the same trick as Sec. 2.1, we can deduce the analytical representation of the inverse semi-analytical FT and the constant factorα and with the coordinate transform vectorβ = β k x ,β k y Up to this point, we have obtained the representation of the inverse semi-analytical FT. Similar to the forward operation, the polynomial of the second orderφ q (κ) in κ-domain is fully analytically handled in the entire process. The sampling effort depends on the two FFTs of the spectrum with less sampling, theÃ (κ) in Eq.(21).

Handling of quadratic phase
The goal of the semi-analytical FT is to improve the computational efficiency of the Fourier transform operation. As we discussed in the previous section, the sampling effort of the (inverse) semi-analytical FT is entirely dependent on the sampling of the residual field U (ρ) or the residual spectrumÃ (κ). One way that works is to minimize the gradient of arg [U (ρ)] or arg Ã (κ) over the given function support. For this purpose, we must find an effective method of determining the optimal factors of the quadratic phase ψ q (ρ) orφ q (κ).
There are two alternatives, analytical or numerical. The most widely recognized analytical approach is the Taylor expansion. For a given phase function γ (ρ), by using Taylor expansion up to the second order, we can get ψ q (ρ) analytically where γ xx , γ xy and γ yy are the entries of Hessian of the phase function γ (ρ) at the reference point ρ 0 = (x 0 , y 0 ). The h(ρ) presents the constant, linear and higher order Taylor polynomials. Indeed, according to the analytical Taylor formula, we can directly calculate the desired quadratic factors without any numerical operations. However, the demand for an analytical phase function γ (ρ) and the fact that the formula has only local validity constrain the usage and the performance of the Taylor expansion. According to Eq. (24), we can say in a neighborhood that around the reference point, the gradient of arg [U (ρ)] would be quite small. But, in a range farther away from the reference point, we cannot conclude this. Besides the Taylor expansion, there are also various analytical approaches, e.g., the Avoort fit [3]. However, all of them have the same limitation that they are only applicable to specific cases.
Because of this kind of shortcoming, for more general situations we decide to use numerical approaches. As an example, we choose the Levenberg-Marquardt method [21] which is devoted to solving minimization problems and is especially suitable for the least squares curve fitting. The fitting model of quadratic terms can be written as Here, arg [U (ρ)] is the deviation between the value of the actual function and that of the fitting result. Typically, a numerical fitting is an iterative process and the merit function is set to the imperfect function value arg [U (ρ)]. However, in this task our goal is to minimize the sampling, which is only related to the gradient of the residual phase. Therefore, we select ∇ (arg [U (ρ)]) as the merit function instead. ψ (x) ψ (a) Quadratic phase fitting of ψ (ρ) (b) Fitting deviation of two methods In order to demonstrate the performance of the numerical and analytical approaches, we select a general phase function as an example. In Fig. 2(a), there is an analytically known spherical phase with coma and astigmatism aberration. Both the Taylor expansion and the Levenberg-Marquardt method are applied for the quadratic phase fitting.
The comparison of the deviation of the Taylor expansion method and the LMM is shown in Fig. 2(b). We can see that the Taylor expansion method only works well locally around the central zone, which is the reference point of the expansion. Both the deviation of the function value and the gradient of the deviation at the edge are much larger than for the other approach. Obviously, compared to the analytical Taylor expansion, the numerical method provides better results. Therefore, we recommend using numerical methods in practical tasks.

Numerical examples
In the previous sections, the background knowledge and theoretical derivation of the semianalytical Fourier transform algorithm have been introduced. This new technique has been implemented in the physical optics modeling and design software VirtualLab Fusion [22]. In the following sections, we will apply this new approach to different numerical experiments to compare its performance with that of the regular FFT. All simulations were done with the optics software VirtualLab Fusion.

Propagation of Laguerre Gaussian 01-Mode in Free Space
An essential problem in physical optics modeling is to simulate the propagation of an optical field in free space (a homogeneous and isotropic medium). Figure 3 illustrates the field tracing diagram for Gaussian beam propagation. initial plane target plane A linearly E x polarized Gaussian field, with a wavelength of 532 nm and a beam waist radius of 1 µm, is employed as the input, Fig. 4(a). Its half-divergence angle is about 19°and its Rayleigh length is z R = 2.95 µm. First of all, a Fourier transform operation is used to tackle the plane-wave decomposition. Since the input field includes only a vertex phase, the regular FFT technique is selected. Then, by multiplying the propagating kernel exp (ik z ∆z) on each plane wave, we obtain the spectrum at the target planeṼ out ⊥ (κ). The last step is to recombine the resulting spectrum via an inverse Fourier transform operation. Because of the phase modulation of the propagating kernel on the spectrum, we can try to use the semi-analytical FT to carry out the inverse Fourier transform operation. Field distributions at distance z = 10 µm and z = 300 µm are respectively presented in Fig. 4(b)-4(c). It should be noted that the semi-analytical FT algorithm provides the phase without the quadratic term, namely the residual phase arg [U (ρ)]. In contrast to the regular inverse FFT, which requires (801 × 801) sampling points for Fig. 4(c), the sampling effort of our approach is much lower, (101 × 101).
Another effect we can observe in Fig. 4 is that as the propagating distance increases, the residual phase arg [U (ρ)] encompasses an ever more complicated high frequency part. This is due to the fact that the semi-analytical FT cannot analytically handle the entire propagating-kernel phase term, only the quadratic phase part will be taken into consideration. The gradient of the residual phase is still proportional to the propagating distance ∆z. To investigate the numerical performance comprehensively, we did more experiments of different propagating distances and recorded the required sampling points for both the FFT and the semi-analytical FT. As a result of the analysis, the curve of the sampling numbers (sampling points in one dimension) versus propagating distance is presented in Fig. 5. Apparently, after getting rid of the quadratic phase, the sampling of U (ρ) becomes much easier than that of the full phase. In the case of considerable propagation distances, the reduction of numerical effort is especially evident. Fig. 6. The schematic representation of an aberrated spherical wave going through its focus in free space, and the corresponding field tracing diagram illustrating the same propagation process.

Ideal/aberrated Spherical Wave Goes Through the Focus
Besides the Gaussian beam, another common propagation issue in optical modeling is the focusing of spherical waves. In practice, one typical application is the calculation of the point spread function (PSF) of an imaging system, i.e., propagating the field from the exit pupil to the focus region. In a similar way as for the first example, a modeling task has been set up as depicted in Fig. 6. The input field used in the simulation is an E x polarized truncated spherical wave with a wavelength of 532 nm. It has an initial diameter of 2.5 mm × 2.5 mm, and the spherical radius is equal to −3 mm. It means that this spherical wave will be converging to a point after 3 mm if we don't consider the diffraction effect. Spherical phase ψ sph = sgn (r) k 0 x 2 + y 2 + r 2 r = −3 mm partly analytical Vertical astigmatism Z −2 2 = x 2 − y 2 10 fully analytical Oblique astigmatism Z 2 2 = 2xy 5 fully analytical Horizontal coma The field distribution at the input plane and at different target planes is shown in Fig. 7. In this example the initial field includes a strong convergent spherical phase so that the forward Fourier transform operation can be a semi-analytical FT. Then, the selection of the inverse Fourier transform technique depends on the characteristics of the resulting spectrum. In the case of , since the field does not include a strong quadratic phase, the FFT is used and the full phase is shown. Fig. 7(b), the target field is located at the focal plane, where the phase component does not include large smooth phase terms. The numerical effort of the Fourier transform mainly depends on the amplitude component |V (ρ)|. Therefore, in such a situation we would use a regular inverse FFT. On the other hand, for the non-focal region cases, e.g. Fig. 7(c), the semi-analytical FT is suitable. By comparing the simulation results on the symmetric planes Fig. 7(a) and Fig. 7(c), we can see the same amplitude distribution but different phase patterns. This is because of the usage of the Levenberg-Marquardt method. The numerical handling of the quadratic phase terms in the two Fourier transform operations present a small difference. It must be emphasized that this kind of minor numerical difference will not lead to a strong effect on the sampling issue. For instance, the required sampling numbers for Fig. 7(a)  Up to now, we have discussed two basic experiments, the Gaussian beam, and the ideal truncated spherical wave. The semi-analytical FT is used to deal with the propagating kernel exp (ik z ∆z) or the spherical phase exp iψ sph (ρ) in either the spatial domain or the spatial-frequency domain. As we have mentioned, the semi-analytical Fourier transform is a rigorous algorithm, and from the physical point of view, there is no restriction. Besides these specific examples, we would like to apply it to some more general but complex situations, e.g., the propagation of an aberrated spherical wave from the exit pupil of a real lens system. To have an intuitive feeling, the aberrated phase information of the field is presented in Table 1 in the form of Zernike polynomials. As described in the table, for different types of phase we have different handling strategies, fully analytical, partly analytical or numerical. Please note that the semi-analytical FT can handle not only the 1D quadratic phase, x 2 or y 2 , but also the cross-talk term xy.
In Fig. 8 and Fig. 9, the simulation results at different planes of the aberrated spherical wave propagating towards, through and away from its focus are shown. In contrast to the propagation of the ideal spherical wave, we can see strong aberration on the phase patterns at different observation planes. Meanwhile, due to these aberrations, there no longer is an Airy pattern in the focal region and the amplitude distributions at symmetric positions with respect to the focal region are different. Finally, we also investigate the numerical effort of the FFT and the semi-analytical FT in the case of the propagation of an aberrated spherical wave. Figure 10 shows the curves of the sampling numbers versus propagating distance. As in the first example, the analysis reveals that the semi-analytical FT is beneficial for the Fourier transform of the field with strong quadratic phase; the numerical effort thereof can be reduced significantly.

Calculation of the E z component for a spherical wave with aberration
The last case under investigation is another plane-wave decomposition-based algorithm: the calculation of the E z component from the E x and E y components. Generally, in practice, instead of storing all six electromagnetic components we will record only two of them.This is allowed by the fact that, according to Maxwell's equations, it is possible to use just two out of the six electromagnetic components to describe the field in free space, as the remaining four components follow from those two via the equations. Additionally, Maxwell's equations in homogeneous media become entirely algebraic in the κ-domain. Naturally, these algebraic equations connect the components.
In terms of the simulation results of the last example, we make an attempt to calculate the E z component at different target planes. In Fig. 11 and Fig. 12, the field distribution of the E z component for the aberrated spherical wave are presented. Analogously to previous examples, the same phenomenon is observed here that for any Fourier transform based algorithm we can take numerical advantage of the semi-analytical FT.

Conclusions
In this paper, we propose a rigorous approach to carry out the Fourier transform efficiently: the semi-analytical Fourier transform. The idea of this method is to extract the quadratic phase from the input field/spectrum with the help of a numerical evaluation technique, e.g., the Amplitude |E z (ρ)| Phase arg [E z (ρ)] Residual Phase arg U E z (ρ)  Levenberg-Marquardt method. Then, the semi-analytical FT can be used to replace the regular FFT of the fully sampled field by two FFTs of complex functions which require significantly fewer sampling points. The sampling issue is entirely dependent on the residual field so that in contrast to the regular FFT, the numerical effort is reduced significantly. The full theoretical derivation of the semi-analytical FT is presented. Additionally, several numerical examples are shown to demonstrate the potential of this approach.