1. INTRODUCTION

The objective of this three-part paper series is to develop a computationally efficient method that can accurately trace the flux lines in the focal region of a lens and to explore the application of this method for understanding lens aberrations and laser spatial modes. The flux lines, which we interchangeably call “nonlinear rays,” are discussed in much more detail in the second paper. Briefly, the concept of wave optics and geometrical ray optics can be related using the eikonal function [1], which can be related to the derivative of the phase of the wavefront. By following paths that are guided by the derivative of the phase of the wavefield in three dimensions, these flux lines can be traced. Indeed, these nonlinear rays and the classical rays are one and the same if the wavefront is slowly varying. A physical interpretation of the flux lines and their relationship to classical ray optics is provided in the next paper as is a detailed discussion of the method used to trace the lines. A necessary framework for this method is the availability of the samples of the converged wavefield in the focal region of the lens over a three-dimensional grid. The computation of this 3D sample distribution in an accurate and timely manner is the subject of this first paper in the three-part series.

More specifically, in the first paper, we focus on numerical methods that can accurately calculate a propagating complex wavefield over a three-dimensional volume, and the associated sampling requirements and computational efficiencies are considered in some detail. We limit this study to focusing fields, for which the most rapid change of flux can be expected. Two distinct numerical algorithms are developed for the case of a thin lens and an ideal lens, respectively, both of which iteratively apply the angular spectrum method over a range of distances, albeit in different manners. In the second paper, these numerical algorithms form the basis of a nonlinear ray tracing algorithm, whereby the 3D wavefield is sampled in a uniform grid in the focal region of the lens, and the “ray” is guided by the derivative of the phase, which is calculated at a sequence of points travelling through this grid. Interestingly, the nonlinear rays appear to propagate around regions of low intensity. In the third paper, the flux tracing algorithm is applied to convergent lenses that contain aberration, several of which cannot readily be interpreted in terms of linear geometrical ray optics but can be interpreted more easily in terms of the nonlinear rays. In summary, the concept of flux and flux-tracing is explored in the second paper, and its application to aberration theory is highlighted in the third paper; in this first paper, we focus only on calculating a focusing wavefield over a uniform 3D grid of sampling points, which is a necessary precursor to flux-tracing.

At the heart of this paper is the concept of numerical propagation. Numerical propagation algorithms that can accurately calculate diffraction patterns at variable distances within an homogeneous isotropic medium are an active area of research, with far reaching applications including in the areas of quantitative phase microscopy [2], computer generated holography [3], confocal microscopy [4], computation of the point spread function in lens design and lithography [5,6], and digital holography [7]. Solutions based on paraxial scalar optics include various discretizations of the Fresnel and linear canonical transforms [8], which continue to find application in quantitative phase imaging. However, algorithms that are founded on nonparaxial scalar optics are, in general, preferable since they offer superior performance in terms of accuracy and can be numerically implemented at least as efficiently as paraxial based algorithms.

These algorithms include the angular spectrum method (ASM) [9,10] and the Rayleigh–Sommerfeld convolution (RSC) method [11–13], which can both provide highly accurate solutions for numerical apertures up to at least 0.6 [5,14]; for higher numerical apertures, algorithms based on the more intricate vectorial model are required [6]. Indeed the first Rayleigh–Sommerfeld diffraction integral and the propagation of the angular spectrum can be shown to provide identical solutions [10]. However, their discrete counterparts, the RSC and ASM, have significantly different sampling considerations [11,15–17]. In this paper, we employ the ASM as the numerical propagation algorithm of choice, making use of recent advances in sampling theory to optimize the computational memory efficiency of the ASM for focused fields [18,19]. The choice of the ASM to calculate the diffraction pattern in the focal region of a lens is based on several factors. First, methods that are based on the non-paraxial (Fresnel) model will not be sufficiently accurate for the numerical apertures considered in these papers. The ASM and RSC methods are computationally fast algorithms that can calculate the field with high accuracy up to a numerical aperture of 0.6–0.7, but the ASM performs better for the spatial and sampling requirements of focused fields resulting from wide lens apertures, as explained in Refs. [11,18]. For an excellent review of the various propagation algorithms and their frequency and sampling considerations, we refer the reader to Ref. [20].

In this tutorial we cover several detailed topics:

The breakdown of this paper is as follows: In Section 2, the ASM is described, and the associated sampling conditions are examined. Particular attention is given to the superposition-tile method [18], which is renamed the “pre-aliasing” method in this paper; this is accompanied by a rigorous proof of the underlying sampling theory in Supplement 1. In Section 3, the two algorithms are defined to calculate the 3D diffraction pattern around the focal point of a thin lens and an ideal lens, respectively. In Sections 4 and 5, the algorithms are extended to account for the spatial mode of the laser. This is followed by a brief discussion in Section 6, which will lead into the second paper.

2. ANGULAR SPECTRUM METHOD

A. Propagation of the Angular Spectrum

Propagation of the angular spectrum [10] is described as follows:

B. Angular Spectrum Method

