1. INTRODUCTION

In the first two parts of this three-paper series, a method is defined to trace the lines of flux in the focal regions of a lens. At the heart of this method exist two numerical propagation algorithms, defined in the first paper, which serve to calculate a three-dimensional grid of samples of the complex scalar wavefield around the focal point of the lens. This grid of samples can then be used as a scaffold through which the flux lines can be traced with high accuracy, which is detailed in the second paper. The two algorithms relate to the case of the thin lens approximation, and the case of the ideal lens, which can be modeled using a Fourier transform. In this paper, we extend both algorithms to include lens aberration, modeled using an additive phase term based on the well-known Zernike polynomials. The manner in which this additive phase is applied is fundamentally different for both algorithms. For the thin lens case the additive phase is applied to the transmission function of the thin lens and, therefore, it describes a modulation of the lens thickness function. For the ideal lens case the additive phase is applied to the back focal plane of the lens. This approach is consistent with the traditional method for modeling aberration in an ideal two-lens imaging system [1,2]. We proceed to examine the lines of flux in the focus of lenses containing a variety of different aberrations, and to relate these results to the classical interpretation of aberration theory, and in doing so we clearly demonstrate the potential of nonlinear ray tracing to gain insight into complex optical phenomena.

Aberration in a lens system is the failure of that system to produce an ideal image; this can be interpreted either in terms of geometrical ray optics, whereby all rays passing through a point in the input plane should also pass through a single point in the image plane, or in terms of wavefronts, whereby the spherical wavefronts emanating from a point in the input plane should be portions of spheres centered at a point in the image plane. Characterization of aberrations can employ finite ray tracing through the optical system (see Chapter 4 Ref. [3]) and wavefront aberration theory (see Chapters 7-8 Ref. [3] and Ref. [4]), which is derived using Hamiltonian optics and is often interpreted in terms of the Seidel aberrations. Both theories can be related using the characteristic function of Hamilton [3]. In this paper, we attempt to use nonlinear ray tracing as an alternative to classical wavefront aberration theory. In the next section we relate the classical aberration theory to nonlinear ray tracing and explain how it has the potential to provide predictions and insights that are more accurate and more consistent with the tenets of Hamiltonian optics. This is followed in Section 3 by a description of how aberrations can be included in the proposed flux tracing method and in Section 4 with several results showing the nonlinear rays for various aberrations.

2. CLASSICAL ABERRATION THEORY

A. Ray Aberration Theory

We begin with a brief discussion on the interpretation of lens aberration using the principles of geometrical optics and classical wavefront aberration theory. Some low order aberrations are most easily understood purely in terms of geometrical optics whereby we can say that the geometrical rays have deviated from the ideal path due to the presence of a lens aberration [3]; the ideal path refers to the path we might expect from an ideal lens system that focuses parallel rays to a single point in three dimensions at the focal plane of the lens. We proceed with a brief description of some of the most commonly encountered low order aberrations that can be approximately described in terms of geometrical optics.

Spherical aberration results in the rays that are converging from the lens coming to focus at different planes along the optical axis as illustrated in Fig. 1(a); for this case the “focal length” for a given ray is related to its position of origin on the lens aperture. Rays that originate farther from the center of the lens will come to focus at a plane some distance farther than the correct focal plane. Spherical aberration can significantly affect resolution and clarity across an image. Spherical aberration is an inherent property of perfectly spherical lenses, such as those described by the thin lens phase transmission functions ${t_1}$ and ${t_2}$ in the first paper. Defocus is a similar aberration in the sense that it also relates to a deviation from the correct focal plane. In this case, however, all of the rays are understood to focus at a single plane, but not the correct image plane. Defocus is illustrated in Fig. 1(b).

 

Fig. 1. Four low order aberrations that can be interpreted in terms of geometrical optics: (a) spherical aberration for which rays focus at different distances depending on their position in the lens aperture; (b) defocus aberration, which causes the rays to focus at an incorrect focal plane; (c) comatic aberrations for which off-axis rays focus at different positions in the focal plane depending on their position in the lens aperture; (d) astigmatism aberration, which results in two or more focal points, most commonly for rays in the vertical and horizontal axes.

Download Full Size | PDF

Comatic aberration or “coma” results from imperfect lens designs that cause off-axis point sources to appear distorted in the image, appearing to trail off like a comet tail. More specifically, coma is defined as a variation in the magnification over the lens pupil. Obliquely parallel rays will focus at different points in the lens focal plane depending on the position of their origin in the lens aperture. The effect of comatic aberration is illustrated in Fig. 1(c). Astigmatism relates to an irregular lens curvature that is elliptical rather than spherical, which results in light coming to focus at two or more different distinct points, with obvious deleterious effects on the imaging properties of the lens. There are various classes of astigmatism depending on the shape of the irregularity. Vertical astigmatism, illustrated in Fig. 1(d), relates to different focal lengths in the $x$ and $y$ directions resulting in rays located along these two axes coming to focus at two different distances. Oblique astigmatism relates to a lens curvature whereby the slow and fast lens meridians are not at 90° or 180° as for the case of vertical astigmatism.

B. Wavefront Aberration Theory

We note that a simple interpretation purely in terms of geometrical optics is limited to only some of the most commonly encountered low order aberrations; many more forms of low order optical aberrations exist [1] that cannot easily be interpreted in terms of geometrical optics, examples of which include trefoil and quadrafoil. Trefoil can be thought of as a three-axis astigmatism but this does result in three distinct focal points; however, triangular focal spots can be seen in the images of lenses with this distortion. To model the effect of these more complicated aberrations, wavefront aberration theory has been developed (see Chapters 7-8 Ref. [3] and Ref. [4] for an excellent review of the topic).

The classical wavefront aberration theory is essentially based on Hamiltonian optics. Hamiltonian optics uses the principles of classical mechanics to describe light propagation. It employs the Hamiltonian function, a scalar function of coordinates and momenta, to describe how light propagates through an optical system. Hamilton’s equations of motion describe how light rays change as they move through an optical system. In the Hamiltonian framework, the optical path length is central. The concept of optical path length is identical to the eikonal function derived through the wave equation in the second paper. Indeed, the eikonal function can also be derived using Hamiltonian optics, independently of the Helmholtz equation [1]. Classical wavefront theory utilizes the Hamiltonian framework to provide an understanding of Seidel aberrations whereby the wavefront can be described by a three-dimensional surface (wavefront of constant phase) at some position in the optical system, and the lines normal to this surface determine the position and direction of geometrical rays propagating from that surface, which can then be traced over a long distance to the image plane. Distortions in the wavefront due to aberration, will result in deviation in the arrival position of the rays in the image plane.

Classical wavefront theory has been demonstrated time and time again to provide valuable insights for a wide range of optical design problems. However, it must be noted that at its core it is based on a simplified form of Hamiltonian optics that is not entirely accurate, whereby geometrical rays are traced linearly over long distances after they exit an optical system. In general, an aberrated wavefront (neither linear nor spherical in nature) can never be described by a set of linear rays propagating over an extended distance. This is a simplification that, while useful for many applications, does not account for the continually evolving nature of the wavefront during propagation, especially in highly diffractive cases. Classical wavefront aberration theory will often calculate the eikonal direction once and never update it despite a long propagation distance, and although this can provide very meaningful results for many examples and indeed aid in the design if compensating for optical elements, it does not fully capture the intricacies of light propagation.

The eikonal equation provides a more comprehensive description of light’s behavior and a more complete representation of Hamiltonian optics; as the wavefield propagates and diffracts, the eikonal will evolve and the “ray” direction will change continuously. Therefore, we believe that the method presented in this paper, which although is much more computationally expensive, can provide a deeper understanding and a richer insight into the effects of complex optical aberrations as well as other forms of complex wavefronts such as a spiral phase. One example is that the simplified Hamiltonian framework often used in wavefront aberration theory could not be used to describe the spiraling rays associated with a Laguerre–Gaussian wavefront.

3. EXTENDING THE THIN LENS ALGORITHM AND 2 TO INCLUDE LENS ABERRATIONS

Here we describe how to augment the Thin Lens Algorithm and 2 developed in the first paper to include the effect of lens aberration in the numerical model such that the lens can be distorted, and the resulting focal pattern can be examined in three dimensions. Rather than use the Seidel aberrations we employ the Zernike polynomials to account for lens aberration for both the thin lens case and the ideal lens case, albeit in different manners. We note that the Seidel aberrations could easily be substituted if desired. The Zernike polynomials are a set of orthogonal two-dimensional polynomials defined over the unit circle. These are particularly useful in describing the wavefront distortions that can exist for optical systems with circular pupils due to the property that they form a complete and orthogonal basis set over the unit circle. Therefore, any given function defined over this unit circle, $W(x,y)$, can be represented as a weighted (infinite) sum of these polynomials as follows:

Table 1. Description of the First 15 Zernike Polynomials, ${Z_n}$, in Cartesian Form^a

View Table

A. Extending the Ideal Lens Algorithm

When wavefront error exists in an imaging system that contains aberration, a classical way to account for it is to describe its effect in terms of a phase delay distortion in the pupil function of the ideal lens system. An ideal diffraction limited optical system can be described as having a frequency response or pupil function defined as the circle function, $P({{{x^\prime_0}}/r\lambda f,{{y^\prime_0}}/r\lambda f})$ as defined in Eqs. (20) and (21) in the first paper where $P({{{x^\prime_0}},{{y^\prime_0}}})$ is the unit circle over $({{{x^\prime_0}},{{y^\prime_0}}})$. An aberrated optical system can be described in terms of the same pupil function with a spatially varying additive phase delay $kW({{{x^\prime_0}},{{y^\prime_0}}})$. Therefore, the pupil function, $P^\prime $, of an aberrated system is defined as follows:

B. Extending the Thin Lens Algorithm

The method above for the case of an ideal lens system, in which aberrations are described as a phase delay in the Fourier plane, or pupil function of the optical system is well established [2]. However, the thin lens case requires a different approach. Here we account for aberrations as direct distortions of the phase delay associated with the transmission function of the thin lens, $t(x,y)$. Therefore, in the description of the Thin Lens Algorithm, which simulates a thin lens system, Eq. (10) in the first paper is replaced with the following definition of the transmission function:

Thus, the phase delay imparted by the lens $\phi$ is added to the phase delay associated with each Zernike polynomial, ${w_n}{Z_n}$. The first step in the Thin Lens Algorithm is replaced with the following.

We note that aberrations have previously been investigated for the thin lens model; in Ref. [6] the Rayleigh–Sommerfeld diffraction formula was used to theoretically investigate the impulse response from a thin planar lens, which is similar in concept to the numerical approach applied here. In that paper, the authors showed that a thin quadratic lens produced an impulse response given by the integral of the pupil function multiplied by a phase function containing the Seidel aberrations. For our RS lens, we can expect no such aberrations (at least for on-axis imaging), since the thickness function of the lens has been modeled on the kernel of the Rayleigh–Sommerfeld integral.

 

Fig. 2. Intensity and phase distributions of focused spot for RS lens with Gaussian laser and with spherical and defocus aberrations, calculated using the Thin Lens Algorithm. See text for more details.

Download Full Size | PDF

C. Results for the Thin Lens Algorithm

We proceed to examine the effect of aberration on the RS thin lens for the Zernike polynomials using the Thin Lens Algorithm. Equivalent results are shown in Supplement 1 for the ideal lens using the Ideal Lens Algorithm with remarkable similarities. For the RS lens with a Gaussian profile that was simulated in the first paper (for which $f = 3.3\,\,{\rm{mm}}$, $R = 2.5\,\,{\rm{mm}}$, $\sigma = R/\sqrt 2$), the Thin Lens Algorithm was implemented once again using an identical set of 3D sampling conditions as described in the first paper. This time, however, the aberration term was included, i.e., ${u_0}({x_0},{y_0}) = {t_{\text{RS}}}({x_0},{y_0})\def\LDeqbreak{}{A_{{{\rm TEM}_{00}}}}({x_0},{y_0})Z({x_0}/R,{y_0}/R)$ and the weights ${w_n}$ were varied in order to examine the effect of specific cases of aberration on the 3D focal volume of light.

1. Spherical Aberration and Defocus

The first case that was examined is that of spherical aberration for which the weight ${w_{11}} = 2 \times {10^{- 7}}$ and ${w_n} = 0\;\forall n \ne 11$. The intensity and phase distributions of the resulting field in the $xz$-plane are shown in Figs. 2(a1) and 2(a2) over a size of $20\; \unicode{x00B5}{\rm m} \times 40\; \unicode{x00B5}{\rm m}$. The effect of diffraction caused by the aberration on both the intensity and phase distributions is evident when compared with Fig. 5 in the first paper. The correct focal plane of the lens (i.e., $z = f$) is indicated in the figure using the symbol “$f$.” In addition to the diffraction there is also a clear shift in the axial position of the focus. The intensity and phase in the $xy$-plane over an area of $10\; \unicode{x00B5}{\rm m} \times 10\; \unicode{x00B5}{\rm m}$ are shown at the brightest point along the $z$-axis in Figs. 2(a3) and 2(a4). The intensity and phase of this focused spot are similar in form to the unaberrated case shown in the first paper. It should be noted, however, that increasing the value of ${w_{11}}$ will result in further distortion of the spot as well as higher levels of diffraction and increased displacement of the focal point relative to the focal plane of the lens.

 

Fig. 3. Intensity and phase distributions of focused spot for the RS lens with Gaussian laser and primary astigmatism, calculated using the Thin Lens Algorithm. See text for more details.

Download Full Size | PDF

Also shown in the same figure is the case of defocus aberration for which the weight ${w_4} = 2 \times {10^{- 6}}$ and ${w_n} = 0\;\forall n \ne 4$. The intensity and phase distributions of the resulting field in the $xz$-plane are shown in Figs. 2(b1) and 2(b2) over the same size as for the previous case. The effect of diffraction caused by the defocus aberration on both the intensity and phase distributions is remarkably similar in form to the previous results for spherical aberration. However, there is a significant difference in the position of the focused spot along the $z$-axis for both cases. Whereas the focused spot with spherical aberration occurs a few micrometers after the correct focal plane, $z = f$, it appears a few tens of micrometers before this plane for the case of defocus. The intensity and phase in the $xy$-plane over an area of $10\; \unicode{x00B5}{\rm m} \times 10\; \unicode{x00B5}{\rm m}$ are shown at the brightest point along the $z$-axis in Figs. 2(b3) and 2(b4). Once again this is similar to the equivalent result for the unaberrated case but larger values of ${w_4}$ will result in further distortion and displacement.

 

Fig. 4. Intensity and phase distributions of focused spot for the RS lens defined with Gaussian laser and vertical coma, calculated using the Thin Lens Algorithm. See text for further details. We note that the phase values in regions of very low or negligible amplitude have inaccurate or indeterminate value and have no physical meaning.

Download Full Size | PDF

2. Astigmatism

The third aberration to be investigated is primary astigmatism, which is commonly perceived as producing two focal points at different distances resulting from different lens curvatures along orthogonal axes. It should be noted that the polynomials ${Z_5}$ (vertical astigmatism) and ${Z_6}$ (oblique astigmatism) are similar, in that the latter is a $\pi /4$ rotation of the former. In Figs. 3(a1)–3(a10) the results are shown for the case of vertical astigmatism with ${w_5} = 4 \times {10^{- 7}}$ and ${w_n} = 0\;\forall n \ne 5$. The intensity of the slice through the center of the $xz$-plane and the $yz$-plane is shown in Figs. 3(a1) and 3(a3), respectively, with the corresponding phase images given in Figs. 3(a2) and 3(a4). The intensity and phase patterns of the $yz$ case appear to be identical to the patterns for the $xz$ case with an inversion of the $z$-axis. It is clear that these two intensity images differ in terms of the distance at which the light is converging to a focus. For the $xz$ case, the light is seen to converge at a distance of $z = f - 3.8\; \unicode{x00B5}{\rm m}$, and the intensity and phase distributions in the $xy$-plane at this distance are shown in Figs. 3(a5) and 3(a6) with the profile of the intensity along the $x$-axis shown in the small inset image; the light is obviously focused along the $x$-axis but has a wider support along the $y$-axis. The opposite is true in the $xy$-plane at a distance of $z = f + 3.8\; \unicode{x00B5}{\rm m}$, as demonstrated by the intensity and phase distributions shown in Figs. 3(a9) and 3(a10), which are the same pattern as in Figs. 3(a5) and (a6) rotated by $\pi /2$. At the correct focal distance of $z = f$, half way between the two “focused spots,” the intensity pattern shown in Fig. 3(a7) takes the shape of a diamond and the phase pattern shown in Fig. 3(a8) has the interesting property of each profile through the center being equal to the inverse of the profile at a rotation of $\pi /2$. A similar set of results is shown in Figs. 3(b1)–3(b10) for the case of oblique astigmatism with ${w_6} = 4 \times {10^{- 7}}$ and ${w_n} = 0\;\forall n \ne 6$. In this case, the intensity and phase distributions for the center slices through the $xz$- and $yz$-planes are identical to those shown in Figs. 3(b1)–3(b4). The reason for this is that the two slices are passing though the elongated patterns at an angle of $\pi /4$. This becomes more apparent by observing the $xy$-plane intensity and phase distributions at $z = f \pm 3.8\; \unicode{x00B5}{\rm m}$ shown in Figs. 3(b5), 3(b6), 3(b9), and 3(b10), which are clearly $\pi /4$ rotations of the equivalent images for the vertical astigmatism case. The intensity and phase distributions at the correct focal plane shown in Figs. 3(b7) and 3(b8) are also a $\pi /4$ rotation of the vertical astigmatism case.

 

Fig. 5. Intensity and phase distributions of focused spot for the case of a Gaussian beam incident on the RS lens with vertical trefoil aberration shown in (a1)–(a6) and vertical quadrafoil aberration shown in (b1)–(b4), calculated using the Thin Lens Algorithm; see text for further details. We note that the phase values in regions of very low or negligible amplitude have inaccurate or indeterminate value and have no physical meaning.

Download Full Size | PDF

3. Comatic Aberration

Comatic aberration was also applied to the RS lens using identical simulation parameters. In Fig. 4 the intensity and phase patterns are shown for the case of vertical coma with ${w_8} = 3 \times {10^{- 7}}$ and ${w_n} = 0\;\forall n \ne 8$. The intensity of the slice through the center of the $xz$-plane and the $yz$-plane is shown in Figs. 4(a1) and 4(a3), respectively, with the corresponding phase images given in Figs. 4(a2) and 4(a4). The curvature of the focused spot along the $z$-axis is evident both in terms of intensity and phase. The intensity and phase distribution in the $xy$-plane at the correct focal plane is shown in Figs. 4(a7) and 4(a8) with the profile of the intensity along the $x$-axis shown in the small inset image; the light is clearly displaced along the $x$-axis with an the characteristic elongated distortion of the spot that is associated with coma. The intensity and phase images in the $xy$ plane are also shown at a distance of $z = f - 5\; \unicode{x00B5}{\rm m}$ in Figs. 4(a5) and (a6). Here, it can be seen that the spot is further displaced from the center of the $x$-axis due to the curvature of the focused light over the $z$-axis.

4. Trefoil and Quadrafoil

