Nonlinear ray tracing in focused fields, part 2. Tracing the flux: tutorial

Qin Yu; Bryan M. Hennelly; Bryan M. Hennelly

doi:10.1364/JOSAA.503925

1. INTRODUCTION

In the first of this three-part paper series, two algorithms were defined that can be used to compute the samples of the complex wavefield on a three-dimensional grid over the focal region of a lens. The differences in these algorithms relate to the physical model used to describe the lens, where one algorithm employs the thin lens approximation and the other models the ideal lens using a Fourier transform. Both algorithms allow for different laser spatial modes to be focused. Either algorithm can be used as a basis for tracing the lines of flux using the methods described in this paper. Understanding how these algorithms could be employed in fast time, and the sampling requirements over the three dimensions governed by the rules of Nyquist-Shannon interpolation, was an important consideration in the first paper since tracing the lines of flux will necessarily require accurate interpolation of the wavefield at arbitrary points within the 3D grid of samples.

In wave optics, the angular spectrum [1,2], which is the basis for the angular spectrum method [2] that is employed at the heart of both of the two algorithms developed in the first paper, represents the decomposition of a wavefront into a superposition of plane waves with different wavevectors (i.e., directions of propagation). Each plane wave in the superposition (which is defined elegantly using the Fourier transform) [3] has a different direction and amplitude, and collectively they form a complete representation of the original wavefront. The angular spectrum is related to the eikonal function [1,2,4], which describes the phase of the wavefront as a function of position and time. The eikonal function can be used to define a “local wavevector” at each point in space, which describes the direction and magnitude of the local wavefront propagation. This local wavevector can be defined in terms of the local gradient of the phase of the wavefront [2]; and it is specifically this quantity that we measure for the wavefield in order to trace the direction of the eikonal function through space. The resultant lines effectively trace the flux, i.e., the movement of the power that is the intensity of the diffraction pattern as it moves fluidly through space, and we use the terms “flux-lines” and “nonlinear ray tracing” interchangeably throughout the text. The use of the term “nonlinear ray tracing” is no doubt controversial and immediately forces a comparison with classical geometrical optics. The concept of the classical “linear-ray” in geometrical optics is also related to the direction of propagation of the local wavevector at each point on the wavefront; the nonlinear ray and the linear ray are one and the same so long as the gradient of the phase is well-described by only the first two terms of its Taylor series expansion; in such case the second derivatives of the phase and higher must be approximately zero [2].

The breakdown of this paper is to first introduce the concept of flux in Section 2. This will begin with an overview of how the angular spectrum can be physically interpreted in terms of plane waves that span all of space that are uniquely defined by a wavevector. This leads to a discussion of the eikonal and the local wavevector, which naturally precipitates a discussion of the ray and the local spatial frequency. This forms a basis for a computational method to trace the flux lines, described in detail in Section 4. Since the method is based on sequentially propagating the ray over small steps, and small errors could quickly accumulate, and an important part of the discussion is an investigation of the accuracy of the proposed methods in terms of tracing along the correct path governed by the eikonal function. This is followed by results for various cases of laser spatial modes.

2. FLUX

A. Physical Interpretation of the Angular Spectrum

In the previous paper we defined the propagation of the angular spectrum in terms of the scalar wavefield ${u_0}({x_0},{y_0})$. Here, we review the physical interpretation of the angular spectrum and optical propagation. This interpretation is based on the concept of the simple plane wave propagating with wave vector $\vec k$, where $\vec k$ has magnitude $2\pi /\lambda$ and has direction cosines $(\alpha ,\beta ,\gamma)$, as illustrated in Fig. 1. Such a plane wave has a complex representation of the form

(1)$$p(x,y,z;t) = \exp \!\left[{j\left({\vec k \cdot \vec r - \frac{{2\pi}}{\lambda}ct} \right)} \right],$$

where $\vec r = x{\hat x} + y{\hat y} + z{\hat z}$ is a position vector (the $\wedge$ symbol signifies a unit vector), while $\vec k = \frac{{2\pi}}{\lambda}({\alpha {\hat x} + \beta {\hat y} + \gamma {\hat z}})$. Dropping the time dependence, the complex phase amplitude of the plane wave across a constant $z$-plane is given by

(2)$$\begin{split}p(x,y,z) &= \exp \!\left[{j\vec k \cdot \vec r} \right] = \exp \!\left[{j\frac{{2\pi}}{\lambda}(\alpha x + \beta y)} \right]\\&\quad\times\exp \!\left[{j\frac{{2\pi}}{\lambda}\gamma z} \right],\end{split}$$

where the direction cosines are interrelated through

(3)$$\gamma = \sqrt {1 - {\alpha ^2} - {\beta ^2}} .$$

Thus, across the plane $z = 0$, a complex-exponential function $\exp [{j2\pi ({f_x}x + {f_y}y)}]$ may be regarded as representing a plane wave propagating with direction cosines,

(4)$$\alpha = \lambda {f_x},\quad \beta = \lambda {f_y},\quad \gamma = \sqrt {1 - {{(\lambda {f_x})}^2} - (\lambda {f_y})^2} .$$

In the Fourier decomposition of ${u_0}$, the complex amplitude of the plane-wave component with spatial frequencies $({f_x},{f_y})$ can be rewritten in terms of the direction cosines as follows:

(5)$$\begin{split}{{U_0}\left({\frac{\alpha}{\lambda},\frac{\beta}{\lambda}} \right) = \iint_{- \infty}^{+ \infty} {u_0}(x,y)\exp \!\left[{- j2\pi \left({\frac{\alpha}{\lambda}x + \frac{\beta}{\lambda}y} \right)} \right]{{\rm d}_x}{{\rm d}_y}}\end{split},$$

which is called the angular spectrum of the disturbance ${u_0}(x,y)$. Optical propagation can be interpreted in terms of this angular spectrum by simply adding a delay (which is a function of the plane wave direction) to each of these angular plane waves as follows:

Fig. 1. Illustration of the wave vector $\vec k$ for a plane wave uniquely defined by the angle cosines $\alpha$, $\beta$, and $\gamma$.

Download Full Size | PDF

(6)$$\begin{split}{{U_z}\left({\frac{\alpha}{\lambda},\frac{\beta}{\lambda}} \right) = {U_0}\left({\frac{\alpha}{\lambda},\frac{\beta}{\lambda}} \right)\exp \!\left({\frac{{j2\pi z\gamma}}{\lambda}} \right)}\end{split}.$$

B. Concept of a Ray as a Flux Trace

A key part of this paper is the concept of a ray in terms of wave optics. The geometrical ray can be related to wave optics through a mathematical function named the eikonal [1]. The eikonal describes a surface in space where the wave has the same phase and can also be interpreted as the path that the ray takes through space such that a constant phase is maintained along that path. The disturbance $u(\vec r)$ (that we have previously described with the notation ${u_z}(x,y)$) can be rewritten as follows:

(7)$$u(\vec r) = a(\vec r)\exp \!\left[{- \frac{{j2\pi}}{\lambda}S(\vec r)} \right],$$

where $a(\vec r)$ is the real-valued amplitude and ${-}2\pi /\lambda S(\vec r)$ is the phase of the wave; the refractive index $n$ of the medium is contained in the definition of $S$, which is called the eikonal function. Surfaces defined by $S(\vec r) =$ constant are called wavefronts of the disturbance. The flux (i.e., the direction of power flow) and the direction of the wave vector $\vec k$ are both normal to these wavefronts at every point $\vec r$ in an isotropic medium. A ray, therefore, is defined as a trajectory or a path through space that begins at any point on the wavefront and moves through space with the wave, always remaining perpendicular to the wavefront at every point on the trajectory. Thus, a ray traces out the flux in an isotropic medium. If the disturbance $u(\vec r)$ is to represent an optical wave, it must satisfy the scalar wave equation,

(8)$${\nabla ^2}u - \frac{{{n^2}}}{{{c^2}}}\frac{{{\partial ^2}u}}{{\partial {t^2}}} = 0,$$

