Generalized Debye integral

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

doi:10.1364/OE.397010

1. Introduction

The study of diffraction theory for imaging systems is one of the most crucial parts of physical-optics modeling and design. From a physical-optics perspective, it is a problem of light propagation in free space (where, by “free space”, we understand any isotropic and homogeneous medium) from the exit pupil to the image plane. Even though there are several rigorous solutions, e.g., the spectrum-of-plane-waves (SPW) analysis [1] and the Rayleigh-Sommerfeld integral [2], solving this problem in a fast and accurate manner is still a big challenge. This is because rigorous approaches suffer from high numerical effort when it comes to propagating arbitrary fields. For instance, the Rayleigh-Sommerfeld integral is a complex-valued integral operator, which presents a complexity of $\mathcal {O} \!\left (N^{2}\right )$ for a function with $N$ samples. The SPW operator involves two Fourier transform integrals but enables fast evaluation by FFT that reduces the computational effort to $\mathcal {O} \!\left (N \log N\right )$. Importantly, for both approaches, in the case of an intensive wavefront phase, the required sampling points $N$ are massive to resolving the “$2\pi$-modulo”.

In 1909 Debye published his work to demonstrate the general features of a diffracted field near its focus [3]. When comparing the Debye and the rigorous Rayleigh-Sommerfeld integrals, we can see how the Debye approximation reduces a double integral to a single one, and in the numerical implementation of this single integral, the required sampling points are far fewer. Thus, the Debye integral provides a considerable reduction in the computational effort [4–6]. However, the Debye approximation breaks down for systems with relatively small Fresnel numbers. And, even though the Debye integral can predict the electromagnetic field distribution around the focal plane, it is only strictly valid for observation planes close to the focus [7–10]. The Fresnel number limitation is linked to the presence of aberrations in the sense that, the smaller the Fresnel number, the less aberration is allowed when using the standard Debye integral, if the accuracy of the operation is to remain constant. Different authors have attempted during the last decades to generalize the scope of application of the Debye integral with various algorithms by investigating special types of aberrations [11–13]. In this work, we generalize the Debye integral in order to include any type of aberration and, at the same time, reduce the limitation of the Fresnel number. The solution provides a fully vectorial method to tackle the fast calculation of the electromagnetic field into the focal region of a general incident convergent field. In Section 2., the theoretical derivation of the generalized Debye integral is presented. We start from the rigorous angular-spectrum-of-plane-waves (SPW) approach and then replace the rigorous forward Fourier transform by a homeomorphic Fourier transform [14–16]. We show that the standard Debye integral follows from the generalized version if the wavefront phase is restricted to adopt exclusively a spherical shape. In Section 3. we present several numerical examples to demonstrate the improved accuracy of the proposed method and its application potential.

2. Theory

2.1 Derivation of the generalized Debye integral

It is well known that the Helmholtz equation follows Maxwell’s equations in the homogeneous media. From the Helmholtz equation, we can conclude that all field components can be propagated independently by the same operator. In other words, the six field components decouple with respect to the propagation operator. Please note that in some focusing scenarios, the resulting field components are coupled [17,18]. It is because, in their modeling, not only the free space propagation but also the lens component are taken into account. For example, the Richards-Wolf integral combines the diffraction integral with an ideal lens model [19]. In contrast, our analysis concentrates on the free space propagation and assumes the input field is given on the plane behind the last lens. Thus, it is possible to formulate a fully vectorial solution in terms of $V_\ell \!\left ( {\boldsymbol {r}} \right )$, with $\ell = 1, \ldots , 6$, where $V_\ell$ acts as shorthand to describe the six components of a harmonic electromagnetic field $\left \{ E_x, E_y, E_z, H_x, H_y, H_z \right \}$. In our notation, ${\boldsymbol {r}} = \left ( x, y, z \right )$ and ${\boldsymbol {\rho }} = {\boldsymbol {r}}_\perp = \left ( x, y \right )$ denote the three-dimensional position vector and its projection onto the transversal plane respectively.

Let us consider a convergent incident light field which is propagating in the positive $z$ direction towards the target plane $z = z^{\prime }$ located in the focal region of the field. The situation is illustrated in Fig. 1. In lens design, $z = z_0$ might be the position of the exit pupil and $\Delta z =f$ the back focal length. As any electromagnetic field component will, in general, be complex-valued, we can write the field at $z_0$, $V_{\ell }\!\left ( {\boldsymbol {\rho }}, z_0 \right )$, in terms of its amplitude and phase,

(1)$$V_{\ell}\!\left( {\boldsymbol{\rho}}, z_0 \right) = \left|V_{\ell} \! \left( {\boldsymbol{\rho}}, z_0 \right) \right| \exp\left\{ \mathrm{i} \arg \left[ V_{\ell} \! \left( {\boldsymbol{\rho}}, z_0 \right) \right] \right\} \textrm{.}$$

In general, the field at $z_0$ possesses a smooth wavefront which is common to all field components. For instance, for an aberration-free focusing system, the wavefront phase would be spherical without any additional residual phase. In general, we can reformulate Eq. (1) to yield

(2)$$V_{\ell}\! \left( {\boldsymbol{\rho}}, z_0 \right) = U_{\ell}\!\left( {\boldsymbol{\rho}} \right) \exp\left[ \mathrm{i}\psi^{\textrm{in}}\!\left( {\boldsymbol{\rho}} \right) \right] \textrm{,}$$

