Generalized far-field integral

Zongzhao Wang; Zongzhao Wang; Olga Baladron-Zorita; Olga Baladron-Zorita; Christian Hellmann; Frank Wyrowski

doi:10.1364/OE.414314

1. Introduction

In physical-optics modeling and design, the purpose is to have access to the expression of the electromagnetic field as it progresses through a system, before and after its interaction with the different optical components, or any type of structured media that is part of the system [1]. The problem of light propagation in free space (where, by “free space”, we understand any isotropic and homogeneous medium), is the most fundamental and widely used modeling scenario. For this subject of interest (namely, the mathematical problem of free-space propagation of electromagnetic fields), rigorous solutions to Maxwell’s equations have been found both in the spatial domain ($x$ domain) and the spatial-frequency domain ($k$ domain): the Rayleigh-Sommerfeld diffraction integral [2,3] and the method of the angular spectrum of plane waves (SPW) [4–6]. In principle, both of these methods can be applied in all kinds of propagation scenarios that match the above-given description of free space. However, from the point of view of their numerical implementation, both of them are integral operators. Due to the complexity of calculus, rigorous approaches always suffer from severe numerical challenges [7,8]. Therefore, several approximate algorithms have been developed to deal with the problem of propagation in free space when certain specific assumptions hold: for instance, the Debye integral is specially designed for the problem of focusing light in lens design [9,10].

The well-known far-field integral is used to calculate the diffraction pattern in the region far from the object plane [11]. This approximate method is applicable to propagation from the focal region to the far-field zone. In contrast to the rigorous SPW operator, the far-field integral reduces two Fourier integral operations to a single one, for which it is allowed to employ the fast Fourier transform (FFT) [12,13] technique in the numerical implementation. The standard far-field integral is essentially a product of the input field’s angular spectrum and a perfect spherical wave under the far-field approximation. Although the far-field integral provides significant numerical superiority compared against the rigorous methods, as soon as the angular spectrum has phase variations (which it has if the input field has aberrations), its scope of application would be constrained. Specifically, the propagation distance must tend to infinity, or the optical system must be aberration-free. In short, it is urgent and vital to overcome the limitations of far-field integral so as to achieve more flexibility in the treatment of physical optics modeling.

The Debye integral has been generalized to include any type of aberration and, at the same time, reduce the limitation of the Fresnel number in [14]. It is possible, following a similar approach, to also generalize the far-field integral to a more general concept and extend its application scope. In the present work, our discussion starts with a rigorous analysis of the propagation via the $k$ domain. With the help of the pointwise Fourier transform, an approximate approach inspired by the stationary phase theory which can reduce the Fourier transform integral to a pointwise operation [15], the generalized far-field integral is concluded. The full derivation process is presented in Section 2. As a particular case, we review and interpret the far-field integral in comparison with the generalized approach in Section 3. The theoretical discussion culminates with a derivation of how to compute the propagation of an electromagnetic field between arbitrarily oriented planes by the generalized far-field integral. The numerical performance of the proposed approach is analyzed along several numerical examples, as presented in Section 5.

2. Generalized far-field integral

All electromagnetic fields have six vectorial components, i.e., $\left \{E_x, E_y, E_z, H_x, H_y, H_z \right \}$. In homogeneous and isotropic media the propagation operator of the electromagnetic field is, simply, the same for all field components. In other words, the six field components are decoupled with respect to the propagation operator, so that it is enough to perform an analysis of a single field component: the results of the analysis can directly be extrapolated to all the other components. Thus, we describe the field with $V_\ell \!\left ({\boldsymbol {r}}\right )$, $\ell = 1, \ldots , 6$, where $V_\ell$ acts as shorthand for the six components of the electromagnetic field, and ${\boldsymbol {r}} = \left (x, y, z\right )$ denotes the three-dimensional position vector in the Cartesian $xyz$-coordinate system. We can assume the light field propagates along the $+z$ axis, $z= z_0$ being the position of the input plane $\textrm {P}^{\textrm {in}}$ where we investigate the field. The field in this plane is known, and will be given as a function of the transversal position vector ${\boldsymbol {\rho }} = {\boldsymbol {r}}_\perp = \left (x, y\right )$, i.e. $V_\ell \! \left ({\boldsymbol {\rho }}, z_0 \right )$. Considering propagation between two parallel planes, the light field at the target plane can be formulated by

(1)$$V_\ell \! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) = {\boldsymbol{\mathcal{P}}} V_\ell \! \left({\boldsymbol{\rho}}, z_0 \right)\textrm{,}$$

where ${\boldsymbol {\mathcal {P}}}$ indicates the free-space propagation operator in the $x$ domain, ${\boldsymbol {\rho }}^{\prime } = \left (x^{\prime }, y^{\prime }\right )$ represents the transversal position vector at the target plane, and $V_\ell \! \left ({\boldsymbol {\rho }}^{\prime }, z^{\prime } \right )$ represents the propagated field component at the target plane. Equation (1) establishes the problem statement. Any analysis of the free-space propagation will deal with the mathematical operations involved in realizing said equation. As mentioned above, there are various approaches, rigorous and approximated, to tackle this problem. In what follows, we derive the generalized far-field integral via a $k$-domain analysis of the propagation, which must also involve the Fourier transform, as it is a necessary tool to move back and forth between the $x$ and $k$ domains. We will, in fact, show how the behavior of the Fourier transform specifically is a key aspect of the whole process.

Any arbitrary field can be described as a superposition of plane waves with different propagation directions (or, equivalently, with different spatial frequencies). With the help of the Fourier transform, it is possible to retrieve this plane-wave description from the spatial expression of the field, and vice versa. Mathematically, the forward Fourier transform of $V \! \left ( {\boldsymbol {\rho }} \right )$ can be formulated by

(2)$$\tilde{V}\! \left( {\boldsymbol{\kappa}} \right) = {\boldsymbol{\mathcal{F}}}_{\kappa} \left[V \! \left( {\boldsymbol{\rho}} \right) \right] = \frac{1}{2\pi} \iint^{+\infty}_{-\infty} V \! \left( {\boldsymbol{\rho}} \right) \exp\left(-\textrm{i} {\boldsymbol{\rho}}\cdot {\boldsymbol{\kappa}}\right) \textrm{d}^{2} \rho \, ,$$

where ${\boldsymbol {\kappa }} = \left (k_x, k_y \right )$ is the spatial frequency vector and $\tilde {V}\! \left ( {\boldsymbol {\kappa }} \right )$ denotes the corresponding complex amplitude of the plane waves in the $k$ domain. The related inverse Fourier transform is given by