The angular spectrum method (ASM) can be used to calculate the samples of the diffracted wavefield, ${u_z}$, using the samples of ${u_0}$. We stipulate that the two continuous wavefields have spatial support $(\Delta {x_0},\Delta {y_0})$ and $(\Delta {x_z},\Delta {y_z})$ in space and both have spatial frequency support of $(\Delta {f_x},\Delta {f_y})$, which is a consequence of the properties of the propagation of the angular spectrum as described in Section 2.A. We note that the spatial supports at distance $0$ and at $z$ can be significantly different for the case of a converging or diverging wavefield. It is well known that a function cannot be truly bounded with finite support in both space and spatial frequency; however, as discussed in Supplement 1, this can be said to be approximately true in many cases. Both continuous functions ${u_0}$ and ${u_z}$ can be recovered (via Shannon interpolation) from their discrete samples, if uniform sampling intervals $({T_x},{T_y})$ have been applied where the Nyquist condition is satisfied, i.e., ${T_x} \le 1/\Delta {f_x}$ and ${T_y} \le 1/\Delta {f_y}$.

The ASM is based on a direct “discretization” of Eqs. (1) and (2) and is implemented using the DFT in place of the Fourier transform operator in Eq. (1) as follows:

We refer the reader to Supplement 1 for a definition of the DFT as well as a detailed derivation of the sampling constraints of the ASM.

C. Pre-aliasing for Convergent Wavefields

Recently, there have been several contributions on the numerical calculation of diffraction fields converging at the focal point of a lens [11,17,18]. Of particular interest in this paper is the method developed by Kanka et al. [19], and further investigated by [18], which exploits redundancy for the case of a small spatial support in the output field for convergent focusing fields; in such cases the input to the ASM can be broken up into smaller tiles, which are superimposed, and this smaller matrix serves as input to the ASM without any loss in accuracy. In this paper we refer to this method as “pre-aliasing,” and a rigorous proof of its validity is provided in Supplement 1 using fundamental sampling theory. The method is summarized below.

Taking the definition of the ASM above, the discrete signals ${u_0}({{n_x},{n_y}})$ and ${u_z}({{n_x},{n_y}})$, which are input and output to the ASM, must have the same support, $({\Delta {x_z},\Delta {y_z}}) = ({\Delta {x_0},\Delta {y_0}}) = ({{N_x}{T_x},{N_y}{T_y}})$. However, for convergent wavefields, the propagated wavefield will have significantly smaller support: $({\Delta {x_z},\Delta {y_z}}) \ll ({\Delta {x_0},\Delta {y_0}})$. In such case, a smaller number of samples ${M_x} \times {M_y}$ is required to fully represent the continuous function ${u_z}$, where $({\Delta {x_z},\Delta {y_z}}) = ({{M_x}{T_x},{M_y}{T_y}})$. It is possible to generate a discrete input function in Step 1 that will produce this reduced set as samples. In short, Step 1 in the algorithm described in Section 2.B is replaced with the following:

It should be noted that the subject of pre-aliasing in the manner described above can be most easily interpreted in terms of phase space description as described in Chapter 10 of Ref. [21], whereby the signal support in both space and frequency can be visualized in phase space, and the effect of propagation is a shearing of this support.

D. Sampling Conditions in Three Dimensions for the ASM

Taking the definition for the propagation of the angular spectrum in Section 2.A, the required sampling conditions can be defined in terms of Nyquist and Shannon sampling theory [22]. First, the input function ${u_0}$ must be Nyquist sampled such that the sampling intervals ${T_x} \le | {1/\Delta {f_x}} |$ and ${T_y} \le | {1/\Delta {f_y}} |$. However, if the support $({\Delta {f_x},\Delta {f_y}})$ is greater than the support of ${H_z}$ defined in Eq. (2), which has radius $1/\lambda$, then this latter support will constrain the sampling intervals. More explicitly, the values for ${T_x}$ and ${T_y}$ should be chosen to satisfy the constraint:

In general, $1/\Delta {f_{x,y}}$ will be the larger of the two values. In such case the spatial frequency support of ${u_z}$ will be the same as that of ${u_0}$; therefore, the same sampling intervals in $x$ and $y$ are applicable in both planes. This case is illustrated in Fig. 1(a). We note that there have been several contributions in the literature that relax the Nyquist constraint for the propagated wavefield, or indeed for any chirped wavefield, which can be related to the shape of the signal’s phase space distribution [21,23,24]. These generalized sampling theorems can tolerate the aliasing of the function in the Fourier domain caused by breaking of the Nyquist condition. However, the ASM cannot generally be applied to such a sub-Nyquist sampled function due to the multiplication in Eq. (4); see Ref. [11] for an excellent analysis of this subject.

 

Fig. 1. Illustration of the bandwidth of the three-dimensional function ${u_z}({x_z},{y_z})$. (a) The bandwidth in $x$ and $y$ is in general determined by the bandwidth of the initial function ${u_0}$, and (b) the bandwidth in the $z$-direction is determined by the size of the surface area on the Ewald sphere over which the function exists.

Download Full Size | PDF

The previous paragraph deals with the sampling conditions for the input/output functions for the ASM for the $x$- and $y$-dimensions. However, if a three-dimensional wavefield was to be sampled, for example, by repeated application of the ASM to ${u_0}$ over a sequence of different distances, ${n_z}{T_z}$, where ${n_z}$ is an integer, then the sampling interval, ${T_z}$, also needs to be considered. Similar to the case of sampling in $x$ and $y$, the relevance of the sampling condition in $x$ relates to the correct reconstruction (using Shannon interpolation) of the continuous field ${u_z}(x,y)$ from a uniform sampling of the field over $z$ using a sampling interval ${T_z}$. In order to define the Nyquist condition for ${T_z}$ we must first investigate the frequency distribution of ${u_z}(x,y)$ in the $z$ dimension.