where we have extracted the smooth wavefront phase $\psi ^{\textrm {in}} \!\left ( {\boldsymbol {\rho }} \right )$ and grouped the residual of the phase alongside the amplitude $\left | V_{\ell } \!\left ( {\boldsymbol {\rho }}, z_0 \right ) \right |$ into $U_{\ell }\!\left ( {\boldsymbol {\rho }} \right )$. It should be noted that the choice of function $\psi ^{\textrm {in}} \!\left ( {\boldsymbol {\rho }} \right )$ is arbitrary and Eq. (2) entails no approximation. It merely constitutes an alternative way of expressing $V_{\ell } \! \left ( {\boldsymbol {\rho }}, z_0 \right )$.

Fig. 1. Schematic illustration of a convergent wave in free space, 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 propagating distance $\Delta z = z^{\prime } - z_0$. Please note that we do not restrict the input field and, therefore, the focus can also be located off-axis, without affecting the theory and the results presented in that paper.

Download Full Size | PDF

Next, we need to consider how to tackle the problem of light propagation in free space. A recognized rigorous method is the analysis of the spectrum of plane waves (SPW). In this approach, firstly, the input field $V_{\ell } \!\left ( {\boldsymbol {\rho }}, z_0 \right )$ is decomposed into its plane-wave components by a Fourier-transform operation. Then, for each plane wave, we know how it transforms as it propagates from the input plane to the target plane. Finally, via an inverse Fourier-transform operation, we recombine all the plane waves and obtain the output field $V_{\ell }\! \left ( {\boldsymbol {\rho }}^{\prime }, z^{\prime } \right )$. Mathematically, we can write this process as follows:

(3)$$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[ \mathrm{i} \check{k}_z \!\left( {\boldsymbol{\kappa}}\right) \Delta z \right] \right\} \textrm{,}$$

where $\check {k}_z \! \left ( {\boldsymbol {\kappa }}\right ) = \sqrt { k_0^{2} \check {n}^{2} - {\boldsymbol {\kappa }}^{2}}$ is the $z$ component of the wave-number in the embedding medium, ${\boldsymbol {\kappa }} = \left (k_x, k_y \right )$ is the transversal wave vector, $\Delta z$ is the distance from the input plane to the target plane along the $z$ axis and ${\boldsymbol {\mathcal {F}}}_{\kappa }$ represents the Fourier-transform operation from ${\boldsymbol {\rho }}$ to ${\boldsymbol {\kappa }}$,

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

where $\tilde {V}_{\ell }\! \left ( {\boldsymbol {\kappa }}, z_0 \right )$ denotes the resulting spectrum on the input plane. This is a well-known result and provides a rigorous and fully vectorial method for the propagation of electromagnetic fields in free space. However, in practice, it suffers from a serious numerical challenge because of the need to sample the term $\exp \left [ \mathrm {i}\psi ^{\textrm {in}}\!\left ( {\boldsymbol {\rho }} \right ) \right ] = \cos [\psi ^{\textrm {in}}\!\left ( {\boldsymbol {\rho }} \right )] + \mathrm {i}\sin [\psi ^{\textrm {in}}\!\left ( {\boldsymbol {\rho }} \right )]$ to enable the calculation of the Fourier transform by the FFT. In this present work, however, we intend to concentrate on systems with low Fresnel numbers in combination with stronger aberrations. These aberrations might appear in a lens design optimization or tolerancing process, or as a result of the intentional use of aberrations for shaping the intensity distribution in the focal region. In such situations the Fourier transform can be replaced by a very fast, approximated, pointwise operation: the homeomorphic Fourier transform [14]. Let us use this method for the first (in order of application in the algorithm, which corresponds to the innermost one) Fourier-transform operator in Eq. (3) and substitute Eq. (2) into Eq. (4) so that we obtain

(5)$$\tilde{V}_{\ell} \!\left( {\boldsymbol{\kappa}}, z_0 \right) = \frac{1}{2 \pi} \iint^{+\infty}_{-\infty} U_{\ell}\!\left( {\boldsymbol{\rho}} \right) \exp\left[ \mathrm{i}\psi^{\textrm{in}}\!\left( {\boldsymbol{\rho}} \right) -\mathrm{i} {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}} \right] \mathrm{d}^{2}\rho \textrm{.}$$

The idea of the homeomorphic Fourier transform is quite straightforward: we consider a situation in which the Fourier integral operator reduces mathematically to a pointwise operation which is governed by the mapping

(6)$${\boldsymbol{\nabla}}_{\!\perp}\psi^{\textrm{in}}\!\left({\boldsymbol{\rho}}\right) = {\boldsymbol{\kappa}}({\boldsymbol{\rho}}) \textrm{.}$$

This behavior of the Fourier transform was discussed by Bryngdahl in the context of geometric transformations in optics [20] and is based on the method of stationary phase [21] applied to the Fourier integral. Using the mapping from Eq. (6) allows us to reformulate the Fourier integral as a pointwise operation. The application characteristics and the validity of the resulting homeomorphic Fourier transform (HFT) have been demonstrated in Wang et al. [14]. In short, to use the HFT, the working field must fulfill two preconditions: the first requires the part of the field which remains after extraction of the wavefront, $U_{\ell }\!\left ( {\boldsymbol {\rho }} \right )$, to be slowly varying. The other condition is a constraint on the extracted wavefront phase – the mapping from the spatial domain to the spatial frequency domain must be a bijective mapping, i.e., ${\boldsymbol {\rho }} \leftrightarrow {\boldsymbol {\kappa }}$ must be a homeomorphism. When both conditions are fulfilled or, in other words, when the HFT is allowed for the given input field, we can employ it to deal with Eq. (5). The resulting spectrum is then written as