(3)$$V \! \left( {\boldsymbol{\rho}} \right) = {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \left[\tilde{V}\! \left( {\boldsymbol{\kappa}} \right)\right] = \frac{1}{2\pi} \iint^{+\infty}_{-\infty} \tilde{V}\! \left( {\boldsymbol{\kappa}} \right) \exp\left(\textrm{i} {\boldsymbol{\rho}}\cdot {\boldsymbol{\kappa}}\right) \textrm{d}^{2} \kappa \, .$$

Each point in the spatial frequency domain, each spatial frequency component, can be understood as an ideal plane wave, whose physical propagation follows the plane-wave ansatz. We can refer to this as the $k$-domain propagation kernel and denote it by $\tilde {{\boldsymbol {\mathcal {P}}}}$. Using the knowledge we have of ideal plane waves, the propagation kernel in the $k$ domain is simply a pointwise operation, i.e., a multiplication with an exponential phase term,

(4)$$\tilde{{\boldsymbol{\mathcal{P}}}} := \exp\! \left[\textrm{i} \check{k}_z \! \left({\boldsymbol{\kappa}}\right) \Delta z \right] = \exp \!\left(\textrm{i} \sqrt{k^{2}_0 \check{n}^{2} - k_x^{2} - k_y^{2}} \Delta z \right) \, ,$$

where $\Delta z = z ^{\prime } - z_0$ denotes the propagation distance along the $z$ axis, $\check {n}$ indicates the complex refractive index of the embedding medium and $k_0 = \frac {2 \pi }{\lambda }$ is the wavenumber of the harmonic field. Combining the forward and backward Fourier transforms with the propagation kernel in the $k$ domain, we can write down the compact mathematical expression of the SPW operator,

(5)$$V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) = {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \! \left\{ {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[V_\ell\! \left({\boldsymbol{\rho}}, z_0 \right)\right] \times \exp \!\left[\textrm{i} \check{k}_z \! \left({\boldsymbol{\kappa}}\right) \Delta z \right] \right\} \, \textrm{.}$$

Up to this point, we have not introduced any approximation: Eq. (5) can, in theory, be used to tackle any free-space propagation problem between two parallel planes. But let us know consider the numerical complexity of the operation. Apart from the pointwise operation of Eq. (4) in the $k$ domain, two complex Fourier integrals are involved in the calculation. Even though using the FFT (instead of a pure discrete Fourier transform, DFT) already reduces the computational effort from $\mathcal {O} \! \left (N^{2}\right )$ to $\mathcal {O} \! \left (N \log N \right )$ for a function with $N$ sampling points, we still face a serious numerical challenge because of the need to sample the $2\pi$-modulo phase, especially light fields presenting an intense wavefront phase.

In this work, we concentrate on the propagation scenario contemplated in the assumptions of the far-field integral, i.e. that the input and output planes are located in the focal region and in the far-field region, respectively. To illustrate this physical process, we present a sketch, as well as the corresponding field-tracing diagram, in Fig. 1. First of all, let us have a look at the forward Fourier-transform operation at the input plane. As the input plane is located in the focal region, the incident field does, per definition, not contain an intense wavefront phase. Therefore, the homeomorphic Fourier transform (HFT) is not allowed, as explained in more depth in a previous publication by the authors of this work [15]. Consequently, we cannot make any approximation on the forward Fourier transform operator, it must be an integral operation. The resulting spectrum can be written in terms of its amplitude and phase,

(6)$$\tilde{V}_{\ell}\!\left( {\boldsymbol{\kappa}}, z_0 \right) = {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[V_\ell\! \left({\boldsymbol{\rho}}, z_0 \right)\right] = \left|\tilde{V}_{\ell}\!\left( {\boldsymbol{\kappa}}, z_0 \right) \right| \exp\left\{ \mathrm{i} \arg \left[ \tilde{V}_{\ell}\!\left( {\boldsymbol{\kappa}}, z_0 \right) \right] \right\} \textrm{.}$$

Without loss of generality, we can reformulate Eq. (6) to yield

(7)$$\tilde{V}_{\ell}\!\left( {\boldsymbol{\kappa}}, z_0 \right) = \tilde{A}_\ell \!\left( {\boldsymbol{\kappa}} \right) \exp \! \left[ \mathrm{i} \tilde{\psi}^{\textrm{in}} \!\left( {\boldsymbol{\kappa}} \right) \right] \, \textrm{,}$$

where we have extracted the smooth wavefront phase $\tilde {\psi }^{\textrm {in}} \!\left ( {\boldsymbol {\kappa }} \right )$ and grouped the residual of the phase with the amplitude $\left |\tilde {V}_{\ell }\!\left ( {\boldsymbol {\kappa }}, z_0 \right ) \right |$ into $\tilde {A}_\ell \!\left ( {\boldsymbol {\kappa }} \right )$. Typically, the residual phase in $\tilde {A}_\ell \!\left ( {\boldsymbol {\kappa }} \right )$ is not unwrappable, e.g., the vortex phase. For a typical system with some aberrations, the extracted wavefront phase corresponds to the mapped aberration phase in the $k$ domain. Thus, considering the phase term of the propagation kernel in Eq. (4), the output wavefront phase in the $k$ domain can be written as

(8)$$\tilde{\psi}^{\textrm{out}} \! \left({\boldsymbol{\kappa}} \right) = \tilde{\psi}^{\textrm{in}} \! \left({\boldsymbol{\kappa}} \right) + \check{k}_z \! \left({\boldsymbol{\kappa}}\right) \Delta z \, \textrm{.}$$

In the case of a large propagation distance, $\check {k}_z \! \left ({\boldsymbol {\kappa }}\right ) \Delta z$ dominates the output wavefront phase, and $\tilde {\psi }^{\textrm {in}} \! \left ({\boldsymbol {\kappa }} \right )$ can be regarded as a correction term for the aberration phases. Then, inserting Eqs. (7)–(8) into Eq. (5) provides

(9)$$\begin{aligned} V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) &= {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \! \left\{ {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[V_\ell\! \left({\boldsymbol{\rho}}, z_0 \right)\right] \times \exp \!\left[\textrm{i} \check{k}_z \! \left({\boldsymbol{\kappa}}\right) \Delta z \right] \right\} \\ & \approx {\boldsymbol{\mathcal{F}}}^{\textrm{h}, -1}_{\kappa} \! \left\{ \tilde{A}_\ell \! \left({\boldsymbol{\kappa}} \right) \exp \! \left[\textrm{i} \tilde{\psi}^{\textrm{out}} \! \left({\boldsymbol{\kappa}} \right) \right] \right\} \, \textrm{.} \end{aligned}$$

Here, we can directly use the conclusion of the inverse homeomorphic Fourier transform, published in [15].

(10)$$V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) = \tilde{a}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}}^{\prime} \right) \right] \tilde{A}_{\ell}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}}^{\prime} \right) \right] \exp\!\left\{ \mathrm{i}\tilde{\psi}^{\textrm{out}}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}}^{\prime} \right)\right] + \mathrm{i}{\boldsymbol{\rho}}^{\prime}\cdot{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}^{\prime}\right)\right\} \textrm{,}$$

With the Jacobian determinant factor

