1. Introduction

There is an increasing demand in the construction of modern optical systems for the inclusion of diffractive optical elements (DOEs), nanostructured components, scatterers, etc. alongside more traditional elements like lenses. In order for a simulation of such systems to account for some of the fundamental physical effects introduced by these components it is necessary to take an electromagnetic physical optics approach. Most traditional algorithms fulfilling these conditions, however, are burdened with impracticable memory loads and simulation times: this is the case of traditional rigorous Maxwell solvers (Finite Differences in Time Domain—FDTD, Fourier Modal Method—FMM, and others), since they require sub-wavelength sampling for the entire optical system to be simulated. The field-tracing approach [1] overcomes this drawback by tearing the system into sub-domains: each sub-domain can then be modeled by a different method, as long as all the methods yield the result in the form of an electromagnetic field. These field results can finally be interconnected into a global solution of the whole system. Depending on the specific characteristics of each element, different assumptions can be made, different approximations can be deemed sufficiently accurate, and, consequently, different techniques may be employed accordingly. For a system containing both gratings and lenses, for instance, we could employ FMM for the gratings and the Local-Plane-Interface Approximation (LPIA) [2] for the lenses. This tearing and interconnecting method allows for increased adaptability of the simulation when it comes to striking a balance between accuracy and computational efficiency.

It is worth noting, furthermore, that most fundamental, rigorous methods used in physical optics are defined in the spatial frequency domain: that is the case of the spectrum-of-plane-waves (SPW) propagation operator and of FMM [3, 4]. But this is not true of all methods: some of them (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]. $\tilde{\mathcal{P}}$ indicates the physical free-space propagation operator in the $\mathit{\kappa }$-domain. $\mathcal{B}$ and $\tilde{\mathcal{B}}$ represent the physical response of an optical element in the $\mathit{\rho }$- or in the $\mathit{\kappa }$- domain respectively. The entire process includes the handling of six electromagnetic field components.

Download Full Size | PDF

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 $\mathcal{O}\left(N^2\right)$ 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–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\pi$-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²) or, at most, separable problems (x², y²), 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.

2. Theory

2.1. Derivation of the semi-analytical Fourier transform

The electromagnetic field at a certain plane is given by