where ${\nabla ^2}$ is the Laplacian operator, $n$ represents the refractive index of the dielectric medium within which light is propagation, and $c$ represents the velocity of light in vacuum. Substituting $u$ as defined in Eq. (7) into the Helmholtz equation yields the following equation that must be satisfied by both $a(\vec r)$ and $S(\vec r)$:

(9)$$\begin{split}{\left({\frac{{2\pi}}{\lambda}} \right)^2}&\left[{{n^2} - |\nabla S{|^2}} \right]a + {\nabla ^2}a\\ &- \left({\frac{{j2\pi}}{\lambda}} \right)\left[{2\nabla S \cdot \nabla a + a{\nabla ^2}S} \right] = 0.\end{split}$$

Examining the real part only reveals the following relationship:

(10)$$|\nabla S{|^2} = {n^2} + {\left({\frac{\lambda}{{2\pi}}} \right)^2}\frac{{{\nabla ^2}a}}{a}.$$

Setting $\lambda \to 0$ produces the so-called eikonal equation, which is a fundamental equation in geometrical optics:

(11)$$|\nabla S(\vec r{)|^2} = {n^2}(\vec r).$$

This equation serves to define the wavefront in terms of the eikonal function, $S$. Once the wavefronts are known, the trajectories defining rays can be determined. The eikonal can be considered as a bridge between wave optics governed by the Helmholtz equation and ray optics.

C. Rays and Local Spatial Frequency

It is also possible to approximately define the ray trajectory in terms of the local spatial frequency of the disturbance $u$, and in some cases, such as for a plane wave, this definition is exact. Since each of the spatial frequency components of the angular spectrum will have an infinite support over $(x,y)$, it is not possible to directly relate a specific point in a disturbance with a single spatial frequency for the general case. However, if the disturbance is a wavefront described by slowly varying amplitude and phase functions $a(x,y)$ and $\phi (x,y)$, each point in the disturbance can be approximately interpreted as being related to a single frequency. Goodman [2] reports that this interpretation is accurate if $\phi (x,y)$ varies sufficiently slowly such that it is well approximated by only three terms of its Taylor series expansion about any point, i.e., a constant term and two first-partial-derivative terms. We begin by rewriting the disturbance using similar notation as in the first paper:

(12)$${u_z}(x,y) = {a_z}(x,y)\exp [{j{\phi _z}(x,y)} ].$$

The local spatial frequency of ${u_z}$ is defined as a frequency pair $({f_{\textit{lx}}},{f_{\textit{ly}}})$ as follows:

(13)$${f_{\textit{lx}}} = \frac{1}{{2\pi}}\frac{\partial}{{\partial x}}{\phi _z}(x,y),\quad {f_{\textit{ly}}} = \frac{1}{{2\pi}}\frac{\partial}{{\partial y}}{\phi _z}(x,y).$$

In Section 2.A we demonstrated that ${u_z}$ could be decomposed by means of a Fourier transform into a collection of plane-wave components traveling in different directions, each defined by unique wave vector with direction cosines $(\alpha ,\beta ,\gamma)$ defined by Eq. (3). The spatial frequencies defined through the Fourier decomposition exist at all points in space and cannot be regarded as being localized; however, for slowly varying disturbances, the definitions of the previously defined local spatial frequencies $({f_{\textit{lx}}},{f_{\textit{ly}}})$ can be interpreted as defining the local direction cosines $({\alpha _l},{\beta _l},{\gamma _l})$ of the wavefront through the relations:

(14)$${\alpha _l} = \lambda {f_{\textit{lx}}},\quad {\beta _l} = \lambda {f_{\textit{ly}}},\quad {\gamma _l} = \sqrt {1 - {{(\lambda {f_{\textit{lx}}})}^2} - (\lambda {f_{\textit{ly}}}{)^2}} .$$

These local direction cosines are the direction cosines of the ray through the $(x,y)$ plane at each point. This definition of a ray is consistent with the definition of a ray that was given earlier in terms of the eikonal function and is the basis for the nonlinear ray tracing method proposed in this paper.

3. LIMITATIONS OF GEOMETRIC OPTICS AND THE EIKONAL APPROXIMATION

A. Relationship between the Eikonal and Ray Optics

It is important to emphasize that the key principle of geometric optics (or ray tracing) is fundamentally based on the eikonal approximation. In the eikonal approximation, the normals to the wavefronts (surfaces of constant phase) effectively behave like rays in geometric optics. This then leads to the concept of light rays traveling in straight lines when the medium is uniform and without significant diffraction or interference effects, which are central assumptions in geometric optics. While reflection and refraction can be accounted for in the geometric approximation, wave effects cannot. One might state that geometric optics (or ray optics) therefore is a further approximation of the eikonal approximation, and in those areas where wave effects such as diffraction or interference cannot be ignored, then a more thorough treatment of the eikonal approximation can be considered.

B. Failure of the Eikonal Approximation

The eikonal approximation is only valid in regions where we can say that the amplitude varies more slowly that the phase. There are several examples of where this approximation very obviously breaks down, including near field effects and propagation through media with rapidly changing refractive indices or dispersive properties, as well for cases of short wavelength for which classical wave theory begins to break down. More nuanced examples include optical speckle [5,6] and caustics [7,8]. Speckle patterns arise from the coherent superposition of multiple scattered waves and are characteristic of rapid changes in amplitude and phase, and consequently the eikonal approximation is inadequate for describing such phenomena. Speckle patterns represent an extreme example that can be described only by the wave model and cannot be interpreted in terms of geometrical ray optics. We can imagine, therefore, that the eikonal, which is in essence a semiclassical approach that goes some way to bridge the gap from geometrical to wave optics, cannot be used to describe speckle formation. Conversely, caustics are a concept that originates from geometric optics, and their translation into wave optics leads to disagreement between the two models. In geometric optics, caustics are well-defined in their origin as the envelopes of light rays that converge after refracting through or reflecting off complex surfaces, leading to regions where rays converge to form bright light in the form of singularities of infinite intensity. Since such singularities cannot exist it is clear that the geometric approximation has failed. In wave optics, caustics can still be described as regions where light waves converge and interfere constructively, leading to high intensity. These intense regions (where geometric optics predicts singularities) can only be described accurately by diffraction effects using the wave model. Similar to the case of speckle, since only the wave model can be used in the region of the caustic, the interpretation of the eikonal cannot be applied in caustic regions.

C. Identifying When the Eikonal Approximation Breaks Down

Mathematically, we can state that the eikonal approximation fails in regions where the amplitude varies quickly relative to the phase. Such regions are typically characterized by rapid changes in phase with high curvature where the field cannot be described as smooth [1]. In terms of physical optics, we can say that the approximation fails where the geometric approximation breaks down. For the case of caustics and speckle described above we can say that both the mathematical and physical interpretations of breakdown are true.

So what then for focal regions of a high numerical aperture lens? This is more nuanced. The focal region of a high NA lens is characterized by significant curvature of the electromagnetic field, and this curvature leads to a more complex phase distribution as one approaches the focal point. However, the field within the focal region can still be considered relatively smooth compared to phenomena like caustics, where the wavefronts experience abrupt changes.

For smooth, non-zero optical fields, the geometric optics approximation is valid at least in the neighborhood of the source plane. For slowly varying source fields, geometric optics may be valid for relatively long propagation distances, and this is the classical application of the approximation used in ray-tracing. On the other hand, if the optical field is rapidly changing or the source plane is near a zero of the field, it may only be valid for a very short propagation distance. Nevertheless, if the source field is smooth and non-zero, then it will always be valid for some propagation distance. This fact is the fundamental mathematical justification for this entire series of articles.

We emphasize that the approximation is valid at least in a neighborhood of the originating plane, so long as within this neighborhood the field exhibits variations in phase and amplitude that are gradual enough not to necessitate wave optics for accurate description. The assumption of smoothness ensures that the light’s phase front does not undergo abrupt directional changes, at least over short propagation distances, such as those used with the angular spectrum method in these papers. In Section 6.B below we highlight a number of methods that can be used to probe the validity of the geometric approximation along the path of a computed eikonal, and we identify markers where the approximation strains.

