Fast focus field calculations

doi:10.1364/OE.14.011277

« Show journal navigation

Fast focus field calculations

Marcel Leutenegger, Ramachandra Rao, Rainer A. Leitgeb, and Theo Lasser »View Author Affiliations

Optics Express, Vol. 14, Issue 23, pp. 11277-11291 (2006)
http://dx.doi.org/10.1364/OE.14.011277

View Full Text Article

Acrobat PDF (531 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. The Debye diffraction integral as Fourier transform
3. Numerical implementation
4. Selected examples
5. Conclusions
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (18)
Cited By
Figures (12)
Metrics

Abstract

We present a fast calculation of the electromagnetic field near the focus of an objective with a high numerical aperture (NA). Instead of direct integration, the vectorial Debye diffraction integral is evaluated with the fast Fourier transform for calculating the electromagnetic field in the entire focal region. We generalize this concept with the chirp z transform for obtaining a flexible sampling grid and an additional gain in computation speed. Under the conditions for the validity of the Debye integral representation, our method yields the amplitude, phase and polarization of the focus field for an arbitrary paraxial input field on the objective. We present two case studies by calculating the focus fields of a 40×1.20 NA water immersion objective for different 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 input fields.

1. 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).

In microscopy this concept is the essence of the Debye approximation and is often used for the calculation of the EM field [1

P. Debye, “Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie,” Ann. Phys. 30, 755–776 (1909). [CrossRef]

, 2

E. Wolf, “Electromagnetic diffraction in optical systems, I. An integral representation of the image field,” Proc. R. Soc. London Ser. A 253, 349–357 (1959). [CrossRef]

, 3

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

] near the focus of high numerical aperture (NA) objectives. Török et al. considerably expanded this concept for studying the focal field distribution and its distortions in stratified media commonly encountered in optical microscopy [5

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

]. For a general and historical review on diffraction theory the reader is referred to Stamnes [6

J.J. Stamnes, Waves in Focal Regions: propagation, diffraction and focusing of light, sound and water waves , Hilger, Bristol UK (1986).

However, for focal field calculations in microscopy, in particular for optical systems with high NA, this classical problem turns into a computational challenge due to the highly oscillatory behavior of the involved functions. In addition, polarization effects cannot be neglected rendering this calculation long and tedious. Recent techniques in microscopy and tomography such as the extended focus field [7

G. Mikula, A. Kolodziejczyk, M. Makowski, C. Prokopowicz, and M. Sypek, “Diffractive elements for imaging with extended depth of focus,” Opt. Eng. 44, 058001 (2005). [CrossRef]

], microscopy beyond the Abbe resolution limit and pointspread function engineering as advanced by S. Hell and his group [8

N. Huse, A. Schönle, and S.W. Hell, “Z-polarized confocal microscopy,” J. Biomed. Opt. 6, 480–484 (2001). [CrossRef]

], or rigorous ab initio calculations for fluorescence fluctuation spectroscopy [9

J. Enderlein, I. Gregor, D. Patra, T. Dertinger, and U.B. Kaupp, “Performance of Fluorescence Correlation Spectroscopy for Measuring Diffusion and Concentration,” Chem. PhysChem. 6, 2324–2336 (2005). [CrossRef]

] amplify the demand for fast focal 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.

The organization of this paper is as follows: Section 2 introduces the Debye approximation, i.e. the general framework and formulae used in the remainder of this work. Section 3 outlines the implementation based on the fast Fourier transform (FFT) and establishes the sampling and border conditions for obtaining accurate numerical results. Finally, section 4 presents selected examples, firstly the calculation of the EMfield for a 40×1.20NA water immersionmicroscope objective, and secondly, for a 100×1.45 NA oil immersion objective taking into account the evanescent field contribution.

2. The Debye diffraction integral as Fourier transform

This section establishes the basic formalism based on the Debye diffraction integral and the formulation of this integral as a Fourier transform. The basic optical layout and the respective coordinate systems are shown in Fig. 1. We assume that this optical setup, i.e. the imaging system obeys Abbe’s sine condition (as usually fulfilled for microscope objectives).

A coherent, monochromatic wave field parallel to the optical axis crosses the aperture stop A, propagates towards the principal plane ℙ₁ and is transferred to the principal plane ℙ₂. At ℙ₂, the wave field is refracted and focused towards the focal point F ₂. The point P lies on the principal plane ℙ₂ and illustrates the focusing of a ray at ℙ₂ towards the focal point F ₂. The spherical surface ℙ₂ is centered at F ₂ and the deflection angle θ at the position P is given by

Fig. 1. Optical setup. The objective is represented by the aperture stop A with radius R, the principal planes ℙ₁ and ℙ₂ with vertex points V ₁ and V ₂, and the foci F ₁ and F ₂. The focal length f is given as

f = \bar{F_{1} V_{1} .}

. The point P is the intersection point of a ray with ℙ₂ and shows the relation of the position r at ℙ₁ of the incident wave E⃗_i to the propagation angle θ at ℙ₂ of the transmitted wave E⃗_t .

\sin θ = \frac{r}{R} \frac{NA}{n_{t}}

(1)

where r is the off-axis coordinate of the incident wave, R the aperture stop radius, NA the numerical aperture of the objective and n_t the index of refraction behind the ℙ₂ surface. In our setup, the aperture A is placed in the back focal plane, which results in a telecentric imaging system.

Instead of the principal planes, pupils are frequently used for modeling the wave propagation through the objective. However, diffraction at the aperture stop inside the objective is not obvious if the incident wave is transferred directly from the entrance pupil to the exit pupil. Within our representation, the wave propagation from the aperture plane A to the principal plane ℙ₁ is easily calculated with the PWS method or in most cases based on classical Fourier optics principles.

The incident field E⃗_i (r,ϕ) at ℙ₁ is decomposed into a radial component (p-polarized) and a tangential component (s-polarized). The unit vectors for p- and s-polarization are

{\vec{e}}_{p} = (\begin{matrix} \cos ϕ \\ \sin ϕ \\ 0 \end{matrix}) and {\vec{e}}_{s} = (\begin{matrix} - \sin ϕ \\ \cos ϕ \\ 0 \end{matrix})

(2)

where ϕ is the azimuth angle around the z-axis. Upon transmission, the unit vector e⃗_p is deflected by θ and becomes

{\vec{e}}_{r} = (\begin{matrix} \cos ϕ \cos θ \\ \sin ϕ \cos θ \\ \sin θ \end{matrix}) .

(3)

Hence, the amplitude, phase and polarization of the transmitted field at ℙ₂ is

{\vec{E}}_{t} (θ, ϕ) = t_{p} ({\vec{E}}_{i} . {\vec{e}}_{p}) {\vec{e}}_{r} + t_{s} ({\vec{E}}_{i} . {\vec{e}}_{s}) {\vec{e}}_{s}

(4)

where t_p (θ,ϕ) and t_s (θ,ϕ) are the transmission coefficients (viz pupil function, apodization) for p- and s-polarization, respectively. Accumulated phase distortions, i.e. aberrations at ℙ₂, as well as attenuations, i.e. amplitude factors, are integrated in the complex parameters 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⃗ near F ₂. Hence, the electric field E⃗ at a point (x, y, z) is obtained by integrating the propagated plane waves, viz

\vec{E} (x, y, z) = - \frac{if}{λ 0} \iint_{Ω} {\vec{E}}_{t} (θ, ϕ) e^{i (k_{z} z - k_{x} x - k_{y} y)} d Ω

= - \frac{if}{λ 0} \int_{0}^{Θ} \sin θ \int_{0}^{2 π} {\vec{E}}_{t} (θ, ϕ) e^{i (k_{z} z - k_{x} x - k_{y} y)} d ϕ d θ .

(5)

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 ℙ₂ is observed at F ₂, i.e. sinΘ=NA/n_t . The wave vector k⃗_t is given in spherical coordinates θ and ϕ by

{\vec{k}}_{t} (θ, ϕ) = k_{0} n_{t} (\begin{matrix} - \cos ϕ \sin θ \\ - \sin ϕ \sin θ \\ \cos θ \end{matrix}) where k_{0} = \frac{2 π}{λ_{0}} .

(6)

The evaluation of Eq. (5) is usually performed with a direct numerical integration taking into account the coordinate transformations, which results in the Richard-Wolf integral representation [2

E. Wolf, “Electromagnetic diffraction in optical systems, I. An integral representation of the image field,” Proc. R. Soc. London Ser. A 253, 349–357 (1959). [CrossRef]

, 3

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

]. Instead of the common ansatz, a (θ,ϕ)-sampling keeping dΩ=sinθdθdϕ constant is obtained by using cosθ_m =1-mΔΘ with m∊ℕ. For m∊{1…M} and n∊{1…N}, the sampling grid is defined by

θ_{m} = \arccos (1 - m \frac{1 - \sqrt{\frac{1 - {NA}^{2}}{n_{t}^{2}}}}{M}) and ϕ_{n} = (n - \frac{1}{2}) \frac{2 π}{N} .

(7)

At θ=0, a sampling point with a weight of dΩ=

{π θ}_{1}^{2}

/4 is added. 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 [4

Typically, a good accuracy is achieved for M≳50 and N≳200 sampling points.

The outlined evaluation of the Debye diffraction integral (5) is quite fast, but still much slower than the conventional computation of a Fraunhofer diffraction integral. However, Eq. (5) 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 ℙ₁ instead of ℙ₂. Using Eq. (1) and (6), the integration step dΩ for a sampling over ℙ₂ is projected onto ℙ₁, which yields

d Ω = {(\frac{NA}{{Rn}_{t}})}^{2} \frac{r d r d ϕ}{\cos θ} = {(\frac{NA}{{Rn}_{t}})}^{2} \frac{d x d y}{\cos θ} = \frac{1}{k_{t}^{2}} \frac{d k_{x} d k_{y}}{\cos θ} .

(8)

Insertion of this sampling step into Eq. (5) results in

\vec{E} (x, y, z) = - \frac{if}{λ_{0} k_{t}^{2}} \iint_{r < R} (\frac{{\vec{E}}_{t} (θ, ϕ) e^{{ik}_{z^{z}}}}{\cos θ}) e^{- i (k_{x} x + k_{y} y)} {d k}_{x} {d k}_{y} .

(9)

Extending now the integration over (k_x , k_y )∊ℝ² 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

\vec{E} (x, y, z) = - \frac{if}{λ_{0} k_{t}^{2}} 𝓕 (\frac{{\vec{E}}_{t} (θ, ϕ) e^{{ik}_{z^{z}}}}{\cos θ}) .

(10)

This is the main result of this work. The Debye integral is now expressed as a Fourier transform of the field distribution in the aperture A. The 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. (10) is equivalent to the Fraunhofer diffraction integral.

3. Numerical implementation

The numerical implementation is straightforward. A fast Fourier transform (FFT) of the weighted field at ℙ₂ is used for the numerical evaluation of Eq. (10). For an equidistant sampling k_x =mΔK and k_y =nΔK with ΔK=k ₀ NA/M, viz M sampling points over the aperture radius, the sampling points on ℙ₂ are

θ_{mn} = \arcsin (\frac{Δ K}{k_{t}} \sqrt{m^{2} + n^{2}}) and ϕ_{mn} = \arctan (\frac{n}{m}) for ∣ m ∣, ∣ n ∣ \leq M .

(11)

Multiplication of the integration step (ΔK)² with the prefactor of Eq. (10) yields the numerical implementation of Eq. (10) as

\vec{E} (x_{kl}, y_{kl}, z) = - \frac{i R^{2}}{λ_{0} f M^{2}} FFT (e^{{ik}_{{zmn}^{z}}} \frac{{\vec{E}}_{t} (θ_{mn}, ϕ_{mn})}{\cos θ_{mn}}) .

(12)

Typically, the FFT is more than 100× faster than the direct integration of Eq. (5) with matrix multiplication. A good accuracy is achieved for 4M ² ≳ 100×100 sampling points over Ω, but care has to be taken in order to avoid artifacts due to sampling and aliasing. Subsequently, the necessary conditions for obtaining accurate results are investigated [10

For simplification, the sample indices kl and mn will be omitted further on.

3.1. Sampling condition

The propagation factor

e^{{ik}_{z^{z}}}

in Eq. (10) has to be calculated with high resolution for accurate results [11

M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6, 786–805 (1989). [CrossRef]

]. This imposes a condition on the phase discretization, i.e. the phase k_zz must not change by more than π between neighboring sampling points in the aperture plane A. With

k_{z} = \sqrt{k_{t}^{2} - k_{xy}^{2}}

, the sampling condition can be expressed as

max ∣ \frac{d (k_{z} z)}{d k_{xy}} ∣ = max ∣ z \frac{k_{xy}}{k_{z}} ∣ = max ∣ z \tan θ ∣ < \frac{π}{Δ K}

(13)

where ΔK=k ₀ NA/M and

max ∣ \tan θ ∣ = \frac{NA}{\sqrt{n_{t}^{2} - {NA}^{2}} .}

. This immediately leads to a condition for the minimum number of sampling points

M > \frac{{2 NA}^{2}}{\sqrt{n_{t}^{2} - {NA}^{2}}} \frac{∣ z ∣}{λ_{0}},

(14)

solely determined by the system parameters. For the numerical evaluation, an oversampling of about 3× is sufficient for improving the accuracy of the result. In addition, a lower limit of M≳50 reveals necessary for an accurate sampling of ϕ. Deviations from these sampling conditions result in granular artifacts as seen in Fig. 3(a). As a typical value for M, we have chosen M=125 for the focus field calculation of a 1.20 NA water immersion objective (see the example 4.1). A high accuracy is obtained for |z|≲25λ ₀, corresponding to ≈12 µm at a wavelength of 488 nm.

Fig. 2. Two-dimensional fast Fourier transform FFT(E⃗_t /cosθ)=E⃗(x, y, 0) limited to the region of interest (dotted square). Left: Field E⃗t, aperture matrix padded with zeros (dotted rectangle). Center: FFT along the first dimension, cropped and padded with zeros. Right: FFT along the second dimension. The arrows indicate the transformed dimension.

3.2. Sampling step

The focus field E⃗ is obtained for the sampling positions (mΔx,nΔy, z). With Δk=2π/NΔr and Δr=fΔK/k_t , the sampling step in the xy-plane is

Δ x = Δ y = f \frac{Δ k}{k_{t}} = \frac{M}{N} \frac{λ_{0}}{NA}

(15)

where N>4M is the number of FFT sampling points per transformed dimension (see also Fig. 2, where the arrows span over 2M+1 samples and the padded dimension over N samples). For optimal FFT performance, it is best to set N=2 ^s with s∈ℕ. Respecting the condition (14), M can be adjusted to fit Δx and Δy. Along the z-direction, the sampling can be chosen arbitrarily by respecting the limits given above.

3.3. Aliasing suppression

Due to the Debye diffraction integral expressed in Eq. (10), the field E⃗_t is the plane wave spectrum of the focus field E⃗. Usually, the smallest area (aperture matrix) containing E⃗_t ≠0 is transformed (see Fig. 2). The spectral product

e^{{ik}_{z^{z}}} \times \frac{{\vec{E}}_{t}}{\cos θ}

in Eq. (10) represents a spatial convolution

\vec{E} = 𝓕 (e^{{ik}_{z^{z}}}) * 𝓕 (\frac{{\vec{E}}_{t}}{\cos θ})

. In general, the result of the convolution is non-zero on an area larger than the aperture size, which may cause aliasing [12

M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Comm. 116, 43–48 (1995). [CrossRef]

]. Therefore, the aperture matrix is enlarged by zero padding to at least twice its dimensions before performing the transform. In a final step, simple cropping of the transform output removes the padding.

Because we are only interested in the field near the focus, typically over a range of several wavelengths, we limit the computation of the FFT to this region of interest (Fig. 2). The transmitted field E⃗_t is padded with zeros along the first dimension. In this dimension, the FFT is calculated and the result cropped. Along the second dimension, the same procedure is applied on the intermediate result. Zero padding simultaneously suppresses aliasing and refines the sampling grid for the focus field. Using two one-dimensional FFTs with intermediate cropping and zero padding minimizes the numerical processing cost.

3.4. Aperture rim smoothing

Fig. 3. Spectrum (logarithmic scale) with binary sampling of the aperture rim (a), respectively with smoothing as given by Eq. (16) (b). Binary sampling leads to discretization errors at the aperture rim, which results in granular artifacts at high frequencies. Therefore, (a) is only accurate at low frequencies over ≲ 20% of the focal field. In (b) these artifacts are almost suppressed for ≳ 70% of the focal field.

Fig. 4. Comparison of cross-sections through the ’sharp’ and ‘smooth’ focal fields.

Figure 3 shows the spectra log|FFT(U)| for a circular aperture with radius R. As already stated, the field U vanishes outside the aperture for r>R, whereas inside the aperture for r<R, the field is given as U=U ₀. This discretization leads to a serrated aperture rim inducing granular artifacts at higher frequencies. Hence, the expected Airy function is only seen at low frequencies (central region in Fig. 3(a), please note the logarithmic scale). A smooth sampling of the aperture rim improves the accuracy of the spectrum [13

P. Luchini, “Two-dimensional numerical integration using a square mesh,” Comp. Phys. Comm. 31, 303–310 (1984). [CrossRef]

]. In Fig. (Fig. 3(b)) the rim was sampled with the hyperbolic tangent as

U (r) = \frac{1}{2} (1 + \tanh (\frac{1.5}{Δ R} (R - r))) U_{0}

(16)

where ΔR=R/30 was the sampling step. The granular artifacts are efficiently reduced and the FFT approximates the Airy function with a good accuracy over a much larger area. Figure 4 shows a comparison of cross-sections of the spectra on the meridian k_y =0. Overall, for values |k_x |>60×2π/256ΔR, the ‘sharp’ spectrum shows granular artifacts, whereas the ‘smooth’ spectrum approximates well the Airy function.

3.5. Generalization based on the chirp z transform

We demonstrated the importance of zero padding while respecting the sampling condition (14). These constraints led to a minimal number of sampling points N=2^s for the FFT (s∈ℕ). The corresponding number of sampling points M over the aperture radius often exceeds the initial guess based Eq. (14). In such cases, the chirp z transform (CZT) is computationally faster than the FFT. In summary, the CZT (a) allows breaking the relationship between M and N, (b) allows an implicit frequency offset, and (c) internalizes the zero padding. Applying this generalization, we adapted the sampling step in the focus field independently of the sampling step in the input field, introduced an additional shift of the region of interest, and finally improved the computational efficiency.

Let z_m ∀m∈[0,M-1] be a discrete representation of a spatial signal z(r=mΔr). The discrete Fourier transform (DFT) at a frequency k=nΔk∀n∈[0,N-1] is then obtained with

F_{n} = \sum_{m = 0}^{M - 1} z_{m^{e^{- imn Δ k}}} .

(17)

The FFT is a particular case of the DFT with Δk=2π/MΔr and N=M. For Δk<2π/MΔr, a zero padding is implicitly contained in Eq. (17). Comparing the DFT with the CZT defined by

Z_{n} = \sum_{m = 0}^{M - 1} z_{m^{a^{- m_{w^{mn}}}}}

(18)

yields a=1 and w=e ^-iΔk for obtaining the DFT as a particular case of the general CZT. Setting

a = e^{{ik}_{0}}

shifts the frequency domain by k ₀ (see above). Furthermore, Eq. (18) can be rewritten as a convolution

Z_{n} = w^{\frac{n^{2}}{2}} \sum_{m = 0}^{M - 1} z_{m} a^{- m} w^{\frac{m^{2}}{2}} . w^{\frac{- {(n - m)}^{2}}{2} = ((z_{m} a^{- m} w^{\frac{m^{2}}{2}}) * (w^{\frac{- m^{2}}{2}})) w^{\frac{n^{2}}{2}}}

(19)

that can be evaluated using two (M+ N-1) point FFTs (a third one can be precomputed) [14

J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt. 41, 4897–4903 (2002). [CrossRef] [PubMed]

Z = {CZT}_{a, w} (z) = w^{\frac{n^{2}}{2}} {FFT}^{- 1} (FFT (z_{m} a^{- m} w^{\frac{m^{2}}{2}}) . FFT (w^{\frac{- m^{2}}{2}}))

(20)

Based on the CZT, our computationmethod can be extended for low NA systems or for focus fields with a large axial span. In such cases, the sampling grid becomes distorted over the focus depth [15

Y. Li and E. Wolf, “Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A 1, 801–808 (1984). [CrossRef]

, 16

W. Hsu and R. Barakat, “Stratton-Chu vectorial diffraction of electromagnetic fields by apertures with application to small-Fresnel-number systems,” J. Opt. Soc. Am. A 11, 623–629 (1994). [CrossRef]

]. But within the framework of the CZT, this distortion can be compensated by a non-linear scaling proportional to the effective NA under which the aperture A is observed at ℙ₂ from a point (0,0,z) on the axis. As a result, the sampling Δk depends upon the axial position z, i.e. Δk(z)=Δk(0)f/(f+z) with Δk(0)=Δk as defined before. Using the CZT, the additional calculations remain restricted to the repeated computation of

FFT (w^{\frac{- m^{2}}{2}})

because w varies now with z.

4. Selected examples

This section presents example calculations for two high NA microscope objectives. In the first example of a 1.20 NA water immersion objective, the variation of different amplitude distributions (apodization) in the aperture A are discussed. For the second example, a 1.45 NA oil immersion objective was chosen as used in total internal reflection microscopy. The refraction at a cover glass-water interface at the focus is added and the effect of different polarization distributions in the aperture plane A are discussed.

Before presenting these specific examples, the transmission coefficients t_p and t_s between the principal planes ℙ₁ and ℙ₂ need to be defined. We present the microscope objective as an optical system of only 2 optical interfaces and a convex interface into the immersion medium n_t . To this end, the three interfaces provide a physical model for deflection angles θ∈[0,π/2). The amplitude transmission efficiency, i.e. apodization, and the polarization are obtained based on the Fresnel equations.

If the glass lens has an index of refraction n_g and the immersion medium n_t , the Fresnel transmission coefficients are calculated for the succession of the air(n_a )-glass(n_g )-air(n_a )-immersion(n_t ) interfaces. The corresponding deflection angle θ_ij at each interface was chosen proportional to the difference of the index of refraction, viz θ_ij ∝|n_i -n_j |. With n_a =1, the Fresnel transmission coefficients are then

t_{p} = \frac{1}{n_{t}} (1 - {(\frac{2 n_{g} - (n_{g}^{2} + 1) \cos θ_{ag}}{(n_{g}^{2} - 1) \cos θ_{ag}})}^{2}) \frac{2 - 2 n_{t} \cos θ_{at}}{2 n_{t} - (n_{t}^{2} + 1) \cos θ_{at}}

(21)

for p-polarization and

t_{s} = (1 - {(\frac{n_{g}^{2} - 2 n_{g} \cos θ_{ag} + 1}{n_{g}^{2} - 1})}^{2}) \frac{2 - 2 n_{t} \cos θ_{at}}{n_{t}^{2} - 2 n_{t} {\cos θ}_{at} + 1}

(22)

for s-polarization, respectively.

4.1. 1.20 NA water immersion objective

Figure 5 shows the focus intensity for a nearly uniform and a Gaussian illumination in the back aperture of a 1.20 NA water immersion objective. For Δx=Δy=20 nm, Δz=50 nm and M=100, a 2.0 GHz Pentium 4 processor computed the field within a volume of 3 µm×3 µm× 5 µm i.e. 150 ×150 ×100 sampling points in less than 40 seconds. Taking the symmetry into account, the volume was further extended to 6 µm×6 µm×10 µm.

In Fig. 5(a), the aperture was overfilled and the resulting focus field shows the well-known symmetry break of vectorial focus fields, for comparison the Airy profile was added. In Fig. 5(b), the aperture was underfilled to about 60% and the field becomes approximately gaussian. Figure 6 shows the electric fields along two major axes through the focus. For an overfilled aperture, the Airy profile (based on a scalar, paraxial approximation) is a good estimation of the electric field along the y-axis. For an underfilled aperture, the diameter of the central lobe is ≈25% larger but the side lobes vanish quickly. In both cases, the polarization leads to a larger x-extension compared to the y-extension.

Fig. 5. Intensity distribution at the focus of a 1.20 NA water immersion objective for a x-polarized laser beam with a wavelength of λ ₀=488 nm. The aperture had a diameter of 6.5 mm and the e ^-2 beam diameter was 10 mm (a) and 4 mm (b), respectively. The iso-intensity surfaces show the surfaces I(x,y,z)=e ^-1…^-4max(I).

Fig. 6. Electric field profiles along the x- and y-axes, respectively, for the 40×1.20 NA water immersion objective with overfilled and underfilled aperture. The full laser beam power was 1 mW. The Airy profile is given for comparison.

Fig. 7. Cross-sections through the focus intensity distribution of Fig. 5(a). The full laser beam power was 1 mW.

Fig. 8. Cross-sections through the focus intensity distribution of Fig. 5(b). The full laser beam power was 1 mW.

Figure 7 and 8 show the intensity on the major planes through the focus. The polarization dependent extensions of the lobes along the major axes x and y creates a transition zone where the fringe contrast is diminished.

4.2. 1.45 NA oil immersion objective

As a second example, we calculate the focus field of an objective designed for total internal reflection fluorescence (TIRF). The objective uses immersion oil with an index of refraction matching the cover slip. Its NA of 1.45 is higher than the index of refraction of the sample (n_s =1.33, aqueous solution). This generates a partially evanescent focus field at the cover slip–sample interface. Depending upon the illumination of the aperture, the focus field can be fully propagating or fully evanescent. A fully propagating field can be calculated easily with the procedure outlined above. However, the evanescent field contribution needs an additional consideration for obtaining the total focus field.

First we determine the plane wave spectrum E⃗_t at the immersion oil-cover slip interface. Next, the refraction at this interface and the cover slip–sample interface is calculated in order to obtain the plane wave spectrum E⃗_s in the sample (water). Finally, applying the Fourier transform on the weighted and propagated spectrum

\frac{e^{{ik}_{z^{z}}} {\vec{E}}_{s}}{\cos θ}

yields the focus field. As before, the angle θ and the weighting factor 1/cosθ are calculated in the immersion oil. But concerning the sampling condition, a specific issue related to the cover slip–sample interface (14) needs to be considered. The highest angles θ result in total internal reflection at the cover slip–sample interface. At the critical angle θ_c =arcsin(n_s /n_t ), k_z vanishes. For higher angles, k_z takes an imaginary value and the sampling condition (14) is relaxed because

e^{{ik}_{z^{z}}}

becomes just an amplitude factor. The problem arises at θ_c where the sampling condition (13) results in a singularity. Let M′ be the number of sampling points over θ<θ_c . For avoiding this singularity at θ_c , the sampling is chosen such that (M′+1/2)ΔK=k_s , i.e. θ_c falls between two sampling points. Inserting M=(M′+1/2)NA/n_s ,

k_{xy} = k_{s} \frac{M'}{M' + \frac{1}{2}} and k_{z} = k_{s} \frac{\sqrt{M' + \frac{1}{4}}}{M' + \frac{1}{2}}

(23)

into Eq. (13) then yields a generalized sampling condition

M ≳ 4 n_{s} NA \frac{z^{2}}{λ_{0}^{2}} .

(24)

A 7× oversampling is used for improving the accuracy of the result, in particular at off-axis points. In addition, a lower limit of M≳100 was used for |z|→0.

Because the field is calculated in the sample space, k⃗_t is replaced by

{\vec{k}}_{s} (θ, ϕ) = k_{0} n_{s} (\begin{matrix} \cos ϕ \sin θ' \\ \sin ϕ \sin θ' \\ \cos θ' \end{matrix}) = k_{0} n_{t} (\begin{matrix} \cos ϕ \sin θ \\ \sin ϕ \sin θ \\ \sqrt{\frac{n_{s}^{2}}{n_{t}^{2}} - \sin^{2} θ} \end{matrix})

(25)

where n_s sinθ′=n_t sinθ. The unit vector e⃗_r for p-polarization becomes

{\vec{e}}_{r} = (\begin{matrix} \cos ϕ Re (\cos θ') \\ \sin ϕ Re (\cos θ') \\ min (1, \sin θ') \end{matrix}) .

(26)

Figure 9 shows the focus field of a 100×1.45NA oil immersion objective. The aperture of the objective was overfilled, resulting in a partially evanescent field at the focus, where the cover slip–sample (water) interface was placed. As for the former example, the central lobe extends less in the y- than the x-direction for linear polarization (Fig. 9(a)). The focal volume is reduced to about 1/7 compared to the former water immersion objective. Selecting a radially polarized input field results in a rotationally symmetric focus field as shown in Fig. 9(b). On the optical axis, the electric field becomes purely z-polarized. For a distance z≲0.3 µm, this z-component is dominant. Further away from the cover slip–sample interface, the xy-components prevail, which results in an annular field distribution.

Fig. 9. Intensity distribution near the focus of a 1.45 NA oil immersion objective for a laser beam with a wavelength of λ ₀=488 nm. The aperture had a diameter of 5.5 mm and the e ^-2 beam diameter was 10 mm. The iso-intensity surfaces show the surfaces I(x,y,z)=e ^-1…^-4 I ₍₀₎ in the sample space.

Fig. 10. Cross-sections through the focus intensity distribution of Fig. 9(a). The full laser beam power was 1 mW.

The fine structure at the interface is due to the evanescent wave contribution with incidence angles above the critical angle. For instance, Fig. 12 shows the weighted field E_s /cosθ for the linear polarization. At the critical angle (NA=1.33), the field amplitude triples in the x-direction and doubles in the y-direction, respectively, hence marking the abrupt transition from propagating to evanescent fields.

Fig. 11. Cross-sections through the focus intensity distribution of Fig. 9(b). The full laser beam power was 1 mW.

Fig. 12. Field E_s /cosθ of Fig. 9(a). For NA<1.33, the field corresponds to a free propagation in the sample space, whereas for NA>1.33 an evanescent field is induced.

5. Conclusions

We showed a fast and simple 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 transformis extremely fast and allows a high flexibility of the input field. The result is accurate under the conditions for the validity of the Debye integral representation of focused fields [17

E. Wolf and Y. Li, “Conditions for the validity of the Debye integral representation of focused fields,” Opt. Comm. 39, 205–210 (1981). [CrossRef]

, 18

P. Török, “Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number,” J. Opt. Soc. Am. A 15, 3009–3015 (1998). [CrossRef]

] and the given sampling conditions. For low NA, it converges quite naturally to a focus field given by the Fraunhofer approximation. With the chirp z transform, we extended our calculations to low NA focus fields requesting a non-linear scaling as shown by Li and Hsu [15

Y. Li and E. Wolf, “Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A 1, 801–808 (1984). [CrossRef]

, 16

]. Table 1 summarizes the performance of the different calculation methods on a personal computer.

Table 1. Performance of different calculation methods.

Method	Input fields, constraints	Integration	Output	Computation time (for 100³ points)
Classic	Analytic functions (rotational symmetry)	Quadrature of Bessel functions	Points	20 min to hours
Direct	Any, high NA (polar sampling)	Matrix product	Lines	≈30 min
FFT	Any, high NA (carteesian sampling)	FFT	xy planes	≈1 min
CZT	Any (carteesian sampling)	CZT	xy planes	≈30 s

In addition, we used a generalized pupil function (apodization) of high NA objectives taking into account amplitude and polarization distributions. The pupil function incorporates wave front aberrations as contained in real objectives as well as Fresnel transmission coefficients. Based on these Fresnel coefficients, it is straightforward to include wave propagation through stratified media.

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.

Acknowledgements

We are grateful to Herbert Gross, Carl Zeiss Oberkochen for many valuable comments and discussions. The support of the Swiss National Science Foundation (SNSF) (contract number 200021-103333) is greatly acknowledged.

References and links

1.	P. Debye, “Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie,” Ann. Phys. 30, 755–776 (1909). [CrossRef]
2.	E. Wolf, “Electromagnetic diffraction in optical systems, I. An integral representation of the image field,” Proc. R. Soc. London Ser. A 253, 349–357 (1959). [CrossRef]
3.	B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system,” Proc. R. Soc. London Ser. A 253, 358–379 (1959). [CrossRef]
4.	Typically, a good accuracy is achieved for M≳50 and N≳200 sampling points.
5.	P. Török and P. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt. 36, 2305–2312 (1997). [CrossRef] [PubMed]
6.	J.J. Stamnes, Waves in Focal Regions: propagation, diffraction and focusing of light, sound and water waves , Hilger, Bristol UK (1986).
7.	G. Mikula, A. Kolodziejczyk, M. Makowski, C. Prokopowicz, and M. Sypek, “Diffractive elements for imaging with extended depth of focus,” Opt. Eng. 44, 058001 (2005). [CrossRef]
8.	N. Huse, A. Schönle, and S.W. Hell, “Z-polarized confocal microscopy,” J. Biomed. Opt. 6, 480–484 (2001). [CrossRef]
9.	J. Enderlein, I. Gregor, D. Patra, T. Dertinger, and U.B. Kaupp, “Performance of Fluorescence Correlation Spectroscopy for Measuring Diffusion and Concentration,” Chem. PhysChem. 6, 2324–2336 (2005). [CrossRef]
10.	For simplification, the sample indices kl and mn will be omitted further on.
11.	M. Mansuripur, “Certain computational aspects of vector diffraction problems,” J. Opt. Soc. Am. A 6, 786–805 (1989). [CrossRef]
12.	M. Sypek, “Light propagation in the Fresnel region. New numerical approach,” Opt. Comm. 116, 43–48 (1995). [CrossRef]
13.	P. Luchini, “Two-dimensional numerical integration using a square mesh,” Comp. Phys. Comm. 31, 303–310 (1984). [CrossRef]
14.	J. L. Bakx, “Efficient computation of optical disk readout by use of the chirp z transform,” Appl. Opt. 41, 4897–4903 (2002). [CrossRef] [PubMed]
15.	Y. Li and E. Wolf, “Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A 1, 801–808 (1984). [CrossRef]
16.	W. Hsu and R. Barakat, “Stratton-Chu vectorial diffraction of electromagnetic fields by apertures with application to small-Fresnel-number systems,” J. Opt. Soc. Am. A 11, 623–629 (1994). [CrossRef]
17.	E. Wolf and Y. Li, “Conditions for the validity of the Debye integral representation of focused fields,” Opt. Comm. 39, 205–210 (1981). [CrossRef]
18.	P. Török, “Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number,” J. Opt. Soc. Am. A 15, 3009–3015 (1998). [CrossRef]

OCIS Codes
(070.2580) Fourier optics and signal processing : Paraxial wave optics
(180.0180) Microscopy : Microscopy
(220.2560) Optical design and fabrication : Propagating methods
(260.1960) Physical optics : Diffraction theory

ToC Category:
Optical Design and Fabrication

History
Original Manuscript: October 3, 2006
Revised Manuscript: October 27, 2006
Manuscript Accepted: October 28, 2006
Published: November 13, 2006

Virtual Issues
Vol. 1, Iss. 12 Virtual Journal for Biomedical Optics

Citation
Marcel Leutenegger, Ramachandra Rao, Rainer A. Leitgeb, and Theo Lasser, "Fast focus field calculations," Opt. Express 14, 11277-11291 (2006)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-23-11277

Sort: Author | Year | Journal | Reset

References

P. Debye, "Das Verhalten von Lichtwellen in der Nähe eines Brennpunktes oder einer Brennlinie," Ann. Phys. 30,755-776 (1909). [CrossRef]
E. Wolf, "Electromagnetic diffraction in optical systems, I. An integral representation of the image field," Proc. R. Soc. London Ser. A 253,349-357 (1959). [CrossRef]
B. Richards and E. Wolf, "Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system," Proc. R. Soc. London Ser. A 253,358-379 (1959). [CrossRef]
Typically, a good accuracy is achieved for M & 50 and N & 200 sampling points.
P. Török and P. Varga, "Electromagnetic diffraction of light focused through a stratified medium," Appl. Opt. 36,2305-2312 (1997). [CrossRef] [PubMed]
J. J. Stamnes, Waves in Focal Regions: propagation, diffraction and focusing of light, sound and water waves, (Hilger, Bristol UK 1986).
G. Mikula, A. Kolodziejczyk, M. Makowski, C. Prokopowicz, M. Sypek, "Diffractive elements for imaging with extended depth of focus," Opt. Eng. 44,058001 (2005). [CrossRef]
N. Huse, A. Schönle, and S. W. Hell, "Z-polarized confocal microscopy," J. Biomed. Opt. 6,480-484 (2001). [CrossRef]
J. Enderlein, I. Gregor, D. Patra, T. Dertinger, and U.B. Kaupp, "Performance of Fluorescence correlation Spectroscopy for measuring diffusion and concentration," Chem. Phys. Chem. 6,2324-2336 (2005). [CrossRef]
For simplification, the sample indices kl and mn will be omitted further on.
M. Mansuripur, "Certain computational aspects of vector diffraction problems," J. Opt. Soc. Am. A 6,786-805 (1989). [CrossRef]
M. Sypek, "Light propagation in the Fresnel region. New numerical approach," Opt. Commun. 116,43-48 (1995). [CrossRef]
P. Luchini, "Two-dimensional numerical integration using a square mesh," Comput. Phys. Commun. 31,303-310 (1984). [CrossRef]
J. L. Bakx, "Efficient computation of optical disk readout by use of the chirp z transform," Appl. Opt. 41,4897-4903 (2002). [CrossRef] [PubMed]
Y. Li and E. Wolf, "Three-dimensional intensity distribution near the focus in systems of different Fresnel numbers," J. Opt. Soc. Am. A 1,801-808 (1984). [CrossRef]
W. Hsu and R. Barakat, "Stratton-Chu vectorial diffraction of electromagnetic fields by apertures with application to small-Fresnel-number systems," J. Opt. Soc. Am. A 11,623-629 (1994). [CrossRef]
E. Wolf and Y. Li, "Conditions for the validity of the Debye integral representation of focused fields," Opt. Commun. 39,205-210 (1981). [CrossRef]
P. Török, "Focusing of electromagnetic waves through a dielectric interface by lenses of finite Fresnel number," J. Opt. Soc. Am. A 15,3009-3015 (1998). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper