We present an algorithm for calculating the field distribution in the focal region of stratified media which is fast and easy to implement. Using this algorithm we study the effect on the electric field distribution of an air gap separating a solid immersion lens and a sample, where we analyse the maximum distance for out-of-contact operation. Also, we study how the presence of a metallic substrate affects the field distribution in the focal region; the interference effects of the reflected field could be used as an alternative for 4Pi-microscopy.
©2004 Optical Society of America
For a sensible interpretation of the image formation and information retrieval in optical far- and near-field microscopy, lithography and optical recording, it is essential to have a good understanding of the field used for illumination of the sample or data structure. Due to the demand for an ever increasing resolution, the wavelength of the used light decreases to the ultraviolet and the numerical aperture of the imaging systems used for illumination and imaging has approached its theoretical limit. The field distribution in the focal region of a homogeneous medium can be obtained with the well-known diffraction integrals given by Ignatowsky or Richards and Wolf [1, 2]. However, in general the sample is placed in a multilayered environment, where one or more medium transitions should be taken into account. This has been done previously for the case of two media in [3–7], and for a multilayer system in [8–11]. Although these theoretical descriptions of the multilayer systems are fully vectorial and valid for a high numerical aperture imaging system, they have minor disadvantages regarding their validity domain or are not easy to interpret. The description discussed in  concentrates on calculation of the field distribution in the second layer and is based on concepts of ray optics using thin-film matrices; that of  is restricted to the calculation of the field distribution in the last layer; and that of  uses a complete set of basis functions to calculate the field distribution in the focal region but is limited to a specific set of aberrations.
In this paper we use a different method to obtain the formulae which describe the field distribution in the focal region including stratified media. Our approach leads to an easy and transparent way to take into account the changes in magnitude and orientation of the electric field vector due to the transition through the lens system and through the assembly of layers. In contrast with earlier analyses including multilayered assemblies, we use cylindrical coordinates both in the exit pupil and the focal region which leads, after some manipulations, to an appreciable simplification of the calculation of the transmitted and reflected field components. In the case of circular symmetry we are left with one single integral to be evaluated numerically; in the case of arbitrary symmetry we can profit from the fact that, in general, the azimuthal integration shows a quicker convergence than the radial one. In this way, the method can also yield a gain in computation time although this aspect has not been quantified. We would like to point out that the results calculated with the algorithm presented in this paper are easy to interpret because of our uniform treatment of s- and p-polarization in a Cartesian basis, in contrast with the use of special matrices required in . In the second section of this paper, the theoretical description of the algorithm is presented, and two applications are discussed in the third section. In the first example, we consider an immersion lens which operates slightly out of contact with the sample. In the second example, we show how the field distribution used to illuminate the sample can be shaped by adding a metallic layer onto its substrate, which could be considered as an alternative for 4Pi-microscopy.
In order to obtain the field distribution in the focal region embedded in a multilayer system, we start our derivation with the field distribution in the focal region of a single homogeneous medium, which formulae were first described by Ignatowsky  and rederived by Richards and Wolf . These diffraction integrals are valid in the Debye approximation, therefore our point of observation should be not too close to the spherical shell over which the integration takes place. In this paper, we consider an aplanatic imaging system, as depicted in Fig. 1. Furthermore, we assume a monochromatic time-harmonic electric field 𝓔(r, t)=Re[E(r)exp(-iωt)] and magnetic field 𝓑(r, t)=Re[B(r)exp(-iωt)]. Also, we consider only nonmagnetic materials, µ=µ 0. The absolute value of the propagation vector in vacuum is defined as =|k 0|2≡µ 0 ε 0 ω 2 with the wavelength in vacuum λ=2π/k 0. The field distribution in the focal region within a single homogeneous medium is given by 
where the integration domain Ω is restricted to the exit pupil of the imaging system. The vectorial field distribution on the exit pupil, a spherical shell with radius R, is proportional to a(kx,ky ). We introduce a set of cylindrical coordinates in the exit pupil k=(kr,kϕ,kz ), and in the focal region r=(r,ϕ,z), which changes Eq.(1) in
For a comparison of our integration variables (kr,kϕ ) and those used in , (θ,ϕ), the relationships kϕ =ϕ, kr /k=sinθ and kz/k=cosθ hold. The longitudinal component of the propagation vector is defined as kz =(k 2- )1/2; in order to make sure that the radiation condition is fulfilled, only exponentially decaying contributions should be considered, i.e., Im(kz )≥0. Since we consider an aplanatic imaging system, the transition from the entrance pupil, a disc defined by Ω′, to the exit pupil, a spherical shell defined by Ω, can be considered as a rotation of the field vector around the angular axis (ϕ̂) in a cylindrical basis (r̂, ,ẑ). This lens operation is described by the propagation matrix M and derived in Appendix A. Furthermore, the field distribution on the exit pupil should satisfy Abbe’s sine condition , yielding a squareroot multiplication factor (kz /k)1/2. Taken all these considerations into account, we rewrite the diffraction integral as
The integration domain Ω′ is restricted to 0≤kϕ <2π and 0≤kr <NAk 0, where NA denotes the numerical aperture of the imaging system. The field distribution in the entrance pupil is denoted as E 0(kr,kϕ ), scalar multiplication with the lens operator yields the field distribution in the exit pupil E 1(kr,kϕ )=M·E 0(kr,kϕ ); if the vector length R and the sine condition are taken into account, we retrieve the original exit pupil distribution a(kr,kϕ )=R(kz /k)1/2 E 1(kr,kϕ ).
Finally, we are ready to include the multilayer character of the formalism. A subscript i is given to each quantity that depends on the medium in which it is located. Since the addition of layers gives rise to reflections at the interface transitions, we have to introduce a set of directionally dependent propagation matrices for the medium of observation i, where the positive or negative sign denotes forward or backward traveling waves, respectively. The resulting electric field distribution in medium i is given by
The propagation matrices include both the lens operation and the geometry of the focal region, and are derived in Appendix A. Note that for a single homogeneous medium, for which the total amount of layers N=1, we find for the propagation matrices =M and =0, resulting in the same expression as in Eq. (3). The propagation matrices are defined as
where and are the effective transmission and reflection coefficients in medium i, respectively, and are derived in Appendix B. For a multilayer system, the cylindrical basis is again a natural choice, since there is a clear separation of the transverse electric (TE) component of the electric field, associated with the angular direction, and the transverse magnetic (TM) component of the electric field, associated with both the radial and longitudinal directions. The coefficients include different contributions for the TE and TM polarisations. The purely transversal nature of electromagnetic radiation is visible from the zeros in the third column of the matrix. The integer n=0,1,2 in the terms cosnkϕ and sinnkϕ relates to the angular dependence for the integration over kϕ . Since the effective transmission and reflection coefficients are independent of kϕ , we can easily decompose the entrance pupil distribution in sineand cosine-functions and perform the integral over the angular coordinate analytically, using the integral representation of the Bessel function of the first kind 
The remaining integral over the radial variable kr can be solved using a numerical integration scheme. The convergence of the numerical integration is an important issue and should be carefully analysed, as was discussed in detail in . Since the derivative of the reflection and transmission coefficients as a function of kr can be discontinuous, these points should be identified and the integration domain split up accordingly. To make sure that the convergence constraint is met, we have used a Gauss-Kronrod integration routine  which splits the integration domain in half when the difference in the sum for n-points and for (n+m)-points is not smaller than a specified accuracy.
In the field of optical microscopy, many different systems and types of illumination have been studied and some effects relevant for our examples, such as the near field contributions , the polarisation-dependent focal shift , multiple reflections , solid immersion systems [18,19], the dependence of the amplitude of the longitudinal component on the medium and the shaping of the field distribution in the focal region  have been described in previous papers. These studies are mostly restricted to the field transmitted to the last medium. In this paper, we present two practical examples; a solid immersion lens and the illumination of a photo-resist layer, where the field distribution is calculated in each individual layer, so that every transition can be analysed with detail including near-field and multiple reflection effects.
3.1. Solid immersion lens operated out-of-contact
As a first example, we study the field distribution in the focal region of a solid immersion lens separated from the sample by an air gap. For a good understanding of the problem, we will begin with a study of two preliminary configurations. Since the TE polarisation is associated with the angular direction () and TM polarisation is associated with both the radial and longitudinal direction (r̂,ẑ), the contributions do not mix and we can study both distributions simultaneously [21–23]. Therefore, we use an equally weighted combination of radial and azimuthal polarisation states with a uniformly distributed amplitude as illumination.
The first configuration we consider concerns a lens with a numerical aperture of NA=0.933 in air and a glass substrate with refractive index n 2=1.5 placed in the focal region. The obtained field distributions in the focal region for the radial, the azimuthal and the longitudinal components are shown in Figs. 2(a), (b) and (c), respectively. Note the different colormap for each component. The medium is located at a distance of two wavelengths in front of the geometrical focus (the origin) and results in a significant defocus and the appearance of interference fringes in the first medium. When we only consider the contribution of the radial and the longitudinal components to the intensity distribution, the transversal width of the spot, located at the maximum value along the optical axis in the second medium, is almost equal to the diffraction limit λ/NA, which is only influenced by the angular dependence of the transmission coefficient. The diffraction limit in the first medium is the same as in the second medium since the numerical aperture remains constant, although the geometrical aperture scales with n 1/n 2. This can be expected from the boundary conditions derived from Maxwell’s equations, which makes sure that the transversal component is continuous over the medium transition; note that the longitudinal component scales with ε 1/ε 2.
The second configuration we consider concerns an immersion lens with numerical aperture NA=1.4, as has been studied in [6, 7]. The transition is now from glass (n 1=1.5) to air, the geometrical aperture NA/n1 is identical to that of the first configuration. The field distribution obtained in the focal region is shown in Fig. 3. Here, the optical focus lies in front of the geometrical focus. When we only consider the contribution of the radial and the longitudinal components to the intensity distribution, the transversal width of the spot, located at the maximum value along the optical axis in the second medium, is slightly bigger than that without the medium transition. The spot width is influenced by three effects. The first effect is caused by the angle dependent transmission coefficients; since the larger angles are more favoured than smaller angles, we would expect a decrease of the spot width. The second effect is caused by the multiplication of the longitudinal component by a factor ε 1/ε 2; since this affects the relative strength of both contributions to the spot width, the spot width should be decreased. However, a third effect dominates, since the maximum numerical aperture in the second medium is limited to NA=1; the higher spatial frequencies do not contribute to the formation of the spot, therefore the width is increased. Again we also observe fringes caused by interference of the incident and reflected field in the first medium.
In the second configuration, the reflected focus in the first medium is more apparent than in the first configuration and there is a clear distinction between the radial and azimuthal components. The focal region of the reflected focus extends over a larger region for the radial component than for the azimuthal component, which is caused by the stronger contribution of the higher spatial frequencies. To simplify the interpretation, we have plotted in Fig. 4 the reflection and transmission coefficients as a function of the angle for the two configurations. In the second configuration, the maximum frequency transmitted to the air medium is, of course, k 0, which corresponds to an angle of incidence of θ=41.8°, known as the critical angle. The critical angle is the angle above which all waves are totally reflected and the components transmitted to the second medium are exponentially decaying evanescent waves.
After the discussion of various important effects in the previous configurations, such as the near field contributions, the polarisation-dependent focal shift, multiple reflections, the dependence of the amplitude of the longitudinal component on the medium and the shaping of the field distribution, we are ready to treat the example of an air gap between a solid immersion lens and a substrate for the case of radially polarised light. In Fig. 5 we show, for an air gap of (a) , (b) , (c) and (d) λ, the square root of the time averaged electric energy distribution in the focal region when the air gap moves from z=λ behind the geometrical focus to z=-2λ in front of the geometrical focus in 4 different movies. Since the effect of the air gap for the transmitted far field, i.e., no evanescent waves are taken into account, is not dependent on its position, the far-field distribution in the last medium remains stationary. The full width at half maximum (FWHM) of the axial maximum in the third medium are 0.52λ, 0.64λ, 0.96λ and 1.16λ for air gap values of , , and λ, respectively. The FWHM of the spot in glass for spatial frequencies limited to NA=1.0 is equal to 1.20λ. A small reduction of the spot size can be attributed to the effect of the angle-dependent transmission through the air gap, since careful analysis shows that the lower frequencies are slightly more attenuated than the higher frequencies. But the main effect which causes a reduction of the spot size should be attributed to a residual amplitude associated to frequencies surpassing the critical angle, which still propagate to the glass substrate for an air gap as large as .
3.2. The use of interference effects in microscopy due to reflecting surfaces
As a second example, we show how interference effects can be used to obtain a smaller lateral and axial spot size in the focal region. We consider radially polarised light to illuminate a lens with NA=0.95; the centre of the lens (NA<0.80) is blocked. The object placed in the focal region consists of a photo-resist layer (n 2=1.5), spun on an aluminum layer (n 3=0.49+4.87i) on top of a glass substrate (n 4=n 2). The metallic layer reflects the light which results in an interference pattern. From the boundary conditions at the interface, it follows that the reflected radial component (parallel to the interface) is 180° out-of-phase with the longitudinal component, which yields a shift of half a wavelength (λ=400nm). Because of the high value of the reflection coefficient, the local maxima in front of the photo-resist layer, see Fig. 6(a), are alternately determined by mostly one component since the other component cancels itself interferometrically at the location of that maximum, e.g. at z=-0.8λ and z=-0.3λ a maximum of the longitudinal and radial component occurs, respectively. If the photoresist layer is thin and placed at a maximum of the longitudinal component, the width of the spot can be reduced significantly and still be symmetric, in contrast to the case were the illumination of the lens is linearly polarised. The ring configuration of the lens is necessary to obtain a longitudinal field which is much larger than the transversal field in the focal region in air. In this way, the longitudinal component is still dominating when transmitted to the photo-resist medium, since the boundary conditions decrease the amplitude of the longitudinal field by a factor of ε 2=2.25 in respect to the transversal component. The developed optical effect in the photo-resist layer depends on the electric energy density in that layer. After averaging over the layer thickness in the longitudinal direction, the exposure profile shows circular symmetry and has a full width at half maximum equal to 0.535λ, see Fig. 6(b). It is possible to tune the field distribution in the focal region by selecting layers with a specific refractive index and thickness, so that a 4Pi-microscope  configuration is approximated.
In this paper we have presented an algorithm to calculate the field distribution in the focal region which is both fast in operation and easy to implement. With this algorithm we have studied the transmission of the near-field components through an air gap of the order of the wavelength. For an air gap as large as a significant decrease of the spot size with respect to the free-space diffraction limit is found for radially polarised light, indicating an appreciable transmission of the high frequency components to the substrate. Also, we have shown that it is possible to obtain a smaller lateral spot size in a photo-resist layer by using the interference effects that arise when a metal substrate is placed just beyond the sample; this produces a 3D size reduction effect that is comparable to the effects observed in a 4Pimicroscopy arrangement. Finally, we believe that with the method presented above, several problems related to far- and near-field optical microscopy, lithography and optical recording can be analysed in detail with our transparent and efficient calculation scheme.
Appendix A: Calculation of the transmission and reflection matrix
In order to calculate the propagation matrix of an imaging system when the focal field is embedded in a multilayer system, we derive a general formalism describing the rotation of the field vector by the lens and the influence of various reflections and transmissions suffered in multilayered media by a traveling wave with a certain oscillation direction (polarisation). The total field can always be decomposed in plane waves by a Fourier transform and each plane-wave component can be described as
Due to the rotational symmetry of the lens, it is logical to decompose the propagation vector into cylindrical coordinates =(kri ,kϕi ,±kzi ). The subscript i denotes the medium in which the vector is defined and superscript ± denotes the sign of the propagation direction. Since the transition of one medium i to another medium i+1 at position z=di should be independent of the coordinates (r,ϕ), it follows from "eA1">Eq. (A1) that kri =k ri+1=kr and k ϕi=k ϕi+1=kϕ . For monochromatic waves, as is the case considered here, we have a fixed length of the propagation vector = + . The sign of the square root which should be taken to acquire the propagation vector in the z-direction kzi follows from energy conservation; i.e. the sign should be chosen such that the wave is exponentially decreasing in the direction of propagation.
To be able to include the rotation effect of the lens on the field vector, we project the Cartesian basis on a cylindrical basis denoted by the operation P, perform a rotation of the field vectors in the cylindrical basis denoted by the operation R, and project the cylindrical basis back onto the Cartesian basis, which is the inverse of the first operation P -1, yielding a lens operator matrix L, as depicted in Fig. 7. The matrices describing the operations P, R and L are given by
where k 1=|k 1| is the amplitude of the propagation vector in the first medium. Next, we define a set of unit vectors (k̂,l̂,m̂) necessary to describe the effect of a medium transition. The unit vector k̂=k/k is parallel to the direction of propagation of the wave, the unit vector l̂=(k̂×ẑ)/|k̂×ẑ| is perpendicular to the propagation vector and the normal of the interface (TE-polarisation), and the unit vector m̂=k̂×l̂ completes the orthonormal basis (TM-polarisation).
By performing a scalar multiplication of the field vector with the three unit vectors as defined in Eqs. (A5–A7), we obtain three components associated with the oscillation direction. For the transition to the following medium we multiply these components with the corresponding unit vector in the second medium; performing both operations at once yields three matrices k̂2 k̂1,l̂2 l̂1 and m̂2 m̂1. Now, for the successive transitions of the media, we have to multiply the matrices in each consecutive medium yielding k̂ i k̂i-1·k̂i-1 k̂i-2·…·k̂2 k̂1=k̂ i k̂1, and similarly for the other two directions, resulting in
where we assume that we know that the field vector in the first medium travels in the positive direction. The unit vector l̂ is independent of both the medium and the direction of propagation, therefore the same holds for matrix l̂l̂. Now, we can obtain the propagation matrix associated with the three components as
where , and denote the effective transmission and reflection coefficients for the longitudinal, TE and TM polarisation components, respectively. Since optical waves traveling in the longitudinal direction are purely transversal, we can set the effective transmission and reflection coefficients for the longitudinal contribution to zero, =0, therefore the contribution of the matrix =0. Finally, we find for the total propagation matrix
with the effective reflection and transmission coefficients given by
where the number n=0,1,2 relates to the angular dependence of the sine- or cosine-function that the coefficient is multiplied with in the propagation matrices for the integration over kϕ . Note that for a single medium N=1, =M and =0, where the lens rotation matrix L would be equal to the propagation matrix M if the last column is taken equal to zero.
Appendix B: Calculation of the transmission and reflection coefficients
To derive the transmission and reflection coefficients for optical waves, we follow the derivation as given in . We assume a time-harmonic electric 𝓔(r, t)=Re[E(r)exp(-iωt)] and magnetic field 𝓑(r, t)=Re[B(r)exp(-iωt)]. The Maxwell’s equations yield the boundary conditions between medium i and medium i+1 given by
We consider only non-magnetic materials (µi =µ0 ), and use the relation B=ω -1 k×E from the Maxwell’s equations, which yields two independent equations for the transverse electric (TE) and for the transverse magnetic (TM) polarised components. For the TE-polarisation (s) the boundary conditions evaluate to
and for the TM polarisation (p)
Solving both sets of equations for and results in the ratio of backward traveling waves over forward traveling waves
where the Fresnel reflection coefficients for the TE and TM polarisations are given by
Since we know that in the last medium N there are no backward traveling waves coming from infinity, the ratio . Now, it is possible to calculate this ratio in each preceding layer until we obtain the ratio in the first layer. The transmission and reflection coefficients in each layer are calculated following the definition =/. In the first medium we find f s/p+1=1 and as the effective transmission and reflection coefficients, respectively. The effective reflection and transmission coefficients pertaining to a general transition i can be calculated by progression through the successive layers by using Eqs. (B2–B5)
with two polarisation-dependent prefactors defined as
The authors would like to acknowledge useful discussion with Peter Török from Imperial College London, UK. Part of this research has been carried out in the frame-work of the European project IST-2000-26479 (SLAM) and A.S. van de Nes acknowledges financial support from Philips Research Laboratories, Eindhoven, The Netherlands.
References and links
1. V.S. Ignatowsky, “Diffraction by a lens of arbitrary aperture,” Tr. Opt. Inst. Petrograd 1(4), 1–36 (1919).
2. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. Roy. Soc. A 253, 358–379 (1959). [CrossRef]
3. H. Ling and S. Lee, “Focusing of electromagnetic waves through a dielectric interface,” J. Opt. Soc. Am. A 1, 965–973 (1984). [CrossRef]
4. P. Török, P. Varga, Z. Laczik, and G.R. Booker, “Electromagnetic diffraction of light focused through a planar interface between materials of mismatched refractive indices: an integral representation,” J. Opt. Soc. Am. A 12, 325–332 (1995). [CrossRef]
5. A. Egner and S.W. Hell, “Equivalence of the Huygens-Fresnel and Debye approach for the calculation of high aperture point-spread functions in the presence of refractive index mismatch,” J. Microsc. 193, 244–249 (1999). [CrossRef]
6. D.P. Biss and T.G. Brown, “Cylindrical vector beam focusing through a dielectric interface,” Opt. Express 9, 490–497 (2001), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-9-10-490 [CrossRef] [PubMed]
7. A.S. van de Nes, P.R.T. Munro, S.F. Pereira, J.J.M. Braat, and P. Török, “Cylindrical vector beam focusing through a dielectric interface: comment,” Opt. Express 12, 967–969 (2004), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-5-967 [CrossRef] [PubMed]
8. D.G. Flagello and T. Milster, “3D Modeling of high numerical aperture imaging in thin films,” in Design, modeling, and control of laser beam optics , Y. Kohanzadeh, G.N. Lawrence, J.G. McCoy, and H. Weichel, eds., Proc. SPIE 1625, 246–276 (1992).
9. D.G. Flagello, T. Milster, and A.E. Rosenbluth, “Theory of high-NA imaging in homogeneous thin films,” J. Opt. Soc. Am. A 13, 53–64 (1996). [CrossRef]
11. R. Kant, “Vector diffraction problem of focussing a category 1 aberrated wavefront though a multilayered lossless medium,” J. Mod. Opt. 51, 343–366 (2004).
12. W. Welford, Aberrations of optical systems (Adam Hilger, Bristol, 1986).
13. G.N. Watson, A treatise on the theory of Bessel functions (Cambridge University Press, Cambridge, 1966).
14. W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical recipes in C: The art of scientific computing, 2nd ed. (Cambridge University Press, Cambridge, 1992).
15. P. Török, C.J.R. Sheppard, and P. Varga, “Study of evanescent waves for transmission near-field microscopy,” J. Mod. Opt. 43, 1167–1183 (1996). [CrossRef]
16. A. Egner, M. Schrader, and S.W. Hell, “Refractive index mismatch induced intensity and phase variations in fluorescence confocal, multiphoton and 4Pi-microscopy,” Opt. Commun. 153, 211–217 (1998). [CrossRef]
17. C.J.R. Sheppard and M. Gu, “Axial imaging through an aberrating layer of water in confocal microscopy,” Opt. Commun. 88, 180–190 (1992). [CrossRef]
18. L.E. Helseth, “Roles of polarization, phase and amplitude in solid immersion lens systems,” Opt. Commun. 191, 161–172 (2001). [CrossRef]
20. S.W. Hell, “Increasing the resolution of far-field fluorescence light microscopy by point-spread-function engineering,” in Topics in fluorescence spectroscopy, Vol. 5, J.R. Lakowicz ed. (Kluwer Academic/Plenum, New York, 1997), pp. 361–426. [CrossRef]
21. K.S. Youngworth and T.G. Brown, “Focusing of high numerical aperture cylindrical-vector beams,” Opt. Express 7, 77–87 (2000), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-7-2-77 [CrossRef] [PubMed]
22. S. Quabis, R. Dorn, M. Eberler, O. Glöckl, and G. Leuchs, “Focusing light to a tighter spot,” Opt. Commun. 179, 1–7 (2000). [CrossRef]
24. M. Paulus, P. Gay-Balmaz, and O.J.F. Martin, “Accurate and efficient computation of the Green’s tensor for stratified media,” Phys. Rev. E 62, 5797–5807 (2000). [CrossRef]