D. Beyond the Eikonal Approximation

The eikonal approximation is an example of a semiclassical method that bridges between the wave nature of light (as described by wave optics) and the ray-like behavior of light (as described in geometric optics). However, this bridge is limited in extent and leans more heavily towards geometrical optics than wave optics. Therefore, in those regions (such as speckle and caustics) where the geometric optics description utterly fails, so too does the eikonal approximation. Other semi-classical models offer a more complete bridging between the two models. One such well-known semiclassical model in optics is the Wentzel-Kramers-Brillouin (WKB) approximation [9]. Like the eikonal, this method is used to find solutions to the wave equation in situations where the wavelength of the light is small compared to the characteristic scale of variation in the system. Unlike the eikonal, however, WKB can account for regions where geometric optics fails, such as for caustics. The WKB method is often applied in quantum mechanics but has its analogues in optical theory [10]. WKB bridges the gap between geometric and wave optics by using a wave-based approach to derive ray-like equations. It treats the phase of the wave in a manner similar to the eikonal approximation but incorporates the varying amplitude of the wave, accounting for some wave effects. WKB is more complex for calculating the eikonal, but it is important to acknowledge the method as being preferable in those regions where the eikonal approximation fails or begins to strain.

4. TRACING THE FLUX

A. Method

The crux of the proposed method is that we iteratively trace the flux between subsequent $z$-planes (separated by a distance ${h_z}$ that is selected by the user) making use of the local spatial frequencies defined in Eq. (13) to guide the ray angle at each plane. This is illustrated in Fig. 2(a) for which case 16 rays, with initial points uniformly arranged in a circle around a focusing Gaussian laser, have been traced through the focal plane. The method is illustrated in more detail in Fig. 2(b), where we highlight a single ray close to the focal plane. It is important to note that field ${u_z}(x,y)$ will have been sampled over a 3D grid with uniform sampling intervals ${T_x}$, ${T_y}$, and ${T_z}$ as described in the first paper; however, the 3D positions at which the ray is sampled will not, in general, be coincident with the points on this grid. Therefore, interpolation is required in order to determine the complex amplitude of ${u_z}(x,y)$ at the desired points. For a given point on the ray $({x_0},{y_0},{z_0})$, the field must be interpolated at four points, two of which are slightly displaced in $x$ ($({x_0} + {\delta _x},{y_0},{z_0})$ and $({x_0} - {\delta _x},{y_0},{z_0})$) and two that are displaced in $y$ ($({x_0},{y_0} + {\delta _y},{z_0})$ and $({x_0},{y_0} + {\delta _y},{z_0})$); this facilitates calculation of the angle cosines ${\alpha _l}$ and ${\beta _l}$ according to Eqs. (13) and (14). The ray then traces a straight line from $({x_0},{y_0},{z_0})$ to its new position $({x_0} + {h_x},{y_0} + {h_y},{z_0} + {h_z})$, where ${h_x}$ and ${h_y}$ are defined in terms of the angle of the wave vector as follows:

(15)$$\begin{split}{h_x} = {h_z}\tan (\arccos ({\alpha _l})),\\{h_y} = {h_z}\tan (\arccos ({\beta _l})),\end{split}$$

where

(16)$$\begin{split}{\alpha _l} = \frac{\lambda}{{2\pi}}\frac{{\phi ({x_0} + {\delta _x},{y_0},{z_0}) - \phi ({x_0} + {\delta _x},{y_0},{z_0})}}{{2{\delta _x}}},\\{\beta _l} = \frac{\lambda}{{2\pi}}\frac{{\phi ({x_0},{y_0} - {\delta _y},{z_0}) - \phi ({x_0},{y_0} + {\delta _y},{z_0})}}{{2{\delta _y}}}.\end{split}$$

Equation (15) can be rewritten more simply in terms of ${\gamma _l}$ as follows:

(17)$$\begin{split}{h_x} = {h_z}\frac{{{\alpha _l}}}{{{\gamma _l}}},\\{h_y} = {h_z}\frac{{{\beta _l}}}{{{\gamma _l}}},\end{split}$$

where ${\gamma _l}$ is calculated using Eqs. (14) and (16). This method is then repeated at the new point resulting in another straight line being tracked to the next plane. The method, therefore, is to trace straight lines between subsequent $z$-planes in a piece-meal fashion. For a general disturbance, this method is accurate so long as the choice of ${h_z}$ is small, which is discussed in more detail in the next section.

Fig. 2. (a) Illustration of 16 rays uniformly arranged in a circle around a focusing Gaussian laser being traced through the focal plane being guided by the local spatial frequency of the field over subsequent $z$-planes; (b) illustration of the method whereby the ray propagates in straight lines between the $z$-planes separated by a user-selected distance ${h_z}$. Note that the field must be interpolated at four points around the ray positions in order to calculate the phase derivatives, which are used to guide the ray to the next $z$-plane.

Download Full Size | PDF

B. Flux Tracing Algorithm

If only a single ray were to be considered, the Rayleigh-Sommerfeld (RS) integral could be calculated with high accuracy using direct numerical integration as described in the first paper. The ray origin $(X,Y)$ is first selected somewhere within the support of ${u_0}$ at distance $z = 0$, and the derivative of the phase of ${u_0}$ at $(X,Y)$ can be calculated by using direct numerical integration to calculate the complex values of ${u_0}$ at four points around $(X,Y)$. For a given step size ${h_z}$, these derivatives are used to determine the new ray positions $X + {h_x}$ and $Y + {h_y}$ as previously described. Direct numerical integration is relatively slow; based on our experiments this method could be employed to trace a single eikonal trace through the focal region over a distance of 20 µm in the order of one minute using a value of $hz \le 10\;{\rm nm} $. It is, therefore, preferable to employ the thin lens algorithm or the ideal lens algorithm when one wishes to trace several eikonals. Our approach can generate ${\gt}{1000}$ traces in the same time frame.

One approach to employ the thin lens algorithm or the ideal lens algorithm would be to use these algorithms to calculate the 3D grid of samples over the entire region of interest as described in the previous paper, using the Nyquist sampling intervals for ${T_x}$, ${T_y}$, and ${T_z}$. It is then possible to identify a starting position for several rays and to trace these rays through the 3D grid by performing 3D interpolation sequentially as the ray progresses in straight lines over steps of length ${h_z}$ along the $z$ axis. This approach is consistent with the graphical representation given in Fig. 2(b). Three-dimensional interpolation can be implemented to calculate the field at these arbitrary 3D points by applying a suitable interpolation scheme, such as a three-dimensional form of linear or cubic spline interpolation to the 3D grid of samples. MATLAB’s interp3 function can perform this operation. However, we have found that 3D interpolation is relatively slow compared with 2D interpolation and requires large 3D matrices to be stored and operated on.

We employ an alternative approach, which is more straightforward and is based on simply selecting the sampling interval to be small such that ${T_z} = {h_z}$ and then applying 2D interpolation in the finely sampled $z$-planes to calculate the field at the four points around each new ray position. In our simulations we employ MATLAB’s interp2 function with spline interpolation. Compared with 3D interpolation this approach requires less memory and time to perform a given interpolation step, it is simple to implement and produces reliable results. The steps of the algorithm are illustrated in Fig. 3 and are defined as follows:

(1) The first step is to select the initial ray origins. A set of $N$ points is defined as $({x_n},{y_n},{z_0})$ for $n{:}1 \to N$, where each $({x_n},{y_n})$ is a unique point in the initial $z$-plane at arbitrary distance ${z_0}$. Values for $\delta x$ and $\delta y$ must also be selected relating to the discussion in Section 4.A. The value of these shifts should be small, but not so small as to result in machine error; we have found the value of $1{e^{- 13}}{\rm m}$ to provide excellent results.
(2) The second step is to apply the thin lens algorithm or the ideal lens algorithm. This includes selecting the properties of ${u_z}(x,y)$ including the aperture radius, $R$, and the Zernike phase term $Z$ and laser mode $A$. This also includes the selection of the sampling intervals ${T_x}$, ${T_y}$, and ${T_z}$ and the number of samples ${N_x}$, ${N_y}$, and ${N_z}$. Note ${T_z}$ is chosen to be very small (equal to ${h_z}$) to ensure accuracy as described in the previous section. The result of Step 2 is that a 3D grid of samples of ${u_z}(x,y)$ is calculated over a region around the focal point of the lens. However, only seven slices are calculated and stored at any one time in order to reduce the memory load; the rays are propagated through these seven slices as described in the next step.
(3) 2D interpolation is applied in the region of each point $({x_n},{y_n},{z_0})$ as described in Section 4.B. For each point, the local cosine angles, ${\alpha _{\textit{ln}}}$ and ${\beta _{\textit{ln}}}$ are calculated according to Eq. (16). Updated ray positions are given by $({x_n},{y_n},{z_0}) \to r({x_n} + {h_x},{y_n} + {h_y},{z_0} + {h_z})$, where a value for ${h_z}$ is defined by the user; in all of the results presented in this paper a value of ${h_z} = 10\,{\rm nm}$ was selected. We proceed to the next $z$-slice, and this step is repeated iteratively seven times until the ray propagates to the last $z$-slice. We then return to Step 2, and another seven $z$-slices are calculated. Steps 2 and 3 are applied repeatedly until the rays have completed their paths over the desired propagation distance ${N_z}{T_z}$.

We note that this choice of seven slices relates to (i) facilitating 3D interpolation and (ii) GPU acceleration. 3D interpolation requires a certain number of slices to be calculated and stored in memory, and this number will depend on the type of interpolation that is being employed, the minimum number being two for linear interpolation. For more accurate methods, such as spline interpolation, it will be necessary to compute and store more than two slices. Although we have not exhaustively examined all interpolation methods, we believe that seven slices chosen with appropriate ${z}$-sampling ${T_z}$ will be sufficient for most 3D interpolation schemes.

Fig. 3. Illustration of the key steps in the flux tracing algorithm. The second step relates to running one of the two algorithms that were defined in the first paper to generate a 3D grid of samples that is stored in memory. The third step is iteratively applied to propagate the rays in straight lines from one $z$-slice to the next using the derivative of the phase to define the ray angle. To reduce the memory load only a small number of $z$-slices are calculated and stored any one time.

Download Full Size | PDF

If only 2D interpolation is applied, then it is necessary to store only a single slice at one time; the ray directions can be computed and traced to the next plane, and then the next slice can be calculated and so on in an iterative manner. Although we apply 2D interpolation, we choose to compute and store seven slices because it is computationally more efficient to so. Both of the algorithms outlined in the first paper can be implemented using a GPU (the FFT operations, and the elementwise matrix multiplication can both be parallelized for SIMD and implemented at high speed on the GPU [11]). However, the interpolation step does not lend itself to SIMD implementation on a GPU; therefore, this step and the subsequent geometrical ray tracing are implemented on the CPU. It is possible, therefore, to run the propagation algorithm over a given block of several slices on the GPU and then transfer the data to the CPU, while simultaneously implementing the interpolation and tracing steps on the previous block on the CPU. In this way we can trace many eikonals without any additional time delay over running the propagation algorithms, which for the CPU and GPU processing units was given by 40 s to compute all of the volumes and the traces that are calculated in the three papers using both algorithms; see the first paper for further details.

5. RESULTS FOR DIFFERENT LASER SPATIAL MODES

Geometrical ray tracing is typically applied without consideration of laser modes. In some cases, such as for example the calculation of optical trapping forces in the ray-optics regime, the rays are given weights according to the amplitude of the electromagnetic wavefield, i.e., relating to their position across the laser mode profile. In that case, the rays in the center of a Laguerre-Gaussian trapping laser will have zero weight in the trapping force calculation, while those at the edge are stronger [12]. However, this approach is based on linear ray paths regardless of the laser’s spatial mode. In this section, we show the results of nonlinear ray tracing applied to focusing lasers with Gaussian and Laguerre-Gaussian spatial modes. These rays follow the lines of flux of the wavefront as it propagates. In Fig. 4 the results are shown for the same ${{\rm TEM}_{00}}$ and ${{\rm TEM}_{01}}$ laser modes focused using the same Rayleigh-Sommerfeld lens investigated in Section 4 and Fig. 5 in the first paper. All parameters used in the thin lens algorithm are identical except that the value of ${T_z}$ is set equal to ${h_z}$ and a reduced number of slices are calculated in a given stack.

Fig. 4. Nonlinear ray tracing for the same ${{\rm TEM}_{00}}$ and ${{\rm TEM}_{01}}$ laser modes focused using the same Rayleigh-Sommerfeld lens with a numerical aperture of 0.6 investigated in Section 4 and Fig. 5 in the first paper: (a1) ${{\rm TEM}_{00}}$ case with log scale of intensity in the background and focal region magnified in (a2); (b1) 3D image of rays spiraling for the ${{\rm TEM}_{01}}$ case. The two inset images show the initial ray positions at ${-}{20}\;{\unicode{x00B5}{\rm m}}$ and at the focal plane; (a2)–(a4) show three different projects of the ray paths. An equivalent figure is shown for the ideal lens (NA 0.6) case computed using the ideal lens algorithm in Supplement 1. Visualization 1 and Visualization 2 show videos of this ray tracing for the ${{\rm TEM}_{00}}$ and ${{\rm TEM}_{01}}$ cases.

Download Full Size | PDF

In Fig. 4(a1) 9 rays are shown for the ${{\rm TEM}_{00}}$ case, and all but one are on the left side of the image. The intensity is shown in the background on a log scale such that the effects of diffraction due to apodization of the lens aperture are more clearly visible. Far from the focus, most of the rays appear to be linear and are similar to the classical geometrical rays. However, close to the focus, which is highlighted in the inset image Fig. 4(a2), the rays bend rapidly before straightening to vertical paths at the focus and then bend outwards again at an opposite angle to arrive at the focus. Along the axis of propagation, some of the strongest effects of diffraction can be observed resulting in intermittent dark and bright intensity regions, and the rays close to this optical axis are seen to propagate around these dark regions. This feature is expected from the idea of flux.

In Fig. 4(b1) 10 rays are shown for the ${{\rm TEM}_{01}}$ case. The positions of these rays are selected to be uniformly separated along the $y$-axis at a distance of ${-}20\,\,{\unicode{x00B5}{\rm m}}$ before the focus as illustrated in the top inset image. The second inset image shows the positions of these rays at the focal plane. In all cases the rays have moved away from the $y$-axis. The 3D paths that these rays have taken between the two planes can be seen in Fig. 4(b1), and the various 2D projections of these 3D paths are shown in Figs. 4(b2)–4(b4). In all cases, the rays are seen to spiral, with the rate of twisting dependent on the distance from the optical axis and the distance to the focal plane. Most noticeably, the rays close to the center twist rapidly as they approach the focal plane.

A deeper investigation of the ray paths for the Laguerre-Gaussian modes is presented in Fig. 5, once again for the same Rayleigh-Sommerfeld lens and with numerical propagation implemented using the thin lens algorithm. All of the simulation parameters are the same as for the previous case except for the choice of $l$ and $p$ in the definition of the LG-modes $L_p^l$ given in the first paper. The $p$ value varies across the columns from 0 to 2, while the $l$ value varies over the rows from 1 to 3. For the first row, corresponding to $p = 0$, the intensity at the focal plane is shown for the three $l$ values, as well as the intensity at a distance of ${-}{20}\;{\unicode{x00B5}{\rm m}}$ before the focus; it can be seen that the size of the doughnut increases as a function of $l$. In each of the three cases, eight ray origins were selected close to the outer rim of the doughnut in the $z = f - 20\,\,{\unicode{x00B5}{\rm m}}$ plane; these are shown as red dots in the intensity images. Another eight ray origins were selected close to the inner rim, which are shown as green dots. These rays were traced over a distance of 40 µm and are shown as red and black lines in the 3D plot. Although the doughnut converges as the beam focuses, it is interesting to note that the relative position of the rays within the doughnut remains unchanged in the intensity images at $z = f - 20\,\,{\unicode{x00B5}{\rm m}}$ and at $z = f\,\,{\unicode{x00B5}{\rm m}}$ for all three cases. It is also interesting to note that the inner rays spiral rapidly, while the outer rays twist only slightly. The rate of angular velocity of the inner rays as a function of $z$ is particularly striking with the rays appearing to undertake a full $2\pi$ revolution over a distance of $z$ that is shorter than the Nyquist distance defined in the first paper; this is particularly noticeable for the $l = 3$ case. The rate of angular velocity of the rays is clearly dependent on the ray position relative to the center of the doughnut together with the factor $l$; focusing on the inner rays, the rate of spin appears to approximately linearly increase as a function of $l$. This is consistent with the concept of angular momentum associated with Laguerre-Gaussian beams, which is directly proportional to the value of $l$; per photon, this momentum is given by $hl$, where $h$ is the reduced Planck constant [13].

Fig. 5. Nonlinear ray tracing for the different Laguerre-Gaussian modes focused using the same Rayleigh-Sommerfeld lens with a numerical aperture of 0.6 investigated in Section 4 and Fig. 5 in the first paper: The values of $l$ and $p$ vary across rows and columns, respectively. Also shown in the figure are the intensity patterns at the focal plane and 20 µm above, which highlight the ray positions in these planes. It is interesting to note that the spin rate of the rays increases linearly as a function of $l$ and also as a function of the position of the ray relative to the center of the doughnut. An equivalent figure is shown in Supplement 1 for the ideal lens (NA 0.6), computed using the ideal lens algorithm.

Download Full Size | PDF

A similar trend is observed for the second row, corresponding to $p = 1$. For all three cases, there exist two distinct doughnuts in the focal plane, although the size and shape of these doughnuts varies significantly as a function of $l$. The diffraction pattern at $z = f - 20\,\,{\unicode{x00B5}{\rm m}}$ is more complex than for the previous case of $p = 0$. It is interesting to observe the movement of the rays as they approach the focal plane; most of the rays diverge towards the outer doughnut, with many requiring a somewhat sudden and significant change in direction as they approach the focal plane in order to do so. Although it appears that the spin rate of the rays is less than for the previous examples for $p = 0$, it must be the case that the aggregate of the spin over all of the rays (weighted according to the wavefield intensity) must be the same given that the angular momentum is related only to $l$ and not to $p$. Similar trends are observed for the $p = 2$ case shown in the bottom row in the figure, although in this case the ray lines are less smooth. Similar results are shown for the ideal lens case using the ideal lens algorithm in Supplement 1.

6. LIMITATIONS AND ACCURACY

There are two fundamental and distinct sources of error to consider in the approach outlined in this paper.

The first and most fundamental source of error is that associated with the eikonal approximation itself. While one may be able to accurately trace the line governed by the phase gradient with low error, we cannot necessarily apply a physical meaning to this “ray” to be following the path of light propagation. We must consider those cases where the eikonal approximation is invalid and make some effort to estimate regions along the traced line where the approximation is questionable.

The second is the error associated with tracing an accurate path defined by the phase gradient. Here, we make no assumption about the validity of the eikonal approximation. We are focused only on accurately tracing the line that is governed by the local phase gradient, and we are not necessarily focused on its physical interpretation. The error here is captured by the expression derived in Section 6.A and can be mitigated by taking (small linear ray) steps over a given segment, or assuming that higher order phase derivatives are negligible. A more intuitive way to interpret this error is in terms of the accuracy of the geometric approximation over a given segment. Two methods to estimate this error are described below in Section 6.B.

A. Accuracy in Tracing the Phase Gradient

In this section we address the error associated with tracing a line governed by the phase gradient in a three-dimensional optical wavefield and make no assumptions on the eikonal approximation. We derive an expression to relate this error to step size and wave function. In the section that follows we outline some simple methods to evaluate this error in our numerical examples.

Following from the definition of the eikonal, the ray must travel a path that is perpendicular to the phase front [1]. Therefore, the phase that it accumulates between two points along its path must be given by $\lambda /2\pi$ times the path length that was travelled between the two points. In general, this path will not be a straight line. However, if ${h_z}$ is sufficiently small, then the path between the two points is in good agreement with linearity, and this is the basis of the method proposed here for tracing the flux. Assuming linearity, we can define the phase accumulation between the two points as follows:

(18)$$\begin{split}&{\phi _{{z_0} \;+\; {h_z}}}({x_0} + {h_x},{y_0} + {h_y}) \\&= {\phi _{{z_0}}}({x_0},{y_0}) + \frac{{2\pi}}{\lambda}\sqrt {{h_x}^2 + {h_y}^2 + {h_z}^2}\end{split}.$$

Substituting the values for ${h_x}$ and ${h_y}$ from Eq. (17) produces the following relationship:

(19)$$\begin{split}{{\phi _{{z_0} \;+\; {h_z}}}({x_0} + {h_x},{y_0} + {h_y}) = {\phi _{{z_0}}}({x_0},{y_0}) + \frac{{2\pi {h_z}}}{{\lambda {\gamma _l}}}}\end{split}.$$

This equation is the model for the method proposed in this paper. If the path of the eikonal function is not linear between the two points, this will result in a phase error given by

(20)$${{\rm error} = {\phi _{{z_0} \;+\; {h_z}}}({x_0} + {h_x},{y_0} + {h_y}) - {\phi _{{z_0}}}({x_0},{y_0}) - \frac{{2\pi {h_z}}}{{\lambda {\gamma _l}}}}.$$

Interestingly, the same two-phase values can be related, without any approximation, by applying a three-dimensional Taylor series expansion to the continuous phase field, ${\phi _z}(x,y)$, as follows:

(21)$$\begin{split}&{\phi _{{z_0} \;+\; {h_z}}}({x_0} + {h_x},{y_0} + {h_y})\\& = {\phi _{{z_0}}}({x_0},{y_0}) + \sum\limits_{n = 1}^\infty \frac{1}{{n!}} \left[{h_x}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta x}}}_{{x = {x_0}\atop{y = {y_0}\atop z = {z_0}}}}} \right.\\&\quad+\left. {h_y}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta y}}}_{{x = {x_0}\atop{y = {y_0}\atop z = {z_0}}}}} +{h_z}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta z}}}_{{x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} \right]^n,\end{split}$$

where the notation has the following meaning: ${[\delta {\phi _z}(x,y)/\delta x]^n} = {\delta ^n}{\phi _z}(x,y)/\delta {x^n}$. Therefore, the error function in Eq. (20) can be rewritten in terms of the partial derivatives of the phase:

(22)$$\begin{split}{\rm error}& = - \frac{{2\pi {h_z}}}{{\lambda {\gamma _l}}} + \sum\limits_{n = 1}^\infty \frac{1}{{n!}}\left({\frac{{{h_z}}}{{{\gamma _l}}}} \right)^n\\&\quad\times \left[{\alpha _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta x}}}_{{x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} + {\beta _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta y}}}_{{ x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} \right.\\&\quad+\left. {\gamma _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta z}}}_{{ x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} \right]^n.\end{split}$$

Using the substitutions from Eq. (14) this further reduces as follows:

(23)$$\begin{split}{\rm error} &= {\sum\limits_{n = 2}^\infty \frac{1}{{n!}}{{\left({\frac{{{h_z}}}{{{\gamma _l}}}} \right)}^n}}\\&\quad\times \left[{\alpha _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta x}}}_{{x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} + {\beta _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta y}}}_{{ x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}}\right.\\&\quad +\left. {\gamma _l}{{\frac{{\delta {\phi _z}(x,y)}}{{\delta z}}}_{{ x = {x_0}\atop{ y = {y_0}\atop z = {z_0}}}}} \right]^n.\end{split}$$

