Fast focus field calculations

We present a fast calculation of the electromagnetic field ne ar the focus of an objective with a high numerical aperture (NA) . Instead of direct integration, the vectorial Debye di ffraction integral is evaluated with the fast Fourier transform for calculating the electromagn etic field in the entire focal region. We generalize this concept with the chi rp z transform for obtaining a flexible sampling grid and an additional gain in computation speed. Under the conditions for the validity of the Debye int gral representation, our method yields the amplitude, phase and polar ization of the focus field for an arbitrary paraxial input field on the object ive. We present two case studies by calculating the focus fields of a 40 × 1.20 NA water immersion objective for di fferent amplitude distributions of the input field, and a 100× 1.45 NA oil immersion objective containing evanescent field contributions for both linearly and radially polarized inp ut fields. © 2006 Optical Society of America OCIS codes: (220.2560) Optical design and fabrication, focus; (260.19 60) Physical optics, diffraction theory; (070.2580) Fourier optics and optical sign al processing, Fourier optics; (180.0180) Microscopy. References and links 1. P. Debye, “Das Verhalten von Lichtwellen in der Nähe eine s Brennpunktes oder einer Brennlinie,” Ann. Phys. 30,755–776 (1909). 2. E. Wolf, “Electromagnetic di ffraction in optical systems, I. An integral representation o f the image field,” Proc. R. Soc. London Ser. A253,349–357 (1959). 3. B. Richards, E. Wolf, “Electromagnetic di ffraction in optical systems, II. Structure of the image field i n an aplanatic system,” Proc. R. Soc. London Ser. A 253,358–379 (1959). 4. Typically, a good accuracy is achieved for M & 50 andN & 200 sampling points. 5. P. Török, P. Varga, “Electromagnetic di ffraction of light focused through a stratified medium,” Appl. Opt. 36, 2305–2312 (1997). 6. J.J. Stamnes, Waves in Focal Regions: propagation, di ffraction and focusing of light, sound and water waves, Hilger, Bristol UK (1986). 7. G. Mikula, A. Kolodziejczyk, M. Makowski, C. Prokopowicz , M. Sypek, “Diffractive elements for imaging with extended depth of focus,” Opt. Eng. 44,058001 (2005). 8. N. Huse, A. Schönle, S.W. Hell, “Z-polarized confocal mi croscopy,” J. Biomed. Opt. 6, 480–484 (2001). 9. J. Enderlein, I. Gregor, D. Patra, T. Dertinger, U.B. Kaup p, “Performance of Fluorescence Correlation Spectroscopy for Measuring Di ffusion and Concentration,” Chem. PhysChem. 6, 2324-2336 (2005). 10. For simplification, the sample indices kl andmnwill be omitted further on. 11. M. Mansuripur, “Certain computational aspects of vecto r diffraction problems,” J. Opt. Soc. Am. A 6, 786–805 (1989). 12. M. Sypek, “Light propagation in the Fresnel region. New n umerical approach,” Opt. Comm. 116,43–48 (1995). 13. P. Luchini, “Two-dimensional numerical integration us ing a square mesh,” Comp. Phys. Comm. 31, 303–310 (1984). 14. J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt.41,4897– 4903 (2002). #75759 $15.00 USD (C) 2006 OSA Received 3 October 2006; accepted 28 October 2006; corrected 26 April 2007 13 November 2006/ Vol. 14, No. 23/ OPTICS EXPRESS 11277 15. Y. Li, E. Wolf, “Three-dimensional intensity distribut ion near the focus in systems of di fferent Fresnel numbers,” J. Opt. Soc. Am. A1, 801–808 (1984). 16. W. Hsu, R. Barakat, “Stratton-Chu vectorial di ffraction of electromagnetic fields by apertures with applica tion to small-Fresnel-number systems,” J. Opt. Soc. Am. A 11,623–629 (1994). 17. E. Wolf, Y. Li, “Conditions for the validity of the Debye i ntegral representation of focused fields,” Opt. Comm. 39,205–210 (1981). 18. P. Török, “Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number,” J. Opt. Soc. Am. A15,3009–3015 (1998).


INTRODUCTION
The plane wave spectrum (PWS) method is a well-known and efficient technique for calculating the propagation and diffraction of electromagnetic (EM) fields.Its efficiency lies in the ability to propagate EM fields from one plane to another using the fast Fourier transform (FFT).This concept is the essence of the Debye approximation and is often used for the calculation of the EM focus field. 1,2 owever, for focal field calculations in microscopy, in particular for optical systems with high numerical aperture (NA), this classical problem turns into a computational challenge due to the highly oscillatory behaviour of the involved functions.In addition, polarization effects cannot be neglected rendering the calculation long and tedious.Recent techniques in microscopy and tomography such as the extended focus field, 3 microscopy beyond the Abbe resolution limit and point-spread function engineering as advanced by S. Hell and his group, 4 or rigorous ab initio calculations for fluorescence fluctuation spectroscopy 5 amplify the demand for fast focus field calculations.
In this paper we revisit the Debye approximation and propose a novel and flexible implementation of the Debye integral incorporating the effects of amplitude, phase and polarization in an overall manner.This new implementation is particularly suited for rapid numerical evaluation and requires substantially less effort for calculating the amplitude, phase and polarization of an EM field distribution generated by a high NA microscope objective.Section 2 introduces the Debye approximation (c.f.Leutenegger et al. for details 6 ).Section 3 presents selected examples, firstly the calculation of a few focus fields obtained with a standard 100 × 1.45NA oil immersion objective, and secondly, the calculation of an extended polychromatic focus field (Bessel beam).
The objective is represented by the aperture stop A with radius R, the principal planes P 1 and P 2 with vertex points V 1 and V 2 , and the foci F 1 and F 2 .The focal length f is given as The point P is the intersection point of a ray with P 2 and shows the relation of the position r at P 1 of the incident wave E i to the propagation angle θ at P 2 of the transmitted wave E t .

THEORY
This section establishes the basic formalism based on the Debye diffraction integral and the formulation of this integral as a Fourier transform. 6The basic optical layout and the respective coordinate systems are shown in Fig. 1.We assume that this imaging system obeys Abbe's sine condition, which is usually fulfilled for microscope objectives.
A coherent, monochromatic, paraxial wave field crosses the aperture stop A, propagates towards the principal plane P 1 and is transferred to the principal plane P 2 .At P 2 , the wave field is refracted and focused towards the focal point F 2 .The point P lies on the principal plane P 2 and illustrates the focusing of a ray at P 2 towards the focal point F 2 .The spherical surface P 2 is centered at F 2 and the deflection angle θ at the position P is given by where r is the off-axis coordinate of the incident wave on P 1 , R the aperture stop radius, NA the numerical aperture of the objective and n t the index of refraction behind the P 2 surface.Because the aperture A is placed in the back focal plane, the imaging system is telecentric.Within our representation, the wave propagation from the aperture plane A to the principal plane P 1 can be calculated in most cases with classical Fourier optics principles.
For calculating the transmitted field E t (θ, φ) at P 2 , the incident field E i (r, φ) at P 1 is typically decomposed into a radial (p-polarized) and a tangential component (s-polarized).Hence, the amplitude, phase and polarization of the transmitted field at P 2 is where t p (θ, φ) and t s (θ, φ) are the transmission coefficients and e p (θ, φ), e r (θ, φ) and e s (θ, φ) the unit vectors for p-and s-polarization, respectively.Accumulated phase distortions (aberrations) at P 2 as well as attenuations (apodization) caused by the objective are integrated in the complex transmission coefficients t p and t s .As we assume the incident field to be paraxial, the axial component E iz is small against the lateral components E ix,y and can be neglected even if the incident phase is not constant.In the Debye approximation, the transmitted field E t is the plane wave spectrum of the focus field E Send correspondence to Marcel Leutenegger: marcel.leutenegger@a3.epfl.chnear F 2 .Therefore, the electric field E at a point (x, y, z) near the focus F 2 is obtained by integrating the propagated plane waves, that is The phase factor e ik z z accounts for the phase accumulation when propagating along the z-axis, whereas the term e −i(k x x+k y y) represents the phase difference of the wave front at off-axis points (x, y, z) with respect to the on-axis point (0, 0, z).The integration extends over the solid angle Ω under which P 2 is observed at F 2 , i.e. sin Θ = NA/n t .The wave vector k t is simply given by where The evaluation of Eq. ( 3) is usually performed with a direct numerical integration taking into account the coordinate transformations, which results in the Richard-Wolf integral representation. 7,8 nstead of the common ansatz, a (θ, φ)sampling with constant dΩ = sin θ dθ dφ was introduced.Besides minimizing the number of sampling points along θ, the calculation of the integrand and its integration can be merged in a single matrix product resulting in a further reduction of the computation time.This evaluation of the Debye diffraction integral (3) is quite fast but still much slower than the conventional computation of a Fraunhofer diffraction integral.However, Eq. ( 3) can be easily rewritten as a Fourier transform by splitting the phase factor into a lateral and an axial term, and by performing the integration over P 1 instead of P 2 .Using Eq. ( 1) and ( 4), the integration step dΩ for a sampling over P 2 is projected onto P 1 , which yields Insertion of this sampling step into Eq.( 3) results in Extending now the integration over (k x , k y ) ∈ R 2 by setting | E t | = 0 for r > R allows to rewrite the Debye diffraction integral as a Fourier transform of the weighted field E t , which finally results in This is the main result of this section.The Debye integral is now expressed as a two-dimensional Fourier transform F x,y of the field distribution in the aperture A projected on P 2 , which can be evaluated with the fast Fourier transform (FFT) or the chirp z transform (CZT). 6The similarity of this expression with the conventional Fraunhofer diffraction integral is obvious.For a low NA imaging system, the weighting factor is approximated by 1/ cos θ ≈ 1 and Eq. ( 7) is equivalent to the Fraunhofer diffraction integral.