The final aberrations that were investigated with the RS lens are vertical trefoil (${w_9} = 3 \times {10^{- 7}}$) and vertical quadrafoil (${w_{14}} = 1 \times {10^{- 6}}$) with resulting intensity and phase patterns shown in Fig. 5. For the case of trefoil, the intensity and phase distributions in the $xz$-plane are shown in Figs. 5(a1) and 5(a2). These patterns are symmetrical in both the $x$-axis and the $z$-axis. However, the intensity and phase in the $yz$-slice, shown in Figs. 5(a1) and 5(a2), are not symmetrical in the $y$-axis. This results from the triangular shape of the distribution in the $xy$-plane, which is shown in Figs. 5(a5) and 5(a6) and the fact that the Zernike polynomial ${Z_9}$ is periodic over rotations of $2pi/3$. We note that the results presented here are for vertical trefoil; see Table 1. Similar results are obtained for oblique trefoil, which is effectively a $\pi /6$ rotation of vertical trefoil. Therefore, for the case of (${w_9} = 3 \times {10^{- 7}}$) we would observe a rotation of the $xy$-distribution by $\pi /2$ or equivalently an inversion of the triangular pattern along the $y$-axis. Also, the $xz$- and $yz$-distributions would transpose and the latter would be inverted along the $x$-axis. The results for quadrafoil are shown in Figs. 5(b1)–(b4). The associated Zernike polynomial, ${Z_{14}}$, is periodic over rotations of $\pi /2$ while the intensity and phase distributions of the focused spot in the $xy$-plane ($z = f$), shown in Figs. 5(b3) and 5(b4), are periodic over rotations of $\pi /4$; in the intensity pattern eight small spots can be seen positioned in an annular pattern around the focused spot. In Figs. 5(b1) and 5(b2) the $xz$-distributions of the intensity and phase are shown. The presence of two of the aforementioned spots can be seen as side lobes around the center spot. We note that the same pattern is observed at angles of $n\pi /4$ around the $xy$-plane, where $n$ is an integer. Not shown here is the opposing pattern at angles of $n\pi /4 + \pi /8$, which contains only the center spot without side lobes. We note that the results presented here are for vertical quadrafoil; see Table 1. Similar results are obtained for oblique quadrafoil, which is effectively a $\pi /8$ rotation of vertical quadrafoil. For that case the $xz$ projection would not contain side lobes.

A similar set of results for these various aberrations is shown in Supplement 1 for the case of the ideal lens with a Gaussian beam, simulated using the Thin Lens Algorithm and using an identical set of parameters as detailed in Sections 3.B and 4 in the first paper. The similarities between the intensity and phase patterns that are produced for the various aberrations are strikingly similar to those for the RS lens calculated using the Thin Lens Algorithm, and using a fundamentally different approach to account for the aberrations, i.e., the latter accounts for the distortion as a variation in the thickness of the transmission function of the thin lens, while the former accounts for these in the pupil function or Fourier domain.

 

Fig. 6. Nonlinear ray tracing applied to different cases of (a1)–(a3) spherical and (b1)–(b3) defocus aberration on the Rayleigh–Sommerfeld lens with numerical aperture of 0.6 and ${{\rm{TEM}}_{00}}$ illumination computed using the Thin Lens Algorithm. This simulation is based on the identical method and parameters described in the first paper. Three different weights, $w$, of spherical and defocus aberration are applied, shown in the figures. Note that the center figures in both rows (a3) and (b3) correspond to the specific results given in Fig. 6 in the first paper.

Download Full Size | PDF

4. TRACING THE FLUX WITH LENS ABERRATION

The results of nonlinear ray tracing applied to the various aberrations are presented in this section, all for the case of the Raleigh–Sommerfeld lens with a numerical aperture of 0.6 with ${\rm TEM}_{00}$ illumination and all computed using the Thin Lens Algorithm in the flux-tracing method described in the second paper using identical parameters. Note that in most cases the same Zernike weights $w$ are applied as given in the previous section, except where stated. The ray origins are selected at positions that best highlight the effects of the aberration on a case by case basis. We note that equivalent results relating to all types of aberration are presented in Supplement 1 for the ideal lens computed using the Ideal Lens Algorithm with remarkable similarity. In all cases a value of ${h_z} = 10\;{\rm{nm}}$ was selected as offering the best compromise of speed and accuracy.

A. Spherical Aberration and Defocus

The results of nonlinear ray tracing applied to spherical aberration are presented in Figs. 6(a1)–6(a3). Note that the third case [Fig. 6(a3)] corresponds to the result given in Fig. 2 above. In each case, the log of the intensity is shown in the background on which 12 rays are overlaid with origins that are distributed along the $x$-axis $20\; \unicode{x00B5}{\rm m}$ before the focal plane. The focal plane, which we define as the plane with the maximum intensity, is seen to shift as a function of $w$, with a significant shift in offset occurring for the case of $w = 5 \times {10^{- 8}}$. For the case with lowest weight shown in Fig. 6(a1) only a slight offset in focus is observed and the results are similar to those for the aberration free case presented in Fig. 4 in the second paper with the innermost rays following paths that avoid regions of low intensity. The outermost rays appear to come to focus slightly before the focal plane. This effect is a little more pronounced for the second case in Fig. 6(a2) for which the weighing has been doubled $w = 1 \times {10^{- 7}}$. The third case, Fig. 6(a3), for which the weight is doubled again to $w = 2 \times {10^{- 7}}$, best exemplifies the effects of spherical aberrations. Here it can be seen that the outer rays focus before the inner rays, and the distance at which they focus depends on the position of the ray along the $x$-axis. This is consistent with the conventional interpretation of spherical aberration in terms of geometrical rays that was discussed earlier.

A similar set of results is shown in Fig. 6(b1)–6(b3) for the case of defocus aberration. For the lower weight of $w = 5 \times {10^{- 7}}$, the result looks similar to case of low spherical aberration in Fig. 6(a1) and the effects of the aberration appear to be minimal. However, the shift in the focal plane is significantly greater, which is consistent with the interpretation of defocus aberration in terms of geometrical optics described earlier. For the second case in Fig. 6(b2) with the value of $w$ doubled, this shift in focus is also doubled. We also see in this case that the outer most rays appear to converge at a plane slightly before the focal plane. This result is comparable to the second case of spherical aberration shown in Fig. 6(a2). The result of a further doubling of $w$ is shown in Fig. 6(b3). As expected, the shift in focus is doubled again with respect to the previous case. However, in this case we can see that the outer rays follow a trend that is similar to the case of high spherical aberration in Fig. 6(a3). We can conclude that high levels of defocus and spherical aberration exhibit similar effects on the flux lines, with defocus resulting in greater shifting of the focal plane. It is also notable that the shift in focus is linearly dependent on $w$ for the case of defocus aberration but this relationship appears to be nonlinear for the case of spherical aberration, although the amount of defocus is much less for the latter.

 

Fig. 7. Nonlinear ray tracing for vertical and oblique astigmatisms for the Rayleigh–Sommerfeld lens; see text for details.

Download Full Size | PDF

B. Astigmatism

The results of nonlinear ray tracing applied to vertical astigmatism aberration are presented in Figs. 7(a1)–7(a12). In the left most column a 3D plot is shown that traces various sets of 10 or 12 rays from $z = - 20\; \unicode{x00B5}{\rm m}$ before the focal plane to $z = 20\; \unicode{x00B5}{\rm m}$ after, and also shows two inset (log of) intensity images that correspond to the intensity at the plane $z = - 20\; \unicode{x00B5}{\rm m}$ before the focal plane, and at the focal plane. The larger inset image has dimensions of $20\; \unicode{x00B5}{\rm m} \times 20\; \unicode{x00B5}{\rm m}$ and the smaller inset image has dimensions of $10\; \unicode{x00B5}{\rm m} \times 10\; \unicode{x00B5}{\rm m}$. The different rows correspond to a different choice in the ray origin positions in the plane at $z = - 20\; \unicode{x00B5}{\rm m}$. The first row of images in Figs. 7(a1)–(a4) illustrate a uniform distribution of 10 ray origins along the $x$-axis. It is clear from the various projections, as well as from the two intensity images, that the rays do not deviate from the plane $y = 0$. It is also clear that the rays come to focus a number of micrometers before the expected focal plane. The second row of images in Figs. 7(a5)–(a8) relate to a uniform distribution of 10 ray origins along the $y$-axis. In this case, the intensity images and the various projections indicate that the rays do not deviate from the plane $x = 0$ and come to focus a number of micrometers after the expected focal plane. These effects are consistent with the conventional interpretation of vertical astigmatism in terms of geometrical rays that was discussed earlier. In the third row of images in Figs. 7(a9)–(a12) the effect of vertical astigmatism on nonaxial ray origins is shown. In this case, 12 ray origins are uniformly distributed around a circle in the plane at $z = - 20\; \unicode{x00B5}{\rm m}$. It can be seen from the intensity images that the ray positions form an ellipse in the focal plane. From the various projections, it can be seen that those rays that are located on the $x$- and $y$-axes remain in the $x$- and $y$-planes, while rays that are not located on these axes appear to twist in three dimensions. Once again the two different focal planes can be observed in the $xz$- and $yz$- projections shown in Figs. 7(a10) and (a11). A more detailed investigation of the effects of astigmatism, specifically relating to non-axial ray positions, can be provided by oblique astigmatism, which is shown Figs. 7(b1)–7(b12). It is important to emphasize that vertical and oblique astigmatisms are related by a simple $\pi /4$ rotation of the wavefront aberration as described in the previous paper. Therefore, sampling ray origins along the $x$- or $y$-axis for the case of oblique astigmatism is equivalent to sampling along a ${\pm}45\,\deg$ line in the $xy$-plane for the case of vertical astigmatism. These two cases are presented in Figs. 7(b1)–(b4) and Figs. 7(b5)–(b8), respectively. It is interesting to note from the two intensity images in Fig. 7(b1) that the rays that are sampled along the $x$-axes in the plane at $z = - 20\; \unicode{x00B5}{\rm m}$ are not found to be along the $x$-axes in the focal plane; rather they are linearly distributed along some diagonal line in the $xy$ plane. This explains why the rays never come to a tight focus in the $xz$-plane in the projection shown in Fig. 7(b2). However, in Fig. 7(b3), it is clear that the rays do focus tightly in the $yz$-plane; first, the rays deviate away from the $x$-plane before turning back towards the focal point. These ray paths are further illuminated by the $xy$ projection shown in Fig. 7(b4), which appears to show that the rays twist around by an amount that depends on the distance of the ray origin from the center. The results of sampling the ray origins along the $y$-axis are shown in Figs. 7(b5)–(b8) and are found to be identical to the previous case but for opposite axes. The final row of images in Figs. 7(b9)–(b12) relates to a circular sampling of 12 ray origins in the plane at $z = - 20\; \unicode{x00B5}{\rm m}$. This result is identical to the corresponding result for vertical astigmatism, with a simple $\pi /4$ rotation of the $xy$ plane. In this case, however, it is a little clearer to see the three-dimensional movement of the ray paths from the $xz$ and $yz$ projections in Figs. 7(b10) and 7(b11).

 