(11)$$\tilde{a}\!\left({\boldsymbol{\kappa}}\right) = \left\{ \begin{array}{ccl} \sqrt{\frac{\mathrm{i}}{\tilde{\psi}^{\textrm{out}}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)}} \sqrt{-\frac{\mathrm{i}\tilde{\psi}^{\textrm{out}}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)}{\left[\tilde{\psi}^{\textrm{out}}_{k_xk_y}\left({\boldsymbol{\kappa}}\right)\right]^{2}-\tilde{\psi}^{\textrm{out}}_{k_xk_x}\left({\boldsymbol{\kappa}}\right)\,\tilde{\psi}^{\textrm{out}}_{k_yk_y}\left({\boldsymbol{\kappa}}\right)}} & , & \tilde{\psi}^{\textrm{out}}_{k_xk_x} \! \left({\boldsymbol{\kappa}}\right) \neq 0 \\ \frac{1}{\left| \tilde{\psi}^{\textrm{out}}_{k_xk_x}\left({\boldsymbol{\kappa}}\right) \right|} & , & \tilde{\psi}^{\textrm{out}}_{k_xk_x} \!\left({\boldsymbol{\kappa}}\right) = 0 \end{array} \right. \, \textrm{.}$$

Fig. 1. Schematic illustration of a divergent field propagating in free space into the far field, and the corresponding field-tracing diagram illustrating the same propagation process according to the analysis of the spectrum of plane waves (SPW). The incident beam is propagated from the input plane located at $z = z_0$ to the target plane, with propagation distance $\Delta z = z^{\prime } - z_0$.

Download Full Size | PDF

Substituting Eqs. (7)–(8) into Eq. (10), the expression can be reformulated to yield

(12)$$V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) = \tilde{a}\!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}}^{\prime} \right) \right] \exp\!\left\{ \mathrm{i} \check{k}_z \!\left[ {\boldsymbol{\kappa}}\!\left( {\boldsymbol{\rho}}^{\prime} \right)\right]\Delta z + \mathrm{i}{\boldsymbol{\rho}}^{\prime}\cdot{\boldsymbol{\kappa}}\!\left({\boldsymbol{\rho}}^{\prime}\right)\right\} {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[V_\ell\! \left({\boldsymbol{\rho}}, z_0 \right)\right]_{{\boldsymbol{\kappa}} \rightarrow {\boldsymbol{\rho}}^{\prime} } \, \textrm{,}$$

where the backward mapping ${\boldsymbol {\kappa }} = {\boldsymbol {\kappa }}\!\left ({\boldsymbol {\rho }}^{\prime }\right )$ is provided by the gradient of the output wavefront phase (following the method of stationary phase [15])

(13)$$\tilde{{\boldsymbol{\nabla}}}_{\perp} \tilde{\psi}^{\textrm{out}}\!\left({\boldsymbol{\kappa}}\right) = -{\boldsymbol{\rho}}^{\prime} \textrm{.}$$

In general, there is no analytical expression for the pointwise mapping relation ${\boldsymbol {\kappa }} = {\boldsymbol {\kappa }}\!\left ({\boldsymbol {\rho }}^{\prime }\right )$, which means that Eqs. (12)–(13) must be handled entirely in a numerical manner. More importantly, in most cases, this coordinate transformation would be a non-uniform mapping, i.e., the target grid of a uniform grid will become non-equidistant. But this should not pose a problem, since all related issues have been fully addressed in the implementation of the homeomorphic Fourier transform. Comparing the generalized far-field integral formula with the rigorous SPW operator of Eq. (5), we can see that using the inverse HFT reduces the computational complexity of the second Fourier integral from $\mathcal {O}\! \left (N \log N \right )$ to $\mathcal {O}\! \left (N \right )$. Additionally, avoiding the sampling of the "$2\pi$-modulo" of the wavefront phase at the target plane leads to a significant reduction in the number of sampling points $N$. Let us roughly summarize the numerical effort of the SPW operator versus the generalized far-field integral: $\mathcal {O}^{\textrm {SPW}} \sim \mathcal {O}\! \left (N^{\textrm {in}} \log N^{\textrm {in}} \right ) + \mathcal {O}\! \left (N^{\textrm {out}} \right ) + \mathcal {O}\! \left (N^{\textrm {out}} \log N^{\textrm {out}} \right )$ and $\mathcal {O} ^{\textrm {far-field}} \sim \mathcal {O}\! \left (N^{\textrm {in}} \log N^{\textrm {in}} \right ) + \mathcal {O}\! \left (N^{\textrm {in}} \right ) + \mathcal {O}\! \left (N^{\textrm {in}} \right )$, where $N^{\textrm {in}}$ and $N^{\textrm {out}}$ respectively denote the required sampling points for the input and output fields, and, $N^{\textrm {out}} \gg N^{\textrm {in}}$. Furthermore, unlike other diffraction integrals, the only approximation used in the derivation process is in the inverse Fourier transform. Thus, the proposed approach has a very flexible application scope. For example, the Fresnel integral makes an approximation in the propagation kernel that restricts its validity to paraxial fields. Also, we would like to emphasize that the inclusion of the aberrant phase $\Delta \tilde {\psi }^{\textrm {in}} \! \left ({\boldsymbol {\kappa }} \right )$ dramatically improves the accuracy of the inverse HFT. In practice, by examining the level of accuracy of the inverse HFT, a decision can be made whether to employ the proposed algorithm.

3. Standard far-field integral as a special case

As we mentioned in the last section, the mapping relation ${\boldsymbol {\kappa }} \rightarrow {\boldsymbol {\rho }}^{\prime }$ is available, in general, only numerically. In order to solve the diffraction integral analytically, we must simplify the experimental model and introduce more approximations. So, we would neglect the step of extracting the aberrant wavefront phase from the input spectrum $\tilde {\psi }^{\textrm {in}} \!\left ({\boldsymbol {\kappa }} \right )$ and assume the wavefront phase of the output spectrum in the $k$ domain is fully contributed by the propagation kernel. Thus, Eq. (8) becomes

(14)$$\tilde{\psi}^{\textrm{out}} \! \left({\boldsymbol{\kappa}}\right) = \check{k}_z \! \left({\boldsymbol{\kappa}}\right) \Delta z \, \textrm{.}$$

Plugging Eq. (14) into Eq. (9) leads to the formula

(15)$$V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) \approx {\boldsymbol{\mathcal{F}}}^{\textrm{h}, -1}_{\kappa} \! \left\{\tilde{V}_\ell\! \left({\boldsymbol{\kappa}}, z_0 \right) \exp \!\left[\textrm{i} \tilde{\psi}^{\textrm{out}}\left({\boldsymbol{\kappa}}\right) \right] \right\} \, \textrm{.}$$

Following Eq. (13), the bijective mapping relation of the $\tilde {\psi }^{\textrm {out}} \! \left ({\boldsymbol {\kappa }}\right )$ can be written in an explicit form,