Taking the Fourier transform over $z$ of the continuous function ${u_z}({x_z},{y_z})$ as defined by Eq. (1), it is straightforward to deduce that the three-dimensional Fourier transform of ${u_z}$ exists only over a surface defined by the Dirac-delta functional $\delta ({{f_z} - \sqrt {1/{\lambda ^2} - f_x^2 - f_y^2}})$, which relates to the well-known Ewald sphere [1], as illustrated in Fig. 1(b). Therefore, the maximum frequency in the $z$-direction is given by ${f_{z,\max}} = 1/\lambda$, which occurs at ${f_x} = {f_y} = 0$; and the minimum frequency in the $z$-direction occurs at ${f_{x,y}} = \Delta {f_{x,y}}/2$ for the general case and is illustrated in Fig. 1. From this illustration it can be seen that the minimum spatial frequency in ${f_z}$ is given by

Defining the minimum sampling interval in ${f_z}$ requires further consideration. As illustrated in Fig. 1(b) it can be seen that in the general case, where ${(\Delta {f_x}/2)^2} + {(\Delta {f_y}/2)^2} \lt 1/{\lambda ^2}$, the three-dimensional Fourier transform of ${u_z}({x_z},{y_z})$ does not exist for the center band of ${f_z}$, and two separate bands exist on each side with a width defined as $\Delta {f_z}/2$. For signals such as this, it is sufficient to apply a sampling interval given by $\Delta {f_z}$ even though the support in ${f_z}$ is distributed over two separate bands. Therefore, the sampling interval in the $z$-direction is given by

 

Fig. 2. Computational and memory savings provided by pre-aliasing the input to the ASM to calculate the 3D distribution of a convergent wavefield where the input has a large number of samples, ${N_x} = {N_y} = 10{,}000$. (a) Computational and memory saving for the 2D case only for a range of values ${P_x}$ and (b) the same result over a continuous range of ${M_x}$. (c) and (d) show the same set of results for the 3D case where ${N_z} = 1000$.

Download Full Size | PDF

Therefore, if the ASM is used to calculate the samples of a wavefield in three dimensions, sampling intervals of ${T_x}$ and ${T_y}$ should meet the constraint set out in Eq. (7), and the sampling interval ${T_z}$ should meet the constraint set out in Eq. (9) such that the continuous three-dimensional function ${u_z}({x_z},{y_z})$ can be recovered via Shannon interpolation, i.e., by convolution with a three-dimensional Sinc convolution as follows:

E. Computational Efficiency of the ASM

The pre-aliasing step for convergent fields described in Section 2.C requires the traditional Nyquist condition to be satisfied. The redundancy that is exploited relates to the small spatial support of ${u_z}$, which places a relaxed sampling condition on its DFT, which in turn translates to pre-aliasing of ${u_0}$. The pre-aliasing step is not required for convergent fields; it is always possible to avoid this step, but the advantage of pre-aliasing relates to computational and memory efficiency benefits. For example, numerical propagation between a lens with large radius and high numerical aperture to its focal plane would necessitate an input field with a small sampling interval applied over a large area; the total number of samples for such a discrete function may be onerous in terms of temporary computer memory. However, following the reduction in the size of the input matrix, ${u^\prime _0}$ defined in Eq. (5), the requirement for memory is significantly relaxed. Although the efficiency of memory afforded by this technique is straightforward to calculate and has been highlighted by other authors [18], there has been little discussion on the numerical efficiency of the technique to date.

It is naive to assume that the pre-aliasing step will translate to a significant computational saving to match the saving in memory. However, the saving in computational load, while still significant, is much less than the saving in terms of memory. The DFT is typically calculated using the fast Fourier transform (FFT) algorithm [22], which will have a number of operations in the order of $O({{M_x}{M_y}\mathop {\log}\nolimits_2 ({{M_x}{M_y}})})$ for radix 2 FFT. Thus, if $({{M_x},{M_y}}) \ll ({{N_x},{N_y}})$, it can be expected that there will be a significant saving in terms of computational time for Steps 2–4 defined in Section 2.C. However, this is partially offset by the additional computational load introduced with the implementation of the pre-aliasing step in the revised Step 1; the greater the saving in Steps 2–4, the greater the additional load in Step 1 due to the larger values of ${P_x}$ and ${P_y}$ in Eq. (6). However, if the purpose is to apply the ASM repeatedly to the same input ${u_0}$ in order to calculate the field over a long sequence of different distances, then far greater computational savings can be made. For this case, Step 1 (the pre-aliasing step) is implemented only once, while Steps 2–4 are repeated iteratively.

 