Fig. 8. Nonlinear ray tracing applied to coma aberration on the Rayleigh–Sommerfeld lens with numerical aperture of 0.6 and ${{\rm{TEM}}_{00}}$ illumination computed using the Thin Lens Algorithm as described in the first paper. The larger inset image has dimensions of $20\; \unicode{x00B5}{\rm m} \times 20\; \unicode{x00B5}{\rm m}$ and the smaller inset image has dimensions of $6\; \unicode{x00B5}{\rm m} \times 6\; \unicode{x00B5}{\rm m}$.

Download Full Size | PDF

C. Comatic Aberration

The results of nonlinear ray tracing applied to coma aberration are presented in Figs. 8(a1)–8(a12). The first row of images in Figs. 8(a1)–(a4) show the sampling of 10 ray origins along the $x$-axis. It can be seen from the various projections and the two intensity images that the rays do not deviate from the $y = 0$ plane. It can also be seen that the rays do not form a tight focus, with the leftmost rays deviating towards the first crescent moon lobe. The second row of images in Figs. 8(a5)–(a8) relates to the sampling of 10 ray origins along the $y$-axis. In this case, the intensity images and the $xz$ and $xy$ projections reveal a subtle protrusion of the rays towards the $x$ direction and the rays come to a much tighter focus that the previous case, as seen in Fig. 8(a7). The final row of images in Figs. 8(b9)–(b12) relates to a circular sampling of 12 ray origins in the plane at $z = - 20\; \unicode{x00B5}{\rm m}$. It is interesting to observe that those rays on the left side all deviate towards the first crescent moon lobe, while the other rays travel paths into the center focal point. Once again, these observations are broadly consistent with the interpretation of coma aberration in geometrical optics, whereby the rays are interpreted as smearing in the direction of the coma.

D. Trefoil and Quadrafoil

While all of the previous results were for those aberrations that are commonly interpreted using classical geometrical optics, the final two sets of results are for the cases of trefoil and quadrafoil, which cannot easily be interpreted in terms of geometrical optics. Here we attempt to provide an interpretation in terms of the nonlinear ray tracing method proposed in this paper. The results of nonlinear ray tracing applied to trefoil aberration are presented in Figs. 9(a1)–9(a12). As for the previous case of astigmatism, two inset (log of) intensity images appear in each row in the figure that correspond to the intensity at the plane $z = f - 20\; \unicode{x00B5}{\rm m}$ before the focal plane, and at the focal plane. The larger inset image has dimensions of $20\; \unicode{x00B5}{\rm m} \times 20\; \unicode{x00B5}{\rm m}$ and the smaller inset image has dimensions of $10\; \unicode{x00B5}{\rm m} \times 10\; \unicode{x00B5}{\rm m}$. The first row of images in Figs. 9(a1)–(a4) relates to the sampling of 10 ray origins along the $x$-axis. It can be seen from the intensity images in Fig. 9(a1) that these ray origins are along a straight line over the triangular shape that is the diffraction pattern at $f - 20\; \unicode{x00B5}{\rm m}$. It is interesting to note that the triangular pattern that appears at the in-focus plane has inverted with respect to the $x$-axis. It is also interesting to note that rays that originate on the right side are close to the flat right side of the triangle at $f - 20\; \unicode{x00B5}{\rm m}$, but these rays end up on the right point of the in-focus triangle. Rays that originated on the left side close to the left point of the triangle at $f - 20\; \unicode{x00B5}{\rm m}$ end up in the left side of the in-focus triangle or in the nearby diffraction lobe. It can also be seen from the various projections and the two intensity images that the rays remain in the $y = 0$ plane and that the rays do not form a tight focus, following various paths to the triangle focus or nearby lobes. The second row of images in Figs. 9(a5)–(a8) relates to the sampling of 10 ray origins along the $y$-axis, which slice through the center of the triangular diffraction pattern at $f - 20\; \unicode{x00B5}{\rm m}$. From the two intensity images in Fig. 9(a1) it can be seen that these rays end up distributed along a vertical arc through the center of the inverted triangular focus pattern. From the various projections, the rays are clearly seen to bend with the degree of bending dependent on their origin position relative to the center. The last row of images in Figs. 9(a9)–(a12) shows the circular sampling of 12 ray origins in the plane at $z = f - 20\; \unicode{x00B5}{\rm m}$. From the two intensity images and the various projections, it can be seen that those rays originating from points close to the flat sides of the triangular diffraction pattern at $z = f - 20\; \unicode{x00B5}{\rm m}$ end up propagating into the three corners of the in-focus triangle. The degree of ray bending is dependent on their position. This is most evident in Fig. 9(a12) where it can be seen that the rays originating at either the center of the flat sides or at the corners of the triangle at $z = f - 20\; \unicode{x00B5}{\rm m}$ propagate through a single plane and do not bend as they propagate through the focus. Rays that are not at these points all bend as they propagate towards the outer lobes or the corners of the in-focus triangle.

 

Fig. 9. Nonlinear ray tracing for trefoil and quadrafoil aberrations for the Rayleigh–Sommerfeld lens; see text for details.

Download Full Size | PDF