(16)$${\boldsymbol{\kappa}}\! \left({\boldsymbol{\rho}}^{\prime} \right) = k_0 \check{n} \frac{{\boldsymbol{\rho}}^{\prime} }{R\!\left({\boldsymbol{\rho}}^{\prime}\right)} \textrm{,}$$

where $R\!\left ({\boldsymbol {\rho }}^{\prime }\right ) = \sqrt {||{\boldsymbol {\rho }}^{\prime }||^{2} + \Delta z^{2}}$. This relation also fits the geometrical analysis of the ideal spherical wave. For this mapping, the scaling factor (Jacobian determinant) can be evaluated by

(17)$$\tilde{a}\left[ {\boldsymbol{\kappa}} \left({\boldsymbol{\rho}}^{\prime} \right) \right] = -\textrm{i} \frac{k_0 \check{n} \Delta z}{R^{2}\!\left({\boldsymbol{\rho}}^{\prime} \right)} \, \textrm{,}$$

and the phase modulation term can be derived,

(18)$$\tilde{\psi}^{\textrm{out}} \! \left[{\boldsymbol{\kappa}}\left({\boldsymbol{\rho}}^{\prime} \right)\right] + {\boldsymbol{\kappa}}\left({\boldsymbol{\rho}}^{\prime} \right)\cdot {\boldsymbol{\rho}}^{\prime} = k_0 \check{n} R\!\left({\boldsymbol{\rho}}^{\prime}\right) \, \textrm{.}$$

According to the inverse homeomorphic Fourier transform formula, we obtain

(19)$$\begin{aligned} V_\ell\! \left({\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) & \approx \tilde{a}\left[ {\boldsymbol{\kappa}} \left({\boldsymbol{\rho}}^{\prime} \right) \right] \tilde{V}_\ell\! \left[{\boldsymbol{\kappa}} \left({\boldsymbol{\rho}}^{\prime}\right), z_0 \right] \exp \!\left\{\textrm{i} \tilde{\psi}^{\textrm{out}}\left[{\boldsymbol{\kappa}}\left({\boldsymbol{\rho}}^{\prime} \right)\right] + \textrm{i} {\boldsymbol{\kappa}}\left({\boldsymbol{\rho}}^{\prime}\right)\cdot {\boldsymbol{\rho}}^{\prime} \right\} \\ & = -\textrm{i} k_0 \check{n} \frac{\Delta z}{R\left({\boldsymbol{\rho}}^{\prime}\right)} \frac{ \exp \left[\textrm{i} k_0 \check{n} R \left({\boldsymbol{\rho}}^{\prime}\right) \right]}{R\left({\boldsymbol{\rho}}^{\prime}\right)} {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[V_\ell\! \left({\boldsymbol{\rho}}, z_0 \right)\right] _{{\boldsymbol{\kappa}} = - k_0 \check{n} \frac{{\boldsymbol{\rho}}^{\prime} }{R\left({\boldsymbol{\rho}}^{\prime}\right)}} \, \textrm{.} \end{aligned}$$

Equation (19) reveals that any light field component in the far-field zone turns into a spherical wave function. This spherical wave function is modulated by the Fourier transform of the input field. The modulation scales are not linear with respect to the transversal coordinates, but instead obey the mapping relation of Eq. (16). In contrast to the generalized far-field integral of Eq. (12), the far-field integral makes the additional approximation on the inverse homeomorphic Fourier transform step: it assumes that the mapping used in the computation of the inverse HFT follows that dictated by a spherical phase function. Therefore, in the case of an aberration-free system, both approaches would provide precisely the same result. But, for a system presenting aberrations, the inclusion of the aberration phase of Eq. (8) has a direct influence on the mapping and would significantly improve the accuracy of the inverse homeomorphic Fourier transform.

Furthermore, we would like to discuss the application scope of the far-field integral and that of the generalized approach. In the literature [16–18], the suitable propagation distances of several diffraction propagation approaches are demonstrated comprehensively. We would interpret the definition of the “far field” differently from the cited works, and discuss said definition and the scope of the proposed method from the perspective of the validity of the homeomorphic Fourier transform. Hence, we introduce the following two new notions:

• Homeomorphic field zone (HFZ): that in which the HFT is allowed and can provide an accurate result when compared with the rigorous FFT. The field presents a general homeomorphic wavefront phase.
• Far-field zone (FFZ): that in which the HFT is allowed and provides an accurate result when compared against the rigorous FFT, even with the additional restriction that the wavefront used in the computation of said HFT (i.e. the one that provides the mapping) be of spherical shape.

Hereby, the far-field zone indicates the application scope of the far-field integral, and the homeomorphic field zone corresponds to the generalized far-field integral. All fields present a spherical wavefront at an infinite distance from the source. It is evident from the above definitions that the far-field zone must be a subset of the homeomorphic field zone. Importantly, to examine the far-field zone and the homeomorphic field zone, we must compare FFT and HFT results in practical simulations. Considering the numerical difficulty, we recommend performing several 1D-cut tests instead of the full 2D-field check. It should be noted that the range of the zones would vary with the customized threshold, i.e., the deviation between FFT and HFT.

4. Propagation to arbitrarily oriented planes

So far, we have derived the mathematical expression of the far-field integral and the generalized far-field integral, for propagation between parallel planes. Equation (19) reveals that when the input spectrum is obtained, we can calculate the field value at any position in the far-field zone by using the inverse mapping as well as the interpolation technique. There is no doubt that the far-field integral can deal with the light field propagation between non-parallel planes. Thus, as a logical generalization of the far-field integral, it is meaningful that we demonstrate how to tackle the propagation to arbitrarily oriented planes by the generalized far-field integral.

Fig. 2. Schematic illustration of the propagation of a divergent wave to an arbitrarily oriented plane, and the corresponding $k$-domain analysis illustrated by the field tracing diagram. The input field ${\boldsymbol {V}}_\bot \! \left ({\boldsymbol {\rho }}, z_0 \right )$ is a divergent wave located in the focal region. The output field ${\boldsymbol {V}}^{\textrm {out}}_\bot \! \left ({\boldsymbol {\rho }}^{\textrm {out}}, z^{\textrm {out}} \right )$ is investigated on a tilted plane.

Download Full Size | PDF

For the problem of free-space propagation between two non-parallel planes, a rigorous physical-optics solution was presented by Zhang et al. [19]. This solution is based on the analysis of the angular spectrum of plane waves; the propagation operator $\tilde {{\boldsymbol {\mathcal {P}}}}$ in the $k$ domain still behaves in a pointwise manner, even between non-parallel planes [20], as shown in Fig. 2. Without showing the derivation process, we simply summarize $\tilde {{\boldsymbol {\mathcal {P}}}}$ by the following equation:

(20)$$\begin{aligned} \tilde{{\boldsymbol{\mathcal{P}}}} \! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] &= \tilde{J} \!\left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] \tilde{{\boldsymbol{\mathcal{B}}}} \! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] \tilde{{\boldsymbol{\mathcal{M}}}} \! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] \\ &= \left[ \begin{array}{ll} \tilde{J} & 0 \\ 0 & \tilde{J} \end{array}\right] \left[ \begin{array}{ll} \tilde{B} & 0 \\ 0 & \tilde{B} \end{array}\right] \left[ \begin{array}{ll} \tilde{M}_{k_x k_x} & \tilde{M}_{k_x k_y} \\ \tilde{M}_{k_y k_x} & \tilde{M}_{k_y k_y} \end{array}\right] \,\textrm{,} \end{aligned}$$

