Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

Chapter

4 Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy Michele Marrocco

Contents

1. Introduction 2. Description of the Problem 2.1. Linear Polarization 2.2. Other States of Polarization 3. Methods of Approximations for Linearly Polarized Laser Beams 3.1. Expansion of the Phasor 3.2. Expansion of Apodization and Pupil Functions 3.3. Multipole Expansion 3.4. Methods of Local Phase and Amplitude Approximation 3.5. Method Based on Transition from Continuous to Discrete Diffraction 3.6. Comparisons Among the Methods for Linearly Polarized Optical Beams 4. Methods of Approximations for Radially Polarized Laser Beams 4.1. Method Based on the Lax Series 4.2. Method Based on Complex-Source-Point Spherical Wave 4.3. Method of Transverse Electric and Transverse Magnetic Decomposition 4.4. Methods of Eigenfunction Representation, Local Phase and Amplitude Approximation, and Discrete Diffraction

132 134 136 139 141 142 145 146 148 150 152 157 159 160 161 162

Italian National Agency for New Technologies, Energy and Sustainable Economic Development (ENEA), via Anguillarese 301, 00123 Santa Maria di Galeria, Rome, Italy Advances in Imaging and Electron Physics, Volume 165, ISSN 1076-5670, DOI: 10.1016/B978-0-12-385861-0.00004-X. c 2011 Elsevier Inc. All rights reserved. Copyright

131

132

Michele Marrocco

4.5. Comparisons Among the Methods for Radially Polarized Optical Beams 5. Conclusions References

164 167 169

1. INTRODUCTION Undoubtedly, one of the most recent scientific challenges in optical imaging is the achievement of very high spatial resolutions that might reach the so-called nanoscopic scale and ultimately go beyond the diffraction limit (Betzig and Trautman, 1992; Domke and Pettinger, 2010; Durant et al., 2006; Dyba and Hell, 2002; Gramotnev and Bozhevolnyi, 2010; Hayazawa et al., 2002; Novotny and Hecht, 2006; Westphal et al., 2008; and many others not listed). The importance of such a challenge is clear in the disciplines of chemical, biological, and physical interest where imaging of details comparable to or even smaller than the optical wavelengths is crucial to the understanding of the phenomena occurring at that spatial level and, in this regard, optical systems combining laser techniques with microscopy realized with lenses of high numerical aperture (NA) can successfully meet the challenge (Cheng et al., 2002; Denk et al., 1990; Dorn et al., 2003; Dyba and Hell, 2002; Freudiger et al., 2008; Kino and Corle, 1997; Masters and So, 2008; Novotny and Hecht, 2006; Pawley, 2006; Van de Nes et al., 2006; Volkmer et al., 2001; Westphal et al., 2008; Zipfel et al., 2003; Zumbusch et al., 1999). To this end, the theoretical effort needed to define the optical problems inherent in applications of extremely focused laser fields is rooted in the electromagnetic diffraction theory (EDT), as developed by Richards and Wolf (1959). In general, EDT provides the basic framework of the vectorial problem of strongly focused optical beams (i.e., not necessarily coming from laser systems) and, for this reason, it has become commonplace in microscopy (Cheng et al., 2002; Denk et al., 1990; Dorn et al., 2003; Freudiger et al., 2008; Kino and Corle, 1997; Masters and So, 2008; Novotny and Hecht, 2006; Pawley, 2006; Van de Nes et al., 2006; Volkmer et al., 2001; Zipfel et al., 2003; Zumbusch et al., 1999), where the spatial structure of the diffracted electromagnetic field seen by the sample placed at the focus of the microscope is essential to the understanding of the interaction between light and matter. Since its importance, EDT provides the main context of the numerical methods reviewed in the first part of the present paper. Alternatively, attempts at the description of strongly focused fields (that are, of course, beyond the paraxial limit) have been made by using the Maxwell equations reduced to the Helmholtz equation for the vector potential. However, these approaches suffer from important limitations and thus are of minor significance compared with the EDT compared with the numerical tools

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

133

developed under the aegis of EDT theory. This aspect is considered in the second part of the review. It is argued that the recourse to EDT is appropriate for NA values greater than 0.5 (Hopkins, 1943; Jipson and Williams, 1983). Below this limit, the scalar theory based on Kirchhoff’s approach is instead preferred (Born and Wolf, 1970). In the former context, on the other hand, the effect of polarization (hence the vectorial character) cannot be neglected and many instructive examples can be found in several branches of spectroscopy and nano-optics, where three-dimensional (3D) spatial resolution and signal collection are pushed to the extreme (Cheng et al., 2002; Denk et al., 1990; Dorn et al., 2003; Freudiger et al., 2008; Kino and Corle, 1997; Masters and So, 2008; Pawley, 2006; Van de Nes et al., 2006; Volkmer et al., 2001; Zipfel et al., 2003; Zumbusch et al., 1999). To restrict the extensive list, it is sufficient to recall some paradigmatic applications in nanoscience such as high-density data storage (Van de Nes et al., 2006) or physical and chemical uses of confocal microscopy (Kino and Corle, 1997; Pawley, 2006), multiphoton fluorescence (Denk et al., 1990; Zipfel et al., 2003), and coherent Raman scattering (Cheng et al., 2002; Freudiger et al., 2008; Volkmer et al., 2001; Zumbusch et al., 1999). In addition to these “classic” examples, EDT is fundamental to novel concepts in optics, such as microscopy with radially polarized laser beams (Dorn et al., 2003), spatial resolution beyond the diffraction barrier (Novotny and Hecht, 2006), or perfect reflection of light by a single oscillating dipole (Dyba and Hell, 2002; Zumofen et al., 2008). In this chapter, we explore the methods aimed at reducing the complexity of the numerical calculation necessary to identify the physical and chemical phenomena explored by means of optical techniques falling within the EDT realm. The main focus is on the role of the so-called diffraction integrals that account for the fundamental dependences on the radial and axial coordinates ρ and z (Novotny and Hecht, 2006; Richards and Wolf, 1959). Such dependences play important roles in the correct interpretation of the light-matter interaction and, unfortunately, the continuous distribution of field modes (allowed by the diffraction process) complicates the analysis to such an extent that for each field point, numerical evaluation of the diffraction integrals is required (the sentence in italics is the quotation extracted from the book by Novotny and Hecht, 2006). This means that precise reconstruction of the focal volume is subjected to numerous calculations needed to cover the significant spatial ranges of the two coordinates ρ and z. This problem is bearable in linear laser microscopy when a few spatial points are sufficient to model the physical or chemical processes. However, as soon as many field points are considered, the calculation understandably becomes lengthier and fast methods of calculation may prove useful in imaging where two-dimensional (2D) or 3D maps could easily count millions of nodes or pixels. In nonlinear

134

Michele Marrocco

laser microscopy (also known as multiphoton laser microscopy), the problem is actually much worse because the measured signal stems from the 3D spatial integration over the volume within which two or more laser fields are simultaneously acting. In other words, the direct numerical simulation (DNS) based on the Debye–Wolf approach is no longer convenient in that the calculation speed decreases significantly according to the nonlinearity involved in the interaction. As an example, DNS elaboration of images in third-order laser microscopy can last days versus a few minutes with approximated methods (Marrocco, 2009a). It is then not surprising that, besides the direct solution of the vectorial problem, researchers have developed various strategies to simplify the entire calculation to realize substantial gains in time and speed. Here we attempt to provide an account of several efforts devoted to this challenge. The review is organized as follows. First, the description of the problem is necessary to give an overview of the general context of the EDT theory with its final results condensed in the expressions known as Debye–Wolf diffraction integrals (Section 2). Having settled the main frame of reference, we begin by reviewing the numerous techniques that have been used in place of DNS calculations to solve EDT for the special case of linear polarization of the fields delivered to the entrance of the microscope (Section 3). The increasingly important case of radial polarization, which is progressively attracting the interest of the optical community, is also considered (Section 4). In addition to EDT models from the simpler case of linear polarization and adapted for the new context of radial polarization, analytical solutions to the Helmholtz equation for the vector potential are examined even though they are constrained by a lack of accuracy at very high NA. Proper comparisons among the methods conclude the discussions of the two different states of polarizations. Finally, we draw some conclusions (Section 5).

2. DESCRIPTION OF THE PROBLEM The scalar version of EDT dates back to Debye (1960), who pictured the focusing through an aperture as the resultant composition of an angular spectrum of plane waves that propagates from points in the aperture toward the focal region. This concept was assumed by Ignatovsky (1919, 1920), who attempted to solve the vectorial problem that later was better formalized by Wolf (1959) and Richards and Wolf (1959) for an uniform illumination of the aperture. Yoshida and Asakura (1975) instead generalized the approach to a Gaussian distribution of the light emerging from the aperture. In this review, we omit the details of the derivation of the Debye–Wolf integral but readers are encouraged to follow the simple reasoning in the book by Novotny and Hecht (2006).

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

135

Lens x Fo

c

p al

lan

θ

e

ϕ r

z

y

FIGURE 1 Geometry of the optical system.

Considering a homogeneous space whose points are given by the displacement vector r whose coordinates are x, y, z and defining the focal plane by z = 0 (Figure 1), the seminal Debye idea of angular spectrum representation of optical fields leads to E(x, y, z) =

ir exp(−ikr) 2π

ZZ



 1 dkx dky , E∞ kx , ky exp i kx x + ky y ± kz z kz

kx2 +ky2 ≤k2

(1)
1/2 where r = x
2 + y2 + z2 , k is the absolute value of the wave vector k = kx , ky , kz , z is the propagation direction, and E∞ (kx , ky ) is the field in the far zone. Note that whenever kz is approximated by the wave vector amplitude k (or kx , ky << kz , paraxial limit), then E∞ and E are related by a Fourier transform at z = 0 (Fourier optics). Equation (1) is known as the Debye–Wolf integral and its application to an aplanatic lens of focal length f with cylindrical symmetry leads to the result ikf exp(−ikf ) E(ρ, ϕ, z) = 2π

θZmaxZ2π

E∞ (θ , φ) exp{ik [ρ sin θ cos(φ − ϕ)

0

0

+ z cos θ ]} sin θdθdφ

(2)

136

Michele Marrocco

where the transverse coordinates x and y of the field point are now x = ρ cos ϕ and y = ρ sin ϕ, and the integration is over the surface of the aperture whose NA is defined by NA = n sin θmax with n the index of refraction of the medium surrounding the lens and θmax limited between zero and the maximum aperture at π/2. The meaning of Eq. (2), or alternatively Eq. (1), is simple: The vectorial properties of the field in the focal region are related to the field in the far zone through the rapidly oscillating phasor (which is often the main cause of numerical concern). In addition, the relationship suggests that it is possible to engineer the spatial structure of the focal field by adequate manipulation of the field in the far zone (Adamson, 2004; Foreman et al., 2008; Kant, 2000; and references therein). Furthermore, the idea lies at the core of the inverse scattering problem where molecular properties can be extracted from the experimental images if proper knowledge of the focused field can be somehow acquired (Novotny et al., 2001; Sick et al., 2000).

2.1. Linear Polarization The typical application of Eq. (2) is for incident fields with linear polarization. The condition emerges naturally by virtue of the most common laser emission, which indeed is characterized by linearly polarized beams. As a consequence, it is important to direct attention to this special case that is recurrent among the techniques currently used in laser microscopy (Cheng et al., 2002; Denk et al., 1990; Dyba and Hell, 2002; Freudiger et al., 2008; Kino and Corle, 1997; Masters and So, 2008; Novotny and Hecht, 2006; Pawley, 2006; Van de Nes et al., 2006; Volkmer et al., 2001; Westphal et al., 2008; Zipfel et al., 2003; Zumbusch et al., 1999). Let us take the more general example of the different index of refraction between the spaces separated by the aplanatic lens of Figure 1, which is assumed to operate under ideal conditions (i.e., no influence of the aberrations). The focused electric field of the lowest Hermite–Gaussian mode (0,0) created by the lens of focal length f is thus described as follows (Novotny and Hecht, 2006; Richards and Wolf, 1959):  I00 + I02 cos 2ϕ ik f (n1 /n2 ) E0 exp(−ik f )  I02 sin 2ϕ  , E(0,0) (ρ, ϕ, z) = 2 −2iI01 cos ϕ 1/2