Fig. 3. Illustration of the two algorithms developed in this paper. (a) For the case of thin lens, the input to the algorithm is ${u_0}({x_0},{y_0})$, which is the product of the complex thin lens transmission function $t$ shown in gray, the laser profile $A$ shown in blue, and the lens pupil function of radius $R$ shown in red. The ASM is the key component in the algorithm and is iteratively applied to generate a 3D matrix of samples that represent the focused wavefield over a volume of interest. (b) The algorithm comprises five steps, which are detailed in the text. Note that Step 3 is a “pre-aliasing” step whereby multiple copies of the input are shifted and added. This step is applied only once and serves to reduce the number of samples fed as input to the iteratively applied ASM. (c) For the case of ideal lens, the input to the algorithm is ${U_0}({x^\prime _0},{y^\prime _0})$, which is the product of the laser profile $A$ and the lens pupil function of radius $R$. A first DFT operation calculates the focused plane, and following this the ASM is applied iteratively to generate a 3D matrix of samples that represent the focused wavefield over a volume of interest around the focused plane. (d) The ideal lens algorithm comprises six steps, which are detailed in the text.

Download Full Size | PDF

The computational efficiency of the ASM with and without the pre-aliasing step can be quantitatively compared. With pre-aliasing, Step 1 will require a number of operations in the order of $O({P_x}{P_y}{M_x}{M_y})$. Therefore, the total number of operations without pre-aliasing is given by $O({N_x}{N_y}\mathop {\log}\nolimits_2 ({{N_x}{N_y}}))$, and with pre-aliasing this becomes $O({P_x}{P_y}{M_x}{M_y} + {M_x}{M_y}\mathop {\log}\nolimits_2 ({{M_x}{M_y}}))$. When the ASM is applied to the same input repeatedly in order to calculate the 3D diffraction distribution, the comparison of these efficiencies changes dramatically; the total number of operations without pre-aliasing is given by $O({N_z}{N_x}{N_y}\mathop {\log}\nolimits_2 ({{N_x}{N_y}}))$, and with pre-aliasing this becomes $O({P_x}{P_y}{M_x}{M_y} + {N_z}{M_x}{M_y}\mathop {\log}\nolimits_2 ({{M_x}{M_y}}))$. To illustrate the practical importance of the pre-aliasing step we examine the case where then input ${u_0}$ is sampled with ${N_x} = {N_y} = 10{,}000$ samples. Taking ${M_x} = {M_y} = 50,100,200,500,1000,2000,5000$ (i.e., ${P_x} = {P_y} = 200,100,50,20,10,5,2$), the memory saving and the computational saving are shown in Fig. 2(a); for the first five cases the computational saving is approximately the same at ${\approx} \times 22 - 26$. The memory saving, which is plotted on a $\mathop {\log}\nolimits_{10}$ scale, is significantly higher and varies over a much wider range from 40,000 for ${P_x} = 200$ down to $4$ for ${P_x} = 2$. These savings are plotted over a continuous range of ${M_x}$ in Fig. 2(b). For the three-dimensional case, significantly greater computational savings are made, while memory saving remains the same. Taking the same parameters as above and now taking ${N_z} = 1000$, the savings are illustrated in Fig. 2(c). In this case the computational saving ranges from $\times 20{,}000$ for ${P_x} = 200$ to $\times 4$ for ${P_x} = 2$. These savings are plotted over a continuous range of ${M_x}$ in Fig. 2(d).

3. NUMERICAL CALCULATION OF FOCUSED WAVEFIELDS

In this section, we outline two methods for calculating the complex amplitude of a focused wavefield in a three-dimensional volume, which are illustrated in Fig. 3. Both methods use the ASM over a sequence of distances in order to generate a point cloud of complex values over a uniform three-dimensional grid. The first method is based on the thin lens approximation, whereby the lens can be described as a simple two-dimensional (2-D) phase only function. The second method is based on modeling the effect of the lens as a simple pupil function in the Fourier domain for the case of coherent imaging.

A. Thin Lens Algorithm

The thin lens approximation [10] is commonly used within the paraxial approximation, to model a lens as a simple 2-D phase transformation, where the phase relates to the spatially variant delay imparted by the lens due to a varying thickness and/or refractive index. This approximation is based on the position of a ray that is incident on the lens in $(x,y)$ being identical to the position of the ray as it leaves the lens; this requires negligible propagation of the ray within the lens. The phase transformation of the lens is given by

Assuming an incident plane wave of complex amplitude $A(x,y)$ is input to the plane of the lens, the resultant wavefield immediately after the lens is given by

The wavefield ${u_0}$ will converge at the focal point of the lens, i.e., at a distance $\approx f$, where $f$ is the focal length of the lens as illustrated in Fig. 3(a). The corresponding algorithm to calculate the three-dimensional complex field in the region of this convergence is illustrated in Fig. 3(b) as a sequence of five steps, which are further defined below.

Some initial results from this algorithm are provided in Fig. 4. Here, the three-dimensional field is calculated for three different thin lens functions; for all cases, “top-hat” illumination is assumed ($A = 1$) such that ${u_0} = t$, and for all cases the aperture of the lens was given by $R = 2.5$ mm, and the wavelength of light was defined to be $\lambda = 1064$ nm. Following the approach taken by Goodman [10], the first lens is defined as follows:

 