where ${\boldsymbol {\kappa }}^{\textrm {out}} = \left (k_x^{\textrm {out}}, k_y^{\textrm {out}}\right )$ indicates the spatial frequency vector on the tilted target plane, and ${\boldsymbol {\kappa }}^{\textrm {out}}\! \left ( {\boldsymbol {\kappa }}\right )$ follows the coordinate transformation denoting the pointwise mapping between the input spatial frequency vector and the output spatial frequency vector, $\tilde {J} \! \left ( {\boldsymbol {\kappa }}^{\textrm {out}} \right )$ is the Jacobian determinant of this mapping, $\tilde {{\boldsymbol {\mathcal {B}}}} \! \left ( {\boldsymbol {\kappa }}^{\textrm {out}} \right )$ is the propagation kernel, $\tilde {{\boldsymbol {\mathcal {M}}}} \! \left ( {\boldsymbol {\kappa }}^{\textrm {out}} \right )$ is the projection matrix of the field components.

It should be noted that because of the field-components projection matrix $\tilde {{\boldsymbol {\mathcal {M}}}}$ in Eq. (20), the vectorial field components are coupled. This means that the output field components cannot be calculated individually, as was previously the case. However, it can be concluded from Maxwell’s equations in homogeneous media that the other four field components can be calculated from $\left (E_x, E_y \right )$ on demand. As a consequence, combining the above with the forward Fourier transform at the input plane and the inverse Fourier transform operation at the target plane, the propagation between two non-parallel planes can be formulated by

(21)$${\boldsymbol{V}}^{\textrm{out}}_\bot ({\boldsymbol{\rho}}^{\textrm{out}}, z^{\textrm{out}} ) = {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa^{\textrm{out}}} \! \left\{ \tilde{{\boldsymbol{\mathcal{P}}}}\! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[{\boldsymbol{V}}_\bot \! \left({\boldsymbol{\rho}}, z_0 \right)\right] \right\} \, ,$$

where ${\boldsymbol {\mathcal {F}}}_{\kappa }$ denotes the Fourier transform operation from ${\boldsymbol {\rho }}$ to ${\boldsymbol {\kappa }}$, and ${\boldsymbol {\mathcal {F}}}^{-1}_{\kappa ^{\textrm {out}}}$ represents the inverse Fourier transform operation from ${\boldsymbol {\kappa }}^{\textrm {out}}$ to ${\boldsymbol {\rho }}^{\textrm {out}}$, and ${\boldsymbol {V}}_\bot = \left (E_x, E_y \right )$ indicates the transverse components of the electric field.

From the definition in Eq. (7), we know that the wavefront phase is independent of the index of the field components, and thus bears the same expression for all of them. Thus, inserting Eq. (7) into Eq. (21) leads to

(22)$${\boldsymbol{V}}^{\textrm{out}}_\bot ({\boldsymbol{\rho}}^{\textrm{out}}, z^{\textrm{out}} ) \approx {\boldsymbol{\mathcal{F}}}^{-1, \textrm{h}}_{\kappa^{\textrm{out}}} \! \left\{ \tilde{{\boldsymbol{A}}}^{\textrm{out}}_\bot \! \left( {\boldsymbol{\kappa}}^{\textrm{out}} \right) \exp \!\left[ \mathrm{i} \tilde{\psi}^{\textrm{out}} \!\left( {\boldsymbol{\kappa}}^{\textrm{out}} \right) \right] \right\} \, ,$$

with

(23)$$\tilde{{\boldsymbol{A}}}^{\textrm{out}}_\bot \! \left( {\boldsymbol{\kappa}}^{\textrm{out}} \right) = \tilde{J} \! \left( {\boldsymbol{\kappa}}^{\textrm{out}}\right) \tilde{{\boldsymbol{\mathcal{M}}}} \! \left( {\boldsymbol{\kappa}}^{\textrm{out}}\right) \tilde{{\boldsymbol{A}}}^{\textrm{in}}_\bot \!\left[ {\boldsymbol{\kappa}} \! \left({\boldsymbol{\kappa}}^{\textrm{out}}\right) \right]$$

and with

(24)$$\exp \!\left[ \mathrm{i} \tilde{\psi}^{\textrm{out}} \!\left( {\boldsymbol{\kappa}}^{\textrm{out}} \right) \right] = \tilde{B} \! \left( {\boldsymbol{\kappa}}^{\textrm{out}}\right) \exp \! \left\{ \mathrm{i} \tilde{\psi}^{\textrm{in}} \!\left[ {\boldsymbol{\kappa}} \! \left({\boldsymbol{\kappa}}^{\textrm{out}}\right) \right] \right\} \, \textrm{.}$$

It is worth mentioning that the wavefront phase of the output spectrum $\tilde {\psi }^{\textrm {out}} \!\left ( {\boldsymbol {\kappa }}^{\textrm {out}} \right )$, in Eq. (24), consists of not only the propagating kernel phase but also the input aberration phases. Afterward, by employing the inverse homeomorphic Fourier transform, we can conclude that

(25)$$\begin{aligned} {\boldsymbol{V}}^{\textrm{out}}_\bot ({\boldsymbol{\rho}}^{\textrm{out}}, z^{\textrm{out}} ) & \approx \tilde{a} \!\left[{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right) \right] \tilde{{\boldsymbol{A}}}^{\textrm{out}}_\bot \! \left[{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right)\right] \\ &\exp\left\{ \mathrm{i}\tilde{\psi}^{\textrm{out}} \!\left[{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right) \right] + \mathrm{i}{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right) \cdot {\boldsymbol{\rho}}^{\textrm{out}} \right\} \end{aligned}$$

where the bijective mapping, i.e., ${\boldsymbol {\kappa }}^{\textrm {out}} \leftrightarrow {\boldsymbol {\rho }}^{\textrm {out}}$, follows

(26)$$\tilde{{\boldsymbol{\nabla}}}_\bot \tilde{\psi}^{\textrm{out}} \!\left({\boldsymbol{\kappa}}^{\textrm{out}} \right) = - {\boldsymbol{\rho}}^{\textrm{out}} \textrm{.}$$

The amplitude scaling factor $\tilde {a} \!\left ({\boldsymbol {\kappa }}^{\textrm {out}} \right )$ is governed by the Jacobian determinant of the mapping from Eq. (26) so that, using the definition $\psi _{x_ix_j} \overset {\textrm {def}}{=} \frac {\partial \psi }{\partial x_i \partial x_j}$:

(27)$$\tilde{a}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right) = \left\{ \begin{array}{ccl} \sqrt{\frac{\mathrm{i}}{\tilde{\psi}^{\textrm{out}}_{k_x k_x}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right)}} \sqrt{-\frac{\mathrm{i}\tilde{\psi}^{\textrm{out}}_{k_x k_x}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right)} {\left[\tilde{\psi}^{\textrm{out}}_{k_x k_y}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right)\right]^{2}-\tilde{\psi}^{\textrm{out}}_{k_x k_x}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right)\,\psi^{\textrm{in}}_{k_y k_y}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right)}} & , & \tilde{\psi}^{\textrm{out}}_{k_x k_x}\!\left( {\boldsymbol{\kappa}}^{\textrm{out}} \right) \neq 0 \\ \frac{1}{\left| \tilde{\psi}^{\textrm{out}}_{k_x k_y}\!\left({\boldsymbol{\kappa}}^{\textrm{out}}\right) \right|} & , & \tilde{\psi}^{\textrm{out}}_{k_x k_x}\!\left( {\boldsymbol{\kappa}}^{\textrm{out}}\right) = 0 \end{array} \right. \textrm{.}$$

Inserting Eq. (20) into Eq. (25), the compact expression can be written as

(28)$$\begin{aligned} {\boldsymbol{V}}^{\textrm{out}}_\bot ({\boldsymbol{\rho}}^{\textrm{out}}, z^{\textrm{out}} ) &\approx \tilde{a} \!\left[{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right) \right] \exp\left[ \mathrm{i}{\boldsymbol{\kappa}}^{\textrm{out}} \left({\boldsymbol{\rho}}^{\textrm{out}}\right) \cdot {\boldsymbol{\rho}}^{\textrm{out}} \right] \\ &\left( \tilde{{\boldsymbol{\mathcal{P}}}}\! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] {\boldsymbol{\mathcal{F}}}_{\kappa} \!\left[ {\boldsymbol{V}}_\bot\! \left({\boldsymbol{\rho}}, z_0 \right)\right] \right)_{{\boldsymbol{\kappa}} \rightarrow {\boldsymbol{\kappa}}^{\textrm{out}} \rightarrow {\boldsymbol{\rho}}^{\textrm{out}} }\textrm{,} \end{aligned}$$

where ${\boldsymbol {\kappa }} \rightarrow {\boldsymbol {\kappa }}^{\textrm {out}} \rightarrow {\boldsymbol {\rho }}^{\textrm {out}}$ represents the pointwise mapping sequence from the input spatial frequency vector to the output spatial frequency vector, then to the position vector on the target plane.

In conclusion, we extend the generalized far-field integral to allow the calculation of the field distribution on an arbitrarily oriented plane. Comparing this to the situation of propagation between parallel planes, the extension does not consume more computational resources, since all additional operations in this process are still pointwise.

5. Numerical examples

In the physical-optics modeling and design software VirtualLab Fusion [21] we implemented the far-field integral, the generalized far-field integral, as well as the SPW operator. In this section, we build up a basic simulation setup with three freely controlled parameters, shown in Fig. 3. There is flexibility in selecting and configuring these three parameters in order to carry out the experiments. Based on such settings, we would make several numerical examples and demonstrate the level of accuracy and performance of the generalized far-field integral.

Fig. 3. Illustration of the simulation task for propagating a light field in free space, from the focal region to an arbitrarily oriented target plane. The test field is initialized on the plane $z = {-40}$ mm. The bottom left panel shows the amplitude distribution of the $E_x$ component. The phase on this plane consists of two parts, the fixed spherical phase $\psi ^{\textrm {sph}} \! \left ({\boldsymbol {\rho }}\right )$ with radius $R = {-40}$ mm and the varying aberration phase $\psi ^{\textrm {zer}} \! \left ({\boldsymbol {\rho }}\right )$. From this information it is of course possible to obtain the field at the input plane $z = {0}$ mm. After setting the three given controlled variables, we can start our investigation.

Download Full Size | PDF

First of all, let us have a look at the specification of the experimental setup. The electromagnetic field, with a wavelength of 32 nm, is initialized at a distance of $z = {-40}$ mm away from the input plane. It is linearly polarized in $E_x$ and its amplitude distribution $\left | V_\ell \! \left ({\boldsymbol {\rho }}, z^{-} \right )\right |$ is modulated as a photo of Carl Zeiss. The phase distribution is specified by

(29)$$\arg \left[ V_\ell \! \left( {\boldsymbol{\rho}} , z^{-}\right) \right] = \psi^{\textrm{sph}} \left({\boldsymbol{\rho}}\right) + \psi^{\textrm{zer}} \left({\boldsymbol{\rho}}\right) \, \textrm{,}$$

with

(30)$$\psi^{\textrm{zer}} \! \left({\boldsymbol{\rho}}\right) = k \sum_{m=0}^{M}\sum_{n=0}^{N}c^{m}_n Z^{m}_n\!\left(r, \theta\right) \textrm{,}$$

where $k=\frac {2\pi }{\lambda }\check {n}$, $\lambda$ being the wavelength, $r=\frac {\left |{\boldsymbol {\rho }}\right |}{\left |{\boldsymbol {\rho }}_\textrm {max}\right |}$ and $\theta = \arctan \left (\frac {y}{x}\right )$; $c^{m}_n$ denotes the coefficients of the corresponding Zernike term. In our configuration, the spherical phase $\psi ^{\textrm {sph}} \left ({\boldsymbol {\rho }}\right ) = -k_0 \check {n} \sqrt {{\boldsymbol {\rho }}^{2} + R^{2}}$ is fixed with a constant spherical radius of $R = {-40}$ mm while we set the aberration phase $\psi ^{\textrm {zer}} \! \left ({\boldsymbol {\rho }}\right )$ as a free variable. By modifying the type and weighting coefficients of the Zernike phase, we can control the aberration in this system. The complete light field should be numerically sampled on the initial $xy$ plane. Here, we employ a hybrid sampling strategy handling the residual field and the wavefront phase independently: the amplitude $\left | V_\ell \! \left ({\boldsymbol {\rho }}, z^{-} \right )\right |$ is sampled on an equidistant grid using sampling points of $\left (115 \times 121 \right )$ with sampling distances of $({86.9}\;\mu\textrm{m} \times {82.6}\;\mu\textrm{m})$; on the other hand, the phase term is sampled by $1091$ points in a non-equidistant manner. By performing some appropriate free-space propagation technique, e.g., the generalized Debye integral, we obtain the input field on the plane $z = {0}$ mm. Before we investigate the propagation from the focal region to the target plane, we need to set another two controlled parameters, the propagation distance $\Delta z$ and the orientation of the target plane $\left (\theta , \phi \right )$. Then, we can employ the generalized far-field integral, standard far-field integral, and the rigorous free-space propagation operator. To have a more in-depth look at the level of accuracy of the proposed approach, we introduce a relative error that expresses the deviation $\sigma$ of the approach under analysis with respect to the rigorous reference. The deviation is mathematically defined by