(3)

where k is the wave vector amplitude, n1 and n2 are the indexes of refraction of the medium before and after the lens, E0 is the incident electric field amplitude, ϕ is the angle of the cylindrical coordinates of E, and finally the

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

137

terms I01 , I02 , and I03 denote the diffraction integrals, given by θZmax

finc (θ )g0m (θ )Jm (kρ sin θ ) exp(ikz cos θ )dθ

I0m =

(4)

0

with m = 0, 1, 2. In Eq. (4), θ is the variable angular aperture with maximum value θmax = sin−1 (NA/n2 ), g0m (θ ) equals (cos θ )1/2 sin θ (1 + cos θ ), (cos θ )1/2 sin2 θ , or (cos θ )1/2 sin θ (1 − cos θ ) for m = 0, 1, 2, and Jm (kρ sin θ ) indicates the Bessel function of the first kind and mth order depending on the radial coordinate ρ. The quantity finc (θ ) appearing in I0m is associated with the spatial filtering by the lens pupil and is often called the apodization function. In the case of Gaussian
beams, the apodization function is given as finc (θ ) = exp −f 2 sin2 θ/w20 with w0 the laser beam waist before
the lens pupil (or finc (θ ) = exp −β 2 sin2 θ/ sin2 θmax if β = fsinθmax /w0 defines the filling factor of the lens pupil). The dependences of Eq. (3) on the diffraction integrals I0m are not only peculiar to electric fields, but magnetic components are also known to behave similarly and it is clear that what is developed here for electric fields could be extended to problems where magnetic fields are invoked (Novotny and Hecht, 2006; Richards and Wolf, 1959). The behaviors of the real and imaginary parts of I00 , I01 , and I02 are shown in Figure 2 for a lens of NA = 1.4 and n2 = 1.5 that focuses a Gaussian laser beam whose waist matches the aperture of the microscope (i.e., β = 1). In comparing the graphs, I00 displays the largest values in dependence on the spatial coordinates. This means that the vectorial structure of E(0,0) (ρ, ϕ, z) in Eq. (3) will be mainly affected by I00 and hence dominated by the component along the direction of the polarization of the incident field. It is then clear that the other components will contribute with much less importance; this explains why the case of linear polarization is often treated as if it were a scalar problem (see, for example, the theoretical assumptions of Cheng et al., 2002; Volkmer et al., 2001; and Zumbusch et al., 1999 regarding coherent anti-Stokes Raman scattering microscopy). The complexity of the diffraction integrals in Eq. (4) is self-evident and research on analytical solutions is surely doomed to failure. The common solution lies in numerical methods that, beyond their clumsiness, are costly in terms of computation times, especially if we recall that optical signals in microscopy are often the resultant of spatial integrals of fields interacting with specific samples of different geometries. This is particularly true in biochemical applications where the target molecule fills the volume of biological specimens placed at the microscope focus (Cheng et al., 2002; Denk et al., 1990; Kino and Corle, 1997; Masters and So, 2008;

138

Michele Marrocco

0.6 −0.6

10 kz

0 kρ

0 10 (a)

(b)

(c)

(d)

(e)

(f)

FIGURE 2 Comparison among the dependences of the real (Re) and imaginary (Im) parts of I00 , I01 , and I02 for a microscope of NA = 1.4 and n2 = 1.5 illuminated by a Gaussian laser beam whose waist is such that β = 1. In particular, (a) ReI00 , (b) ImI00 , (c) ReI01 , (d) ImI01 , (e) ReI02 , and (f) ImI02 . The spatial range is 0 ≤ kρ ≤ 10 and 0 ≤ kz ≤ 10 in each graph. The vertical axis is identical and between −0.6 and 0.6. These limits are indicated for the plot of ReI00 only and it is demonstrated that the most relevant diffraction integral for the case of linearly polarized laser beams is I00 .

Pawley, 2006; Volkmer et al., 2001; Zipfel et al., 2003; Zumbusch et al., 1999). Furthermore, the spectroscopic techniques based on nonlinear optical phenomena introduce the additional complication of multiple laser fields interacting with the matter at the focus. This means that the calculation speed rescales with the nonlinearity of the problem and methods to simplify the simulation are welcome.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

139

For the sake of generality, diffraction integrals of the type in Eq. (4) also appear in the presence of other profiles of incident laser beams. For instance, the other lowest Hermite–Gaussian modes (1,0) and (0,1) are responsible for fields of the following kind (Novotny and Hecht, 2006):   iI cos ϕ + iI14 cos 3ϕ ik f 2 (n1 /n2 )1/2 E0 exp(−ikf )  11 −iI12 sin ϕ + iI14 sin 3ϕ  , E(1,0) (ρ, ϕ, z) = 2w0 −2I + 2I cos 2ϕ 10

13

(5)   i(I + I12 ) sin ϕ + iI14 sin 3ϕ ik f 2 (n1 /n2 )1/2 E0 exp(−ik f ) 11 −iI12 cos ϕ − iI14 cos 3ϕ , E(0,1) (ρ, ϕ, z) = 2w0 2I sin 2ϕ 13

(6) with the new diffraction integrals given by the expression in Eq. (4) except for the function g0m (θ ), which is now replaced by the following functions: g10 (θ ) = (cos θ )1/2 sin3 θ , g11 (θ ) = (cos θ )1/2 sin2 θ (1 + 3 cos θ ), g12 (θ ) = (cos θ )1/2 sin2 θ (1 − cos θ ), g13 (θ ) = g10 (θ ), and g14 (θ ) = g12 (θ ).

2.2. Other States of Polarization Beyond the ordinary use of Eq. (2) for linearly polarized fields, other states of polarization are increasingly attracting the interest of the community (Bomzon et al., 2002; Dorn et al., 2003; Lerman and Levy, 2008; Novotny and Hecht, 2006; Quabis et al., 2000; Shoham et al., 2006; Yew and Sheppard, 2007a; Youngworth and Brown, 2000). For instance, radially polarized laser beams have been recently introduced in optics as one of the most ingenious ways to overcome the diffraction limit that characterizes ordinary laser microscopy (Dorn et al., 2003; Lerman and Levy, 2008; Quabis et al., 2000; Yew and Sheppard, 2007a). However, azimuthally polarized laser beams also pose interesting questions regarding their injection into microscopes with high NA (Bomzon et al., 2002; Novotny and Hecht, 2006; Shoham et al., 2006; Youngworth and Brown, 2000). The subject can be further enriched by research connected to spatial profiles that differ from the mere Gaussian distribution. Examples are the so-called doughnut modes that can be created by combining the Hermite–Gaussian modes described previously (Novotny and Hecht, 2006) and, in general, the associated fields present dependences that are still related to the contribution of the diffraction integrals as for linear polarization. The reference work for the definition of radially and azimuthally polarized tightly focused laser beams was written by Youngworth and Brown

140

Michele Marrocco

(2000). In particular, considering a radially polarized beam that passes through an annular pupil with angular apertures between α1 and α2 , the focused field with wave vector k can be reconstructed through the projections onto the radial plane (x, y) and the axial direction z according to Eρ (ρ, z) = A

Zα2

finc (θ ) cos1/2 θ sin 2θ J1 (kρ sin θ ) exp(ikz cos θ )dθ

(7)

α1

Ez (ρ, z) = 2iA

Zα2

finc (θ ) cos1/2 θ sin2 θ J0 (kρ sin θ ) exp(ikz cos θ )dθ ,

(8)

α1

where k is the wave vector amplitude, ρ = (x2 + y2 )1/2 is the radial coordinate, A is the field amplitude, and finc (θ ) indicates the pupil function as discussed previously for linearly polarized laser beams. The physical meaning of Eqs. (7) and (8) is manifest: The field components Eρ and Ez are the result of the interplay between the propagation function (i.e., the plane wave with the phase term) and the corresponding radial function (i.e., the Bessel function) that describes the spreading of the light within planes lying orthogonal to the optical axis (or axial direction coincident with the z-axis). The product of these two functions is modulated by the spatial structure of the light beam at the pupil, plus the residual effect of the focusing lens mediated by the sine and cosine functions. It is finally understood that the vectorial structure of the field is reconstructed through the relationship E(ρ, z) = Eρ (ρ, z)vˆ ρ + Ez (ρ, z)vˆ z with vˆ ρ and vˆ z unit vectors belonging, respectively, to the radial plane and the axial direction. An example of the behavior of diffraction integrals appearing in Eqs. (7) and (8) is shown in plots (a)–(d) of Figure 3. Here, the real and imaginary parts of Eρ and Ez are shown for the same optical conditions as in Figure 2 and it is seen that radially polarized beams create a strong contribution of the longitudinal component (i.e., Ez ), which lies along the propagation direction of the beam itself. This is remarkably different from the plots of Figure 2 and constitutes one of the main properties of radially polarized beams. For azimuthally polarized beams, the focused fields instead become (Youngworth and Brown, 2000) E8 (ρ, z) = 2A

Zα2

finc (θ ) cos1/2 θ sin θ J1 (kρ sin θ ) exp(ikz cos θ )dθ ,

(9)

α1

An example of the real and imaginary components of E8 is given in plots (e) and (f) of Figure 3. These were obtained under the same optical conditions of plots (a)–(d) drawn for radially polarized beams.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

141

0.4 −0.4

10 0

kρ (a)

10

0

kz (b)

(c)

(d)

(e)

(f)

FIGURE 3 Comparison among the dependences of the real (Re) and imaginary (Im) parts of the diffraction integrals in Eρ , Ez , and E8 given in Eqs. (7), (8), and (9). The optical parameters of the microscope are those of Figure 2. The angular extension is chosen such that α1 = 0 and α2 = θmax . The plots specifies the dependences for (a) ReEρ , (b) ImEρ , (c) ReEz , (d) ImEz , (e) ReE8 , and (f) ImE8 . The spatial range is 0 ≤ kρ ≤ 10 and 0 ≤ kz ≤ 10 in each graph. The vertical axis is identical and between −0.4 and 0.4. These limits are indicated for the plot of ReEρ only.

In complete analogy with the equations characterizing linear polarization, it is immediately concluded that Eqs. (7)–(9) can only be treated numerically. This means that precise reconstruction of the focal region is possible only if a detailed map is made of the spatial dependences on the two coordinates ρ and z. In other words, an accurate description of the focused fields is possible by means of a 2D grid where the spatial nodes are sufficiently close to each other and each node is associated with numerical solutions of the diffraction integrals seen in Eqs. (7)–(9). For this reason, the entire commentary in Section 2.1 about the advantages of computing routines faster than DNS holds true for other states of polarization of the incident fields.

3. METHODS OF APPROXIMATIONS FOR LINEARLY POLARIZED LASER BEAMS Having outlined the importance of the diffraction integrals with reference to the description of the fields at the focus of a high-NA microscope, we

142

Michele Marrocco

now turn to the central argument of our review—namely, the numerical methods that have been conceived over the years to shorten the calculation of diffraction integrals or, in turn, focused fields. Understandably, we start with the methods developed to ensure detailed handling of linearly polarized laser beams. Before doing this, let us clarify a very general comment. The methods presented in the following text are not exact in that approximated solutions are used. However, the approximations are studied in such a way that they introduce specific simplifications without meaningful loss of numerical accuracy, the evaluation of which takes as reference DNS calculation. This is the declared intention of the authors of the proposals reported in the following and we verify that only some attempts are successful, whereas others are very difficult treatments or, even worse, unproductive. This neat response will be reached in Section 3.6 after appropriate and instructive comparisons among the methods. As for linearly polarized lasers, many are the strategies reported in the literature. They are alternatively concentrated on various kinds of ¨ ok, ¨ expansions: phasor . . . (Agrawal and Pattanayak, 1979; Sherif and Tor ¨ ok ¨ et al., 1997), apodization and pupil func2005; Sherif et al., 2008; Tor ¨ ok, ¨ tions (Kant, 1993), multipole (Asatryan et al., 2004; Sheppard and Tor 1997), whole kernel of the diffraction integrals (Marrocco, 2009a), but other authors have suggested simplifications in the functional dependences appearing in the kernel of the integrals (Gravelsaeter and Stamnes, 1982; Hopkins, 1957; Hopkins and Yzuel, 1970; Stamnes et al., 1983; Stamnes, 1986). We explain the different methods and compare them under several interesting conditions pertaining to laser microscopy with high NAs.