Fig. 4. The intensity profiles of convergent light from four different lenses with the same numerical aperture and with the same top-hat uniform illumination, $A = 1$. For the case of the thin lens defined in Eq. (14), (a1) and (a2) are slices through the computer 3D volume for $y = 0$ (i.e., the $xz$-plane) and $z = f - 49\;{\rm{\unicode{x00B5}{\rm m}}}$, which was found to be the $xy$-plane with the brightest spot. Also shown in (a2) is the inset image of the profile of this spot. In (b1) and (b2) the same figures are shown for the thin lens with purely quadratic phase given in Eq. (17). The profiles for these two cases are almost identical owing to the high accuracy of the approximation of obtaining ${t_2}$ from ${t_1}$. For the case of the third thin lens defined in Eq. (18), (c1) and (c2) are slices through the volume for $y = 0$ and $z = f$; for this case the light converges to a much smaller spot due to the profile of the lens function. Finally, the same slices are shown in (d1) and (d2) for the case of the ideal lens. Remarkably, the intensity profiles obtained using the RS lens and the ideal lens are very similar even though these were computed using different algorithms and definitions. Note that the scale used for images (c1), (c2), (d1), and (d2) is 0.05 times that used for the other four images due to the significantly higher intensity values; (e1) shows a comparison of the amplitude profile over the line indicated in (c3). The result from the thin lens algorithm (TLA) is plotted with the corresponding result using direct numerical integration, which is described in the text. The difference in amplitude is shown on a log scale in (e2); equivalent results are shown for the ideal lens case in (e3) and (e4).

Download Full Size | PDF

We note that ${t_2}$ can be obtained from ${t_1}$ based on the paraxial approximation. The same values for $f$ and $r$ are used as for ${t_1}$ in the simulation. The results for ${t_1}$ and ${t_2}$ are almost identical, highlighting the validity of the paraxial approximation in modelling this thin lens.

It is interesting to note that the thin lens algorithm predicts that these lenses produce their bright spots at a distance that is slightly defocused. The quadratic phase lens function that defines these lenses is essentially a paraxial model, by which we mean that such a phase function will produce accurate focusing when the NA is small (resulting from a small lens aperture and/or a long focal length). However, a lens with an NA of 0.6 cannot be considered paraxial. If we had employed the Fresnel transform (or rather a numerical equivalent), which is of course itself based on the paraxial approximation, then this lens would indeed produce the brightest spot at the focal plane as expected. However, here we use the AS method, which can be related to the RS diffraction integral and is more accurate in modelling light propagation for numerical apertures of 0.6. Therefore, we find that the quadratic phase lens is imperfect for an NA of 0.6 and produces the brightest spot at a distance that is slightly defocused.

A third thin-lens function was investigated, where

This definition is borrowed from the Rayleigh–Sommerfeld (RS) diffraction integral, which relates a field ${u_0}({x_0},{y_0})$ to the diffracted field at distance $z$ as follows:

This integral can be rewritten as a convolution as follows:

The complex term in square brackets has an approximately constant phase angle of $\pi /2$. Thus, a point source of light, represented by a Dirac-delta functional (see Supplement 1), will produce a diffraction pattern at distance $f$ with a phase function given by the conjugate of ${t_{\text{RS}}}$. Such a lens function, which we refer to as a Rayleigh–Sommerfeld (RS) lens going forward, can therefore be expected to produce a highly focused spot after propagation a distance $z = f$. Selecting $f = 3.3\;{\rm{mm}}$ and choosing all of the same sampling parameters as for the previous two cases, the resulting intensity is shown in Figs. 4(c1) and 4(c2). The spot in this case is significantly more tightly focused and as a consequence is more intense than for the previous cases, and it is scaled $\times 0.05$ for comparison. The wide range of intensity allows only a small region around the focus to be visualized using a linear scale; hence, an inset image shows the logarithm of the intensity to reveal the light pattern in the region of the focus in more detail. The inset image in Fig. 4(c2) shows the profile of the focused spot, which is significantly narrower and $\times 20$ more intense that the corresponding profiles shown in Figs. 4(a2) and 4(b2).

In order to evaluate the accuracy of the thin lens algorithm a comparison is given with direct numerical integration in Figs. 4(e1) and 4(e2). Direct integration is implemented using a Gauss–Kronrod quadrature algorithm [25], which has previously been applied to evaluate the angular spectrum method [11]. This algorithm is used to numerically calculate the Rayleigh–Sommerfeld diffraction integral using the RS lens as the input function, and it produces a value of the complex amplitude of the propagated wavefield at a single point in space. While this algorithm is comparatively slow for the large number of points that describe our RS lens, we can be reasonably confident that it has high accuracy, and we use it as a standard for evaluating the accuracy of the thin lens algorithm, and later using the ideal lens algorithm. The profile of the amplitude over a line close to the focal plane [indicated in Fig. 4(c3)] is calculated using both methods; the results are shown to be in close agreement in Fig. 4(e1), and the difference in amplitude is given in Fig. 4(e2). Although this difference is relatively small, it is important to point out that small errors in calculating the complex wavefield could potentially lead to erroneous results when tracing the eikonal in the two papers that follow. Nevertheless, this comparison provides a validation that the thin lens algorithm is operating as expected within a small margin of error.

B. Ideal Lens Algorithm

Assuming an incident plane wave of complex amplitude $A(x,y)$, we define the effect of an ideal lens with focal length $f$ and numerical aperture NA in terms of the optical Fourier transform [10]:

Setting $A = 1$ for the case of an ideal infinite plane wave, the result is in the form of a Bessel function of the first kind [10]: $\pi {r^2}{J_1}(2\pi r\rho)/(\pi r\rho)$, where $\rho$ is the radial coordinate. The intensity of this pattern is known as the Airy pattern and the width of the central lobe given by $1.22\lambda /{\rm NA}$, which is the well-known limit of optical resolution for coherent imaging. Therefore, two ideal microscope objectives with the same $NA$ but different focal lengths, $f$, will produce identical Airy patterns at their respective focal length for the same uniform plane wave illumination but would have different input physical pupil radius given by $R = fNA$. Since the pupil plane of the microscope objective is typically close to the back aperture of the objective, we can take this radius to be equivalent to the physical radius of the lens back aperture. The wavefield ${u_0}$ will converge at the focal point of the lens, and we proceed to define an algorithm to calculate the three-dimensional volume of light in this region. This algorithm is illustrated in Fig. 3 as a sequence of six steps, which are further defined below.

An initial result from this algorithm is shown in Figs. 4(d1)–4(d3). Here, for Step 1 the following selections were made: As for the previous cases, “top-hat” illumination was used, i.e., $A(x,y) = 1$, and no aberration is included, such that ${U_0}({x^\prime _0},{y^\prime _0}) = P({{{x^\prime_0}}/fNA,{{y^\prime_0}}/fNA})$, where $f$ and NA were chosen to be $3.3\;{\rm{mm}}$ and $0.6\;{\rm{mm}}$, respectively, resulting in an aperture size of $R = 2.5\;{\rm{mm}}$. As for the previous cases, the wavelength of light was defined to be $\lambda = 1064\;{\rm{nm}}$. The parameters selected in Step 2 of the algorithm are as follows: ${M_x} = {M_y} = {M_z} = 2000$, ${T_x} = {T_y} = {T_z} = 0.1\;{\rm{\unicode{x00B5}{\rm m}}}$. As for the previous examples, to save memory only one $xz$-slice and one $xy$ image were stored; cropped areas of both are shown in Figs. 4(d1) and 4(d2) with $\times 10$ interpolation. The profile of the focused spot is shown in the inset image in Fig. 4(d1). The spot in this case is similar to that for the case of the RS lens. The magnified image in Fig. 4(d3) shows the logarithm of the intensity of the $xz -$ image, in which it can be seen that there are striking similarities to the same image for the thin (RS) lens.

As for the previous case, we evaluate the accuracy of the ideal lens algorithm by comparison with direct numerical integration using a Gauss–Kronrod quadrature algorithm. In this case this algorithm is used to numerically calculate the Rayleigh–Sommerfeld diffraction integral using a Bessel function as input, which is taken to be the analytical result of the first part of the ideal lens algorithm (i.e., DFT of a circ function). The profile of the amplitude over a line close to the focal plane [indicated in Fig. 4(d3)] is calculated using both methods, and the results are shown to be in close agreement in Fig. 4(e3); the difference in amplitude is given in Fig. 4(e3). As for the previous case we can state that this comparison provides a validation that the ideal lens algorithm is operating as expected within a small margin of error, with the caveat that even small errors could lead to slightly inaccurate tracing of the eikonal in the papers that follow.

4. INCLUDING THE LASER SPATIAL MODE

In the discussion thus far, top hat illumination has been assumed for which $A = 1$, and the laser power, $P$, is given no consideration. In practice, the laser will have a particular spatial mode that defines its complex spatial profile and a total power, both of which should be taken into account. Many applications employ lasers with a ${{\rm TEM}_{00}}$ or ${{\rm TEM}_{01}}$ mode that are characterized by a Gaussian or, more generally, by a Laguerre–Gaussian (LG) distribution. While “top-hat” illumination would produce the tightest laser spot, most applications will require maximal use of the laser power, and as such the laser cannot be widely expanded to achieve a uniform “top-hat” profile across the lens aperture; it is common for the laser to be expanded such that the lens aperture apodizes the Gaussian/LG distribution at a particular radius in the distribution, which is deemed to provide the most acceptable compromise between the loss in laser power and the enlargement of the spot size. Examples that apply this approach include optical trapping, confocal microscopy, and laser cutting. We examine two commonly used laser spatial modes here: ${{\rm TEM}_{00}}$ or ${{\rm TEM}_{01}}$. The former is the most common laser mode, and in its simplest form for which the laser is assumed to have negligible expansion, the electric field is defined by a Gaussian distribution as follows:

 

Fig. 5. Phase and intensity of the focused spot for the case of the Rayleigh–Sommerfeld lens defined in Eq. (18) for the case of both Gaussian and ${{\rm TEM}_{01}}$ laser spatial modes. All images were calculated using the thin lens algorithm; (a1) and (b1) are the intensity and phase of the center XZ slice for $y = 0$ for the case of a Gaussian beam; Visualization 1 shows a video of the $xz$ intensity/phase slice as a function of different values of $y$. (c1) shows the intensity of the $xy$-plane at five different distances indicated in the lower left of the figure; (d1) shows the profile of the focused spot in the five planes. (e1) and (e2) show the corresponding phase image and profiles. Visualization 2 shows a video of the $xy$ intensity/phase slice as a function of different values of $z$. A similar set of results is shown in (a2)–(f2) for the case of a ${{\rm TEM}_{01}}$ laser mode. The characteristic doughnut shape of the focused spot is seen to expand with defocus, and the spiral phase pattern is seen to twist. Remarkably, the $xz$-slice of the phase distribution is identical to that for the Gaussian case, except for a constant phase shift on the right side of the image. Visualization 3 and Visualization 4 show videos of the variation in the $xz$- and $xy$-planes as for the Laguerre Gaussian case. The same results for the ideal lens calculated using the ideal lens algorithm are provided in Visualization 5, Visualization 6, Visualization 7, and Visualization 8.