(31)$$\sigma := \frac{\sum_{x,y} \!\left| V^{\textrm{ref}} \!\left({\boldsymbol{\rho}}\right) - V^{\textrm{test}} \!\left({\boldsymbol{\rho}} \right) \right|^{2}}{\sum_{x,y} \!\left|V^{\textrm{ref}} \!\left({\boldsymbol{\rho}}\right) \right|^{2}} \, \textrm{,}$$

where $V^{\textrm {test}} \!\left ({\boldsymbol {\rho }}\right )$ denotes the result obtained by the generalized or the standard far-field integral respectively, and the rigorous reference $V^{\textrm {ref}} \!\left ({\boldsymbol {\rho }}\right )$ is typically provided by the SPW operator.

5.1 Accuracy vs. the single aberrant phase

In the first set of experiments, we would like to investigate the influence of the aberration phase on the level of accuracy for both the generalized and the standard far-field integral. Firstly, let us clarify the configuration of the experimental setup for this simulation task: (1) only consider propagation between parallel planes, i.e., the orientation of the target plane is set to $\left (\theta , \varphi \right ) = \left (0^{\circ }, 0^{\circ }\right )$; (2) fix the propagation distance as $\Delta z = {60}$ mm (the target plane is not located in the far-field region); (3) select some single Zernike phase to perform the comparison.

We selected three types of single Zernike aberrant phase and individually superimposed them onto the light field. The corresponding mathematical expression and scaling coefficients are described in Table 1. In the following tests, we firstly compared the generalized far-field integral result with the rigorous SPW solution. For all cases, the deviation between them is smaller than $0.01\%$. It means that the target plane is located in the so-called homeomorphic field zone, where the HFT is accurate enough. So, in the subsequent step, we concentrate on the deviation between the generalized and the standard far-field integral. In Fig. 4, the curves of how the standard deviation between the two approaches varies with the scaling coefficients of the aberrant phase are presented. We can see that all three types of aberrant phases reveal the same tendency: (1) For the aberration-free or those cases with very weak aberration, the standard far-field integral can also provide an accurate result. (2) As the scaling factor of the aberration phase grows, the deviation between the two approaches increases rapidly. These experimental behaviors reflect sufficiently the approximation used in the derivation of the two approaches. In the generalized far-field integral, not only the primary spherical phase but also all aberration phases are taken into account in the mapping of the inverse homeomorphic Fourier transform. But the standard far-field integral considers just the primary spherical phase, so the inverse homeomorphic Fourier transform fails fast in the situation of a strong aberrant phase.

Fig. 4. Comparison of accuracy between the generalized and the standard far-field integral for the simulation task presented in Section. 5.1. The abscissa records the scaling factor of the aberrant phase $\psi ^{\textrm {zer}} \! \left ({\boldsymbol {\rho }}\right )$, while the ordinate marks the deviation between the two approaches, with the generalized far-field integral as reference. It should be noted that in this experiment, the deviation between the generalized far-field integral result and the rigorous SPW technique is lower than 0.01%.

Download Full Size | PDF

Table 1. Simulation parameters of the aberrant phase $\psi ^{\textrm {zer}} \! \left ({\boldsymbol {\rho }}\right )$ for the numerical example presented in Section. 5.1.

View Table | View all tables in this article

5.2 Accuracy vs. propagation distance

In the previous example, we have investigated the influence of three single Zernike phases on the level of accuracy. In this second set of experiments, we would like to turn our attention to another parameter, i.e., the propagation distance. Similarly, let us clarify the configuration of the experimental setup for this simulation task: (1) only consider the propagation between parallel planes, i.e., the orientation of the target plane is set to $\left (\theta , \varphi \right ) = \left (0^{\circ }, 0^{\circ }\right )$; (2) specify the aberrant phase $\psi ^{\textrm {zer}} \! \left ({\boldsymbol {\rho }}\right )$ and fix it as a constant; (3) customize the value of propagation distance $\Delta z$.

The simulation parameters for the phase of the light field $\textrm {arg}\! \left | V_\ell \!\left ({\boldsymbol {{\rho }}}, z^{-}\right ) \right |$ on the initial plane $z = {-40}$ mm are presented in Table 2. The propagation distance is set in the range of $\left [{100}\;\textrm{mm}, {1000}\;\textrm{mm}\right ]$. The same as in the first experiment, for all propagation distances in the definition domain, the generalized far-field integral can predict the same result as the reference SPW operator. In contrast, the standard far-field integral results, are more or less distinct from the reference. To compare the accuracy of the generalized and the standard far-field integral, we calculated the deviation between the two approaches at different propagation distances. And we plotted the deviation curve as it varies with the propagation distance in Fig. 5. The experimental result coincides very well with our expectation that the deviation drops as the propagation distance increases.

Fig. 5. Comparison of the accuracy of the generalized far-field integral and the standard far-field integral for the numerical example presented in Section 5.2. The abscissa records the propagation distance along the $z$ axis, while the ordinate marks the deviation of the two methods, with the generalized far-field integral as reference. It should be noted that in this set of experiments, the deviation between the generalized far-field integral result and the rigorous SPW technique is always lower than 0.01%.

Download Full Size | PDF

Table 2. Simulation parameters of the primary spherical phase and the aberrant phase for the numerical example presented in Section. 5.2.

View Table | View all tables in this article

The physical reason for this result is that the propagation kernel will attenuate the wavefront phase $\tilde {\psi } \left ({\boldsymbol {\kappa }}\right )$ in Eq. (24) as the propagation distance becomes larger and larger. However, we cannot, in this case, neglect the influence of the aberrant phase in the homeomorphic Fourier transform. Thus, there is still about a $1\%$ difference in the situation of $\Delta = {1}\;\textrm{m}$.

5.3 Propagation to inclined planes

Fig. 6. Simulation results for the task of propagating a light field in free-space from the focal region to an arbitrarily oriented target plane. Panels (a)-(c) present the obtained amplitude distribution of the $E_x$ component at three tilted target planes along the $z$-axis.

Download Full Size | PDF

Last but not least, we would like to employ the proposed approach to deal with propagation between non-parallel planes. With the same phase parameters used in the second example, the light field on the input plane is configured. Then, we select three different titled target planes along the optical axis. All three positions are located in the homeomorphic field zone. To examine the accuracy of the generalized far-field integral, we also perform a rigorous propagation based on the angular spectrum analysis. The obtained simulation results are presented in Fig. 6. The deviation errors between two approaches are always lower than $0.01\%$.

6. Conclusion

In this work, we generalize the far-field integral to improve its accuracy and application scope. The mathematical derivation starts from the SPW analysis. By replacing the inverse Fourier integral with the pointwise inverse Fourier transform, the wavefront phase aberrations are included and handled numerically via a mapping concept. In the derivation process, the propagation to arbitrarily oriented planes is also taken into account. Furthermore, we compare the proposed method’s validity condition with the standard far-field integral according to the Fourier transform property and conclude the criteria to specify the far-field zone. Along with a series of simulation examples, the proposed method’s accuracy and computational speed are verified.