3.1. Expansion of the Phasor Among the methods developed to simplify the calculation, the most recurrent approach is focused on the expansion of the phasor. The objective of this approach is to remove the complications arising from the rapidly oscillating term exp(i k z cos θ ) in order to reduce the integral to a series of integrals of purely real functions. The idea was initially explored by Agrawal and Pattanayak (A-P; 1979), who used the following expansion: exp(ikz cos θ ) = (1)



π kz 2

1/2 X ∞  s=0

kz 2

s

sin2s θ (1) Hs−1/2 (kz), s!

(10)

where Hs−1/2 (kz) represents the Hankel function of the first kind and order s − 1/2 (Abramowitz and Stegun, 1965). In this way, the diffraction

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

143

integrals of Eq. (4) become  I0m =

πkz 2

1/2 X ∞ 

kz 2

s

1 (1) (kz)Qs,m (ρ), H s! s−1/2

(11)

finc (θ )gm (θ )Jm (kρ sin θ ) sin2s θ dθ .

(12)

s=0

where the function Qs,m (ρ) is given by θZmax

Qs,m (ρ) =

0

The difficulties in the calculation of the diffraction integral have now been transferred to the calculation of a series of terms containing integrals of real functions. Despite the simplification, it must be noted that the calculation of Eq. (11) is not faster than the DNS evaluation of Eq. (4) (as a matter of fact, methods of numerical integration revolve around summations). This is explained in greater detail in Section 3.6; here, it suffices to observe that the series in the A-P approach must be stopped at a maximum index smax that guarantees sufficient accuracy and, at the same time, saves time otherwise spent for the calculation of unnecessary (i.e., negligible) terms of Eq. (11). After the attempt by Agrawal and Pattanayak, other authors insisted on ¨ ok ¨ et al. (1997), inspired by the A-P expanding the phasor. In particular, Tor method, tried alternative expansions. One is from Watson (1952), where the phasor is written as exp(ikz cos θ ) =



2π kz

1/2 X ∞  s=0

s+

 1 s i Js+1/2 (kz)Ps (cos θ ) 2

(13)

with Ps (cos θ ) the Legendre polynomial of order s. Hence the diffraction integral becomes  I0m =

2π kz

1/2 X ∞  s=0

 1 s s+ i Js+1/2 (kz)Qs,m (ρ) 2

(14)

with the functions Qs,m (ρ) given by Qs,m (ρ) =

θZmax

finc (θ )gm (θ )Jm (kρ sin θ )Ps (cos θ )dθ .

0

(15)

144

Michele Marrocco

Of course, the final comment for the A-P approach applies again to Eq. (14), where the series replaces the integral of Eq. (4). This means that a maximum index smax must be found that satisfies the criteria of numerical accuracy and rapid evaluation of I0m . ¨ ok ¨ A further approach based on the same strategy is proposed by Tor et al. (1997). The expansion is from Gradshteyn and Ryzhik (G-R; 1980); that is, ∞ X


exp ikz cos θ = J0 kz + 2 is Js kz cos(sθ ),

(16)

s=1

and the diffraction integral becomes ∞ X

I0m = J0 kz Q0,m (ρ) + 2 is Js kz Qs,m (ρ)

(17)

s=1

with the functions Qs,m (ρ) written as Qs,m (ρ) =

θZmax


finc (θ ) gm (θ )Jm kρ sin θ cos (sθ )dθ .

(18)

0

Again, in Eq. (17), summarizing the G-R approach, a series appears and the calculation might be reasonably advantageous if a cutoff index smax is determined to eliminate the terms that do not contribute meaningfully to the result. To complete this introduction to the methods based on phasor expan¨ ok ¨ et al. sions, it is intriguing to emphasize that the short review by Tor (1997) contains the earliest comparison between the approaches seen so far (A-P, Watson, G-R) and another method conceived by Kant (1993). The latter approach is discussed later in this review, but it is worth recalling that ¨ ok ¨ and co-workers concluded that DNS remains preferable in terms of Tor computation speed. The reason is soon explained. All the series are composed of terms that do not ensure a rapid convergence on an outcome that nears the DNS integration, whose results are thus obtained considerably faster. The attempts focused on the phasor expansion are not limited to the A-P, Watson, and G-R approaches. Indeed, more sophisticated methods have been studied. They rely on elaborate concepts beyond the simple mathematical expansions of the previous methods and, although the handling of the approximated diffraction integrals (or fields) is rather difficult, it is fitting to include a concise presentation of these alternative views on the central problem of this paper. ¨ ok ¨ (2005) and Sherif et al. (2008) In two recent works, Sherif and Tor have designed a new route to calculation of the diffraction integrals. The approach is termed eigenfunction representation or eigenfunction expansion

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

145

and the terminology refers to the eigenfunctions of the finite Hankel transform (Karoui and Moumni, 2009; Moore and Cada, 2004) that has attracted some recent interest in various applicative fields (Karoui and Moumni, 2009). The starting point is again the expansion of the phasor according to the formula (Arfken, 1985) ∞ X

exp ikz cos θ = Js kz exp[is (π/2 − θ)] ,

(19)

s=−∞

and the eigenfunctions φm,n (c, ρ) are defined by the integral equation Zρ0

p φm,n (c, ρ)Jm (ωρ) ρdρ = (−1)n (ρ0 / ) λm,n φm,n (c, ρ0 ω/ )

(20)

0

with ρ0 and Ω the cutoffs of the radial and frequency coordinates ρ and ω. The functions φm,n are invariant in Eq. (20) and are termed circular prolate spheroidal functions, whose eigenvalues are given by λm,n . The derivation of φm,n and λm,n is not easy and, concerning the determination of reliable numerical evaluations of the eigenfunctions as well as their eigenvalues, the debate is still open (Karoui and Moumni, 2009; Moore and Cada, 2004). ¨ ok ¨ have achieved the result that, after suitNonetheless, Sherif and Tor able manipulation of the exponential in Eq. (19), the diffraction integrals become " ∞ # ∞ p
X
sin θmax X I0m = Js kz bm,s,n (−1)n λm,n φm,n kρ0 , ρ/ρ0 sin θmax , kρ s=−∞ n=0

(21) where the coefficients bm,s,n are given by

bm,s,n =

1 λm,n

sin Zθmax

am,s (µ)φm,n (c, µ) µdµ

(22)

0

with the functions am,s (µ) indicating proper combinations of Chebyshev polynomials (Abramowitz and Stegun, 1965). As is apparent, the entire mechanism leading to the approximated calculation of I0m is quite involved. Plus, Eq. (21) contains two series that slow the convergence to the required accurate value of the diffraction integral.

3.2. Expansion of Apodization and Pupil Functions If the majority of the works are concentrated on the treatment of the rapid oscillations appearing with the phase term, Kant (1993) tried a different

146

Michele Marrocco

method based on the expansion of apodization and pupil functions. This idea could be regarded as an extension of the technique based on Legendre polynomials appearing naturally in harmonic analysis and, to put the method into perspective, the essentials of Kant’s proposal can be formulated as follows The apodization and pupil functions are defined on a finite interval of the azimuth angle and can be rewritten as suitable periodic functions called f (θ ) and qm (θ ) (the same notation used by Kant is used here). Their product can be expanded as the Fourier-Gegenbauer series f (θ )qm (θ ) =

∞ X

(m+1/2)

am,s Cs

(cos θ ),

(23)

s=0 (m+1/2)

where Cs (cos θ ) are the Gegenbauer polynomials of the first kind (Abramowitz and Stegun, 1965) and the coefficients are calculated according to Z am,s = Nm,s

f (θ )qm (θ ) 1 − cos2 θ


m

(m+1/2)

Cs

(cos θ ) sin θ dθ

(24)

with the normalization factor given by Nm,s =

22m−1 (2m + 2s + 1) 0 2 (m + 1/2) (s + 1) , π 0 (s + 2m + 1)

(25)

Under these circumstances, the diffraction integrals take on the following appearance: I0m = 2 sinm (ψ)

∞ X

(m+1/2)

am,s is Cs

(cos ψ)jm+s (ω),

(26)

s=1

where the functions jm+s (ω) are the spherical Bessel functions of the first kind and the variable ψ is such that kρ = ωsinψ. Equation (26) is the final result of the Kant approach and suffers from the same problem of other series expansions previously mentioned: The series replaces the integral with no real advantage. This was already noted in the review work of ¨ ok ¨ et al. (1997) but more insight is provided in Section 3.6. Tor

3.3. Multipole Expansion ¨ ok ¨ (1997) suggested another viewpoint of the problem Sheppard and Tor (different from the preceding approaches) inspired by Kant’s work (1993) based on Gegenbauer polynomials (see Section 3.2). These polynomials can be connected to the associated Legendre functions Pls (cos θ) of degree

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

147

l (with the condition that −s ≤ l ≤ s). The latter are, in turn, related to the spherical harmonics, which appear in multipole expansions of the electromagnetic fields. In effect, the solution of Maxwell equations can be written as the superposition of multipole fields that show dependences on the product of spherical Bessel functions of the third kind and spherical harmonics (Messiah, 1962; Panofsky and Phillips, 1962). This fact can be used to equate the far field of such a superposition to the field ¨ ok, ¨ after proper algeat the lens focus. In this way, Sheppard and Tor bra, achieved the following result given in terms of spherical coordinates r, θ , and φ: Er =

 ∞ X (−1)s qs s=1

Eθ =

3 2s + 1

 ∞ X (−1)s qs

1/2



 js−1 (kr) + js+1 (kr) P1s (cos θ ) cos φ

(27)

1/2 

 js−1 (kr) js+1 (kr) d 1 − P (cos θ ) s s + 1 dθ s s=1  i 2s + 1 1 + js (kr)Ps (cos θ ) cos φ (28) sin θ s(s + 1)

∞ X

3 2s + 1



! js−1 kr js+1 kr 1 Eφ = − P1s (cos θ ) − (−1) qs sin θ s s+1 s=1 
d 1 2s + 1 +i P (cos θ) sin φ, (29) js kr s (s + 1) dθ s s



3 2s + 1

1/2 "

where the j functions are the spherical Bessel functions as defined previously and the coefficients qs are given by Zπ qs = 2π

a (θ ) F1s (sinθ ) sinθ dθ

(30)

0

with a (θ ) a function related to the strengths of the multipoles and the function F depending on the associated Legendre function and its derivative. ¨ ok ¨ (which In addition to the interesting idea of Sheppard and Tor avoids the calculation of diffraction integrals at the expense of simplicity in the expressions of the fields), other authors proposed the use of spherical harmonics to reach a treatable result. As an example, we consider the work of Asatryan et al. (2004). Although this work cannot be qualified as a

148

Michele Marrocco

¨ ok ¨ pure multipole expansion in the sense explained by Sheppard and Tor et al. (1997), nevertheless it is closely related. As a matter of fact, spherical harmonics and spherical Bessel functions are commonly used in the multipole expansion of plane waves (Jackson, 1998), and the final fields can be rendered in terms of expressions containing coefficients that depend on the associated Legendre polynomials.

