## Abstract

We show that the volumetric field distribution in the focal region of a high numerical aperture focusing system can be efficiently calculated with a three-dimensional Fourier transform. In addition to focusing in a single medium, the method is able to calculate the more complex case of focusing through a planar interface between two media of mismatched refractive indices. The use of the chirp z-transform in our numerical implementation of the method allows us to perform fast calculations of the three-dimensional focused field distribution with good accuracy.

©2012 Optical Society of America

## 1. Introduction

The field distribution in the focal region of a focusing lens has been studied extensively due to the presence of this element in a large number of optical systems [1]. The Huygens-Fresnel principle, with a few approximations to simplify the diffraction integral, is often applied in the analysis of this particular diffraction problem. For instance, the Fresnel diffraction integral can be used to calculate the focused field for a lens with low numerical aperture (NA) [2], whereas the scalar Debye integral is a good approximation for a focusing system with higher NA [3]. It is well known, however, that the focused field experiences a change in its polarization state under the high-NA condition [4]. Therefore, a vectorial calculation has to be adopted in order to evaluate all the polarization components at the focus. The vectorial Debye integral, usually known as the Debye-Wolf integral, established by Richards and Wolf [4], is the most widely used vectorial diffraction method. For some simple cases, this method can give a closed-form expression for the three-dimensional (3D) field distribution around the focus. Nevertheless, numerical evaluation of the integral is necessary for a more general analysis with an arbitrary field incident on the lens. Since the integrand of the Debye-Wolf integral involves highly oscillatory functions, extra caution and a long computation time are often needed to obtain accurate numerical results. Recently, Leutenegger *et al.* demonstrated that the 3D vectorial field distribution at the focus can be computed plane by plane under a proper transformation of the Debye-Wolf integral [5]. The field in a plane at a given axial position is obtained by a two-dimensional Fourier transform, which can be implemented efficiently on a computer with the help of algorithms such as the fast Fourier transform (FFT) and the chirp z-transform (CZT) [6].

In an alternative approach to the evaluation of the diffraction integral, McCutchen showed that the scalar Debye integral can be written as a three-dimensional Fourier transform (3D-FT) [7]. Furthermore, McCutchen argued that the method is applicable to the vectorial case if it is applied to each component of the electromagnetic (EM) field separately [8]. The suitability of McCutchen’s method to model focusing of complex incident polarization states has been previously analyzed [9]. In this work, we show the equivalence between McCutchen’s method and the Debye-Wolf integral. In other words, we show that the Debye-Wolf integral can be written as a 3D-FT, simplifying the calculation of the complete volumetric three-dimensional field distribution in the focal region of a high-NA system. The use of the CZT algorithm in the implementation of our method allows for fast and accurate calculations when compared with the results obtained from direct evaluation of the Debye-Wolf integral. We present a number of examples showing the suitability of the 3D-FT method in different scenarios.

## 2. Light focused into a homogenous medium

The change of polarization in the focal region of a high-NA focusing system comes from refraction at the system, which can be represented by a transformation of the electric field components at the reference surface representing a perfect converging wave. For an aplanatic system, which is the class of systems considered in this work, the reference surface is a spherical shell limited by the NA (Fig. 1 ).

Consider a non-paraxial incident monochromatic wave having three orthogonal local polarization components along the directions depicted by unit vectors ${e}_{\rho}$, ${e}_{\varphi}$, and ${e}_{z}$ shown in Fig. 1. At the reference surface, vectors ${e}_{\rho}$ and ${e}_{z}$ transform into ${e}_{\theta}$ and ${e}_{r}$, respectively, upon refraction, while ${e}_{\varphi}$ remains the same. Thus, for incident linearly polarized light, the transformation results in a space-variant polarization state that generates polarization components other than the original ones. This transformation can be described as follows

The factor ${A}_{0}\sqrt{\mathrm{cos}\theta},$ with ${A}_{0}$ a constant, is the apodization function to account for energy conservation [4]. The transformation given by ${M}_{1}$ in Eq. (2) is a matrix representation of the transformation required to calculate the strength factor in the Debye-Wolf integral. The propagation after refraction of the three orthogonally polarized components (in *x-*, *y-*, and *z-*direction) is governed by the individual scalar Debye integral