(7)$$\begin{aligned} \tilde{V}_{\ell}\! \left( {\boldsymbol{\kappa}}, z_0 \right) & \approx a\!\left[ {\boldsymbol{\rho}} \left( {\boldsymbol{\kappa}} \right) \right] U_{\!\ell}\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] \exp\left\{ \mathrm{i} \psi^{\textrm{in}}\! \left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] - \mathrm{i} {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}}\!\left( {\boldsymbol{\kappa}} \right) \right\} \\ & = \tilde{A}\! \left({\boldsymbol{\kappa}}\right) \exp\left[ \mathrm{i} \tilde{\psi}^{\textrm{in}}\! \left( {\boldsymbol{\kappa}} \right) \right] \end{aligned} \textrm{.}$$

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

(8)$$a\!\left({\boldsymbol{\rho}}\right) = \left\{ \begin{array}{ccl} \sqrt{\frac{\mathrm{i}}{\psi^{\textrm{in}}_{xx}\!\left({\boldsymbol{\rho}}\right)}} \sqrt{-\frac{\mathrm{i}\psi^{\textrm{in}}_{xx}\!\left({\boldsymbol{\rho}}\right)} {\left[\psi^{\textrm{in}}_{xy}\!\left({\boldsymbol{\rho}}\right)\right]^{2}-\psi^{\textrm{in}}_{xx}\!\left({\boldsymbol{\rho}}\right)\,\psi^{\textrm{in}}_{yy}\!\left({\boldsymbol{\rho}}\right)}} & , & \psi^{\textrm{in}}_{xx}\!\left( {\boldsymbol{\rho}} \right) \neq 0 \\ \frac{1}{\left| \psi^{\textrm{in}}_{xy}\!\left({\boldsymbol{\rho}}\right) \right|} & , & \psi^{\textrm{in}}_{xx}\!\left( {\boldsymbol{\rho}} \right) = 0 \end{array} \right. \textrm{.}$$

In general, the mapping $\rho = \rho ({\boldsymbol{\kappa}})$ is known only numerically, and although this requires handling some gridless data in the algorithm, it should not pose a problem when performing the HFT in practice. By introducing appropriate error estimations within the algorithm, we can control the level of accuracy of replacing the FFT by the HFT. This is part of the mathematical concept of the HFT and does not rely on decisions of the physical model. To finalize the derivation of the generalized Debye integral, plugging Eq. (7) into Eq. (3) leads to the formula

(9)$$\begin{aligned} V_{\ell}\! \left( {\boldsymbol{\rho}}^{\prime}, z^{\prime} \right) &= {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \left\{ \tilde{V}_{\ell} \!\left( {\boldsymbol{\kappa}}, z_0 \right) \exp \left[ \mathrm{i} \check{k}_z \!\left( {\boldsymbol{\kappa}}\right) \Delta z \right] \right\} \\ &= {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa}\! \left\{ \tilde{A}\! \left({\boldsymbol{\kappa}}\right) \exp\left[ \mathrm{i} \tilde{\psi}^{\textrm{in}}\! \left( {\boldsymbol{\kappa}} \right) + \mathrm{i} \check{k}_z\! \left( {\boldsymbol{\kappa}}\right) \Delta z \right] \right\} \\ &= {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \!\left( a\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] U_{\!\ell}\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] \exp\!\left\{ \mathrm{i} \psi^{\textrm{in}}\! \left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] - \mathrm{i} {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}}\!\left( {\boldsymbol{\kappa}} \right) + \mathrm{i} \check{k}_z \!\left( {\boldsymbol{\kappa}}\right) \Delta z \right\} \right) \textrm{,} \end{aligned}$$

where ${\boldsymbol {\mathcal {F}}}^{-1}_{\kappa }$ indicates the inverse Fourier-transform operation from ${\boldsymbol {\kappa }}$ to ${\boldsymbol {\rho }}^{\prime }$ at the target plane. Equation (9) constitutes the major result of this work and expresses the generalized Debye integral with the use of the Fourier integral operator ${\boldsymbol {\mathcal {F}}}$ for the sake of compactness. As mentioned before, the mapping relation $\rho ({\boldsymbol{\kappa}})$ is in general available only numerically. Consequently, it is not possible to resolve the integral analytically.

Comparing Eq. (9) with the initial SPW formula of Eq. (3), we can see that two FT integrals are reduced to one FT operation. More importantly, according to the homeomorphic Fourier transform, the wavefront phase $\tilde {\psi }^{\textrm {in}}\! \left ( {\boldsymbol {\kappa }} \right )$ can be treated in an isolated manner during the process. This means that up until the last inverse Fourier-transform operation, it is possible to use different approaches to sample the wavefront phase $\tilde {\psi }^{\textrm {in}}\! \left ( {\boldsymbol {\kappa }} \right )$, instead of in the form of a “$2\pi$-modulo” phase. In the end, since the high-frequency phase component of $\tilde {\psi }^{\textrm {in}}\! \left ( {\boldsymbol {\kappa }} \right )$ would be partly compensated by the propagating kernel $k_z({\boldsymbol {\kappa }})\Delta z$, we can select the inverse FFT technique to carry out the inverse Fourier transform with low sampling effort.