3.4. Methods of Local Phase and Amplitude Approximation The calculation of tightly focused fields is within the EDT realm, but the task could be equally fulfilled by adjusting the numerical techniques usually used to evaluate diffraction integrals known in Kirchhoff diffraction (Born and Wolf, 1970). To this end, we briefly describe a family of approaches whose validity is very general and, for this reason, their application cannot be confined to scalar theory only. These approaches are termed methods of local phase and amplitude approximation (LPA) (Gravelsaeter and Stamnes, 1982; Hopkins, 1957; Hopkins and Yzuel, 1970; Stamnes et al., 1983; Stamnes, 1986) and are often considered in optical microscopy. Unlike the methods previously mentioned where the key role is played by expansions with respect to radial and axial coordinates, the goal of LPA is to simplify the dependences on the azimuth angle such that the functions appearing in the amplitude and the argument of the phasor become treatable and, in turn, the integral of Eq. (4) becomes solvable. It is clear that the simplest way to realize the objective is the introduction of linear phase function and a constant amplitude as first explained by Hopkins (1957) more than fifty years ago. More sophisticated methods were later proposed, especially those by Stamnes and co-workers (Stamnes et al., 1983; Stamnes, 1986), and their applications have been demonstrated over the past two decades (Barbero and Marcos, 2008; D’Arcio et al., 1994; Dhayalan and Stamnes, 1997; Jain et al., 2009; Lotsberg et al., 2005; and many others not cited here). We provide a brief account of the method that best adheres to the calculation of Eq. (4). A more detailed review is given by Stamnes (1986). We consider linear approximations in the functions appearing in the kernel of the following general integral: ZθU I0m =

A(θ ) exp[iP(θ )]dθ ,

(31)

θL

where A(θ ) is the amplitude, P(θ ) the phase, and θ extends from the lower extreme θL to the upper extreme θU . The key idea is to divide the integral into many segments of extension from θL,n to θU,n and centered on the

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

149

average angle, θ¯n = (θU,n + θL,n )/2, so that

I0m

θ N ZU,n X = An (θ ) exp[i Pn (θ )] dθ ,

(32)

n=1θ

L,n

where, in each segment, it becomes reasonable to make linearizations of the amplitude and phase, that is, An (θ ) ≈ An,0 + An,1 (θ − θ¯n )/1θn and Pn (θ ) ≈ Pn,0 + Pn,1 (θ − θ¯n )/1θn , with 1θn = (θU,n − θL,n )/2. There are different choices for the coefficients appearing in the linear expansions of A(θ ) and P(θ ). For instance, the best choice for the phase factor is


1 2 Pn θL,n + Pn θU,n Pn θ¯n + 3 6

 1 = Pn θU,n − Pn θL,n 2

Pn,0 =

(33)

Pn,1

(34)

so that the phase error is minimized (Gravelsaeter and Stamnes, 1982; Stamnes et al., 1983). An analogous reasoning applies for the coefficients of the amplitude An (θ ). The solution of the integrals in Eq. (32) can be obtained as a simple equation in the form of ZθU,n An (θ) exp[i Pn (θ )] dθ = Xn + iYn ,

(35)

θL,n

where Xn = 21θn

Yn = 21θn

! Pn,1 cos Pn,1 − sin Pn,1 sin Pn,1 An,0 cos Pn,0 + An,1 sin Pn,0 Pn,1 P2n,1 (36) ! sin Pn,1 Pn,1 cos Pn,1 − sin Pn,1 An,0 sin Pn,0 − An,1 cos Pn,0 Pn,1 P2n,1 (37)

At this point, the whole integral I0m is now reduced to a simple summation P of Eqs. (36) and (37), I0m = N n=1 (Xn + iYn ). It is self-evident that the LPA method is straightforward and holds promise of fast execution. Indeed, the novelty is that local approximations of the type discussed herein do not involve special functions that

150

Michele Marrocco

are always fundamental to the methods considered previously. However, the only unanswered question regards the number of intervals required to achieve a sufficient accuracy for the local approximation to work well. On the other hand, this is equivalent to the problem arising with the truncated series in the methods based on expansions. These aspects are addressed in Section 3.6. To conclude, it must be noted that the reasoning leading to Eq. (35) can be reformulated by assuming more-complex functional dependences for the amplitude and the phase. In effect, parabolic functions represent the next level of complication and, fortunately, solutions as simple as Eq. (35) can still be found (Stamnes et al., 1983). This upgraded version of the LPA approach has inspired different works. For instance, D’Arcio et al. (1994) presented an attempt where the local approximation of the diffraction integrals is realized by using polynomials for the amplitudes and parabola for the phase. The combination of such simple functions was successfully applied in the calculation of diffraction from apertures of complicated shapes that were simulated by patching together rectangular and arched elements.

3.5. Method Based on Transition from Continuous to Discrete Diffraction The integral of Eq. (4) shows an important property of diffraction: its continuous nature. By contrast, any course on optics teaches from the very beginning that the boundary conditions on the electric and magnetic fields determine the quantization of electromagnetic modes (and hence quantization of wave vectors) (Jackson, 1998; Loudon, 2000). Although, the classic quantization does not apply to our problem, we can nonetheless ask if it is possible to think of a method where the discrete characteristic of classic quantization can somehow be transferred to diffraction. where the wave vectors are continuously distributed between the extremes allowed by the angular aperture (or, more in general, by the shape of the microscope aperture). The advantage is clear: The removal of the integral without touching the formal dependences in its kernel would greatly simplify the calculation. The idea was recently developed by Marrocco (2009a) to simulate the focused intensity of a high-NA microscope and a scattering problem of nonlinear microscopy. This review summarizes the main result as follows. Diffracted wave vectors depend on the diffraction angle θ and, in the case of a clear pupil of circular shape, their continuity is guaranteed by the variable angular aperture lying between the axial direction (θ = 0) and a maximum value (θ = ϑmax ) corresponding to the NA. As a consequence, the transition to a discrete set of wave vectors is established in agreement

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

151

with the replacement of the integral with respect to the azimuth angle θ with a summation running over a number N of θ values. The strategy to achieve this goal is based on the simple physical observation that the most important spatial region is around the focal point (ρ, z) = (0, 0). This suggests a perturbative approach where Taylor expansions can be used to replace the radial and axial dependences in the diffraction integrals of Eq. (4). This means that, in marked contrast with other methods based on specific expansions, we replace both the Bessel function and its phasor with ∞
X
2l+m h 2l+m
i Jm kρ sin θ = / 2 l! m + l ! (−1)l kρ sin θ

(38)

l=0

∞
X
t exp i k z cos θ = i k z cos θ /t!.

(39)

l=0

With the further Taylor expansion of the apodization function and applying the theorem for the integration of uniformly convergent series (Apostol, 1967), we find after some algebra that we can decompose I0m into a limited number N of elementary terms, that is, N  X     I0m = exp −β 2 / sin2 θmax am,p gm,p Jm kρ sin θm,p exp i k z cos θm,p , p=1

(40) where gm,p = 1 for m = 0, 2, but g1,p = sin θ1,p for m = 1 and the angles θm,p are found together with the coefficients am,p on the basis of the finite summation N
X cm q = am,p cosq θm,p ,

(41)

p=1

where the left-hand side is constrained to some values obtained by means of the parameters belonging to the optical system ( ∞
X (β/ sin θmax )2s 1 − (cos θmax )3/2+q+2s cm q = s! 3/2 + q + 2s s=0 )

δm,2 1 − (cos θmax )5/2+q+2s + 1 − δm,1 −1 (42) 5/2 + q + 2s

152

Michele Marrocco