It has been shown that the scalar Debye integral can be rewritten in the form of a scaled 3D-FT [7,10]. That is, Eq. (3) can be rewritten as

In Eq. (4) the pupil function represents an infinitesimally thin spherical shell defined as $S(Q)=\delta (r-1),$ with $\delta (\cdot )$ the Dirac delta distribution and $r=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}}.$

Now, as stated by Richards and Wolf [4], the Debye-Wolf integral for the electric field can be written as

Considering a perfect lens, i.e. with no aberrations, and substituting $s$ and ${r}_{p}$ from Eq. (6) in Eq. (5), we find that

Noting that the strength factor, $a(\theta ,\varphi )$, is only defined on the surface of a sphere of radius $k$, since $s$ is a unitary vector, it is apparent that Eq. (7) can be written as the following triple integral

Implementing the 3D-FT in either Eq. (4) or Eq. (8), after transformation to Cartesian coordinates, using the FFT algorithm would result in an output space with the same transverse dimension and sampling size as the aperture. However, the region of interest (ROI) in the high-NA case is only a small portion at the center (~2μm × 2μm) as compared with the aperture size of the lens (typically a few mm in radius). To sufficiently resolve details of the ROI, zero padding is often used to increase the output sampling size. Since the ROI only accounts for a small part of the whole output, most of the additional output sampling points brought in by the padding will fall out of the ROI. Furthermore, the zero padding inevitably increases the size of the 3D input matrix by tens of times, increasing the demand for computational resources. Therefore, we opted for the CZT algorithm [6] to evaluate the 3D-FT. At first, the CZT algorithm appears slower than the FFT in dealing with an input with size of a power of two. Nevertheless, this algorithm allows us to avoid the use of zero padding and select only the ROI as the output and, hence, turns out to be much faster than the FFT algorithm in our calculations. For a size of the 3D output space of 128×128×128 pixels, the calculation time for our method is typically a few seconds on a personal computer, whereas it takes tens of minutes for the direct integration of the Debye-Wolf integral. The gain in computation time does not reduce the accuracy of the results. Figure 3
is a comparison between different cross sections of the *E _{x}* component in the focal region of an aplanatic system, with NA=0.95, for incident light linearly polarized in the

*x*-direction, λ

_{0}=532nm, as obtained by direct evaluation of the Debye-Wolf integral and the 3D-FT method. The relative error in the focal volume between the two methods was calculated as [11],

*z*. This figure shows that the accuracy of the 3D-FT method is best at the focal plane, reduces as we move away in the

*z*-direction, and slightly improves at the edge of the volume analyzed. As can be expected, the accuracy is worse at planes close to the minima in the focused field axial distribution. This behavior reflects the role of computational precision rather than a limitation of the 3D-FT method.

To complete our comparison between the two methods, Figs. 4b and 4c show a cross-section along the *x*- and *y*-direction in the focal plane for $\Delta $, the magnitude of the difference of the electric field obtained with the two methods, defined as

^{−3}, which is enough for a large number of applications.

An appealing characteristic of the 3D-FT method presented in this work is that it can handle an incident field with arbitrary polarization and complex amplitude distribution. For instance, we can use the method to calculate the field distribution at the focal region of a high-NA lens for an incident vortex beam [12]. A vortex beam linearly polarized in the *x*-direction, with unity topological charge, has a space-variant phase distribution, in the transverse plane, and can be expressed as${E}^{in}=\left\{{E}_{0}\mathrm{exp}(i\varphi ),0,0\right\},$ where ${E}_{0}$ is a constant and $\varphi $ the beam’s phase. The intensity of the resultant 3D vectorial focal field distribution, for the same set of parameters used in the previous example, is shown in Fig. 5
. As a consequence of the non-negligible longitudinal component of the focused field (Fig. 5c), the total intensity distribution in the focal plane (Fig. 5d) is no longer zero at the center, as it was the case for the scalar calculation [13].

Another interesting example of polarization distributions is that of the so-called cylindrical beams. A typical example of these beams is a radially polarized beam [14], for which the electric field distribution can be written as ${E}^{in}={E}_{0}\left\{\mathrm{cos}\varphi ,\mathrm{sin}\varphi ,0\right\}.$ Fig. 6 shows the different components of the field distribution in the focal region for a radially polarized beam focused by a high-NA system. The parameters used in the calculation are the same as in the previous examples. The radial component in Fig. 6 was obtained as ${E}_{r}={E}_{0}\mathrm{cos}\varphi +{E}_{0}\mathrm{sin}\varphi ,$ whereas the azimuthal component is not shown because it is negligible.