It is worthy of mentioning that the only approximation used in the derivation is on the Fourier transform. The proposed approach has proven to be very powerful for systems with low Fresnel numbers. However, if the system becomes even more paraxial so that the HFT is not accurate enough, we recommend replacing the HFT by the semi-analytical Fourier transform (SFT) [22]. It is a rigorous Fourier transform technique that can also significantly reduce the sampling effort of the wavefront phase.

2.2 Standard Debye integral as a special case

For a spherical wavefront phase $\psi ^{\textrm {in}}$ with radius of curvature $R$, i.e.

(10)$$\psi^{\textrm{in}}(\rho) = \psi^{\textrm{sph}}(\rho) = \textrm{sign}\left(R\right) k_0 n \sqrt{{\boldsymbol{\rho}}^{2} + R^{2}} \textrm{,}$$

the mapping in Eq. (6) and the factor $a(\rho({\boldsymbol{\kappa}}))$ of Eq. (8) can be evaluated analytically. Here, $n$ is the real part of the refractive index. And, in the considered cases, $\psi ^{\textrm {in}}$ is a convergent spherical phase, so the radius of curvature is negative, i.e., $R < 0$. From Eq. (6) follow for the $\psi ^{\textrm {sph}}$ of Eq. (10) the mapping equations

(11)$$\begin{cases} \, {\boldsymbol{\kappa}}({\boldsymbol{\rho}}) = - k_0 n \frac{{\boldsymbol{\rho}}}{\sqrt{{\boldsymbol{\rho}}^{2} + R^{2}}} \\ \, {\boldsymbol{\rho}}({\boldsymbol{\kappa}}) = R \frac{{\boldsymbol{\kappa}}}{ k_z \!\left( {\boldsymbol{\kappa}}\right)} \end{cases}$$

in explicit form. Here, $k_z \!\left ( {\boldsymbol {\kappa }}\right )$ is defined as $k_z \!\left ( {\boldsymbol {\kappa }}\right ) = \sqrt {k_0^{2} n^{2} - {\boldsymbol {\kappa }}^{2}}$. The Jacobian determinant can be evaluated for this mapping and we obtain the scaling factor

(12)$$a^{\textrm{sph}}\big(\rho({\boldsymbol{\kappa}})\big) = \textrm{i} \frac{k_0 n R}{ k^{2}_z({\boldsymbol{\kappa}})}$$

in a compact and analytical form. The expression $\psi ^{\textrm {in}}\big (\rho({\boldsymbol{\kappa}})\big ) - {\boldsymbol {\kappa }} \cdot {\boldsymbol {\rho }}\!\left ( {\boldsymbol {\kappa }} \right ) + \check {k}_z({\boldsymbol{\kappa}})\Delta z$ in Eq. (9) reduces to $\check {k}_z({\boldsymbol{\kappa}})\Delta z + k_z({\boldsymbol{\kappa}}) R$ for the spherical wavefront phase by using the mapping relations of Eq. (12). Inserting these results into Eq. (9) provides

(13)$$V_{\ell}(\rho,z) = \mathrm{i} {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa} \!\left\{ \frac{k_0 n R}{ k^{2}_z \!\left( {\boldsymbol{\kappa}}\right)} U_{\ell}\big({\boldsymbol{\rho}}({\boldsymbol{\kappa}})\big) \exp \left[ \mathrm{i} \check{k}_z\! \left({\boldsymbol{\kappa}}\right) \Delta z + k_z\! \left({\boldsymbol{\kappa}}\right) R \right] \right\} \, \textrm{,}$$

with $\rho({\boldsymbol{\kappa}})$ from Eq. (11). This integral formula constitutes the special case of Eq. (9) which emerges when we select a spherical wavefront phase. The complex amplitude $U_\ell$ comprises magnitude and phase variations, including aberrations.

To recognize that Eq. (13) is identical with the standard Debye integral, we express $U$ as the magnitude of a spherical wave component:

(14)$$U_{\ell}\!\left( {\boldsymbol{\rho}} \right) = \frac{T_{\ell}(\rho)}{r} \textrm{,}$$

with $r=\sqrt {{\boldsymbol {\rho }}^{2} + R^{2}}$. The magnitude $|T_{\ell }(\rho)|$ expresses apodization and vignetting effects in lens design, whereas $\arg [T_{\ell }(\rho)]$ comprises the aberrations. With the mapping from Eq. (11), $r=\sqrt {{\boldsymbol {\rho }}^{2} + f^{2}}$ can be reformulated to yield $1/r=k_z/(|R| k_0n)$. Inserting $U_\ell (\rho)$ of Eq. (14) into Eq. (13) results in

(15)$$V_{\ell}\! \left( \rho, z \right) = - \mathrm{i} {\boldsymbol{\mathcal{F}}}^{-1}_{\kappa}\! \left\{ \frac{T_{\ell}\big(\rho({\boldsymbol{\kappa}})\big)}{k_z({\boldsymbol{\kappa}})} \exp \left[ \mathrm{i} \check{k}_z({\boldsymbol{\kappa}}) \Delta z + k_z({\boldsymbol{\kappa}}) R\right] \right\} \, \textrm{,}$$

which takes us to the standard Debye integral formula in a compact notation with the Fourier integral operator. In the case of a diffraction-limited field in the exit pupil the wavefront phase is purely spherical, and the extended and the standard Debye integral would therefore provide the same result. The difference between the two becomes obvious when aberrations are present. The standard Debye integral deals with them in $\arg \!\big (T_\ell (\rho)\big )$, whereas in the generalized formula of Eq. (9) the aberrations are included in $\psi ^\textrm {in}$, a term which has direct influence on the mapping and which consequently provides more accurate results without incurring higher numerical effort. In Section 3. we present some numerical examples in which we compare the generalized formula in Eq. (9) with the standard one from Eqs. (15) and (13).