with q taking on integer numbers and δm,j the usual notation for the Kronecker delta. Equation (42) is calculated once and for all in dependence on the experimental conditions and represents the sole numerical calculation leading to the simple solution of I0m . In other words, having established the initial conditions (i.e., filling factor β and maximum aperture θmax ), the coefficients cm (q) are calculated, and this makes it possible to determine [through Eq. (41)] the parameters am,p and θm,p needed for the decomposition of Eq. (40). Regarding the latter, its physical meaning is now manifest. Among the various projections of wave vectors onto the optical axis and the radial plane, the entire procedure isolates those that are sufficiently informative about the diffraction pattern. In this manner, the original structure of the focused field—that is, the Bessel function in the radial plane multiplying the propagation function along the optical axis [see Eq. (41)]–is reproduced in the linear combination on the righthand side of Eq. (40). The last question that is still unanswered is about the definition of N. Different alternatives can be chosen. Here, we prefer to set N according to the criterion that, away from the focal zone (that is, at fixed values of ρ and z far from the focus), the deviation between numerical and analytical calculation of the EDT intensity is negligible (for instance, less than 0.1% as in Marrocco (2009a). Following this specific criterion and taking k z = ±10 and kρ = 20 for a microscope of NA = 1.4, n2 = 1.5, and β = 1, it is found that the number of projections is surprisingly small (N = 10) compared with the complex structure of a focus of such a high-NA microscope.

3.6. Comparisons Among the Methods for Linearly Polarized Optical Beams The availability of different approaches to the calculation of the diffraction integrals of Eq. (4) poses a consequent question about the real advantage of working with one of the alternative numerical methods instead of the direct calculation of I0m . As anticipated in Section 3.1, the question ¨ ok ¨ et al. (1997), who examined was partially addressed in the past by Tor four approximations of I0m under uniform illumination of the microscope lens. Two of them were obtained by introducing known mathematical expressions for the phasor appearing in Eq. (4) and the results are here summarized in Eqs. (14) and (17). One equation represents the method based on Watson (1952). The other describes the diffraction integrals calculated according to G-R (1980). The strategy behind these two approximations of I0m does not differ from the method conceived by A-P (1979) except for the use of different special functions in the series expansion of the phasor–namely, Hankel functions in the A-P method versus a combination of Bessel functions and Legendre polynomials in the Watson ¨ ok ¨ approximation or Bessel functions and cosines in the G-R approach. Tor

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

153

and co-workers not only compared these three analogous methods with each other but also scrutinized the idea suggested by Kant (1993), who had introduced the use of Gegenbauer polynomials combined with spherical Bessel functions. Kant’s result is concisely illustrated in Eq. (26). ¨ ok ¨ et al. (1997) is in the The importance of this earlier work by Tor demonstration of the futility of alternative solutions to DNS calculation of I0m . In the 1997 paper, it is clearly stated that “the expansions do not speed up numerical evaluation of the diffraction integrals” and thus these alternative methods are pointless as long as the aim is an improvement in the computation efficiency accompanied by reliable determination of I0m . Indeed, the other theme connected to the use of numerical methods is related to their accuracy, which should also be evaluated. The demonstrated inadequacy of the above-mentioned numerical “recipes” has a simple explanation. In all the proposals, series of special functions make their appearance. Since these more sophisticated functions are computed by means of basic functions, it is reasonable to expect that the series of such special functions are not faster than the DNS solution of the integral of Eq. (4). The same conclusion holds for the approach sug¨ ok, ¨ 2005; Sherif et al., 2008) gested by Sherif and others (Sherif and Tor discussed in Section 3.1. This last method uses the finite Hankel transform that leads to the strenuous use of circular prolate spheroidal functions. Eventually, the diffraction integrals are calculated by means of the product of the two series of Eq. (21) (note, two series instead of one as for the most methods considered in this review). The first runs over the order of Bessel functions of the first kind, the second runs over an index that relates the spheroidal functions to some numerical coefficients, among which the one in Eq. (22) is obtained from the integration of the spheroidal functions multiplied by the functions aN,s (µ) containing proper combinations of Chebyshev polynomials. Unquestionably, since the declared inefficiency ¨ ok ¨ et al., of the cited one-series methods (Watson, G-R, A-P, and Kant) (Tor 1997), it is possible to conclude that, beyond the unavoidable complexity inherent in the handling of eigenfunctions of the finite Hankel transform, the heavy load of two series and different special functions in Eq. (21) impedes the fast calculation of diffraction integrals. For this reason, the approach of Sherif et al. is not included in the following comparisons. The suggestions for the numerical calculation of Eq. (4) are completed by other proposals that could be compared with the DNS, Watson, G-R, A-P, and Kant methods. One of these is the linear version of the LPA approach (Gravelsaeter and Stamnes, 1982; Stamnes et al., 1983) discussed in Section 3.4. The other is suggested by Marrocco (2009a); its brief account is given in Section 3.5. ¨ ok ¨ et al. (1997). The comparisons are made in a fashion similar to Tor Although we might reach an identical conclusion as to the DNS, Watson, G-R, A-P, and Kant methods, the attempt is noteworthy for a number of

154

Michele Marrocco

reasons. First, the calculation speed of modern computers is one or two orders of magnitude higher than what was available more than ten years ¨ ok ¨ et al. used a 90-MHz PC, the following data are obtained ago. If Tor by means of an ordinary PC with an internal clock of 2.66-GHz. In addition, new versions of mathematical tools, such as MatLab (MathWorks, Inc.) or Mathematica (Wolfram Research, Inc.), simplify the effort of treating the series and the integration of special functions. For the purpose of this study, Mathematica 6.0 is used to elaborate the diffraction integrals according to the DNS, Watson, G-R, A-P, Kant, LPA and Marrocco methods. The comparisons are arranged for two ideal optical systems. One is identified by an oil-immersed lens with NA = 1.4 and index of refraction of 1.5. The other system is for a microscope of NA = 0.9 operated in air. The laser beam is supposed with a Gaussian profile and the filling factor is chosen equal to one. The results are reported in tables for the most important diffraction
integral I00
(see Figure 1) calculated at the spatial points kρ, kz = (0, 0) and kρ, kz = (10, 10). These represent the two extremes of the focal point and the region away from the focus, respectively. The tables are organized as follows. First, the series appearing in some methods (i.e., Watson, G-R, A-P, Kant) are truncated to avoid very long computation runs. The number N of terms that are meaningful for the comparative use of the truncated series is given in the first column. The second column indicates the run time T of the calculation evaluated in units of ms. The ratio R = T0 /T between the run time T0 taken by the traditional DNS solution of I00 and T is shown in the third column. Obviously, the condition R > 1 identifies the methods that are convenient in terms of their speed. The last three columns show the relative deviation between the DNS value of I00 and its approximated results obtained for each alternative method. Since I00 is usually a complex function of radial and axial coordinates, the three columns are for the real (Re) and imaginary (Im) parts plus the absolute value (Abs) of the diffraction integral. The truncated series calculated for the Watson, G-R, A-P, and Kant methods have been constructed in such a way that relative deviations less than or on the order of 1% could be guaranteed. A similar criterion was adopted for the LPA method. The last method (Marrocco) does not rely on series expansion and, based on preliminary considerations (Marrocco, 2009a), it was known that about ten terms are sufficient to ensure a good selection of wave vectors characterizing diffraction of the chosen optical system. As a consequence, nine terms are used in the tables. Let us analyze the results in Table 1 relative to the most important point of the focal region (i.e., kρ = 0, kz = 0). As expected, the Watson, G-R, A-P, and Kant methods are very well represented by just one single term, the first of their respective truncated series. Despite this advantage, their speed does not exceed that of DNS. By contrast, the LPA and Marrocco

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

155

TABLE 1 Comparison at the Point kρ = 0, kz = 0∗ Method

N

T (ms)

R

Err Re (%)

Err Im (%)

Err Abs (%)

Watson† G-R A-P† Kant† LPA∗ Marrocco

1 1 1 1 3 9

7.1 7.2 7.2 7.3 0.89 0.19

1.02 1 1 1 8.1 37.4

0 0 4.44 × 10−14 0 −0.04 4.49 × 10−7

No No No No 0.42 No

0 0 4.44 × 10−14 0 −0.039 4.49 × 10−7

∗ For a microscope defined by the following parameter values NA = 1.4, n = 1.5, β = 1. The reference values 2

are I00 = 0.560 and T0 = 7.2 ms. Note that the methods indicated with a dagger (†) diverge at kz = 0 and a point in the very proximity of the origin (kz = 10−20 ) has been chosen. N, Number of terms; T, run time. The ratio R = T0 /T between the run time T0 taken by the traditional DNS solution of I00 and T is shown in the third column. The condition R > 1 identifies the methods that are convenient in terms of their speed. The last three columns show the relative deviation between the DNS value of I00 and its approximated results obtained for each alternative method. Since I00 is usually a complex function of radial and axial coordinates, the three columns are for the real and imaginary parts plus the absolute value of the diffraction integral.

TABLE 2

Comparison at the Point kρ = 10, kz = 10∗

Method

N

Watson G-R A-P Kant LPA Marrocco

15 15 18 13 30 9

T (ms)

235 181 134 136 8.8 0.71

R

Err Re (%)

Err Im (%)

Err Abs (%)

0.05 0.07 0.09 0.09 1.46 18.2

0.012 −0.12 −2.25 × 10−7 −0.76 0.71 0.28

−1.64 4.13 2.74 −10.6 3.8 0.71

9.7 × 10−4 −0.09 0.018 −0.82 0.011 0.28

∗ For a microscope defined by the following parameter values NA = 1.4, n = 1.5, β = 1. The reference values 2

are I00 = 0.056 + i0.0045 and T0 = 12.9 ms. N, Number of terms; T, run time. See notes to Table 1.

methods are, respectively, about 8 and 40 times faster than DNS. Remarkably, a certain number of terms is needed in these methods, whereas the Watson, G-R, A-P, and Kant methods are extremely accurate with one single term. The results in Table 2 are in agreement with the conclusion of Table 1 except for the absolute values. Now the Watson, G-R, A-P, and Kant methods are sensitively slower than DNS, which agrees with the outcome of ¨ ok ¨ et al. (1997). In addition, the LPA method is less the analysis by Tor effective but its run time is still shorter than the DNS calculation. A reduction in efficiency is observed for the Marrocco method as well. However, the Marrocco result is found about 20 times faster than the DNS result. Moreover, with the exclusion of the Marrocco approach, all the remaining methods required a much larger number of terms to safeguard a reasonable convergence to the chosen accuracy.

156

Michele Marrocco

TABLE 3

Comparison at the Point kρ = 0, kz = 0∗

Method

N

T (ms)

R

Err Re (%)

Err Im (%)

Err Abs (%)

Watson† G-R A-P† Kant† LPA† Marrocco

1 1 1 1 2 9

7.1 7.2 7.3 7.4 0.89 0.19

1.01 1 0.99 0.97 8.1 37.9

0 0 2.22 × 10−14 2.22 × 10−14 −0.2 1.90 × 10−7

No No No No 0.28 No

0 0 2.22 × 10−14 2.22 × 10−14 −0.2 1.90 × 10−7

∗ For a microscope defined by the following parameter values NA = 0.9, n = 1, β = 1 . The reference values 2 are I00 = 0.519 and T0 = 7.2 ms. Note that the methods indicated with a dagger (†) diverge at kz = 0 and the point kz = 10−20 has been chosen. N, Number of terms; T, run time. See notes to Table 1.

TABLE 4

Comparison at the Point kρ = 10, kz = 10∗

Method

N

Watson G-R A-P Kant LPA Marrocco

16 16 16 16 35 9

T (ms)

244 189 120 167 9.7 0.7

R

Err Re (%)

Err Im (%)

Err Abs (%)

0.049 0.062 0.1 0.071 1.23 17.0

0.027 −0.11 −2.36 × 10−6 −0.98 0.44 0.043

0.056 −1.16 1.06 −1.52 7.7 1.89

0.24 −0.11 7.79 × 10−4 −0.98 0.45 0.044

∗ For a microscope defined by the following parameter values NA = 0.9, n = 1, β = 1. The reference values 2 are I00 = 0.051 + i0.0014 and T0 = 11.9 ms. N, Number of terms; T, run time. See notes to Table 1.

The other tables (Tables 3 and 4) confirm the results of Tables 1 and 2. This fact implies that the change of the optical system is not a fundamental variable of the problem. For this reason, plots of the speed ratio R as a function of the dimensionless spatial coordinates kρ and k z are shown in Figures 4 and 5 for the optical system with higher NA. In particular, Figure 4 displays the behavior of R in dependence on the radial coordinate taken within the plane at the focus. In contrast, Figure 5 is informative about the speed ratio along the axial axis. Other plots arranged by changing simultaneously both the spatial coordinates kρ and k z looked very similar to Figure 4 and are not shown. Generally, in any of the points of the plane (kρ, kz), it appears clear that the Watson, G-R, A-P, and Kant methods always have ratios R smaller than or comparable to unity. Thus, in ¨ ok ¨ et al. (1997), we should conclude that despite modagreement with Tor ern mathematical tools (here, Mathematica 6.0) that promote and facilitate numerical elaboration, such methods are definitely slower than the direct numerical integration of the diffraction integrals. On the other hand, the Marrocco method provides ratios R that are significantly larger than 1.

157

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

40 Speed ratio along the radial axis (kz = 0) NA = 1.4; n2 = 1.5; β = 1

30 20 10 R 1.00

W 0.75

GR AP

0.50

K LPA

0.25

M

0.00 0

2

4

6

8

10

kρ

FIGURE 4 Comparison among various methods of Section 3 based on the speed ratio R in the focal plane kz = 0 plotted as a function of the dimensionless radial coordinate kρ. Key to methods: W, Watson; GR, Gradshteyn–Ryzhik; AP, Agrawal and Pattanayak; K, Kant; LPA, local phase and amplitude approximation; M, Marrocco.

More precisely, the Marrocco method is about 40 times faster than any other method based on expansions if the calculation is limited to points very close to the focus. As soon as the calculation is made in the region near the focus, the efficiency of the Marrocco method is even more striking in comparison to the other approaches. Among them, the LPA method shows interesting values of R but its results are not as rapidly attained as for the Marrocco approach. In conclusion, the method of discrete diffraction seems to be capable of faster calculation speed and holds promise of very efficient response in applications (i.e., imaging and nonlinear microscopy) where higher calculation speed is recommended.

4. METHODS OF APPROXIMATIONS FOR RADIALLY POLARIZED LASER BEAMS Alternative solutions to the evaluation of diffraction for high-NA microscopes are not exclusive of linearly polarized beams and, among the possible choices, we limit discussion in this section to the radial polarization,

158

Michele Marrocco

40 Speed ratio along the axial axis (kρ = 0) NA = 1.4; n2 = 1.5; β = 1

30 20 10 R 1.00

W GR

0.75

AP K

0.50

LPA M

0.25 0.00 0

2

6

4

8

10

kz

FIGURE 5 Comparison among various methods in Section 3 based on the speed ratio R in the axial axis kρ = 0 plotted as a function of the dimensionless radial coordinate kz. Key to methods: W, Watson; GR, Gradshteyn–Ryzhik; AP, Agrawal and Pattanayak; K, Kant; LPA, local phase and amplitude approximation; M, Marrocco.

which has attracted much interest lately (Bomzon et al., 2002; Dorn et al., 2003; Lerman and Levy, 2008; Novotny and Hecht, 2006; Quabis et al., 2000; Shoham et al., 2006; Yew and Sheppard, 2007a; Youngworth and Brown, 2000). The reason for such interest is found in the potential for beating the diffraction limit of ordinary microscopy. This was demonstrated a few years ago (Dorn et al., 2003) and since then, studies on radial polarization of electromagnetic waves have continued to grow steadily; now applications in various fields of physics have become feasible. Recent examples include the acceleration of massive particles (Salamin, 2007), the design of a plasmonic lens (Yanai and Levy, 2009), measurements with terahertz radiation (Grosjean et al., 2008), tip-enhanced Raman microscopy (Steidtner and Pettinger, 2008), and second harmonic generation microscopy (Yew and Sheppard, 2007b); the list could continue with many other examples in the specialized literature. In high-resolution microscopy, radially polarized laser beams are specifically used to reduce the spot size near focus. This is in line with the original proposal (Dorn et al., 2003; Quabis et al., 2000), and studies concerning the optical conditions leading to spot sizes well below the diffraction limit have been reported (Lerman and Levy, 2008; Yew and Sheppard, 2007a). These works are based on the diffraction theory developed by

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

159

Youngworth and Brown (2000) (which is, in turn, derived from the pioneering work of Richards and Wolf, 1959). Following this theory, the vectorial structure of the focused electric field E(ρ, z) of radially polarized beams resembles what is seen for linear polarization. A short summary is given in Section 2.2 where, considering an annular pupil with angular apertures between α1 and α2 , Eqs. (7) and (8) decompose the total electric field into its radial and longitudinal contributions, respectively. In general, the most common choice for the determination of Eρ (ρ, z) and Ez (ρ, z) is given by the DNS calculation, but other attempts that avoid the difficulty of purely numerical approaches to the diffraction integrals are reported in the literature. In effect, research is under way to look for analytical expressions of the fields. For instance, the Lax series (i.e., the solution of the wave equation written as polynomial expansion with respect to the diffraction angle) is used by some authors (Luo et al., 2007; Salamin, 2006a). Others suggest an analytical approach based on the so-called complex-source-point spherical wave description (Yan and Yao, 2008). However, the vector angular spectrum mentioned for the EDT theory is also adopted to reach analytical results in a fashion different from the original EDT version (Deng and Guo, 2007). The method of the ¨ ok, ¨ 2005), eigenfunctions of the finite Hankel transform (Sherif and Tor the LPA method (Gravelsaeter and Stamnes, 1982; Stamnes et al., 1983; Stamnes, 1986), and the approach of discrete diffraction (Marrocco, 2009b) can also be extended from the case of linear polarization (examined in Section 3) to radially polarized beams. The following text details these approaches and compares them so that the declared advantages of these methods can be put to the test. The results show that important limitations discourage the use of the solutions to the wave equation at the high NA usually used in research on radial polarization so that EDT methods should be preferred.

4.1. Method Based on the Lax Series In past years, considerable attention has been directed to the approach built on the work of Lax et al. (1975) and Davis (1979). Despite the number of contributors to the approach, the methods inspired by these seminal papers are commonly referred to as the Lax series. The idea here is that the breakdown of the paraxial approximation can be judiciously described by perturbative corrections that stem from the expansion of the general fields satisfying the exact Maxwell equations. The initial results were obtained for Gaussian beams with linear polarization and corrections up to the fifth order can be found in the literature (Barton and Alexander, 1989). Although the method itself does not lead to diffraction integrals of the kind shown in Eqs. (7) and (8), the final objective is in common with the work of Youngworth and Brown (2000)—the analytical description of the focused fields.

160

Michele Marrocco

The research on the Lax series for high-resolution laser microscopy with radial polarization started very recently (Luo et al., 2007; Salamin, 2006a). In particular, the series is useful to calculate the expansion of the radial and axial electric fields with respect to the diffraction angle ε (defined as the ratio between the Gaussian beam waist w0 and the Rayleigh length zR = kw20 /2) (Luo et al., 2007; Salamin, 2006a); the result obtained initially by Salamin (2006a) shows the dependence on ε up to the fifth order. Defining the dimensionless coordinates zˆ = z/zR and ρˆ = ρ/w0 , the approximated fields are ! "     ρˆ 5 f 5 ρf ˆ 3 3 4 2 2 3 + ρˆ f − Er = E0 exp −f ρˆ + iη ε ρf ˆ +ε − 2 2 !# 3ρf ˆ 4 3ρˆ 3 f 5 17ρˆ 5 f 6 3ρˆ 7 f 6 ρˆ 9 f 8 + ε5 − − + − + (43) 8 8 16 8 32 " !#     3 ρˆ 2 f 4 5ρˆ 4 f 5 ρˆ 6 f 6 2 2 2 2 3 4 f Ez = iE0 exp −f ρˆ + iη ε f − ρˆ f +ε + − + , 2 2 4 4 (44) where E0 is the field amplitude, f = i/(i + zˆ ) is related to the Gouy phase, and η = ωt − kz is the phase of the plane wave. Immediately after the publication of Eqs. (43) and (44), the same author extended the calculation to the 15th order in the diffraction angle ε (Salamin, 2006b, 2008). These new results are not reported here, but it is interesting to mention that they were completed by the work of Luo et al. (2007), who managed to find expressions corrected to any order, but their coefficients were to be found from recursion relations. Another point of interest about the Lax series is raised by Salamin (2009), who compares his work with the representation of a Gaussian beam according to another technique known as the complex-source-point spherical wave description. This second method is the subject of the next section.

4.2. Method Based on Complex-Source-Point Spherical Wave This method originates from an idea of Couture and Belanger (1982). They were inspired by the work of Lax and co-workers, but they succeeded in solving the Helmholtz equation for the vector potential by using the socalled complex-source-point spherical wave (CSPSW) mentioned above. Indeed, the perturbative solution to the scalar Helmholtz equation can be elaborated in such a way that it takes on the CSPSW form written as exp(−ikRC )/RC with RC = [ρ 2 + (z + izR )2 ]1/2 (Eq. (20) in Couture and Belanger, 1982).

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

161

Using this framework with its relationships among vector potential and electric fields and, finally, adjusting it to radially polarized beams, Yan and Yao (2008) obtained the following result: !
iωA0 exp ikRc (z + izR ) 3 3ik k2 ρ + 4 − 5 (45) Er = k2 R3C RC RC !
iωA0 exp ikRc 3ik (z + izR )2 1 ik k2 (z + izR )2 3 (z + izR )2 Ez = − + 3 − 2+ k2 RC R3C R5C RC R4C
exp ikRC (46) − iω , RC

where A0 is the amplitude of the vector potential. These expressions of the electric fields are very compact and in contrast with the Lax series derived by Salamin (2006a, 2006b). As a matter of fact, Yan and Yao (2008) concluded their work with comparisons with the fifth-order fields of Eqs. (43) and (44) (the comparison is in Table 1 of Yan and Yao). They found that the CSPSW fields better satisfy the Maxwell equations with errors that were, on average, always on the order of 10−4 %, whereas the Lax series was incapable of small errors for values of the diffraction angle ε greater than 0.4. The CSPSW equations were also reformulated by Salamin in a new work (Salamin, 2009), where comparison with the Lax series is arranged such that dependences of the radial and longitudinal fields on the radial coordinate in the focal plane are illustrated. The new equations are quite similar to Eqs. (45) and (46) herein except for some algebraic signs that, if neglected, can culminate in strong disagreement with those of Yan and Yao (2008) and Salamin (2009). In more detail, we found that one important difference is in the exponential factor exp(ikRc ) (other differences do not seem crucial). This appears as exp(−ikRc ) in Salamin’s derivation. This fundamental difference has a consequence that cannot be overlooked because Rc is a complex variable and a plus or minus sign can cause contrasting behaviors. Such a difference is explained in the redefinition of this parameter by Yan and Yao, who rearranged it using the paraxial expansions. On the contrary, Salamin used the original definition of Couture and Belanger (1982). To avoid any confusion in the comparison of Section 4.5, we use Salamin’s formulas for the CSPSW equations (not reported here for the sake of brevity).

4.3. Method of Transverse Electric and Transverse Magnetic Decomposition The starting point of the previous methods conceived for radially polarized beams was the Helmholtz equation that was solved with a

162

Michele Marrocco

perturbative criterion (Lax series) or in closed form (CSPSW). Deng and Guo (2007), however, advance the idea that the angular spectrum representation of the fields in the transverse electric (TE) and transverse magnetic (TM) decomposition can be useful in finding an analytical result. This was obtained for Laguerre–Gaussian beams whose mathematical description in the focal plane at z = 0 can be adequately treated in the angular spectrum representation to determine the field away from the focus. After Fourier transformation of the field components and some additional algebra, Deng and Guo report the following equation for the TE field: ETE

√ (−1)n+1 E0 2zR

   z 1 2 exp[ikr−c(r)] P(r) L xˆ e + yˆ e −ρ /zˆ e = [Q(r)] x y z n 2σ r r2 (1 − izR /r)2 # n X (−1)m−n (n + 1) ! (47) Lm+1 [c(r)] , + iw0 eˆz (n − m) !m! (1 − izR /r)m m=0


where c(r) = s(r)/[2(1 − izR /r)] with s (r) = ρ 2/ 2r2 σ 2 , σ = 1/k w0 , P (r) =
[(1 + izR /r) / (1 − izR /r)]n , Q (r) = s (r) / 1 + iz2R /r2 , and L1n (.) is the generalized Laguerre polynomial of the radial mode number n and the angular mode number 1 and Lm+1 (.) is the Laguerre polynomial of (m + 1)th order. Equation (47) was used to obtained plots of the energy density at distances of z = 10λ and z = 100λ. Shorter axial distances cannot be studied with this method because the angular spectrum representation used here connects the field at the focus to the field away from it and, based on such a premise, the whole work by Deng and Guo cannot be applied to points near focus where z ∼ = 0.

4.4. Methods of Eigenfunction Representation, Local Phase and Amplitude Approximation, and Discrete Diffraction The approaches of Sections 4.2 and 4.3 were born as attempts to overcome the paraxial approximation of radial polarization of Gaussian beams. However, among the methods reviewed in Section 3, there are some that are not conditioned uniquely to linear polarization. The method based on eigenfunctions of the Hankel transform plus the general approximations of the LPA method and discrete diffraction are potentially applicable to the current subject of radially polarized beams. (We do not repeat the general concepts behind these methods; they can be recalled by reading the corresponding text of Section 3.) Nevertheless, a brief review of the relevant conclusions valid for radial polarization is essential to the development of our argument.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

163

With regard to the eigenfunction representation, the expressions of the focused fields can be obtained through the diffraction integrals of Eq. (4) ¨ ok, ¨ 2005) solved for the current case. The result is (Sherif and Tor " ∞ # ∞ X p sin θmax X Eρ = 2E0 Js (kz) b1,s,n (−1)n λ1,n φ1,n kρmax s=−∞ n=0
× kρ sin θmax , ρ/ρmax sin θmax (48) sin θmax Ez = 2iE0 kρmax

"

∞ X s=−∞

# Js (kz)

∞ X

b0,s,n (−1)n

p

λ0,n φ0,n

n=0

× kρ sin θmax , ρ/ρmax sin θmax




(49)

with the parameters of the eigenfunctions already defined in Eq. (22). The general comment to these expressions of the fields is similar to that stated in Section 3.1. The approach is very involved and does not lend itself to a reasonable application. Furthermore, beyond the undeniable difficulty related to the treatment of the various special functions hidden in the parameters of Eqs. (48) and (49), the use of two series is time consuming in the elaboration of the fields in the focal region. Unlike the eigenfunction representation of diffraction integrals, the technique established in some works regarding the LPA method (Gravelsaeter and Stamnes, 1982; Stamnes et al., 1983; Stamnes, 1986) is simple and its general validity allows application to the diffraction integrals of Eqs. (7) and (8). For instance, the reasoning for the purely linear version of the LPA method proceeding from Eq. (31) to Eq. (37) can now be translated in terms of radially polarized beams in agreement with the work of Youngworth and Brown (2000). The application is particularly stimulating because Tables 1 through 4 and Figures 4 and 5 demonstrate the advantage of quick LPA elaboration of the integrals for linear polarization of the fields. Similar to the LPA methods, discrete diffraction (see Section 3.5) is also extendable to radial polarization (Marrocco, 2009b). Repeating the same reasoning for linealy polarized beams and taking as reference the fields of Eqs. (7) and (8), it possible to determine the following expressions: Eρ(ρ, z) = A

Nρ X



aρ,h J1 kρ sin θρ,h exp ikz cos θρ,h

(50)

h=1

Ez (ρ, z) = iA

Nz X h=1



az,h J0 kρ sin θz,h exp ikz cos αz,h ,

(51)

164

Michele Marrocco

where Nρ and Nz indicate the total numbers of terms necessary to approximate Eqs. (7) and (8) with Eqs. (50) and (51) within a certain accuracy and aρ,h , az,h are the coefficients of the two linear combinations, cρ (q) =

Nρ X


q

(52)


q az,h sin θz,h ,

(53)

aρ,h sin θρ,h

h=1

cz (q) =

Nz X h=1

whose left-hand members cρ and cz depend exclusively on the optical parameters characterizing the experimental setup. In other terms, cρ and cz are calculated once and for all when α1 , α2 , and β are given. The operation is realized according to

cρ (q) =

2 sin Z α2

h i exp −(β/ sin α2 )2 ζ (1 − ζ )1/4 ζ (q+1)/2 dζ

(54)

sin2 α1

cz (q) =

2 sin Z α2

h i exp −(β/ sin α2 )2 ζ (1 − ζ )−1/4 ζ (q+1)/2 dζ .

(55)

sin2 α1

This is the only numerical calculation needed in this procedure and the values of cρ (q) and cz (q) can be tabulated and stored in view of the analytical elaboration of the fields.

4.5. Comparisons Among the Methods for Radially Polarized Optical Beams After the foregoing excursion into methods for the calculation of tightly focused electric fields we can categorize them according to their general properties. First, we considered theoretical approaches based on the solution of the Helmholtz equation: Lax series and CSPSW. We briefly comment on them. The Lax series follows a perturbative criterion and, by its definition, it can be applicable as long as the perturbative conditions are valid. These are summarized by the constraint that the diffraction angle must be less than 1: ε < 1 or, recalling that ε = w0 /zR , the constraint is w0 > λ/π. This means that the Gaussian beam for which the Lax series can be used must have a beam waist at the focus adequately greater than λ/π. In numerical terms, if the laser beam injected into the microscope is generically within

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

165

the infrared portion of the electromagnetic spectrum (this is a common situation in biochemical studies), the waist radius must be greater than 300 nm. The natural consequence is that the Lax series cannot be aimed at the description of lenses with very high NA that are capable of beam sizes at the focus on the order of 100 nm or less. In this regard, the restriction on the Lax series is in manifest contrast with one of the most distinguished characteristics of radially polarized beams as a means to reduce the size of the focal spot below the diffraction limit (Dorn et al., 2003). If the Lax series shows this impediment to the broad application of radial polarization to high-resolution microscopy, the CSPSW model does not make an exception. Apparently, the neat analytical formulation of the focused fields in Eqs. (45) and (46) (or the Salamin, 2009, version) is ideal for fast elaboration. The main drawback is that the fields are no longer representative of the optical system whenever the spatial coordinates are such that one of the singularities is reached. This happens at the points where Rc → 0 and, for example, in the focal plane, the singularity appears at ρ = zR . In terms of normalized coordinate ρ 0 = ρ/w0 , the condition of validity of the CSPSW approach reads ρ 0 < 1/ε and, for one of the critical diffraction angles of the Lax series, say ε = 0.5, then the condition becomes ρ 0 < 2 or ρ < 2w0 . The further consequence is that large diffraction angles limit more severely the range of application of Eqs. (45) and (46) and, in this respect, the CSPSW formulation presents the same limitation as the Lax series when high NA values are requested. Common to both formalisms is another restriction implied by the initial hypotheses of the two methods. Indeed, Gaussian beams at the focus (defined by the parameters w0 and zR ) make up the premise of the theoretical analyses leading to the fields according to either the Lax series or the CSPSW technique. On the other hand, the role of an aperture is not contemplated here. On the contrary, the aperture can be mathematically incorporated in works based on the integral representation of Eqs. (7) and (8). This difference draws attention to the second group of the methods reviewed in the preceding sections. The second group is formed by the TE and TM decomposition (Section 4.3) and some of the techniques taken among those valid for linear polarization—namely, eigenfunction representation, LPA, and discrete diffraction (Section 4.4). Starting with the results of Deng and Guo (2007), it was already understood that the main problem of the TE and TM decomposition is its restriction to fields whose axial coordinate is much greater than the beam waist (z  w0 ). Since we are interested in the physics at the focus (z < w0 ) or its nearby regions, it is understandable that the technique would be of scarce interest in practical applications. This is why we do not consider the work by Deng and Guo further. Regarding the remaining methods of eigenfunction representation, LPA and discrete diffraction, the remarks in Section 3 can be extended in the present context

166

Michele Marrocco

of radially polarized beams. Specifically, the method based on the eigenfunctions of the Hankel transform has a very intricate structure and shows a final formula with two series that defies fast numerical description of the fields (for instance, an incredible large number of terms: 29 × 9 = 261, is ¨ ok, ¨ used to truncate the series of the numerical example of Sherif and Tor 2005). These are more efficiently calculated according to the LPA elaboration or the method of discrete diffraction. Between the two techniques, Tables 1 and 4 and Figures 4 and 5 show that discrete diffraction ensures better numerical efficiency; this is replicated in the case of radial polarization by virtue of the formal equivalence between Eqs. (7) and (8) and Eq. (4) (the kernel of these integrals has always the same structure of trigonometric functions multiplying a phasor and a Bessel function). Everything considered, it seems that the most reasonable comparison can be arranged among the Lax series, the CSPSW theory, and discrete diffraction. Note, however, the comparison can only be established for moderate NA values that do not require the high spatial resolution imposed by the requirements of nano-optics (Novotny and Hecht, 2006). At higher NA values, we instead demonstrate that EDT methods (specifically, DNS simulation and discrete diffraction) are by far more accurate. To proceed with the comparison, let us consider the case of ε = 0.4 that was examined in Salamin’s work (2009). Plots of the squared absolute values of the radial and axial fields calculated according to the two methods conceived on the basis of the Helmholtz equation are shown in Figure 6 relative to the focal plane. As expected, the two approaches are quite similar as to the general prediction of the formation of the maxima and decay of the fields (actually, in Salamin’s work, differences in absolute amplitudes are claimed, but these are tolerated here in that we are more concerned with the relative spatial dependences). The situation changes dramatically at larger diffraction angles ε. Let us take the next example of ε = 0.8 (Figure 7) that still corresponds to one of the optical cases discussed by Salamin. If the Lax series seems to yield reasonable results, it is undeniable that the singularities of the CSPSW method alter the fields that become divergent at about kρ ∼ = 2.9. The skyrocketing behavior is an unphysical phenomenon and it is possible to conclude that the CSPSW serves no purpose for high-resolution studies. Apparently, the Lax series succeeds in the description of the general feature of a vanishing radial component for kρ → 0 and, at the same position, it predicts that the longitudinal component of the field reaches the maximum intensity. These effects are typical of radially polarized beams. Nonetheless, the plots are troubled by multiple oscillations for kρ > 2. These concentrated oscillations remain somewhat of a mystery because standard diffraction theory for the focusing of the fundamental Gaussian beam does not predict them (Novotny and Hecht, 2006; Youngworth

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

167

ε = 0.4

|Eρ|2

ε = 0.4

|Ez|2

0

2

4

6

8

10

kρ

FIGURE 6 Plot of the squared absolute values of the radial and longitudinal components of the fields obtained in agreement with the Lax series (lighter line) and the CSPSW approach (darker line) at the chosen diffraction angle ε = 0.4. The fields are normalized to their own maximum and are calculated in the focal plane (z = 0) so that the dependence on the dimensionless coordinate kρ is considered.

and Brown, 2000). To verify this point, the numerical solutions of Eqs. (7) and (8) are reported as continuous lines in Figure 8 for an NA value corresponding to the diffraction angle of Figure 7. Although the abovementioned partitioning of light between radial and axial components near the optical axis (kρ = 0) is replicated here, the remarkable oscillations of Figure 7 do not occur, as expected from the original work of Youngworth and Brown (2000). More importantly, the analytical treatment based on Eqs. (50) and (51) correctly captures the behavior of diffraction theory (see the points in Figure 8) and at the same time, the speed ratio R does not differ from the results of linear polarization (Tables 1 through 4, Figures 4 and 5).

5. CONCLUSIONS Numerous methods to deal with the modeling of tightly focused optical beams typical of highly resolved laser microscopy have been reviewed in reference to the two important cases of linear and radial polarization.

168

Michele Marrocco

ε = 0.8

|Eρ|2

ε = 0.8

|Ez|2

0

2

4

6

8

10

kρ

FIGURE 7 Plot of the squared absolute values of the radial and longitudinal components of the fields obtained in agreement with the Lax series (lighter line) and the CSPSW approach (darker line) at the chosen diffraction angle ε = 0.8. The Lax fields are normalized to their own maximum, whereas the diverging CSPSW fields are normalized to a maximum value chosen for graphical purposes. All the calculations refer to the focal plane (z = 0) so that the dependence on the dimensionless coordinate kρ is shown.

The former is rather common in practical experiments inasmuch as laser beams commercially available are ordinarily linearly polarized. The latter is attracting much interest for its promising characteristic of strong longitudinal field components that tend to decrease the focal spot size. The conclusions of this comparative research strengthen the role played by EDT as means to study the interaction between laser and matter at high NA. In this regard, the numerical technique known as discrete diffraction seems to outrun the other methods of numerical calculation and DNS simulations of linearly polarized beams. From a different perspective, analytical solutions of the Helmholtz equation for the vector potential have also been introduced by some authors who intended to simulate the profiles of radially polarized beams under tight focusing. However, it appears that these attempts are chiefly conditioned to optical microscopes of low NA. Furthermore, the accomplishments of EDT methods exceed the simplest case of linear polarization and more complex elaborations are exemplified through the correct reproduction of the profiles of radially polarized beams.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

169

Diffraction theory Eq. (14)

|Eρ|2

(a)

Diffraction theory Eq. (15)

|2

|Ez

0

10

20

30

kρ (b)

FIGURE 8 Plot of the squared absolute values of the radial and longitudinal components of the fields according to the standard diffraction theory (line) for the NA value equivalent to the diffraction angle of Figure 7. The squares show the results based on the analytical approach of Eqs. (50) and (51). As expected from the work of Youngworth and Brown (2000), the oscillations seen in Figure 7 are absent.

REFERENCES Abramowitz, M., & Stegun, I. A. (1965). Handbook of mathematical functions. New York: Dover. Adamson, P. (2004). High-aperture focusing systems: Control of light concentration in focal region by pupil filtering. Journal of Modern Optics, 51, 65–74. Agrawal, P., & Pattanayak, D. N. (1979). Gaussian beam propagation beyond the paraxial approximation. Journal of the Optical Society of America, 69, 575–578. Apostol, T. M. (1967). Calculus. New York: John Wiley & Sons. Arfken, G. (1985). Mathematical methods for physicists. San Diego, CA: Academic Press. Asatryan, A. A., Sheppard, C. J. R., & de Sterke, C. M. (2004). Vector treatment of secondharmonic generation produced by tightly focused vignetted Gaussian beams. Journal of the Optical Society of America B: Optical Physics, 21, 2206–2212. Barbero, S., & Marcos, S. (2008). Analysis of the optical field on the human retina from wavefront aberration data. Journal of the Optical Society of America A: Optics, Image Science, and Vision, 25, 2280–2285. Barton, J. P., & Alexander, D. R. (1989). Fifth-order corrected electromagnetic field components for a fundamental Gaussian beam. Journal of Applied Physics, 66, 2800–2802.

170

Michele Marrocco

Betzig, E., & Trautman, J. K. (1992). Near-field optics: Microscopy, spectroscopy, and surface modification beyond the diffraction limit. Science, 257, 189–195. Bomzon, Z., Biener, G., Kleiner, V., & Hasman, E. (2002). Radially and azimuthally polarized beams generated by space-variant dielectric subwavelength gratings. Optics Letters, 27, 285–287. Born, M., & Wolf, E. (1970). Principles of optics. New York: Pergamon. Cheng, J.-X., Volkmer, A., & Xie, X. S. (2002). Theoretical and experimental characterization of coherent anti-Stokes Raman scattering microscopy. Journal of the Optical Society of America B: Optical Physics, 19, 1363–1375. Couture, M., & Belanger, P.-A. (1982). From Gaussian beam to complex-source-point spherical wave. Physical Review A, 24, 355–359. D’Arcio, L., Braat, J. J. M., & Frankena, H. (1994). Numerical evaluation of diffraction integrals for apertures of complicated shape. Journal of the Optical Society of America A: Optics, Image Science, and Vision, 11, 2664–2674. Davism, L. W. (1979). Theory of electromagnetic beams. Physical Review A, 19, 1177–1179. Debye, P. (1909). Das Verhalten von Lichtwellen in der N¨ahe eines Brennpuktes oder iener Brennlinie. Annals of Physics, 30, 755–776. Deng, D., & Guo, Q. (2007). Analytical vectorial structure of radially polarized light beams. Optics Letters, 32, 2711–2713. Denk, W., Strickler, J. H., & Webb, W. W. (1990). Two-photon laser scanning fluorescence microscopy. Science, 248, 73–76. Dhayalan, V., & Stamnes, J. J. (1997). Focusing of electric-dipole waves in the Debye and Kirchhoff approximations. Pure and Applied Optics, 6, 347–372. Domke, K. F., & Pettinger, B. (2010). Studying surface chemistry beyond the diffraction limit: 10 years of TERS. ChemPhysChem, 11, 1365–1373. Dorn, R., Quabis, S., & Leuchs, G. (2003). Sharper focus for a radially polarized light beam. Physical Review Letters, 91, 233901 [4 pages]. Durant, S., Liu, Z., Steele, J. M., & Zhang, X. (2006). Theory of the transmission properties of an optical far-field superlens for imaging beyond the diffraction limit. Journal of the Optical Society of America B: Optical Physics, 23, 2383–2392. Dyba, M., & Hell, S. W. (2002). Focal spots of size λ/23 open up far-field florescence microscopy at 33 nm axial resolution. Physical Review Letters, 88, 163901 [4 pages]. ¨ ok, ¨ P. (2008). Inversion of the Debye-Wolf Foreman, M. R., Sherif, S. S., Munro, P. R. T., & Tor diffraction integral using an eigenfunction representation of the electric fields in the focal region. Optics Express, 16, 4901–4917. Freudiger, C. W., Min, W., Saar, B. G., Lu, S., Holtom, G. R., He, C., et al. (2008). Label-free biomedical imaging with high sensitivity by stimulated Raman scattering microscopy. Science, 322, 1857–1861. Gradshteyn, I. S., & Ryzhik, I. M. (1980). Table of integrals, series, and products. London: Academic Press. Gramotnev, D. K., & Bozhevolnyi, S. I. (2010). Plasmonics beyond the diffraction limit. Nature Photonics, 4, 83–91. Gravelsaeter, T., & Stamnes, J. J. (1982). Diffraction by circular apertures—1. Method of linear phase and amplitude approximation. Applied Optics, 21, 3644–3651. Grosjean, T., Baida, F., Adam, R., Guillet, J.-P., Billot, L., Nouvel, P., et al. (2008). Linear to radial polarization conversion in the THz domain using a passive system. Optics Express, 16, 18895–18909. Hayazawa, N., Inouye, Y., Sekkat, Z., & Kawata, S. (2002). Near-field Raman imaging of organic molecules by an apertureless metallic probe scanning optical microscope. Journal of Chemical Physics, 117, 1296–1301. Hopkins, H. H. (1943). The Airy disc formula for systems of high relative aperture. Proceedings of the Physical Society of London, 55, 116–128.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

171

Hopkins, H. H. (1957). The numerical evaluation of the frequency response of optical systems. Proceedings of the Physical Society B, 70, 1002–1005. Hopkins, H. H., & Yzuel, M. J. (1970). The computation of diffraction patterns in the presence of aberrations. Optica Acta: International Journal of Optics, 17, 157–182. Ignatovsky, V. S. (1919). Diffraction by a lens having arbitrary opening. Transactions of the Optical Institute of Petrograd, 1, paper IV. Ignatovsky, V. S. (1920). Diffraction by a parabolic mirror having arbitrary opening. Transactions of the Optical Institute of Petrograd, 1, paper V. Jackson, J. D. (1998). Classical electrodynamics. New York: Wiley. Jain, M., Lotsberg, J. K., Stamnes, J. J., Frette, ø., Velauthapillai, D., Jiang, D., et al. (2009). Numerical and experimental results for focusing of three-dimensional electromagnetic waves into uniaxial crystals. Journal of the Optical Society of America A: Optics, Image Science, and Vision, 26, 691–698. Jipson, V. B., & Williams, C. C. (1983). Two-dimensional modelling of an optical disk readout. Applied Optics, 22, 2202–2209. Kant, R. (1993). An analytical solution of vector diffraction for focusing optical systems. Journal of Modern Optics, 40, 337–347. Kant, R. (2000). Superresolution and increased depth of focus: An inverse problem of vector diffraction. Journal of Modern Optics, 47, 905–916. Karoui, A., & Moumni, T. (2009). Spectral analysis of the finite Hankel transform and circular prolate spheroidal wave functions. Journal of Computational and Applied Mathematics, 233, 315–333. Kino, G., & Corle, T. (1997). Confocal scanning optical microscopy and related imaging systems. San Diego, CA: Academic Press. Lax, M., Louisell, W. H., & McKnight, W. B. (1975). From Maxwell to paraxial wave optics. Physical Review A, 11, 1365–1370. Lerman, G. M., & Levy, U. (2008). Effect of radial polarization and apodization on spot size under tight focusing conditions. Optics Express, 16, 4567–4581. Lotsberg, J. K., Zhao, X., Jain, M., Dhayalan, V., Sithambaranathan, G. S., Stamnes, J. J., et al. (2005). Focusing of electromagnetic waves into a biaxial crystal, experimental results. Optics Communications, 250, 231–240. Loudon, R. (2000). The quantum theory of light. Oxford: Oxford University Press. Luo, H., Liu, S., Lin, Z., & Chan, C. T. (2007). Method for accurate description of a radially polarized Gaussian laser beam beyond the paraxial approximation. Optics Letters, 32, 1692–1694. Marrocco, M. (2009a). High-resolution microscopy with transition from continuous to discrete diffraction. Optics Communications, 282, 3869–3872. Marrocco, M. (2009b). Discrete diffraction for analytical approach to tightly focused electric fields with radial polarization. Optics Communications, 282, 3862–3868. Masters, B. R., & So, P. T. C. (2008). Handbook of biomedical nonlinear optical microscopy. Oxford: Oxford University Press. Messiah, A. (1962). Quantum mechanics. Amsterdam: North-Holland. Moore, I. C., & Cada, M. (2004). Prolate spheroidal wave functions, an introduction to the Slepian series and its properties. Applied and Computational and Harmonic Analysis, 16, 208–230. Novotny, L., Beversluis, M. R., Youngworth, K. S., & Brown, T. G. (2001). Longitudinal field modes probed by single molecules. Physical Review Letters, 86, 5251–5254. Novotny, L., & Hecht, B. (2006). Principles of nano-optics. Cambridge: Cambridge University Press. Panofsky, W. K. H., & Phillips, M. (1962). Classical electricity and magnetism. New York: Wiley. Pawley, J. B. (2006). Handbook of biological confocal microscopy. New York: Springer.

172

Michele Marrocco

¨ Quabis, S., Dorn, R., Eberler, M., Glockl, O., & Leuchs, G. (2000). Focusing light to a tighter spot. Optics Communications, 179, 1–7. Richards, B., & Wolf, E. (1959). Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system. Proceedings of the Royal Society of London Series A-Mathematical Physical and Engineering Sciences, 253, 358–379. Salamin, Y. I. (2006a). Fields of a radially polarized Gaussian laser beam beyond the par axial approximation. Optics Letters, 31, 2619–2621. Salamin, Y. I. (2006b). Accurate fields of a radially polarized Gaussian laser beam. New Journal of Physics, 8, 133 [16 pages]. Salamin, Y. I. (2007). Acceleration in vacuum of bare nuclei by tightly focused radially polarized laser light. Optics Letters, 32, 3462–3464. Salamin, Y. I. (2008). Corrigendum of Salamin (2006b). New Journal of Physics, 10, 069801 [2 pages]. Salamin, Y. I. (2009). Fields of a tightly focused radially polarized laser beam: The truncated series versus the complex-source-point spherical wave representation. New Journal of Physics, 11, 033009 [8 pages]. ¨ ok, ¨ P. (1997). Efficient calculation of electromagnetic diffraction in Sheppard, C. J. R., & Tor optical systems using a multipole expansion. Journal of Modern Optics, 44, 803–818. ¨ ok, ¨ P. (2005). Eigenfunction representation of the integrals of the DebyeSherif, S. S., & Tor Wolf diffraction formula. Journal of Modern Optics, 52, 857–876. ¨ ok, ¨ P. (2008). Eigenfunction expansion of the electric fields Sherif, S. S., Foreman, M. R., & Tor in the focal region of a high numerical aperture focusing system. Optics Express, 16, 3397– 3407. Shoham, A., Vander, R., & Lipson, S. G. (2006). Production of radially and azimuthally polarized polychromatic beams. Optics Letters, 31, 3405–3407. Sick, B., Hecht, B., & Novotny, L. (2000). Orientational imaging of single molecules by annular illumination. Physical Review Letters, 85, 4482–4485. Stamnes, J. J., Spjelkavik, B., & Pedersen, H. M. (1983). Evaluation of diffraction integrals using local phase and amplitude approximations. Optica Acta, 30, 207–222. Stamnes, J. J. (1986). Waves in focal regions. Bristol, UK: Adam Hilger. Steidtner, J., & Pettinger, B. (2008). Tip-enhanced Raman spectroscopy and microscopy on single dye molecules with 15 nm resolution. Physical Review Letters, 100, 236101 [4 pages]. ¨ ok, ¨ P., Hewlett, S. J., & Vargas, P. (1997). On the series expansion of high-aperture, vectorial Tor diffraction integrals. Journal of Modern Optics, 44, 493–503. Van de Nes, A. S., Braat, J. J. M., & Pereira, S. F. (2006). High-density optical data storage. Reports on Progress in Physics, 69, 2323–2363. Volkmer, A., Cheng, J.-X., & Xie, X. S. (2001). Vibrational imaging with high sensitivity via epidetected coherent anti-Stokes Raman scattering microscopy. Physical Review Letters, 87, 023901 [4 pages]. Watson, G. N. (1952). A treatise on the theory of bessel functions. Cambridge: Cambridge University Press. Westphal, V., Rizzoli, S. O., Lauterbach, M. A., Kamin, D., Jahn, R., & Hell, S. W. (2008). Video-rate far-field optical nanoscopy dissects synaptic vesicle movement. Science, 320, 246–249. Wolf, E. (1959). Electromagnetic diffraction in optical systems. I. An integral representation of the image field. Proceedings of the Royal Society of London Series A-Mathematical Physical and Engineering Sciences, 253, 349–357. Yan, S., & Yao, B. (2008). Accurate description of a radially polarized Gaussian beam. Physical Review A, 77, 023827 [4 pages]. Yanai, A., & Levy, U. (2009). Plasmonic focusing with a coaxial structure illuminated by radially polarized light. Optics Express, 17, 924–932.

Methods for Vectorial Analysis and Imaging in High-Resolution Laser Microscopy

173

Yew, E. Y. S., & Sheppard, C. J. R. (2007a). Tight focusing of radially polarized Gaussian and Bessel-Gauss beams. Optics Letters, 32, 3417–3419. Yew, E. Y. S., & Sheppard, C. J. R. (2007b). Second harmonic generation polarization microscopy with tightly focused linearly and radially polarized beams. Optics Communications, 275, 453–457. Yoshida, A., & Asakura, T. (1974). Electromagnetic field in the focal plane of a coherent beam from a wide angular annular aperture system. Optik, 40, 322–331. Yoshida, A., & Asakura, T. (1975). Electromagnetic field near the focus of a Gaussian beam. Optik, 41, 281–292. Youngworth, K. S., & Brown, T. G. (2000). Focusing of high numerical aperture cylindricalvector beams. Optics Express, 7, 77–87. Zipfel, W. R., Williams, R. M., & Webb, W. W. (2003). Nonlinear magic: Multiphoton microscopy in the biosciences. Nature Biotechnology, 21, 1369–1377. Zumbusch, A., Holtom, G. R., & Xie, X. S. (1999). Three-dimensional vibrational imaging by coherent anti-Stokes Raman scattering. Physical Review Letters, 82, 4142–4145. Zumofen, G., Mojarad, N. M., Sandoghdar, V., & Agio, M. (2008). Perfect reflection of light by an oscillating dipole. Physical Review Letters, 101, 180404 [4 pages].