Let us next turn our attention to the phase term $\exp \left(\mathrm{i}{\gamma }_{\mathcal{l}}\left(\rho \right)\right)$ in Eq. (1). The relationship between the complex exponential function and the trigonometric functions, namely $\exp \left(\mathrm{i}{\gamma }_{\mathcal{l}}\left(\rho \right)\right)=\cos {\gamma }_{\mathcal{l}}\left(\rho \right)+\mathrm{i}\sin {\gamma }_{\mathcal{l}}\left(\rho \right)$, is well-known. The periodic nature of the trigonometric functions constrains all phase points to the range $\left(-\pi ,\pi \right]$, also referred to as a “$2\pi$-wrapped phase’’. It is important to note that a phase function that is “smooth’’ (${\mathcal{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 {\gamma }_{\mathcal{l}}\left(\rho \right)$ and $\sin {\gamma }_{\mathcal{l}}\left(\rho \right)$; the periods are inversely proportional to the slope of the unwrapped phase function, i.e., to the gradient of ${\gamma }_{\mathcal{l}}\left(\rho \right)$. The largest slope in the phase function restricts the periodwhich 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 ${\gamma }_{\mathcal{l}}\left(\rho \right)$. 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

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_{\mathcal{l}}\left(\rho \right)$. Although $U_{\mathcal{l}}\left(\rho \right)$ must still be sampled in the numerical operation, $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ has a less pronounced local gradient than the full phase function $\text{arg}\left[V_{\mathcal{l}}\left(\rho \right)\right]={\gamma }_{\mathcal{l}}\left(\rho \right)$. 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 \left(\mathrm{i}{\psi }_{\mathrm{q}}\left(\rho \right)\right)$. It means that we must handle this information analytically in the entire process.

Next, we attempt to calculate the Fourier transform of $V_{\mathcal{l}}\left(\rho \right)$ without sampling the quadratic phase term defined in Eq. (5).

Plugging Eq. (2) into the Fourier transform equation, we obtain

From the convolution theorem, we know that

In principle, the term ${\mathcal{F}}_{\kappa }\left[U_{\mathcal{l}}\left(\rho \right)\right]$ must be treated numerically. On the other hand, we find that the Fourier transform of the polynomial of second order ${\mathcal{F}}_{\kappa }\left[\exp \left(\mathrm{i}{\psi }_{\mathrm{q}}\left(\rho \right)\right)\right]$ can be solved with the integral formula

With the help of Eq. (8), we deduce the Fourier transform of the polynomial of second order

It should be remarked that when $C^2-4D_x{D}_y=0$, the two-dimensional phase function ${\psi }_{\mathrm{q}}\left(\rho \right)$ 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 ${\mathcal{F}}_{\kappa }\left[\exp \left(\mathrm{i}{\psi }_{\mathrm{q}}\left(\rho \right)\right)\right]$ in Eq. (7) and expanding the convolution, we get

And by defining the coordinate transform vector $\mathbf{\beta }=\left({\beta }_x,{\beta }_y\right)$

The interpretation of Eq. (15) as the product of a Fourier transform and a polynomial of the second order

2.2. 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

Using the same trick as Sec. 2.1, we can deduce the analytical representation of the inverse semi-analytical FT

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 ${\tilde{\varphi }}_{\mathrm{q}}\left(\kappa \right)$ in $\kappa$-domain is fully analytically handled in the entire process. The sampling effort depends on the two FFTs of the spectrum with less sampling, the ${\tilde{A}}_{\mathcal{l}}\left(\kappa \right)$ in Eq.(21).

2.3. 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 theresidual field $U_{\mathcal{l}}\left(\rho \right)$ or the residual spectrum ${\tilde{A}}_{\mathcal{l}}\left(\kappa \right)$. One way that works is to minimize the gradient of $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ or $\text{arg}\left[{\tilde{A}}_{\mathcal{l}}\left(\kappa \right)\right]$ over the given function support. For this purpose, we must find an effective method of determining the optimal factors of the quadratic phase ${\psi }_{\mathrm{q}}\left(\rho \right)$ or ${\tilde{\varphi }}_{\mathrm{q}}\left(\kappa \right)$.

There are two alternatives, analytical or numerical. The most widely recognized analytical approach is the Taylor expansion. For a given phase function ${\gamma }_{\mathcal{l}}\left(\rho \right)$, by using Taylor expansion up to the second order, we can get ${\psi }_{\mathrm{q}}\left(\rho \right)$ analytically

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 ${\gamma }_{\mathcal{l}}\left(\rho \right)$ 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 $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ 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, $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ 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 imperfectfunction value $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$. 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 $\nabla \left(\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]\right)$ as the merit function instead.

 

Fig. 2 Optimization of the quadratic phase ψ_q (ρ) for the general phase function ψ(ρ) by the Taylor expansion and Levenberg-Marquardt methods. (a) Using the two methods in the case of an aberrated phase function (spherical phase with coma and astigmatism aberration) in spatial domain. (b) Comparison of the fitting performance of the two approaches. In the range far from the central area which is the reference point of the expansion, both the function value and the gradient of the residual phase are higher with the Taylor approach than with the LMM.

Download Full Size | PDF

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.

3. Numerical examples

In the previous sections, the background knowledge and theoretical derivation of the semi-analytical 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.

3.1. 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.

 

Fig. 3 Sketch of the propagation of the Laguerre Gaussian 01-Mode in free space, and the corresponding field tracing diagram illustrating the same propagation process.

Download Full Size | PDF

 

Fig. 4 Field distributions at different planes for the simulation task: the propagation of the Laguerre Gaussian 01-Mode in free space. (a) Gaussian waist z = 0; (b) z = 10 μm ≈ 3.5 z_R; (c) z = 300 μm ≈ 100 z_R where z_R = 2.95 μm is the Rayleigh length of the Laguerre Gaussian mode. The residual phase $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ is the phase excluding the quadratic terms which is obtained from the semi-analytical FT.

Download Full Size | PDF

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 \left(\mathrm{i}{k}_z\mathrm{\Delta }z\right)$ on each plane wave, we obtain the spectrum at the target plane ${\tilde{V}}_{\perp }^{\text{out}}\left(\kappa \right)$. 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 $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$. 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).

 

Fig. 5 Analysis of the numerical effort: sampling numbers of FFT and semi-analytical FT versus propagation distance for the simulation task of the propagation of a Laguerre Gaussian 01-Mode. The ordinate tracks the required Nyquist sampling points of the two approaches in one dimension.

Download Full Size | PDF

Another effect we can observe in Fig. 4 is that as the propagating distance increases, the residual phase $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ 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 $\mathrm{\Delta }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_{\mathcal{l}}\left(\rho \right)$ 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.

Download Full Size | PDF

3.2. 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.

 

Fig. 7 Simulation results of an ideal spherical wave with spherical radius r = −3 mm going through its focus: amplitude, residual phase or complete phase on the plane (a) z = 0; (b) focal plane z = 3 mm; (c) z = 6 mm. Here, the residual phase $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ is the phase excluding the quadratic terms.

Download Full Size | PDF

Table 1. Aberration phase parameters of the test field.

View Table

 

Fig. 8 Field distribution at different planes of a convergent aberrated spherical wave propagating to its focal region: amplitude and residual phase or complete phase on the target plane (a) input plane z = 0; (b) z = 2.8 mm; (c) focal region z = 2.95 mm. The residual phase $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ is the phase excluding the quadratic terms. In the focal region, panel (c), since the field does not include a strong quadratic phase, the FFT is used and the full phase is shown.

Download Full Size | PDF

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 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 $\left|V_{\mathcal{l}}\left(\rho \right)\right|$. 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) and Fig. 7(c) are (221 × 221) and (251 × 251).

 