2.3 Extension of the Debye integral: propagation to inclined planes

Because of the general mapping concept, the generalized Debye integral of Eq. (9) can be further extended to encompass the evaluation of electromagnetic fields on tilted planes in the focal region. A rigorous physical-optics approach to solve the problem of field propagation between non-parallel planes has been presented by Zhang et al. [23]. It is based on an extended version of Eq. (3), which allows us to combine it with the replacement of the first Fourier transform by a homeomorphic one. The only difference in the derivation is the replacement of the propagation step in the $k$ domain, which is now denoted by $\tilde {{\boldsymbol {\mathcal {P}}}}$ as shown in Fig. 2. This step can be completely described in a pointwise manner. Without showing all derivations in detail, the connection between input and output field values can be summarized as

(16)$$\tilde{V}_\perp^{\textrm{out}} \! \left( {\boldsymbol{\kappa}}^{\textrm{out}}, z^{\textrm{out}} \right) = \tilde{{\boldsymbol{\mathcal{P}}}}\! \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] \tilde{V}_\perp\! \left( {\boldsymbol{\kappa}}, z_0 \right)\,\textrm{,}$$

with

(17)$$\begin{array}{lll} \tilde{{\boldsymbol{\mathcal{P}}}} \left[ {\boldsymbol{\kappa}}^{\textrm{out}}\! \left( {\boldsymbol{\kappa}}\right)\right] & = & \tilde{a} \!\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{a} & 0 \\ 0 & \tilde{a} \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] \end{array}\,\textrm{,}$$

where ${\boldsymbol {\kappa }}^{\textrm {out}}\! \left ( {\boldsymbol {\kappa }}\right )$ denotes the mapping relation between the input and output spatial frequency vectors because of the tilted plane, $\tilde {a} \! \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 field- components projection matrix, and $\tilde {V}_\perp = \left (E_x, E_y \right )$ indicates the transversal components of the electric field. Replacing $\exp \big (\textrm {i} \check {k}_z({\boldsymbol{\kappa}})\Delta z\big )$ in Eq. (9) by $\tilde {{\boldsymbol {\mathcal {P}}}} \left [ {\boldsymbol {\kappa }}^{\textrm {out}}\! \left ( {\boldsymbol {\kappa }}\right )\right ]$ we obtain

(18)$$\begin{array}{lll} V^{\textrm{out}}_\perp ({\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] a\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] U_{\!\perp}\!\left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] \exp\!\left\{ \mathrm{i} \psi^{\textrm{in}}\! \left[ {\boldsymbol{\rho}} \!\left( {\boldsymbol{\kappa}} \right) \right] - \mathrm{i} {\boldsymbol{\kappa}} \cdot {\boldsymbol{\rho}}\!\left( {\boldsymbol{\kappa}} \right) \right\} \right) \textrm{.} \end{array}$$

As a result, Eq. (18) provides a way to calculate the field distribution on a tilted plane. Analogously to the propagation between parallel planes, the extended approach also reduces two Fourier-transform operations to one inverse Fourier transform. The extension of the integral does not require significantly more computation time, since it entails just another pointwise operation which is based, together with the homeomorphic Fourier transform, on the gridless data concept in combination with mapping operations. Because of the tilt and the change of the coordinate system for the electric field vector Eq. (18) is in vectorial form. The other four field components can be calculated from $(E_x,E_y)$ on demand [24].

Fig. 2. Propagation of a convergent wave to a tilted plane in the focal region, and the corresponding field-tracing diagram illustrating the same propagation process according to a spectrum-of-plane-wave (SPW) analysis. There we have two Cartesian coordinate systems, ${\Upsilon }^{\textrm {in}}$ and ${\Upsilon }^{\textrm {out}}$, in which the input field $V_\ell \!\left ({\boldsymbol {\rho }}, z_0 \right )$ and the output field $V^{\textrm {out}}_\ell ({\boldsymbol {\rho }}^{\textrm {out}}, z^{\textrm {out}} )$ are respectively defined. The propagation distance $\Delta z = z^{\prime } - z_0$ is given in the coordinate system ${\Upsilon }^{\textrm {in}}$. $\tilde {{\boldsymbol {\mathcal {P}}}}$ indicates the propagating operator in the spatial frequency domain.

Download Full Size | PDF

3. Numerical examples

The numerical implementation for the standard and generalized Debye integrals and the extension to tilted panes has been done in the physical-optics modeling and design software VirtualLab Fusion [25]. In what follows, we present several examples to demonstrate the accuracy and some essential applications of the generalized Debye integrals.

To have a clear-cut discussion, we would like to emphasize that, for each simulation example, we compare the results of the generalized Debye integral approach, the standard Debye integral, and the rigorous SPW analysis. The deviation between the reference and the testing approach is provided by the following expression:

(19)$$\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 the reference $V^{\textrm {ref}} \!\left ({\boldsymbol {\rho }}\right )$ is typically provided by the rigorous SPW technique and $V^{\textrm {test}} \!\left ({\boldsymbol {\rho }}\right )$ denotes the result obtained by the generalized or standard Debye integral, respectively.

3.1 Fields with aberrant phase in focal plane