Figure 7 is a comparison between different cross sections of the focused field components for incident radial polarization as obtained by direct evaluation of the Debye-Wolf integral and the 3D-FT method. A good agreement between both methods is observed from the figure. The relative error in this case, as obtained from Eq. (9), is $\epsilon =1.508\times {10}^{-3}$. As in the case for incident light linearly polarized, in Fig. 8a we present the axial distribution of $\epsilon $ for planes parallel to the focal plane in order to determine where, within the focal volume, the accuracy of the 3D-FT method is best. The axial distribution of $\epsilon $ shows that the accuracy is, again, best at the focal plane and worst at planes near the minima of the focused field intensity axial distribution (Fig. 7d). The accuracy in the transverse direction, within the focal plane, is given by the distribution of Δ given in Fig. 8b. The total intensity distribution for the radial polarization is rotationally symmetric. Therefore, the distribution of Δ is given in terms of the radial distance from the focus. As expected from the previous discussion about the computational precision, the accuracy of the 3D-FT method is worse –Δ is larger– at distances that correspond to regions near the minima in the intensity distribution.

## 3. Light focused through an interface

So far, we have only considered focusing in a single homogeneous medium. However, it is possible to have in the focal region an interface between two media of mismatch refractive indices. For instance, when imaging biological specimens in an aqueous environment through an oil-immersion microscope objective. The refraction occurred at the interface between the two media is known to cause spherical aberration [15], and a shift in the position of the focus, as shown in Fig. 9 . The optical field right after the interface can be expressed as

*s*-polarized and

*p*-polarized light, respectively [3]. The additional aberration function, due to the interface between the two media, is $\Phi ({\theta}_{2})=\mathrm{exp}\left[i{k}_{0}d\left({n}_{2}\mathrm{cos}{\theta}_{2}-{n}_{1}\mathrm{cos}{\theta}_{1}\right)\right]$, with ${k}_{0}=2\pi /{\lambda}_{0}$ and $d$ the distance between the focus and the interface. Angles ${\theta}_{1}$ and ${\theta}_{2}$ are related through Snell’s law. The vectorial focal field in the second medium is given by

*x*-polarized uniform incident wave (λ

_{0}=532nm, n

_{1}=1.52, n

_{2}=1.33, and NA=1.2) focused through an interface located 5μm in front of the geometrical focus, the focal shift brought by the index mismatch is shown in Fig. 10a . Good agreement is observed when the result is compared with that obtained by direct evaluation of the integral for the axial distribution given in [17] (Fig. 10b).

## 4. Conclusion

We have shown the equivalence between McCutchen’s method and the Debye-Wolf integral by rewriting the latter as a three-dimensional Fourier transform. Thus, the three-dimensional vectorial field distribution in the focal region of a high-NA system predicted by Debye-Wolf integral can be accurately calculated with the 3D-FT method. This method is able to deal with an arbitrary input field, including vortex and cylindrical beams. We have also shown that the method is applicable to focusing through an interface between two media of mismatch refractive indices by introducing the spherical aberration function due to the interface. The method can be extended to a stratified medium [18] and other scenarios where the Debye-Wolf integral is applicable. In this case, the transformation matrix has to be modified accordingly. Note that, however, this modification is also required in the calculation of the strength factor for the direct integration of the Debye-Wolf integral.

The proposed method facilitates the fast vectorial calculation of a focal field. One important advantage is that the calculation time of the 3D-FT method does not depend on the complexity of the incident field, whereas the calculation time to evaluate the Debye-Wolf integral is largely determined by the input. The integration over the azimuthal angle in the Debye-Wolf integral can be done analytically for uniform input polarization (e.g. linearly polarized light). This leaves only one integral to be performed numerically, for each component of the field (*x*, *y*, and *z*), which can speed up the calculation. However, even in these favorable conditions for direct integration of the Debye-Wolf integral, the 3D-FT method is significantly faster. We compared the speed using a volume of 4λ×4λ×4λ around the focal point, for cubic pixels with side λ/30 (121×121×121 points) using MATLAB on a PC with an Intel Core2 Duo 2.53 GHz processor and 4 GB of RAM. The typical calculation time for direct evaluation of the Debye-Wolf integral is around 23 minutes. The same result obtained with the 3D-FT method takes approximately 2 seconds, and the accuracy achieved, measured as the relative error, is in the order of 10^{−3}. The maximum benefit of the method presented here is obtained when the field in the whole 3D focal volume is desired. If the focused field is required only in a small number of points (e.g. a single plane at the focus) direct integration of Debye-Wolf integral may be sufficiently fast.