This error function defines the phase error between two subsequent points that have computed. If the second point were truly along the flux line defined by the eikonal function, this error function would reduce to zero. If this error function is not close to zero, then the second point on the ray has not been selected correctly using the method defined in Eq. (17). In this case, it could be expected that further errors will follow with subsequent ray traces, and that the error would propagate resulting in an erroneous flux line.

Some interesting observations can be made from Eq. (23). First, for the case of a plane wave, such as that defined in Eq. (2), the error is zero regardless of the value of ${h_z}$. This is due to the fact that all partial derivatives are zero for $n \gt 1$. This makes perfect sense since the plane wave is unidirectional; therefore, any step size ${h_z}$ will produce an error free value of ${h_x}$ and ${h_y}$. A second observation is the dependence of the error on the term $h_z^n$, where the smallest value of this term is $n = 2$, which will be applied to the second derivatives only. The error can be expected to be small if (i) the partial derivatives are small for $n \gt 1$ and/or (ii) ${h_z}$ is small. In the results presented in the remainder of this paper, we make no assumptions about the value of the higher order partial derivatives; instead, we make sure to use very small values of ${h_z} \ll 1/\Delta {f_z}$ (that are significantly smaller than the Nyquist interval) in order to guarantee accurate flux lines.

B. Estimating the Accuracy of the Geometric Approximation over a Computed Eikonal Path

In this section, we conceive methods to provide insight into the accuracy of a computed eikonal path. A critical consideration is the applicability of the geometrical optics approximation between subsequent planes.

1. WKB

Given the difficulties in obtaining exact solutions for comparison, one approach to validate the geometric optics approximation is to apply the WKB series as a solution to the Helmholtz equation in the form $u(x,y,z) \approx A(x,y,z)\exp \def\LDeqbreak{}(jk({S_0}(x,y,z)\; +\; \frac{1}{k}{S_1}(x,y,z) \;+\; \frac{1}{{{k^2}}}{S_2}(x,y,z) \;\;+\;\; \cdots))$, where ${S_0}$ and ${S_1}$ represent the phase and amplitude. Solving for the higher order term ${S_2}$ and ensuring $(jk){S_0} \ll {S_1} \ll {S_2}{(jk)^ {-1}} \ll 1$ would provide a rigorous verification of the geometric approximation at a given point $(x,y,z)$ [9]. Since the amplitude and phase are known, in so much accuracy as can be provided by the ASM methods described in the first paper, it should be possible to provide this verification at every point along the path at which we recompute the ray direction. Although this approach is rigorous, it is often difficult to solve for ${S_2}$ and higher orders using numerical methods, especially given the large number of points in a given eikonal path. Furthermore, it will likely not yield a closed form expression to guide the user.

2. Symmetry

Fortunately, there are other means of evaluating the accuracy of geometric approximation over the various steps in a path. The simplest approach is to investigate the symmetry of the computed eikonal paths for the Gaussian case for which we can theoretically predict precise symmetry before and after the focal point. Figure 6(a) shows the positions that were selected as the origins of six rays across the cross-section of a Gaussian laser 10 µ before the focal plane. One ray was selected close to the optical axis at position $x = - 1\,\,\unicode{x00B5}{\rm m}$, one at the extreme end of the wavefields support at $x = - 7\,\,\unicode{x00B5}{\rm m}$, and one in the middle of these two positions at $x = - 4\,\,\unicode{x00B5}{\rm m}$. If the eikonal paths have been computed correctly we can expect the eikonal to pass through these same positions at a propagation distance of 10 µm after the focal plane. The parameters of the simulations are identical to those used in the first paper (see Fig. 5 in the first paper). The eikonal paths for these six rays are computed using five different values of ${h_z}$ ranging from 1 µm to 0.1 nm; these results are presented in Figs. 6(b)–6(f), and these results are overlaid as illustrated in Figs. 6(g) and 6(h). It is clear that the largest step size of ${h_z}$ gives obviously erroneous results, with the other four in closer agreement. Zoomed in regions around the arrival points for the three rays are shown in Figs. 6(i)–6(k). For the three rays, starting from the on-axis one outwards, the error in arrival position using ${h_z} = 10\;{\rm nm} $ is given by 5.5, 22, and 42 nm. For the next smallest step size of ${h_z} = 1\;{\rm nm} $ this reduces to 1, 5, and 28 nm. For the smallest case ${h_z} = 0.1\;{\rm nm} $ this further reduces to 0.02, 3, and 26 nm. For the outer ray, which follows the most non-linear path, it is more difficult to reduce the error by reducing ${h_z}$. Positional error may original from several sources. Although inappropriate choice of ${h_z}$ is a primary source of error, other sources exist including the error associated with the ASM, interpolation error, and error in calculating the phase gradient, and given the numerical precision of the method, machine error is also a possibility. Therefore, reducing the value of ${h_z}$ can only be expected to go so far to improve the results, and this appears to be the case for the outer ray. One must also consider the total computation time for a given value of ${h_z}$. For the case of ${h_z} = 10\;{\rm nm} $, the time taken to calculate the stack of complex amplitudes to generate these results was 40 s using a modern CPU and GPU, and this value will increase proportionally as ${h_z}$ is reduced. The value of ${h_z} = 10\;{\rm nm} $ was selected for the main results in the paper as the best compromise of accuracy and computation time.

Fig. 6. Illustration of symmetry of the rays which can be related to the accuracy of the method. For the Gaussian case we can predict symmetry before and after the focal point. Therefore, a ray originating at a particular $(x,y)$ position at 10 µm before the focal point should arrive at the same $(x,y)$ position at 10 µm after the focal point. (a) This is the amplitude of the Gaussian laser at 10 µm before the focal point and six ray positions are selected; (b)–(f) the six rays are traced using different values of ${h_z}$ ranging from 1 µm to 0.1 nm; (g) these different results are stacked, and the projection is shown in (h); (i)–(k) show the zoomed in regions around the expected arrival positions of ${-}{1}\;{\unicode{x00B5}{\rm m}}$, ${-}{4}\;{\unicode{x00B5}{\rm m}}$, and ${-}{7}\;{\unicode{x00B5}{\rm m}}$, respectively.

Download Full Size | PDF

3. Testing the Geometrical Approximation Using the Intensity Law

The symmetry method described above offers insights into the accuracy of the computed paths with respect to tracing the path governed by the phase gradient and by inferring the accuracy of the geometric approximation over the step size ${h_z}$. However, this approach is limited to the Gaussian case for which no lens aberration exists and, therefore, cannot be applied generally. Furthermore, it does not give insight into the validity of the eikonal approximation over a computed path. To gain a deeper appreciation, and one that can be computed for the general case, other approaches must be considered. Since the “actual” optical field is independently calculated at each $z$-plane in the grid using the ASM, the ASM solution for the amplitude at each end of a ray trace between two planes can be compared to the geometric optics solution; the geometric solution for the amplitude is determined by this time examining the imaginary part of Eq. (9), which reveals the following relationship:

(24)$$\nabla S \cdot \nabla a + \frac{a}{2}{\nabla ^2}S = 0.$$

This result implies the “intensity law” of geometric optics; see Ref. [1] Ch. 3.1.2, or what Ref. [14] refers to as conservation of energy. This implies that if an optical system causes the cross-sectional area of a tube or rays or “ray bundle” to decrease, as in focusing by a lens, then the solid angle the light subtends must increase correspondingly, and vice versa. This relationship ensures that the energy carried by the light beam remains constant, assuming no light is lost. The “intensity law” ensures that the light’s intensity (power per unit area) changes in such a way that the total energy in the light beam is conserved as it progresses through the optical system. This concept underpins many calculations and principles in geometric optics, particularly in the design and analysis of lenses and imaging systems.