The initial field used in the first example is linearly $E_x$-polarized, with a wavelength of $532\; \textrm{nm}$. The amplitude is distributed uniformly and truncated by an aperture with a circular shape and a diameter of $6\; \textrm{mm}$. A convergent spherical phase with a radius of curvature $R = -100\; \textrm{mm}$, is superimposed on this field. We can calculate the corresponding numerical aperture (NA) according to these parameters; $\textrm {NA} = n \sin \theta \approx 0.03$.

Based on the initial configuration, we continue by superimposing different types of Zernike aberrant phases on it. The variable being tested would be the type of Zernike aberration and its weighting factor. We selected three types of aberrant phase, which are described in the form of Zernike polynomials,

(20)$$\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 $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 Table 1, we present the mathematical expression of different Zernike aberrant phases and the range of their scaling factor.

Table 1. Simulation parameters of the aberrant phase for the example presented in Section 3.1.

View Table | View all tables in this article

For the following results we have tested the accuracy of the generalized Debye integral by comparing it with the SPW solution. In all cases we found values of $\sigma < 0.02{\%}$, which demonstrates the high accuracy of the generalized formula. In the subsequent part of Sec. 3.1 we compare the standard and the generalized Debye integral via $\sigma$ as defined in Eq. (19).

In Fig. 3, simulation results for secondary trefoil X with the scaling factor $c^{-2}_{4} = 5 \lambda$ are shown. Comparing, via a naked-eye observation, the obtained amplitude distributions of the $E_x$ component in sub-figures (a) and (d), it is not possible to detect the difference, with a mathematical calculation revealing a standard deviation lower than 0.05 %. However, when we turn our attention to the full complex amplitude (including phase), the standard deviation ascends to 25 %. The differences in this case, we must conclude, are mainly located in the phase. This effect can be observed by the comparison of the real part of the $E_x$ component. As shown in Figs. 3(b) and (e), where the same scaling is used, the generalized Debye integral result exhibits a higher contrast. The high accuracy of the generalized integral can be explained by something we already covered when we presented the two approaches used here: the generalized Debye integral uses a more general extracted wavefront phase than the standard one. In Figs. 3(c) and (f), the Fourier transform of the output field on the focal plane for both approaches is given. In contrast to the generalized Debye integral, the resulting spectrum of the standard Debye integral presents a uniform amplitude distribution. This is due to the fact that, considering we are dealing with a relatively low-NA situation, the pointwise mapping relation generated by the primary spherical phase employed in the standard Debye integral, in Eq. (11), approximately becomes a linear mapping. The experimental results coincide well with our theoretical expectations.

Fig. 3. Comparison of accuracy between the generalized Debye integral and the standard Debye integral for the focusing of a convergent field with a secondary trefoil X aberrant phase (SPW results are not shown, since the deviation between the generalized Debye integral result and the rigorous SPW technique is just 0.02 %). Simulation results for both approaches are respectively shown in the upper and lower rows. Panels (a) and (d) present the amplitude distribution of the $E_x$ component at the focal plane; panels (b) and (e) present the real part of the $E_x$ component; panels (c) and (f) show the amplitude distribution of the $\tilde {E}_x$ component at the focal plane (i.e. the Fourier transform of the output field).

Download Full Size | PDF

In a further investigation we did a full series of comparisons with different aberrant phases. The curves corresponding to the standard deviation between the two approaches are presented in Fig. 4. Three types of aberrant phase are superimposed, individually and in turns, onto the given input field. Even though the variation rate of the standard deviation error is different depending on the type of aberrant phase, the resulting curves reveal the same tendency. At the weak-aberration end of the scale or, in other words, for the ideal focusing system, both approaches provide an identical result. Then, as the scaling factor of the aberrant phase increases, the deviation between the two approaches grows fast. That is, when the aberration of the focusing system is too strong to be neglected, the standard Debye integral cannot predict a correct result at the focal plane. On the other hand, the high accuracy and the validity of the generalized Debye integral was tested in advance, as mentioned before ($\sigma < 0.02\;{\%}$).

Fig. 4. Comparison of accuracy between the generalized Debye integral and the standard Debye integral for the simulation task of calculating the focal field for a convergent field with different single types of aberrant phase. The abscissa records the scaling factor – i.e., the factor $C^{m}_n$ in Tab. 1 – of the corresponding Zernike term, while the ordinate marks the deviation between the two approaches, with the generalized Debye integral as reference. It should be noted that in this experiment, the deviation between the generalized Debye integral result and the rigorous SPW technique is just 0.02 %. Three curves are plotted in this graph: the solid line for secondary astigmatism Y, the dotted line for coma x and, finally, the dashed line for secondary trefoil x.

Download Full Size | PDF

3.2 Evaluation of fields in focal region

Besides its application precisely at the focal plane, both the standard and generalized Debye integrals allow us to calculate the full electromagnetic field at different positions within the focal region. Next, we would like to demonstrate the application potential of the generalized Debye integral in this regard, by applying it, in simulation, to the investigation of the focal region for a given focused field.

From the first example, we know that the proposed approach is suitable for dealing with systems with strong aberrations. On the other hand, under such circumstances, the standard Debye integral suffers from severe inaccuracy. Thus, in the second example, we configure an input field with some strongly aberrant phases and only compare the generalized Debye integral with the SPW reference. Specifically, the fundamental parameters of the system are the same as in the first case. In addition, more complicated aberrant phases are superimposed onto the original convergent field to produce the input field. The parameters for the aberrations used in the simulation are listed in Table 2.

Table 2. Simulation parameters of the aberrant phase for the example presented in Section 3.2.

View Table | View all tables in this article