In conclusion, the 3D-FT method, and its implementation using the CZT algorithm, allows for fast, efficient, and accurate calculation of the volumetric focused field distribution when compared with direct integration of the Debye-Wolf integral.

## Acknowledgments

J. Lin acknowledges the fellowship support from the Agency for Science, Technology and Research (A*STAR), Singapore. O. G. Rodríguez-Herrera, F. Kenny, and J. C. Dainty are grateful to Science Foundation Ireland for its support through grant 07/IN.1/I906.

## References

**1. **J. J. Stamnes, *Waves in Focal Regions* (Hilger, 1986).

**2. **J. W. Goodman, *Introduction to Fourier Optics* (McGraw-Hill, 1996).

**3. **M. Born and E. Wolf, *Principles of Optics*, 7th ed. (Cambridge University Press, 1999).

**4. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems. II. Structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A Math. Phys. Sci. **253**(1274), 358–379 (1959). [CrossRef]

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

**6. **A. V. Oppenheim, R. W. Schafer, and J. R. Buck, *Discrete-time signal processing*, 2nd ed. (Prentice Hall, 1999).

**7. **C. W. McCutchen, “Generalized aperture and the three-dimensional diffraction image,” J. Opt. Soc. Am. **54**(2), 240–244 (1964). [CrossRef]

**8. **C. W. McCutchen, “Generalized aperture and the three-dimensional diffraction image: erratum,” J. Opt. Soc. Am. A **19**(8), 1721–1721 (2002). [CrossRef]

**9. **I. Iglesias and B. Vohnsen, “Polarization structuring for focal volume shaping in high-resolution microscopy,” Opt. Commun. **271**(1), 40–47 (2007). [CrossRef]

**10. **J. Lin, X.-C. Yuan, S. S. Kou, C. J. R. Sheppard, O. G. Rodríguez-Herrera, and J. C. Dainty, “Direct calculation of a three-dimensional diffracted field,” Opt. Lett. **36**(8), 1341–1343 (2011). [CrossRef] [PubMed]

**11. **P. Török, P. R. T. Munro, and E. E. Kriezis, “Rigorous near- to far-field transformation for vectorial diffraction calculations and its numerical implementation,” J. Opt. Soc. Am. A **23**(3), 713–722 (2006). [CrossRef] [PubMed]

**12. **D. Ganic, X. Gan, and M. Gu, “Focusing of doughnut laser beams by a high numerical-aperture objective in free space,” Opt. Express **11**(21), 2747–2752 (2003). [CrossRef] [PubMed]

**13. **J. Lin, X.-C. Yuan, S. H. Tao, and R. E. Burge, “Variable-radius focused optical vortex with suppressed sidelobes,” Opt. Lett. **31**(11), 1600–1602 (2006). [CrossRef] [PubMed]

**14. **Q. Zhan, “Cylindrical vector beams: from mathematical concepts to applications,” Adv. Opt. Photon. **1**(1), 1–57 (2009). [CrossRef]

**15. **C. J. R. Sheppard and M. Gu, “Axial imaging through an aberration layer of water in confocal microscopy,” Opt. Commun. **88**(2-3), 180–190 (1992). [CrossRef]

**16. **P. Török, P. Varga, Z. Laczik, and G. R. Booker, “Electromagnetic diffraction of light focused through a planar interface between materials of mismatched refractive indices: an integral representation,” J. Opt. Soc. Am. A **12**(2), 325–332 (1995). [CrossRef]

**17. **S. H. Wiersma, P. Török, T. D. Visser, and P. Varga, “Comparison of different theories for focusing through a plane interface,” J. Opt. Soc. Am. A **14**(7), 1482–1490 (1997). [CrossRef]

**18. **P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt. **36**(11), 2305–2312 (1997). [CrossRef] [PubMed]