APPLICATIONS
In this section, we present several case studies by calculating the focus fields achieved with a standard high NA oil immersion objective.The focus fields in the vicinity of the interface into the aqueous sample are presented for three cases: (1) a simple glass-water interface, (2) a glass-supported gold film-water interface and (3) a glass-supported silver filmwater interface.Last but not least, the calculation of an extended polychromatic focus field generated by a Bessel beam is presented.This extended focus field is of particular interest for Fourier domain optical coherence tomography (FDOCT) because it preserves a lateral resolution of a few micrometers over an axial distance in the millimeter range.

Focus fields with oil immersion objectives
Current state of the art oil immersion objectives provide a very high lateral and axial resolution for imaging in the vicinity of a cover slide-sample interface.These objectives allow for instance total internal reflection excitation in the evanescent field penetrating only a fraction of a wavelength into the sample, which provides an outstanding axial confinement for highly surface-sensitive measurements.The following calculations outline the performance ideally achieved by a common oil immersion objective in case of confocal illumination.The calculations assume a 100 × 1.45NA oil immersion objective with an aperture diameter Ø A = 6mm and an x polarized laser beam with 633nm wavelength and a Gaussian e −2 beam waist of w A = 5mm at the aperture A.  As a first example, the focus field at a cover glass-water interface is calculated.The intensity distribution of the transmitted field E t is shown in Fig. 2. At supercritical angles, i.e. for the outer aperture region where NA > 1.33, the evanescent field is easily identified by its characteristic intensity increase.As the incident polarization is linear along the x axis, the x component (red) is dominant except for a strong z component (blue) in the evanescent field.Figure 3 shows the corresponding focus field in the aqueous sample.The focus field is outlined with iso-intensity surfaces of I = e −1..−4 I max .In a confocal laser scanning microscope, a lateral resolution of ∆x ≈ 220nm and ∆y ≈ 180nm would be achieved.
As second example, the cover glass is covered with a 38nm thin gold film and the focus field at the gold-water interface is calculated.A relative dielectric constant Au = −9.39+ 1.15i of gold is assumed.The film thickness of 38nm optimizes surface plasmon coupling near the aperture edge as shown in Fig. 4. Due to the dominant z polarization and the angular filtering by the plasmon coupling, the focus field in Fig. 5 shows a very particular profile similar to a tunnel along the y axis.
As third example, the cover glass is covered with a 55nm thin silver film and the focus field at the silver-water interface is calculated.A relative dielectric constant Ag = −18.28+ 0.46i of silver is assumed.With a film thickness of 55nm, the optimal surface plasmon coupling occurs for an NA of 1.40, which is well inside the aperture field shown in Fig. 6. Figure 7 shows the focus field, which inherits its rotational symmetry from the incident field with radial polarization.In the central peak, the z polarization is dominant whereas xy polarization prevails in the first side lobe.In comparison with the previous example, the plasmon coupling through the silver is much stronger than through gold because the optimal film thickness and the optimal incidence angle are both achievable with silver.