Funding

European Social Fund (2017 SDP 0049).

Acknowledgments

We gratefully acknowledge financial support by LightTrans GmbH and the funding from the European Social Fund (ESF) (Thüringen-Stipendium Plus: 2017 SDP 0049).

Disclosures

The authors declare no conflicts of interest.

References

1. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58(5-6), 449–466 (2011). [CrossRef]

2. A. Sommerfeld, “Lectures on theoretical physics,” Vol. I, (Academic, New York, 1964).

3. J. C. Heurtley, “Scalar Rayleigh-Sommerfeld and Kirchhoff diffraction integrals: a comparison of exact evaluations for axial points,” J. Opt. Soc. Am. 63(8), 1003–1008 (1973). [CrossRef]

4. É. Lalor, “Conditions for the validity of the angular spectrum of plane waves,” J. Opt. Soc. Am. 58(9), 1235–1237 (1968). [CrossRef]

5. P. Liu and B. Lü, “The vectorial angular-spectrum representation and Rayleigh Sommerfeld diffraction formulae,” Opt. Laser Technol. 39(4), 741–744 (2007). [CrossRef]

6. W. Zhang, H. Zhang, and G. Jin, “Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range,” Opt. Lett. 45(6), 1543–1546 (2020). [CrossRef]

7. L. Mandel and E. Wolf, Optical coherence and quantum optics (Cambridge university press, 1995).

8. A. Papoulis, Systems and transforms with applications in optics (McGraw-Hill, 1968).

9. P. Debye, “Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie,” Ann. Phys. 335(14), 755–776 (1909). [CrossRef]

10. C. J. Sheppard, “Limitations of the paraxial Debye approximation,” Opt. Lett. 38(7), 1074–1076 (2013). [CrossRef]

11. M. Born and E. Wolf, Principles of optics: electromagnetic theory of propagation, interference and diffraction of light (Elsevier, 2013).

12. J. W. Goodman, Introduction to Fourier optics (Roberts and Company Publishers, 2005).

13. J. Brüning and V. W. Guillemin, Mathematics Past and Present: Fourier Integral Operators (Springer, 1994).

14. Z. Wang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Generalized Debye integral,” Opt. Express 28(17), 24459–24470 (2020). [CrossRef]

15. Z. Wang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Theory and algorithm of the homeomorphic Fourier transform for optical simulations,” Opt. Express 28(7), 10552–10571 (2020). [CrossRef]

16. D. G. Voelz and M. C. Roggemann, “Digital simulation of scalar optical diffraction: revisiting chirp function sampling criteria and consequences,” Appl. Opt. 48(32), 6132–6142 (2009). [CrossRef]

17. K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” Opt. Express 17(22), 19662–19673 (2009). [CrossRef]

18. W. Zhang, H. Zhang, C. J. Sheppard, and G. Jin, “Analysis of numerical diffraction calculation methods: from the perspective of phase space optics and the sampling theorem,” J. Opt. Soc. Am. 37(11), 1748–1766 (2020). [CrossRef]

19. S. Zhang, D. Asoubar, C. Hellmann, and F. Wyrowski, “Propagation of electromagnetic fields between non-parallel planes: a fully vectorial formulation and an efficient implementation,” Appl. Opt. 55(3), 529–538 (2016). [CrossRef]

20. K. Matsushima, H. Schimmel, and F. Wyrowski, “Fast calculation method for optical diffraction on tilted planes by use of the angular spectrum of plane waves,” J. Opt. Soc. Am. 20(9), 1755–1762 (2003). [CrossRef]

21. Physical optics simulation and design software, “Wyrowski VirtualLab Fusion”, developed by Wyrowski Photonics GmbH, distributed by LightTrans International UG, Jena, Germany.

Name & Type	Mathematical Expression	Value
Spherical aberration	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = 0 \sim 10 λ$
Tertiary astigmatism y	$Z_{6}^{- 2} = \sqrt{14} (15 r^{6} - 20 r^{4} + 6 r^{2}) \sin 2 θ$	$c_{2}^{- 2} = 0 \sim 10 λ$
Coma X	$Z_{3}^{1} = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 0 \sim 10 λ$

Name & Type	Mathematical Expression	Value
Spherical phase	$ψ^{sph} (ρ) = - k_{0} \overset{ˇ}{n} \sqrt{ρ^{2} + R^{2}}$	$R = - 40$ mm
Astigmatism Y	$Z_{2}^{- 2} (r, θ) = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{- 2} = 5 λ$
Spherical aberration	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = - 6 λ$
Trefoil Y	$Z_{3}^{- 3} (r, θ) = \sqrt{8} r^{3} \sin 3 θ$	$c_{3}^{- 3} = 8 λ$
Tetrafoil X	$Z_{4}^{4} (r, θ) = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{4} = - 8 λ$
Tertiary Coma X	$Z_{7}^{1} (r, θ) = 4 (35 r^{7} - 60 r^{5} + 30 r^{3} - 4 r) \cos θ$	$c_{7}^{1} = 10 λ$

Name & Type	Mathematical Expression	Value
Spherical aberration	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = 0 \sim 10 λ$
Tertiary astigmatism y	$Z_{6}^{- 2} = \sqrt{14} (15 r^{6} - 20 r^{4} + 6 r^{2}) \sin 2 θ$	$c_{2}^{- 2} = 0 \sim 10 λ$
Coma X	$Z_{3}^{1} = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 0 \sim 10 λ$

Name & Type	Mathematical Expression	Value
Spherical phase	$ψ^{sph} (ρ) = - k_{0} \overset{ˇ}{n} \sqrt{ρ^{2} + R^{2}}$	$R = - 40$ mm
Astigmatism Y	$Z_{2}^{- 2} (r, θ) = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{- 2} = 5 λ$
Spherical aberration	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = - 6 λ$
Trefoil Y	$Z_{3}^{- 3} (r, θ) = \sqrt{8} r^{3} \sin 3 θ$	$c_{3}^{- 3} = 8 λ$
Tetrafoil X	$Z_{4}^{4} (r, θ) = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{4} = - 8 λ$
Tertiary Coma X	$Z_{7}^{1} (r, θ) = 4 (35 r^{7} - 60 r^{5} + 30 r^{3} - 4 r) \cos θ$	$c_{7}^{1} = 10 λ$

Generalized far-field integral

Abstract

1. Introduction

2. Generalized far-field integral

3. Standard far-field integral as a special case

4. Propagation to arbitrarily oriented planes

5. Numerical examples

5.1 Accuracy vs. the single aberrant phase

5.2 Accuracy vs. propagation distance

5.3 Propagation to inclined planes

6. Conclusion

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (6)

Tables (2)

Equations (31)

Optics Express