The results of nonlinear ray tracing applied to Quadrafoil aberration are presented in Figs. 9(b1)–9(b12). Here we use $w = 5 \times {10^{- 6}}$, which was found to produce excellent results for nonlinear ray tracing. The first row of images in Figs. 9(b1)–(b4) relates to the sampling of 10 ray origins along the $x$-axis. It can be seen that all of these rays remain in the $xz$-plane and most propagate to a tight focus. Those rays at the extreme edges of the diffraction pattern $z = f - 20\,\,\unicode{x00B5}{\rm m}$ end up propagating into the left and right diffraction spots. The same effect is observed for the rays sampled over the $y$-axis as shown in Figs. 9(a5)–(a8). The final row of images in Figs. 9(b9)–(b12) relates to a circular sampling of two sets of 12 ray origins in the plane at $z = f - 20\,\,\unicode{x00B5}{\rm m}$ shown as red and green spots. It is very interesting to see that those rays on the $x = 0$ and $y = 0$ planes are not highly focused and all follow paths into the outer diffraction spots around the center focus spot. Most inner rays (red spots) focus around the side of the center in-focus spot, while the outer rays (green spots) all follow paths into the surrounding diffraction spots.

5. CONCLUSION

In this paper, we have extended the flux tracing method that was described in the first two papers to include lens aberrations. Although the aberrations can be modeled using the well-known Zernike polynomials for the two underlying algorithms for numerically propagating light at the focal region of a thin lens and ideal lens, respectively, the application of these polynomials has different physical interpretation for both cases.

Aberrations are a particular area of interest for nonlinear ray tracing. Our results for spherical, defocus, astigmatism, and coma were broadly consistent with the classical geometrical interpretation. Similarities between spherical and defocus aberration were drawn, with the latter exhibiting properties similar to the former for high levels of defocus. For the case of spherical aberration a small amount of defocus was also observed. Astigmatism was particularly interesting, and the results revealed the different focal planes for rays that were on the $x = 0$ and $y = 0$ planes. The results were also interesting for rays that originate at positions that are not on these plans. These rays appear to twist into one of the two focal points.

The results for coma, trefoil, and quadrafoil provided insights that cannot be gleaned from classical geometrical optics. The rays could be traced moving into the characteristic lobes and spots surrounding the center spot. The rays could even be traced into the corners and sides of the triangular center spot for the case of trefoil. These results clearly demonstrate the power of the approach as an extension of classical geometrical ray optics.

A fundamental concept in optical design is the wavefront aberration theory developed by many authors including Welford [3] and Hopkins [7]. One important interpretation in this theory is that the wavefront can be used to define a set of ray trajectories based on the wavefront normal, and these can be traced to the image plane to understand the complex patterns of light distribution resulting from optical aberrations. This approach is based on a simplified form of Hamiltonian optics that cannot fully account for the effects of diffraction. Classical wavefront aberration theory is conceptually similar to the approach taken in this paper where we essentially trace the eikonal nonlinearly through the continuum of wavefronts around the image plane, rather than using a single wavefront to project the rays linearly. For this reason, we believe that the method may be more accurate than the classical wavefront aberration theory described above.

Although classical wavefront aberration theory is primarily used in educational settings it remains a valid theory for understanding how aberrations affect image quality. Furthermore, it remains a basic principle in modern lens design, which heavily relies on computational software like Zemax, Code V, or LightTools. These tools integrate ray tracing and physical optics methods (including wavefront analysis) in an advanced and automated way, especially in high precision optical applications such as advanced microscopy, astronomical imaging, and lithography. For such applications, these tools automatically perform wavefront analysis based on ray data, providing detailed insights into aberrations without manually constructing rays normal to the wavefront as in Welford’s method. For this reason, it is not immediately clear if nonlinear ray tracing methods developed in our three papers can guide in lens design, and if so how they might do so. At the very least tracing the flux may lead to a deeper understanding of optical aberration and its manifestation in diffraction patterns.

A final consideration relates to areas of low amplitude. In the first paper it was emphasized that in practical calculations such as the Thin Lens Algorithm and the Ideal Lens Algorithm, it is common to encounter zones where the amplitude is very small or even negligible, resulting either from destructive interference at these points or the area being outside the support of the field. We noted that in these areas, interpreting the phase should be avoided. First, the phase value in such regions is susceptible to low accuracy caused by small error in the ASM or in machine error. Second, and more importantly, the phase itself loses any physical significance. For these reasons, we should make no attempt to be overly concerned with the accuracy of the phase in areas of very low amplitude, and we should make no attempt to examine the eikonal in these areas. Certainly, we should avoid originating a ray at a point where the amplitude is very low, and so long as the eikonal trace is accurate we can be reasonably confident that it will never pass through an area of low amplitude and indeterminate (or inaccurate) phase that might send it on an erroneous path. Indeed we have seen in the results presented in this paper that the eikonal will seek to avoid areas of low amplitude, always bending away from or around areas of destructive interference in the diffraction pattern, and towards areas of higher amplitude, which is of course in agreement with the concept of flux.

In order to validate the geometrical approximation in the eikonal traces for each of the results in this paper, values are provided in Supplement 1 for the mean and maximum values of the geometric failure function, $F$, over a central eikonal trace as described in the second paper. The code used to implement the two propagation algorithms and the eikonal traces in this paper is provided in Code 1, Ref. [8].