Extended depth of focus for Fourier domain optical coherence tomography
Optical coherence tomography (OCT) is a well-known, non-invasive, three-dimensional optical imaging method.Using a broadband light source, axial sectioning is achieved by interfering the back-reflected light from the sample with unaltered light from the reference arm of similar optical path length.In OCT, the axial resolution ∆z is therefore (only) limited by the (detected) bandwidth ∆λ of the light source.As with conventional optical imaging, its lateral resolution ∆r scales with λ 0 /NA, where λ 0 is the free space wavelength and NA the numerical aperture of the imaging optics.Using a coherence gating in OCT, it is therefore possible to decouple axial and lateral resolution.In addition, due to the coherent amplification, OCT can achieve very high sensitivity for imaging the structure of transparent biological samples, i.e. small changes in refraction index already yield sufficient contrast for quantitative measurements.
Instead of a sequential depth scanning with the reference arm, Fourier domain OCT (FDOCT) simultaneously reads the full depth profile encoded in the measured interferogram spectrum. 9,10 he sample structure is then obtained by Fouriertransforming these measured spectra.The main advantages of FDOCT are (a) a much faster acquisition as only lateral sample scanning is required, 11 (b) an improved phase stability due to fixed optical paths 12 and (c) an even better sensitivity (Fourier advantage).However, as the depth profile is recorded at once, a sharp image of lateral features is only obtained within the depth of focus (DOF).As the DOF ∝ NA −2 and ∆r ∝ NA −1 , enlarging the DOF usually means sacrificing lateral resolution.This limitation was overcome by using a Bessel beam for the sample illumination. 13The Bessel beam provides  a high lateral confinement (∆r in the order of λ 0 ) over a very long axial distance z, i.e.DOF ∆r. Figure 8 shows the optical setup for creating the Bessel beam.The collimated broadband laser beam with a Gaussian beam profile suffers a conical wavefront distortion by passing through a linear axicon.This converging conical beam creates the desired Bessel pattern in the overlap region behind the axicon.A telescope consisting of a lens and the microscope objective demagnify this Bessel pattern and image it into the sample.Figure 9 shows the spectrum of the illuminating Bessel beam. Figure 10 shows the calculated Bessel beam consisting of an incoherent superposition of focus fields for λ ∈ [720nm, 870nm] with 1nm steps.The illuminating needle of 550µm FWHM length has a central spot diameter of 2∆r = 2.5µm only!Note also that the initial linear polarization of the laser beam is well conserved because of the moderate NA.

CONCLUSIONS
We showed a fast and accurate implementation of the vectorial Debye integral for calculating the focus field of high NA objectives for arbitrary amplitude, phase and polarization distributions of the input field.The numerical evaluation with the fast Fourier transform is very efficient and allows a high flexibility of the input field.With the general chirp z transform, we extended our calculations to low NA focus fields requesting a non-linear scaling as shown by Li and Hsu 14,15 (c.f.Sec.3.2 for example).For low NA, the method converges quite naturally to a focus field given by the Fraunhofer approximation.In addition, we used a generalized pupil function (apodization) of high NA objectives taking into account amplitude and polarization distributions.The pupil function incorporates the Fresnel transmission coefficients and can contain wave front aberrations as induced by real objectives.Based on these Fresnel coefficients, it is straightforward to include wave propagation through stratified media as shown with example (2) and (3) in Sec.3.1.
In summary, our method allows fast and accurate calculations of the focus field in the entire focal region, which opens the path to fast simulations for point spread function engineering and image deconvolution in three-dimensional light microscopy.

Figure 2 .
Figure 2. Transmitted field components I xyz ∝ |E t,xyz | 2 .The x component is shown in red, y in green and z in blue, respectively.

Figure 3 .
Figure 3. Focus field at the cover glass-water interface obtained with a 100×1.45NAoil immersion objective.The incident laser beam has a Gaussian e −2 waist of 5mm to overfill the objective aperture A of 6mm diameter.The laser beam is x polarized and has a wavelength of 633nm.

Figure 5 .
Figure 5. Focus field at the cover glass-38nm gold film-water interface obtained with a 100 × 1.45NA oil immersion objective.

Figure 7 .
Figure 7. Focus field at the cover glass-55nm silver film-water interface obtained with a 100×1.45NAoil immersion objective and radial polarization of the incident field.

Figure 8 .Figure 9 .
Figure 8. Sample illumination with a linear axicon and a telescope.The incident laser beam has a Gaussian waist of 3.4mm at λ 0 = 785nm wavelength.The power spectrum of the fs-pulsed Ti:Sapphire laser is shown in Fig. 9.The linear axicon (20mm clear aperture, 178f ull cone angle, BK7 glass) induces a converging conical wavefront.The lens ( f l = 450mm focal distance, 22mm clear aperture) and the 10 × 0.30NA objective ( f o = 16.4mm,Ø A = 10mm aperture) image the Bessel pattern into the sample with a 28× demagnification.

Figure 10 .
Figure 10.Illuminating Bessel beam for an extended DOF.Note the different scaling of the z axis.