The evolution of the electromagnetic field along the optical axis is shown in Fig. 5. The amplitude distribution of the $E_x$ component on the horizontal and vertical planes is displayed in sub-figures (b) and (c) respectively, in the range of $\left ( 95\; \textrm{mm}, 110\;{mm}\right )$. We can observe that the asymmetric aberrations cause an entirely different evolution of the focal field along the vertical and horizontal axes. And, in line with our expectation, the distance of $z = 100\; \textrm{mm}$ is not the position that corresponds with the minimal value of the beam diameter. The deviation between the result of the generalized Debye integral and the SPW reference is in all cases less than 0.04 %.

Fig. 5. Scan of the focal region for the system with strong aberrations using the generalized Debye integral (SPW results are not shown, since the deviation between the generalized Debye integral result and that of the rigorous SPW technique is just 0.04 %). Panel (a) shows the amplitude distribution of the $E_x$ component on the transversal $x,y$ plane at the position $z = 100\; \textrm{mm}$; (b) and (c) show the amplitude distribution of the $E_x$ component on the horizontal and vertical planes respectively, around the focal region, from $z = 95\;\textrm{mm}$ to $z = 110\;\textrm{mm}$.

Download Full Size | PDF

3.3 Fields on inclined plane in focal region

In Sec. 2.3 we provided an extension of the generalized Debye integral which enables its use to obtain fields on tilted planes in the focal region. For its demonstration we use the same input field as in Sec. 3.2. However, now we place the screen at different positions and with different inclination angles along the optical axis. The specification of the screen parameters and the simulation results are presented in Fig. 6. In all cases we verified the accuracy of the generalized Debye integral by comparing it with the SPW method from Eq. (3) and found $\sigma \approx 0.05\;{\%}$.

Fig. 6. Field distribution computed by simulation at three different $z$ positions and for different inclination angles of the plane of observation. Panel (a) shows the amplitude of the $E_x$ component where the tilted detector is located in front of the focal plane; panel (b) shows the field distribution for a tilted detector located around the focal plane; panel (c) shows the result for a tilted detector located behind the focal plane.

Download Full Size | PDF

4. Conclusion

The Debye integral can be understood as a special case of a more general integral formula, which takes the wavefront phase aberrations explicitly into account via a mapping concept. The result follows directly from the spectrum-of-plane-wave (SPW) approach by replacing the first Fourier integral by the pointwise homeomorphic Fourier transform, which was recently discussed in Wang et al. [14]. By further exploiting the pointwise mapping-type approach, the result can be extended to include propagation onto tilted planes in the focal region. The theoretical results are demonstrated and verified by several examples. They show that the generalized Debye integral formulas are more flexible in their application, and typically almost as accurate as the SPW approach for convergent fields while significantly more efficient numerically. The new techniques have been implemented in the software VirtualLab Fusion [25]. The typical calculation time on a standard laptop for the examples presented in this paper are several hundreds of ms for the Debye integral, compared to several tens of s for the SPW technique.

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. L. Mandel and E. Wolf, Optical coherence and quantum optics (Cambridge university press, 1995).

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

3. P. Debye, “Das verhalten von lichtwellen in der nähe eines brennpunktes oder einer brennlinie,” Ann. Phys. 335(14), 755–776 (1909). [CrossRef]

4. V. Dhayalan, “Focusing of electromagnetic waves,” Ph.D. thesis, University of Bergen Norway (1996).

5. J. J. Stamnes, Waves in focal regions: propagation, diffraction and focusing of light, sound and water waves (Routledge, 2017).

6. M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus field calculations,” Opt. Express 14(23), 11277–11291 (2006). [CrossRef]

7. J. J. Stamnes, “Focusing of two-dimensional waves,” J. Opt. Soc. Am. 71(1), 15–31 (1981). [CrossRef]

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

9. G. C. Sherman and W. C. Chew, “Aperture and far-field distributions expressed by the debye integral representation of focused fields,” J. Opt. Soc. Am. 72(8), 1076–1083 (1982). [CrossRef]

10. M. R. Foreman, S. S. Sherif, P. R. Munro, and P. Török, “Inversion of the debye-wolf diffraction integral using an eigenfunction representation of the electric fields in the focal region,” Opt. Express 16(7), 4901–4917 (2008). [CrossRef]

11. C. J. Zapata-Rodríguez, “Debye representation of dispersive focused waves,” J. Opt. Soc. Am. 24(3), 675–686 (2007). [CrossRef]

12. D. Panneton, G. St-Onge, M. Piché, and S. Thibault, “Exact vectorial model for nonparaxial focusing by arbitrary axisymmetric surfaces,” J. Opt. Soc. Am. 33(5), 801–810 (2016). [CrossRef]

13. J. Kim, Y. Wang, and X. Zhang, “Calculation of vectorial diffraction in optical systems,” J. Opt. Soc. Am. 35(4), 526–535 (2018). [CrossRef]

14. 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]

15. O. Baladron-Zorita, Z. Wang, C. Hellmann, and F. Wyrowski, “Isolating the gouy phase shift in a full physical-optics solution to the propagation problem,” J. Opt. Soc. Am. 36(9), 1551–1558 (2019). [CrossRef]

16. F. Wyrowski, “Unification of the geometric and diffractive theories of electromagnetic fields,” in Proc. DGaO, (2017), p. A36.

17. R. Shi and F. Wyrowski, “Comparison of aplanatic and real lens focused spots in the framework of the local plane interface approximation,” J. Opt. Soc. Am. A 36(10), 1801–1809 (2019). [CrossRef]

18. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems, ii. structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A 253(1274), 358–379 (1959). [CrossRef]

19. E. Wolf, “Electromagnetic diffraction in optical systems - i. an integral representation of the image field,” Proc. R. Soc. Lond. A 253(1274), 349–357 (1959). [CrossRef]

20. O. Bryngdahl, “Geometrical transformations in optics,” J. Opt. Soc. Am. 64(8), 1092–1099 (1974). [CrossRef]

21. J. J. Stamnes, “Waves, rays, and the method of stationary phase,” Opt. Express 10(16), 740–751 (2002). [CrossRef]

22. Z. Wang, S. Zhang, O. Baladron-Zorita, C. Hellmann, and F. Wyrowski, “Application of the semi-analytical fourier transform to electromagnetic modeling,” Opt. Express 27(11), 15335–15350 (2019). [CrossRef]

23. 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]

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

25. 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	Scaling Factor
Coma X	$Z_{3}^{1} (r, θ) = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 0 \sim 5 λ$
Secondary Astigmatism Y	$Z_{4}^{- 2} (r, θ) = \sqrt{10} (4 r^{4} - 3 r^{2}) \sin 2 θ$	$c_{4}^{- 2} = 0 \sim 5 λ$
Secondary Trefoil X	$Z_{6}^{4} (r, θ) = \sqrt{14} (6 r^{6} - 5 r^{4}) \cos 4 θ$	$c_{6}^{4} = 0 \sim 5 λ$

Name & Type	Mathematical Expression	Scaling Factor
Tilt Y	$Z_{1}^{- 1} (r, θ) = 2 r \sin θ$	$c_{1}^{- 1} = 1.5 λ$
Astigmatism Y	$Z_{2}^{- 2} (r, θ) = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{- 2} = 2 λ$
Spherical	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = - 2 λ$
Trefoil Y	$Z_{3}^{- 3} (r, θ) = \sqrt{8} r^{3} \sin 3 θ$	$c_{3}^{- 3} = 2.5 λ$
Secondary Coma X	$Z_{5}^{1} (r, θ) = \sqrt{12} (10 r^{5} - 12 r^{3} + 3 r) \cos θ$	$c_{5}^{1} = 0.5 λ$
Secondary Spherical	$Z_{6}^{0} (r, θ) = \sqrt{7} (20 r^{6} - 30 r^{4} + 12 r^{2} - 1)$	$c_{6}^{0} = 1 λ$
Tetrafoil Y	$Z_{4}^{- 4} (r, θ) = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{- 4} = - 1.5 λ$
Tertiary Astigmatism X	$Z_{6}^{2} (r, θ) = \sqrt{14} (15 r^{6} - 20 r^{4} + 6 r^{2}) \cos 2 θ$	$c_{6}^{2} = 2 λ$

Name & Type	Mathematical Expression	Scaling Factor
Coma X	$Z_{3}^{1} (r, θ) = \sqrt{8} (3 r^{3} - 2 r) \cos θ$	$c_{3}^{1} = 0 \sim 5 λ$
Secondary Astigmatism Y	$Z_{4}^{- 2} (r, θ) = \sqrt{10} (4 r^{4} - 3 r^{2}) \sin 2 θ$	$c_{4}^{- 2} = 0 \sim 5 λ$
Secondary Trefoil X	$Z_{6}^{4} (r, θ) = \sqrt{14} (6 r^{6} - 5 r^{4}) \cos 4 θ$	$c_{6}^{4} = 0 \sim 5 λ$

Name & Type	Mathematical Expression	Scaling Factor
Tilt Y	$Z_{1}^{- 1} (r, θ) = 2 r \sin θ$	$c_{1}^{- 1} = 1.5 λ$
Astigmatism Y	$Z_{2}^{- 2} (r, θ) = \sqrt{6} r^{2} \sin 2 θ$	$c_{2}^{- 2} = 2 λ$
Spherical	$Z_{4}^{0} (r, θ) = \sqrt{5} (6 r^{4} - 6 r^{2} + 1)$	$c_{4}^{0} = - 2 λ$
Trefoil Y	$Z_{3}^{- 3} (r, θ) = \sqrt{8} r^{3} \sin 3 θ$	$c_{3}^{- 3} = 2.5 λ$
Secondary Coma X	$Z_{5}^{1} (r, θ) = \sqrt{12} (10 r^{5} - 12 r^{3} + 3 r) \cos θ$	$c_{5}^{1} = 0.5 λ$
Secondary Spherical	$Z_{6}^{0} (r, θ) = \sqrt{7} (20 r^{6} - 30 r^{4} + 12 r^{2} - 1)$	$c_{6}^{0} = 1 λ$
Tetrafoil Y	$Z_{4}^{- 4} (r, θ) = \sqrt{10} r^{4} \sin 4 θ$	$c_{4}^{- 4} = - 1.5 λ$
Tertiary Astigmatism X	$Z_{6}^{2} (r, θ) = \sqrt{14} (15 r^{6} - 20 r^{4} + 6 r^{2}) \cos 2 θ$	$c_{6}^{2} = 2 λ$

Generalized Debye integral

Abstract

1. Introduction

2. Theory

2.1 Derivation of the generalized Debye integral

2.2 Standard Debye integral as a special case

2.3 Extension of the Debye integral: propagation to inclined planes

3. Numerical examples

3.1 Fields with aberrant phase in focal plane

3.2 Evaluation of fields in focal region

3.3 Fields on inclined plane in focal region

4. Conclusion

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (6)

Tables (2)

Equations (20)

Optics Express