The Debye-Wolf electromagnetic diffraction integral is now routinely used to describe focusing by high numerical (NA) lenses. We obtain an eigenfunction expansion of the electric vector field in the focal region in terms of Bessel and generalized prolate spheroidal functions. Our representation has many optimal and desirable properties which offer considerable simplification to the evaluation and analysis of the Debye-Wolf integral. It is potentially also useful in implementing two-dimensional apodization techniques to synthesize electromagnetic field distributions in the focal region of a high NA lenses. Our work is applicable to many areas, such as optical microscopy, optical data storage and lithography.
© 2008 Optical Society of America
The calculation of the field distribution in the focal region of a high numerical aperture (NA) lens is an important problem because of the many applications that use tightly focused beams, such as optical microscopy, optical data storage and lithography. When the NA of a lens exceeds approximately 0.5 the field distribution in the focal region can no longer be described in an accurate manner by the scalar diffraction theory so the use of an electromagnetic diffraction theory is necessary. The electromagnetic field distribution at any point in the focal region of a diffraction-limited, high NA incident field of arbitrary uniform polarization is typically obtained by the Debye-Wolf vector integral [1–3].
To simplify the Debye-Wolf diffraction integral the formulae of Watson , Gradshteyn and Ryzhik  and Agrawal and Pattanayak , were used by Török et al. to express both the infocus and defocus terms by means of a series of analytic functions . Kant reported a series expansion of the diffraction integrals using Gegenbauer polynomials . Sherif and Török  reported an eigenfunction representation of the so-called I integrals of Richards and Wolf. In reality these expansions do not reveal anything about the physical nature of the problem – they merely provide a simplified way to calculate the diffraction integrals when computational time is deemed to be of significance.
Braat et al. obtained the field components in the focal region as a series using Nijboer-Zernike functions , while Sheppard and Török obtained the field components as a multipole expansion . These two expansions are more physical than those listed above because the Nijboer-Zernike expansion aims at obtaining formulae where the incident and focused fields are represented in terms of aberration functions. This representation has immediate significance in the study of realistic focusing systems with aberrations present. Nevertheless, the multipole expansion may be regarded as the one with the most physical content because the focused field is represented in terms of physically existing multipole functions.
In this paper, we obtain an expansion of the electric field components in the focal region of a high NA lens in terms of Bessel functions and the generalized prolate spheroidal functions which are eigenfunctions of the two-dimensional finite Hankel transform. As discussed in Section 5, our expansion has many optimal, desirable and physical properties. Its optimal properties include being expressed in its simplest possible operator form (diagonal) and having maximum energy packing properties. Its desirable properties include separability of its component functions in cylindrical coordinates, and fast convergence in the azimuthal, radial and axial directions. Physically the dominant modes are closely connected with the resolution of the optical system, whilst the generalized prolate spheroidal functions are themselves the in-focus eigen field distributions for focusing by a lens. Such functions are hence important since they are unchanged by the focusing operation (to within a scale factor). Our expansion therefore provides a simple and natural way to carry out both forward and inverse analysis of high NA focusing systems. In contrast to our earlier work  which allowed only for one-dimensional (1-D) apodization techniques , our current expansion could be used to implement two-dimensional (2-D) apodization and masking techniques to synthesize fields at the focal region of a high NA focusing system .
2. Focusing of high numerical aperture beams
The electric field distribution at any point, p(xp,yp,zp), in the focal region of a diffraction-limited, high numerical aperture lens with an arbitrary polarized uniform incident illumination as shown in Fig. 1, is given by the Debye-Wolf integral
where k=2πn/λ is the wavenumber of the illumination, a(sx, sy) is the strength vector of a ray at the Gaussian reference sphere converging to the focal point, sx, sy and sz are the Cartesian components of the unit vector s along this ray and Θ is the solid angle associated with all the rays which reach the image space through the exit pupil of the lens [1–3].
In the case of aplanatic focusing with a linearly polarized incident illumination, the diffraction integral, Eq. (1), can be re-written as follows. For convenience, spherical coordinates r, θ, ϕ(r>0, 0≤θ<π, 0≤ϕ<2π) with the polar axis θ=0 in the z direction and with the azimuth ϕ=0 containing the linearly polarized incident vector are introduced. The unit vector s along a ray and the point p(xp,yp, zp) in the focal region can be expressed as
where A is a constant, α=θ max and ρp=rpsinθp, zp=rp cosθp. We note that the dummy variables of the double integrals in Eq. (4) are given in spherical coordinates while the field distributions are given in cylindrical coordinates.
3. Eigenfunction representation of the fields of the Debye-Wolf focusing integral
To avoid repetition of mathematics we show the details of obtaining the eigenfunction representation for only one field component, Ex in Eq. (4a). Afterwards we show the final result which we can obtain in a similar manner for Ey in Eq. (4b) and Ez in Eq. (4c). We start by noting that the defocus term, exp[ikzpcosθ], in Eq. (4) could be written as 
We could simplify Eq. (6) through a coordinate transformation u=sinθ which results in and dθ=du/cosθ.
Substituting this transformation into Eq. (6), we obtain
where α′=sinα. Since exp(iθ)=cosθ+i sinθ, cosmθ=Tm (cosθ) and sin[(m+1)θ]=sinθUm (cosθ), where Tm(·) and Um(·) are Chebyshev polynomials of the first and second types respectively , we can write
The function axm(u,ϕ) is space-limited, so we could expand it in terms of generalized circular prolate spheroidal functions, which are eigenfunctions of the two-dimensional finite Hankel transform [16, 17]
and c is a parameter equal to or larger than the radial space-bandwidth product of axm(u,ϕ) for all values of 0<ϕ≤2π .
Using the identity exp[ikρpucosϕ]=∑∞ q=-∞iqJq(kρpu)exp(iqϕ)  and the fact that cos(ϕ-ϕp)=cos(ϕp-ϕ) and ∫2π 0 exp(i(N-q)ϕ)dϕ=δNq where δNq is the Kronecker delta function, we obtain
We also further use the fact that for N<0 we could write JN(x)=(-1)N J |N| (x) and iN(-1)N=i-N=i|N| to obtain
Using the eigen equation
where Aym,N,n and Azm,N,n are the expansion coefficients of
in terms of generalized prolate spheroidal functions. Equation (16) is our eigenfunction expansion of the electric fields in the focal region of a high NA focusing system with linearly polarized incident beam and hence it constitutes the main result of this paper.
4. Arbitrary polarization
This result can be easily extended to describe focusing of beams with an arbitrary polarization distribution through amplitude and phase masking optics. Incorporating these additional factors merely modifies the components of the strength vector from which the expansion coefficients are calculated using Eq. (11) and its y and z counterparts. The general form for the strength vector takes the form
where g(u,ϕ) and Ψ(u,ϕ) describe the amplitude and phase variation introduced to the incident field distribution
by an apodization and phase mask respectively and A is the matrix given by 
which describes the action of the lens and maps the field to the Gaussian reference sphere.
We note that this general formalism can readily be used to account for polarization apodization which is becoming increasingly popular in the literature.A subset of these is the focusing of vortex beams. The generalized expressions for these were obtained by Török and Munro .
5. Physical properties of the eigenfunction expansion
On examining our expansion of the Debye-Wolf integral, Eq. (16), we notice many desirable, optimal and physical properties. Firstly, the component functions Jm (kzp), exp(iNϕp), Φ|N|,n (α′ρp/ρ pmax) are separable in cylindrical coordinates which could simplify analysis involving fields in the focal region of a high NA focusing system.
Secondly, the double integral in Eq. (12) represents a finite Hankel transform of the generalized prolate spheroidal functions. As a linear operator, this transform takes its simplest possible form, i.e., diagonal, through its eigen representation given in Eq. (16a).
Thirdly, one of the defining properties of the generalized prolate spheroidal functions is that Φ|N|,0(r,c) maximize the fractional energy within a circular region of radius r 0 over the class of all bandlimited functions . Thus our summation over N has maximum energy packing properties in both radial and azimuthal directions when n=0. More generally the eigenvalues, λ |N|,n, describe the fractional energy within the circular region for each generalized prolate spheroidal function.
Fourthly, we note, as shown in Fig. 2(a), that the eigenvalues of the generalized prolate spheroidal functions decrease monotonically to very small values, compared to their initial values, after certain orders |N|≥|N 0|, n≥n 0. We also note, as shown in Fig. 2(b), that Bessel functions of the same argument, but of increasing orders, decrease to very small values, say 10-4, after some order m≥m 0. Therefore our expansion demonstrates fast convergence in the azimuthal, radial and axial directions.
Our eigen representation of the field in the focal region of a focusing system requires separate expansions for each field component the basis functions of which have physical significance since they are unaltered by the focusing operation. In Fig. 3 we plot the in-focus (zp=0) eigenfunctions which reduce to the generalized prolate spheroidal functions, whilst defocused eigenfunctions are further modulated by a Bessel function dependent on the axial coordinate. For low numerical aperture systems the polarization properties of light become less important often allowing a scalar treatment to be used. Under such circumstances we can interpret the field distributions shown in Fig. 3 as the true eigenfunctions of the focusing operation. At higher numerical apertures the distributions shown are not strictly eigenfunctions of the Debye-Wolf integral since they are scalar functions, however they do remain eigenfunctions on a component-wise basis i.e. if x-polarized light with amplitude distribution given by the (N,n)th order were focused, the x component of the output field would also have the same (although scaled) distribution.
Furthermore, from the shape of our eigenfunctions we see that higher order functions contain a larger fraction of energy in the sidelobe structure, which itself becomes progressively more complicated as the order increases. We also see that only the N=0 modes have a central focal spot. Put another way the higher order functions contain higher spatial frequencies however energy is not concentrated as efficiently into the central region. Consequently more complicated masking optics will, in general, require more terms to be calculated in the eigenfunction expansion to accurately determine the field in the focal region.
Finally in Fig. 4 we show the variation of the eigenvalues, λ |N|,n, with the numerical aperture of the focusing system. It can be seen that for a given order (N,n) the eigenvalue decreases as the NA decreases. This dependence means that higher orders are energetically less significant in the focused distribution and hence the dominant orders are those with lower spatial frequencies or equivalently are less tightly packed. The dominant modes in the focal region hence provide a further means to determine the resolution of the optical system.
6. Numerical examples
To verify the validity of our eigenfunction expansion, we use Eq. (16) to obtain the field distributions in the Gaussian focal plane and a defocused plane, and compare them to the corresponding distributions obtained by the direct evaluation of the Debye-Wolf integral. Assuming the availability of tabulated values or computer routines to evaluate Bessel functions, Chebyshev polynomials  and circular prolate spheroidal functions [21, 22], the main task to evaluate our eigenfunction expansion is to determine the space-bandwidth parameter c and suitable finite limits for the three summations in Eq. (16). The parameter c has to be equal to or larger than the radial space-bandwidth product of the functions a (·) m (u,ϕ) in Eq. (9) and Eq. (17) for all values 0<ϕ≤2π. In addition, c=α′kρ pmax, so for a given value of the NA, c determines the radial field of view, ρ pmax, in the plane of interest. For the following numerical examples, we found that c=20 satisfied these two requirements.
From Fig. 2(a), we note that for c=20, |N 0|=23 and n 0=8 are appropriate limits for the summations with respect to N and n, respectively. From Fig. 2(b), we note that for a defocus distance z=λ, m 0=14 is an appropriate limit for the summation with respect to m.
In Fig. 5 and Fig. 6 we show the pointwise relative error between the optical distributions due to a linearly polarized incident beam, at the Gaussian focal plane and at a defocused plane, z=λ, respectively, obtained by evaluating our eigenfunction expansion, Eq. (16), with finite summation limits and by direct integration. The actual calculated optical distribution is also shown in the insets.
To confirm that these optical distributions are indeed equal to the ones obtained by direct evaluation of the Debye-Wolf integral, we define an overall percentage error factor ε given by
where I direct k,l and I expansion k,l are the intensities at point (k, l) obtained by evaluating the Debye-Wolf integral and our eigenfunction expansion, respectively, and (K,L) are the maximum number of points in the plane of interest. On applying Eq. (21) to the above numerical examples, we found ε equal to 0.0253% and 0.065% at the Gaussian focal plane (Fig. 5) and at a defocus distance z=λ (Fig. 6), respectively. This overall error is very small, bearing in mind the relatively small number of terms used to evaluate Eq. (16). Furthermore it is found the variation of the total percentage error as the NA is increased from 0 to 0.966 is of order 0.003% and is hence negligible.
In this paper, we obtained an eigenfunction expansion of the electric fields in the focal region of a high NA system in terms of Bessel and generalized circular prolate spheroidal functions. We demonstrated, by summing a relatively small number of terms, the fast and accurate convergence of this expansion. Our expansion has many optimal, desirable and physical properties. Its optimal properties include having the simplest possible (diagonal) form and maximum energy packing, in both the radial and azimuthal directions. Its desirable properties include separability of its component functions in cylindrical coordinates, and fast convergence in the azimuthal, radial and axial directions. Thus it provides a simple and natural way to carry out both forward and inverse analysis of high NA focusing systems. Physically the generalized prolate spheroidal functions have also been shown to be the eigenfunctions of the focusing operation. This expansion could also be used to implement 2-D apodization techniques to synthesize arbitrary fields at the focal region of a high NA focusing system in which the dominant modes present give insights into the resolution of the system. By arbitrary we mean bandlimited fields that satisfy the wave equation.
We would like to thank Paul Abbott and Peter Falloon from the Physics Department of the University of Western Australia for supplying the Mathematica code to compute Slepian’s Generalized Spheroidal functions. S. Sherif thanks Costel Flueraru, Shoude Chang and Linda Mao for their support through the Optics Group at the Institute for Microstructural Sciences, National Research Council Canada.
References and links
1. B. Richards “Diffraction in systems of high relative aperture,” in Astronomical Optics and Related Subjects, Z. Kopal, ed. (North Holland Publishing Company, 1955), pp. 352–359.
2. E. Wolf, “Electromagnetic diffraction in optical systems I. An integral representation of the image field,” Proc. Roy. Soc. London A , 253, 349–357 (1959). [CrossRef]
3. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. Roy. Soc. London A , 253, 358–379 (1959). [CrossRef]
4. G. N. Watson, A Treatise on the Theory of Bessel Functions, 2nd Ed. (Cambridge University Press, 1952). [PubMed]
5. I. S. Gradshtyen and I. M. Ryzhik, Table of Integrals, Series and Products (Academic Press, 1980).
6. G. P. Agrawal and D. N. Pattanayak, “Gaussian beam propagation beyond the paraxial approximation,” J. Opt. Soc. Am. 69, 575–578 (1979). [CrossRef]
7. R. Kant, “An analytical solution of vector diffraction for focusing optical systems with seidel aberrations” J. Mod. Opt. 40, 2293–2310 (1993). [CrossRef]
8. P. Török, S.J. Hewlett, and P. Varga, “On the series expansion of high aperture, vectorial diffraction integrals,” J. Mod. Opt. 44, 493–503 (1997). [CrossRef]
9. S. S. Sherif and P. Török, “Eigenfunction representation of the integrals of the Debye-Wolf diffraction formula,” J. Mod. Opt. 52, 857–876 (2005). [CrossRef]
10. J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, and A. S. van de Nes, “Extended Nijboer-Zernike representation of the vector field in the focal region of an aberrated high-aperture optical system,” J. Opt. Soc. Am. A 20, 2281–2292 (2003). [CrossRef]
11. C. J. R. Sheppard and P. Török, “Efficient calculation of electromagnetic diffraction in optical systems using a multipole expansion,” J. Mod. Opt. 44, 803–818 (1997). [CrossRef]
12. S. S. Sherif and P. Török, “Pupil plane masks for super-resolution in high-numerical-aperture focusing,” J. Mod. Opt. 51, 2007–2019 (2004).
13. M. R. Foreman, Department of Physics, Imperial College London, Prince Consort Rd. London SW7 2BZ, United Kingdom, S. S. Sherif, P.R.T. Munro, and P. Török are preparing a manuscript to be called “Inversion of the Debye-Wolf diffraction integral using an eigenfunction representation of the electric fields in the focal region”
14. G. Arfken, Mathematical Methods for Physicists, 3rd Edition (Academic Press, 1985).
15. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions (Dover, 1970).
16. D. Slepian, “Prolate spheroidal wave functions, Fourier analysis and uncertainty- IV: Extensions to many dimensions; generalized prolate spheroidal functions,” Bell Sys. Tech. J. 43, 3009–3057 (1964).
17. J. C. Heurtley, “Hyperspheroidal functions-optical resonators with circular mirrors,” in Proceedings of Symposium on Quasi-Optics, J. Fox, ed., (Polytechnic Press, 1964), pp. 367–371.
18. B. R. Frieden, “Evaluation, design and extrapolation methods for optical signals, based on the prolate functions,” Progress in Optics 9, E. Wolf, ed., (Pergamon, 1971), pp. 311–407.
19. P. Török, P.D. Higdon, and T. Wilson, “On the general properties of polarised light conventional and confocal microscopes,” Opt. Commun. 148300–315 (1998) [CrossRef]
20. S. Zhang and J. Jin, Computation of Special Functions (Wiley, 1996).
21. P. E. Falloon, Hybrid Computation of the Spheroidal Harmonics and Application to the Generalized Hydrogen Molecular Ion Problem (University of Western Australia, 2001). [PubMed]