Fig. 7. (a) The field amplitude 10 µm before the focal plane. Three regions are highlighted for inspection at different distances from the center; (b1)–(b3) show these three regions in more detail and highlight the three points taken in each case. A minimum separation of 100 nm is taken. We attempt to choose these points in an area of uniform amplitude. The various (c) sub-figures show the results of comparing the amplitude ${a_2}$ along the respective eikonal paths computed by the ASM, and the value given by Eq. (25), and (d) shows the corresponding root term that relates the areas of the triangles in the equation.

Download Full Size | PDF

The intensity law can also be interpreted as follows: Given a tube of rays, the integral of the intensity over any cross-sectional area along the path is constant. A simple algorithm for computing the amplitude along a ray is as follows [14]: Let ${a_1}$ be the amplitude at point ${p_1}$ in the source plane. Take three points very close to ${p_1}$ and let $d{\sigma _1}$ be the area of the triangle formed by these points. Trace the three rays to the target plane and let $d{\sigma _2}$ be the area of the triangle in the target plane. Then the amplitude ${a_2}$ at the intersection of the original ray at the target plane is given by

(25)$${a_2} \approx {a_1}\sqrt {\frac{{d{\sigma _1}}}{{d{\sigma _2}}}} .$$

Of course, energy is not in fact conserved within a tube of rays due to diffraction, and this is what will account for any discrepancy between the ASM computed field and the geometric optics solution provided by Eq. (25). The results of this method are shown for the Gaussian case in Fig. 7. Three triangular regions are selected for inspection at various distances from the center, similar to the approach taken for the symmetry test. Right-angled triangles are used with a side of 100 nm, which will reduce to $\approx 10\;{\rm nm} $ at the focus as shown in Figs. 7(a) and 7(b1)–7(b3). For the innermost ray, a comparison is shown in Fig. 7(c1) of the amplitude ${a_2}$ along the eikonal path computed by the ASM, and the value given by Eq. (25) and Fig. 7(d1) shows the corresponding root term that relates the areas of the triangles in the equation. Interestingly we see some discrepancy between the two results, with a sharp deviation occurring at ${-}3\,\,{\unicode{x00B5}{\rm m}}$ before the focus and another sharp change in area occurring at ${-}2\,\,{\unicode{x00B5}{\rm m}}$, apparently due to diffraction. The overall trend is symmetrical, and the geometric approximation appears to recover after the focus. In the following section, we will demonstrate how the field curvature can be measured along the same eikonal paths and that it is at these same points along the path that the curvature is found to suddenly spike, which is indicative of positions where the geometric approximation is straining. Corresponding results are shown in Figs. 7(c2) and 7(d2) and Figs. 7(c3) and 7(d3) for the middle and outer rays, and in both cases we can see good agreement with the geometric approximation.

Equivalent results are shown for the Laguerre-Gaussian case in Fig. 8. In this case, there is close agreement with the computed amplitude along the eikonal path and the intensity law. There is some discrepancy close to the focus where we can see the ratio of the triangle area undergo a rapid oscillation, presumably due to the effects of diffraction. Once again, in comparison with the calculation of curvature in the next section we will see that it is at these positions that the curvature is highest and the geometric approximation is under some strain.

Fig. 8. (a) The field amplitude 10 µm before the focal plane for the case of the Laguerre-Gaussian spatial mode. A single region is highlighted for inspection, which is shown in more detail in (b), where the three points are selected in an area of approximately uniform amplitude. (c) shows the result of comparing the amplitude ${a_2}$ along the eikonal path computed by the ASM, and the value given by Eq. (25), and (d) shows the corresponding root term that relates the areas of the triangles in the equation.

Download Full Size | PDF

It is possible to design a simple metric based on this approach that can be applied to later results in order to gauge the validity of the geometric approximation. We define the failure function $F$ as follows:

(26)$$F(z) = \frac{{\sqrt {{{\left({{a_2} - {a_1}\sqrt {\frac{{d{\sigma _1}}}{{d{\sigma _2}}}}} \right)}^2}}}}{{\max({a_2})}}.$$

The denominator $\max({a_2})$ is the maximum value of the computed amplitude over a trace in $z$ and is used to normalize the difference of the computed ${a_2}$ with respect to that predicted by geometric optics. The mean and maximum values of $F$ over the trace can be used to represent the accuracy of the geometric approximation over the trace, with a value close to zero being desirable.

4. Measuring the Curvature along a Computed Eikonal Path

Another method to inspect the geometric approximation along a trace is to estimate the field curvature at each point. We can conceive a simple method to do this by returning to Eq. (24) and rewriting it as follows:

(27)$$\left| {{\nabla ^2}S} \right| = \left| {\frac{2}{a}\nabla S \cdot \nabla a} \right|.$$

On the left of this equation is the Laplacian of the phase, which is related to the field curvature, and it is this quantity we set out to calculate at each point along the computed eikonal path. The term on the right of the equation $\nabla S \cdot \nabla a$ is the derivative of the amplitude along the direction of the eikonal. Therefore, we arrive at a straightforward method to calculate the Laplacian of the phase (and, therefore, to evaluate the accuracy of the geometric approximation) at every point along the eikonal. Interpolation can be applied to obtain the ASM amplitudes at each end of a ray: ${a_1}$ and ${a_2}$. The difference between these values $da = {a_2} - {a_1}$ can be interpreted as the derivative of the amplitude of the field in the direction of the ray (in units of normalized amplitude/ray path length). Taking the value for the gradient of the amplitude in the direction of the eikonal we can simply divide by the local amplitude, and the Laplacian of the phase can be calculated at every point along a computed path. This gives us an instantaneous measurement of the phase curvature and is broadly independent of step size and is perhaps the best indicator of the validity of the eikonal approximation.

This result is given in Fig. 9 for both the Gaussian and LG cases. It is clear that the phase curvature varies along each point in the ray paths, with larger values close to the focal plane. It is noticeable that the outer ray encounters the highest phase curvature as it nears the focal plane for both laser modes. The values of the Laplacian of the phase never exceed ${\lambda ^{- 1}}$, and although the point of curvature at which the eikonal approximation begins to break down is not strictly defined, it is commonly accepted that the approximation is valid when $| {{\nabla ^2}S} | \ll {\lambda ^{- 2}}$. In such case, we can see that this condition is valid for the three rays. This method of evaluating the eikonal approximation is simple to calculate as we compute an eikonal path using the method presented in this paper. We note that those points along the computed eikonal paths that were detected in the previous section, where the geometric approximation was found to strain, match with those points where we see high values in the Laplacian of the phase.

Fig. 9. The Laplacian of the phase can easily be computed along any computed eikonal trace. This value is related to the phase curvature and must be $\ll {\lambda ^{- 2}}$ to be in agreement with the eikonal approximation. The Laplacian of the phase is shown for the same three rays for both the Gaussian and Laguerre-Gaussian beams in Fig. 6.

Download Full Size | PDF

5. Crossing of Rays

An important feature of caustics is the crossing of rays. In those cases where we compute two rays to be crossing we can clearly identify these regions at worst as invalidating the eikonal approximation or at best being erroneously calculated due to the sources of error already discussed. We have encountered this phenomenon in generating the results for these papers using longer propagation distances but not for values of ${h_z} \lt 10\;{\rm nm} $ when using relatively smooth fields. Crossing can be expected if the ray origins are selected to be close to each other, and the paths are close to regions of high phase curvature, and/or large values of ${h_z}$ are employed. If an exact crossing point happens to occur in a given trace, the two paths must merge into one since the phase gradient at that point must provide a single unique direction for the two intersecting points.

7. CONCLUSION

In this paper, we have proposed and demonstrated a new algorithm for tracing the flux lines through a focusing optical wavefield. This method builds on two algorithms that were developed in the first paper, which can be used to generate a 3D grid of samples that represent the discretized wavefield over three dimensions in the focal region of the lens. The difference between those two algorithms relates to the type of lens that is being simulated; the first algorithm simulates the effect of a thin lens and involves numerical propagation from the plane of the thin lens to the region around the focal plane; the second algorithm deals with the ideal lens, which can be modelled using an optical Fourier transform. These two algorithms can include the laser spatial mode, as well as the lens aberrations, which is explored in the third paper.