Fig. 9 Field distribution at different planes of a convergent aberrated spherical wave going through its focal region: amplitude and residual phase or complete phase on the target plane (d) z = 3.05 mm; (e) z = 3.5 mm; (f) z = 6 mm. The residual phase $\text{arg}\left[U_{\mathcal{l}}\left(\rho \right)\right]$ is the phase excluding the quadratic terms. In the focal region, panel (d), since the field does not include a strong quadratic phase, the FFT is used and the full phase is shown.

Download Full Size | PDF

 

Fig. 10 Numerical effort analysis: sampling numbers of the FFT and the semi-analytical FT versus propagation distance for the simulation task, in which an aberrated truncated spherical wave goes through its focus. The test fields are well-sampled in two dimensions, (N_x, N_y). The “Sampling Numbers” correspond to the larger of these two values.

Download Full Size | PDF

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 \left(\mathrm{i}{k}_z\mathrm{\Delta }z\right)$ 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² or y², 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.

3.3. 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 recordonly 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 $\kappa$-domain. Naturally, these algebraic equations connect the components.

 

Fig. 11 Simulation results for the experiment E_z calculation: amplitude |E_z (ρ)| and residual phase $\text{arg}\left[U_{E_z}\left(\rho \right)\right]$ or complete phase arg[E_z (ρ)] of the E_z component, at the target planes (a) z = 0 ; (b) z = 2.8 mm ; (c) z = 2.95 mm.

Download Full Size | PDF

 

Fig. 12 Simulation results of the experiment E_z calculation: Amplitude |E_z (ρ)| and residual phase $\text{arg}\left[U_{E_z}\left(\rho \right)\right]$ or complete phase arg[E_z (ρ)] of the E_z component, at the target planes (d) z = 3.05 mm; (e) z = 3.5 mm; (f) z = 6 mm.

Download Full Size | PDF

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.

4. 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 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.