Download Full Size | PDF

The Gaussian distribution defined in Eq. (24) is ideal and does not take into account expansion of the laser beam or the change in the phase distribution as a function of propagation. The ${{\rm TEM}_{00}}$ laser mode is more accurately described as a special case of the more general Laguerre–Gaussian mode. The amplitude distribution of the Laguerre–Gaussian mode is uniquely defined by the beam waist width $w(0)$ and the Rayleigh range, ${z_R}$, as well as integers $l$ (the azimuthal index) and $p$ (the number of radial nodes), as follows:

The argument of the last exponential term in Eq. (26) is the Gouy phase. The term $L_p^{| l |}$ is the associated Laguerre–Gaussian polynomial obtained from Laguerre polynomials:

We note that for $l = p = 0$ and for the case of negligible beam expansion ${z_R} \to \infty$, the above expression reduces to the form of the Gaussian mode defined in Eq. (24), where $w(0) = \sqrt 2 \sigma$. In order to normalize the amplitude ${A_{{{\rm LG}_{\text{pl}}}}}$ such that the integral of the intensity distribution is equal to $P$, the laser power, the amplitude distribution can be scaled by a parameter $\alpha$, where

Once again, using the same RS lens, a LG profile was applied at the input, i.e., ${u_0}({x_0},{y_0}) = {t_{\text{RS}}}({x_0},{y_0})\alpha {A_{{{\rm LG}_{01}}}}({x_0},{y_0})$, where $w(0)$ was set equal to $R$ and ${z_R} = w_0^2k/2$. Once again the parameters of the simulation using the thin lens algorithm are the same as those defined in Section 3.A. The intensity and phase distributions of the field in the $xz$-plane is shown in Figs. 5(a2) and 5(b2) over a size of $20\;{\rm{\unicode{x00B5}{\rm m}}} \times 20\;{\rm{\unicode{x00B5}{\rm m}}}$ and a clear discontinuity in the phase is observed at $({x_z},{y_z}) = (0,0)$. The intensity and phase in the $xy$-plane over an area of $10\;{\rm{\unicode{x00B5}{\rm m}}} \times 10\;{\rm{\unicode{x00B5}{\rm m}}}$ are shown at five different distances in Figs. 5(c2) and 5(e2) with profiles of these distributions given in Figs. 5(d2) and 5(f2). The spiral phase shape of the distribution is clearly evident and appears to rotate with propagation distance. The video files cited in the Fig. 5 caption also enable visualization of this rotation.

5. DIFFRACTIVE OPTICAL ELEMENTS

An interesting question arises when we consider diffractive optical elements that are designed as lenses: Is such a thin lens function applicable for the thin lens algorithm? If we are dealing with a diffractive optical element that contains diffraction orders, such as a binary phase plate or a spatial light modulator, in order to fully capture the propagating wavefield we would need to sample this input function with a rapid sampling interval that could encompass the spatial frequency support of all of the diffraction orders up to the point at which their energy is sufficiently small to be considered negligible. This densely sampled function could then be input to a numerical propagation, but the algorithm would have to be chosen such that it will produce an output field over the true spatial and spatial frequency supports of such a field. Certainly the ASM would be a poor choice since the support would never expand over that of the input support and zero-padding or chirp-transform would be necessary.

There is a short-cut one could consider, though, that may validate the use of the thin lens algorithm. If one were confident that in the region of interest, for example, the focal plane of the DOE (assuming it is a focusing phase function), the diffraction orders were well separated with negligible overlapping energy, then the thin lens algorithm (or ASM with pre-aliasing) becomes applicable for each order in isolation. This assumption is likely to be valid if the DOE is well pixelated. Taking the SLM as our example, presuming the central order is the only order of real interest, the approach might be as follows. Taking one sample value to represent the phase for each pixel (i.e., sampling interval equals pixel pitch), perform a DFT on this 2D discrete function. This DFT can then be multiplied by a 2D “sinc” function (since the SLM can be described by the convolution of our samples with a pixel “rect” function). Performing an inverse DFT will produce a discrete function could be used as input to the thin lens algorithm. Here we have isolated only the zero order for propagation in the knowledge that it will be indeed be isolated in the focal region.

6. INTERPRETATION OF PHASE

When calculating diffraction patterns via the two algorithms presented in this paper or by any other method, interpreting the phase at points where the amplitude of the wave is negligible presents a unique challenge. The phase of the wavefield at a given point is a critical parameter that describes the position within the wave’s oscillatory cycle at that specific point in space and time and is crucial for determining interference. However, the physical significance of the phase becomes less pronounced at points where the amplitude, or the wave’s strength, is close to zero. In these regions of negligible amplitude, the phase can become mathematically indeterminate or ill-defined. This indeterminacy arises because the phase is derived from the argument of a complex number representing the wave, and when the modulus of this complex number (the amplitude) approaches zero, small fluctuations or uncertainties in the wavefield caused by small errors inherent to the numerical method, or due to imperfect sampling conditions, or even due to machine error, can lead to disproportionately large and undefined variations in the phase value.