Regardless of which algorithm is employed, the resulting 3D grid of samples serves as the starting point for the flux tracing method proposed in this paper. Section 2.A provides an overview of the eikonal function, which can be used to relate the complex electromagnetic wavefront and the classical interpretation of geometrical ray optics. The eikonal function can be used to define the flux lines, which trace the movement of the power over three dimensions as the wavefront propagates and diffracts. In the same section we also reviewed the concept of the derivative of the phase and explain how this is related to the concept of a ray in the classical geometrical optics sense. We show that the phase derivative can be used to define the angle at which the flux lines are moving at any position. We proceed to develop an algorithm that traces the flux through a sequence of different propagation distances separated by a value ${h_z}$; more specifically, the phase derivative of the wavefield at the position of the ray is recalculated at each plane to direct the ray linearly and trace between successive planes, thereby finding the new ray positions at each plane over a sequence of steps. Using the concept of the eikonal, we prove that the accuracy of this algorithm is high so long as the step size ${h_z}$ is small and the field is smooth.

The independent computation of the 3D optical field is an important feature that gives the scheme a strong chance of working, as long as the $z$-planes are spaced sufficiently close together. The geometric optics solution is valid at least in a neighborhood of each plane. There is no guarantee that even very high-order WKB expansions can be propagated indefinitely without rays intersecting. Furthermore, these equations are nonlinear, and integrating them stably is non-trivial. The approach presented in this paper avoids these issues.

Having validated the method, the algorithm is applied to various cases of focusing different laser spatial modes. The two different algorithms developed in the first paper are applied to compute the 3D wavefield, relating to the thin lens approximation and the ideal lens, and we find the results are very similar for both. The most striking results relate to the Laguerre-Gaussian laser spatial modes for which the rays appear to spiral rapidly with a spin rate dependent on the order of the LG mode and the ray position relative to the optical axis. Those rays closest to the center spin most rapidly, and there appears to be a linear dependence with the Laguerre-Gaussian order, which is consistent with the existing interpretation of orbital angular momentum [13]. It is very interesting to note the movement of the rays close to the focus “doughnut,” with many rays undertaking a rapid change of direction, and some forming a bubble like shape around the focus as they move from the inner to the outer doughnut patterns.

Other than the study of optical aberration and its potential usage in lens design, we envisage the application of non-linear ray tracing as both an educational and a research tool. Educationally, it can provide insights into the movement of light energy through an optical system at high computational speed. As well as tracing the eikonal though the light field it can also produce information on phase curvature at all points along this trace providing an indication on the validity of geometrical optics and the eikonal approximation in these regions. Wavefront sensing methods are commercially available to characterize optical surfaces and are used for a variety of applications. One could envisage the use of our non-linear ray tracing using a recorded wavefront as a starting point to further investigate these applications. Another potential application is the study of Laguerre-Gaussian beams in the context of orbital angular momentum and the transfer of angular momentum to optically trapped particles. The results presented in this paper highlight the very different spin rates of the flux lines at different positions in the focal region, with the fastest spin closest to the center. It is conceivable that further study might reveal a model to relate the eikonals spin to the transfer of angular momentum to a trapped particle. The code to implement the method proposed in this paper is provided in Code 1, Ref. [15]. In the third paper we will highlight the potential of the flux tracing method presented here to investigate the nonlinear rays that pass through the focal region of a lens that contains aberration. This will require only a simple extension of the two numerical propagation algorithms and provides extraordinary insights into the effects of various complex aberrations.

Funding

Irish Research eLibrary; Science Foundation Ireland (15/CDA/3667, 16/RI/3399, 19/FFP/7025).

Acknowledgment

The authors would like to acknowledge the support of Science Foundation Ireland. The authors would also like to thank an anonymous reviewer for their patience and insights and for lending us their wonderful expertise. This paper is much better because of their contribution. Open access funding provided by Irish Research eLibrary. This publication has emanated from research conducted with the financial support of Science Foundation Ireland under Grant numbers 15/CDA/3667 and 19/FFP/7025.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper have been produced using the code provided in Code 1, Ref. [15].

Supplemental document

See Supplement 1 for supporting content.

REFERENCES

1. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Elsevier, 2013).

2. J. W. Goodman, Introduction to Fourier Optics (Roberts and Company, 2005).

3. R. N. Bracewell and R. N. Bracewell, The Fourier Transform and Its Applications (McGraw-Hill, 1986), Vol. 31999.

4. W. T. Welford, Aberrations of Optical Systems (Routledge, 2017).

5. J. C. Dainty, Laser Speckle and Related Phenomena (Springer, 2013), Vol. 9.

6. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company, 2007).

7. O. N. Stavroudis and R. C. Fronczek, “Caustic surfaces and the structure of the geometrical image,” J. Opt. Soc. Am. 66, 795–800 (1976). [CrossRef]

8. O. Stavroudis, The Optics of Rays, Wavefronts, and Caustics (Elsevier, 2012), Vol. 38.

9. C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory (Springer, 1999).

10. M. Berry and P. Shukla, “Semiclassical superoscillations: interference, evanescence, post-WKB,” New J. Phys. 23, 113014 (2021). [CrossRef]

11. L. Ahrenberg, A. J. Page, B. M. Hennelly, et al., “Using commodity graphics hardware for real-time digital hologram view-reconstruction,” J. Disp. Technol. 5, 111–119 (2009). [CrossRef]

12. A. Ashkin, “Forces of a single-beam gradient laser trap on a dielectric sphere in the ray optics regime,” Biophys. J. 61, 569–582 (1992). [CrossRef]

13. A. M. Yao and M. J. Padgett, “Orbital angular momentum: origins, behavior and applications,” Adv. Opt. Photonics 3, 161–204 (2011). [CrossRef]

14. B. D. Seckler and J. B. Keller, “Geometrical theory of diffraction in inhomogeneous media,” J. Acoust. Soc. Am. 31, 192–205 (1959). [CrossRef]

15. Q. Yu and B. M. Hennelly, “Code for ‘Nonlinear ray tracing in focused fields, part 2. Tracing the flux: tutorial’,” figshare, 2024, https://doi.org/10.6084/m9.figshare.25066103.

Name	Description
Code 1	This code can be used to implement the nonlinear ray tracing method outlined in the main paper.
Supplement 1	This document show the results for various LG modes focused using the ideal lens algorithm.
Visualization 1	This video shows the nonlinear rays spinning around the focal point of a focused TEM01 laser.
Visualization 2	This video shows the nonlinear rays spinning around the focal point of a focused TEM01 laser.

Nonlinear ray tracing in focused fields, part 2. Tracing the flux: tutorial

Abstract

1. INTRODUCTION

2. FLUX

A. Physical Interpretation of the Angular Spectrum

B. Concept of a Ray as a Flux Trace

C. Rays and Local Spatial Frequency

3. LIMITATIONS OF GEOMETRIC OPTICS AND THE EIKONAL APPROXIMATION

A. Relationship between the Eikonal and Ray Optics

B. Failure of the Eikonal Approximation

C. Identifying When the Eikonal Approximation Breaks Down

D. Beyond the Eikonal Approximation

4. TRACING THE FLUX

A. Method

B. Flux Tracing Algorithm

5. RESULTS FOR DIFFERENT LASER SPATIAL MODES

6. LIMITATIONS AND ACCURACY

A. Accuracy in Tracing the Phase Gradient

B. Estimating the Accuracy of the Geometric Approximation over a Computed Eikonal Path

1. WKB

2. Symmetry

3. Testing the Geometrical Approximation Using the Intensity Law

4. Measuring the Curvature along a Computed Eikonal Path

5. Crossing of Rays

7. CONCLUSION

Funding

Acknowledgment

Disclosures

Data availability

Supplemental document

REFERENCES

Supplementary Material (4)

Data availability

Cited By

Figures (9)

Equations (27)

Journal of the Optical Society of America A