In practical calculations, it is common to encounter zones where the amplitude is very small. Physically, these areas with negligible amplitude often correspond to points of destructive interference or locations outside the main propagation path of the wave, where its influence is minimal. It is clear in Fig. 5 (and Fig. 2 in Supplement 1) for example that the amplitude of the field is very low away from the center, and yet the phase values in this region appear to be characterized by an organized pattern. From a computational perspective, one might argue that handling these regions of negligible amplitude requires careful consideration since these phase values must be poorly defined based on the arguments above. We would argue the opposite; these phase values can be ignored since these low-amplitude regions must be less significant in terms of the wave’s physical impact. The phase value in these weak amplitude regions is less critical to the overall behavior of the wave, particularly regarding observable quantities like intensity, which is primarily dependent on the amplitude. For regions where the amplitude falls below a certain threshold, the phase might be considered undefined or treated as less reliable, and we should make no attempt to be overly concerned with its accuracy. Therefore, we should make no attempt to examine the eikonal in these areas in subsequent papers. Indeed we will see in those papers that the eikonal will seek to avoid these areas, always bending away from them towards areas of higher amplitude, which is of course in agreement with the concept of flux.

7. CONCLUSION

In this paper, we have developed two algorithms for the numerical calculation of the three-dimensional complex amplitude of coherent wavefields that are focused by a lens. The first algorithm deals with the case of the thin lens approximation, while the second deals with an ideal lens system. In both cases, the angular spectrum method is applied in a loop in order to generate a stack of complex images over a regular grid around the focal region of the lens. However, the manner in which the ASM is applied differs considerably for both cases; for the thin lens, the ASM is applied directly to the lens transmission function, but the difference in spatial support between the input plane (i.e., the lens aperture) and the output plane (the focused spot) requires a “pre-aliasing” step, recently proposed by other authors [11,19] and further developed in this paper, which although numerically intensive, needs only to be applied once before the ASM is applied iteratively; for the ideal lens case, a single discrete Fourier transform is used to transform between the circular pupil function (in the Fourier domain) of the lens system to the focal plane of the lens, and then the ASM is applied iteratively to generate the stack of complex images immediately before and after this plane. A detailed list of steps is provided for both algorithms in Section 3 that includes the selection of the key lens properties as well as the optimal sampling conditions over the three spatial dimensions.

The thin lens functions that were initially selected for analysis included the case of the spherical lens for which the transmission function is directly related to the thickness of the glass. This was compared with the three-dimensional focus of an equivalent lens defined with a simple quadratic phase, and the results were found to be the same, as expected. A third case was also investigated whereby a transmission function defined by the kernel of the Rayleigh–Sommerfeld diffraction integral was proposed. This “lens” was simulated to produce a focused spot that was much smaller than the equivalent results for the spherical lenses. Indeed, the results were shown to be equivalent to those from an ideal lens with similar numerical aperture using the second algorithm.

In Section 4 both algorithms were augmented to account for the laser spatial mode for the case when the input to the lens is a temporally coherent laser input. In particular, ${{\rm TEM}_{00}}$ (Gaussian) and ${{\rm TEM}_{01}}$ (Laguerre–Gaussian) were investigated. Furthermore, in the third paper both algorithms will be further augmented to account for aberrations in the lens, and for both cases the well-known Zernike polynomials will be used to described the various forms of lens aberration, albeit in different ways for both algorithms.

The computational time for the two algorithms presented in this paper is an important consideration. Disregarding the time-costly initial pre-aliasing step in the thin lens algorithm, which is implemented once and can be stored in memory, both propagation algorithms essentially comprise a single element-wise matrix multiplication and a single DFT operation, both of which can be implemented using parallelized SIMD instruction on a graphics processing unit (GPU) [27]. This process involves transferring the input from the CPU to the GPU over a PCI bus (this can be ignored since it happens only once), implementing the GPU operations and transferring the result back to the CPU. Using an Intel i7-9700K CPU Processor and an Nvidia GeForce GTX 1070 GPU, we were able to implement this overall process in 20 ms for a single propagation of size $2000 \times 2000$ samples. To calculate the volume of $2000 \times 2000 \times 2000$ takes 40 s for both algorithms. This is the time taken to calculate all of the volumes used in all of the results across the three papers. Further speed up could be achieved using higher performance GPUs. The code to implement the two algorithms is provided in Code 1, Ref. [28].

This paper is a useful reference for those wish to simulate focused light over three dimensions with lenses that have numerical aperture ${\lt}0.6$ and can therefore be modeled using the paraxial approximation. It is particularly useful for those wishing to simulated focused laser beams by lenses with aberrations. Although this paper is useful in a broader context, it is the first part in a three-part series; the second part relates to calculating and visualizing the flux of focused light fields, which has been investigated by several authors as an extension of the concept of geometrical ray optics. Calculating the flux involves calculating the derivative of the phase of the wavefield as it propagates. A necessary first step, as described in this paper, is to calculate the complex amplitude of the light over three dimensions with sufficiently fine-grained sampling intervals such that the flux rays can be tracked accurately. The third part highlights the power of the method to provide insights into lens aberrations.