Image formation in a multilayer using the extended Nijboer-Zernike theory

J. J. M. Braat Optics Research Group, Faculty of Applied Sciences, Delft University of Technology, Lorentzweg 1, NL-2628 CJ Delft, The Netherlands S. van Haver s.vanhaver@tudelft.nl Optics Research Group, Faculty of Applied Sciences, Delft University of Technology, Lorentzweg 1, NL-2628 CJ Delft, The Netherlands A. J. E. M. Janssen Philips Research Europe, HTC-36, NL 5656 AE Eindhoven, The Netherlands


INTRODUCTION
Numerical methods based on rigorous calculation of the electric field around the focal region of a focusing lens system are a powerful tool to understand and interpret image formation and information retrieval, since the 3D point-spread function is directly related to the resolution and frequency response of the optical system.For lens systems where the Debye approximation is valid, a very common method to calculate the field in focus is to consider the diffraction integral according to the formulation first described by Ignatowsky [1] and re-derived by Wolf [2] and Richards and Wolf [3].Using this general description, one can obtain the field distribution in a homogeneous medium, but in most cases the propagating field, between the exit pupil of the lens system and the focal region, is subjected to one or even several transitions between layers of a stratified focal region.Thus, in order to obtain an accurate description of the focused field at the region of interest, one must take into account the effects of these layers, since the transmission and reflection at the transitions depend on the angle of incidence and state of polarization of an incident plane wave.The total effect on the imaging by the system is obtained by integrating over the plane wave components that are present within the angular aperture of the lens system.The procedure has been carried out in previous work on the transition between two media like in microscopy (air/sample or immersion liquid/sample) [4]- [8], and regarding general stratified media, with applications in optical lithography [9]- [11], optical recording [12]- [15], and confocal microscopy [16].As an alternative to numerical methods to solve directly the diffraction integral, it has been shown that the Extended Nijboer-Zernike (ENZ) theory can be used to compute the through-focus behaviour of the optical image [17,18].The ENZ method is based on first constructing the Zernike expansion of the wave field in the entrance pupil.This expansion is then multiplied by the complex lens transmission function to include the amplitude and phase changes of the optical field accumulated during the transition through the imaging system.Using the resulting Zernike field expansion in the exit pupil, ENZ theory provides the through-focus electric and magnetic field in image space.The advantages of using the ENZ method as compared to direct numerical computations of the diffraction integral is that it is semi-analytical, arbitrarily accurate and that it can deal with aberrations in a straightforward way.The separation in the ENZ-analysis of the aberrational effects and the diffraction effects in the focal region allows for an important reduction of the computational effort.The computation beforehand and the subsequent storage in look-up tables of basic spread-functions greatly enhances the numerical speed once a real aberrated system needs to be analyzed with respect to its imaging properties in the focal volume.In previous papers, we have shown how the ENZ diffraction theory can also be used to obtain an efficient and accurate imaging algorithm for extended objects with general shape and for arbitrary illumination conditions in the case of a homogeneous medium between the exit pupil and the focal region [19,20].In the present paper, we extend the ENZ

EXTENDED NIJBOER-ZERNIKE IMAGING IN A MULTILAYER
In an earlier publication [20] we have shown that the Extended Nijboer-Zernike (ENZ) theory can be used to compute the through-focus aerial image of an extended object.This work already included the possibility to have a different medium in both object and image space, but did not account for a non-uniform or layered configuration in image space.Such configurations are of particular interest as they are often encountered in advanced imaging applications such as advanced lithography.In lithography, the image is created in a resist layer that is enclosed by several other layers such as the wafer stack and protective capping layers.
In this work we will analyze the implications that arise when a multilayer configuration in the focal region of an imaging system is considered within the frame-work of the ENZ formalism.It will be shown that, although the layered configuration gives rise to many light reflections and thus strongly influences the image formation, it still remains possible to apply the main results provided by ENZ theory.
In Figure 1, a schematic presentation is given for the general optical system considered in this paper.It is assumed that the object is very small compared to the transverse dimension of the entrance pupil and that it is illuminated using a K öhler illumination scheme.The optical system, represented by its complex transmission function, T I , and magnification, M, transforms the field captured by the entrance pupil into a field distribution in the exit pupil.So far, the approach remains the same as in [20], but now, instead of constructing an image in an uniform image space we are faced with image formation in a focal region that consists of several layers with different refractive index n h .The problem with such a configuration is that, at the interface between two layers having a different refractive index, the light will be partly refracted and reflected, and thus gives rise in every layer to both forward and backward travelling electric field components, t h and r h .In this paper, we will track all transmitted and reflected contributions in a certain layer in a systematic way, and will show that this will allow us to do simulations of image formation in a layered stack using the ENZ formalism.

ENZ THEORY IN CASE OF A MULTILAYER IN THE FOCAL REGION
To keep track of the forward and backward propagating waves in the multilayer part of the image space, it is convenient to decompose the electric (and magnetic) field components in orthogonal sand p-polarization components.In the refraction or reflection process at the optical surfaces in the optical system and at the interfaces in the multilayer, the sand p-components suffer different changes in amplitude and phase that follow from the complex reflection and transmission coefficients at the interfaces.In this section, we first develop the expressions for the sand p-polarized plane wave components in the entrance pupil of the imaging system.The next step is to follow these components through the optical system and to find their direction and complex amplitudes on the exit pupil sphere in the homogeneous image space with index n 1 .In the last step, we calculate the modified direction and amplitude of the transmitted and reflected plane wave components as they are found in the layer with index h (refractive index n h ).These plane wave components are then used in the integrand of the Debye integral to obtain the field in the specific layer with index n h .The integration is carried out over the exit pupil where the truncated transmitted and reflected plane wave spectra originate.

Field components in the entrance pupil
We start with the Cartesian field components of the incident field distribution in the object plane directly after the object structure.These components are obtained by some other means, for instance, a rigorous diffraction calculation of the coherent field transmitted through the object structure.This object field is then propagated to the entrance pupil sphere of the optical system using a plane wave expansion of the Cartesian object field components.As usually, it is sufficient to consider two components out of three, for instance E x and E y , the third component E z following from the orthogonality of the electric components with respect to the propagation direction of the particular plane wave.Referring to Figure 2 where the object structure has been simplified to a fictitious point source, we project the Cartesian field components in the entrance pupil onto a spherical coordinate base that is given in a general point Q 0 (ρ, θ) by 2 Definition of the local basis for a general point Q 0 on the entrance pupil sphere with an axial cross-section (left-hand graph) and a cross-section perpendicular to the z-axis (right-hand graph).
The normalized radial coordinate ρ with 0 ≤ ρ ≤ 1 is defined by α 0 = arcsin(ρ/R).In general, we present unit vectors by a bold character with a hat on top and the unit vectors p0 , ŝ0 and k0 in the figure form a right-handed coordinate system along the sand p-polarization directions and the local propagation direction of the field.
From the components E x and E y in the entrance pupil, we can find E z by exploiting the orthogonality of the field components with respect to the propagation unit vector k0 yielding The field components on the new basis are then given by E 0,k (ρ, θ) was forced to be zero because of the Fraunhofer approximation.This approximation is valid in the case of diffraction towards an entrance pupil that is very far away from the object in terms of the wavelength of the light.

Field components in the exit pupil
The transmission of the plane wave spectrum through the optical system is treated according to the laws of geometrical optics.Each plane wave component with unit propagation vector k0 at Q 0 in the entrance pupil experiences amplitude and phase changes following from the trajectory of the corresponding geometrical ray through the optical system.The intersection point Q 1 of this optical ray with the exit pupil in the homogeneous part of the image space with index n 1 follows from some ray mapping condition between object and image space.In the case of a high-quality large-field imaging system, the mapping condition is given by the Abbe sine condition.In terms of the unit wave propagation vectors k0 and k1 in object and image space we have with M the lateral magnification of the imaging system.The position of the point Q 1 then follows from the plane wave propagation vector in the homogeneous part of the image space with index n 1 and is given by with R 1 the axial distance of the exit pupil towards the image plane.At this point we also introduce some abbreviations for goniometric quantities that will frequently occur in the expressions for the field components on the entrance and exit pupil sphere and in the multilayer region.We write From the sand p-polarized components of the electric field on the entrance pupil sphere, we derive the sand p-polarized components on the exit pupil sphere for a forward propagating plane wave with unit wave vector k1 according to 3 Definition of the sand p-polarization directions in an azimuthal plane defined by the angles θ and θ + π.The transmitted propagation wave vector in the thin layer with refractive index n h (azimuthal coordinate θ) is derived from an effective wave vector in image space with unit vector kh = (k x,h , k y,h , k z,h ).The sand p-components of the oppositely directed reflected wave with unit wave vector kh,r = − kh are obtained from the diametrically opposed direction of incidence with unit wave vector (−k h,x , −k h,y , k h,z ) and azimuth θ + π.The p-directions of the transmitted and reflected waves are parallel, their s-components have opposite directions.
Here we have used the same steps as in [20], Eqs. ( 7)-( 15) with R p given by the axial distance from the exit pupil to the paraxial focal point of the imaging system.See also Figure 3 for the direction of the sand p-components in image space.In these expressions we have The ρ-dependent amplitude factor T R follows from the chosen lens mapping factor, in our case the Abbe sine condition.
As customarily, T I (ρ, θ) is the complex lens transmission and aberration function, f 1 is the focal distance, n 0 is the refractive index of the object space and n 1 the refractive index of the homogeneous part of the image space before entering the multilayer structure.

Field components in layer h of the multilayer system
When we introduce a multilayer in the focal region the expressions for the sand p-components of the electric field change due to various reflections and refractions at the multilayer interfaces.Snell's law applied to the various interfaces of the mulitlayer yields for the unit propagation vector in the h-th layer The expressions for the sand p-components of the electric field in the h-th layer in the focal region are the sum of those pertaining to the forward propagating wave with unit propagation vector kh,t = (k h,x , k h,y , k h,z ) and complex transmission factor t h and those associated with the reflected wave in the h-th layer with unit propagation vector kh,r = (−k h,x , −k h,y , −k h,z ) and complex reflection coefficient r h .The reflected wave originates from a forward propagating wave with unit propagation vector kh = (−k h,x , −k h,y , +k h,z ); in the far field, its polar coordinates are given by (ρ, θ + π).We have to define the sign convention for the reflection coefficient r in the case of sand p-polarization.Like in most textbooks, we choose the s-direction for the reflected wave in the same direction as for the incident wave; the p-direction of the reflected wave is chosen opposite to the p-direction of the incident wave.With this convention, we find the following expressions for the sand p-field components in the h-th layer: The minus sign in front of r h,s in Eq. ( 14) follows from the opposite sign convention for s-polarization of the transmitted wave and the reflected wave, see Figure 3.The Cartesian components in the h-th layer are given by The substitution of E h,s and E h,p in Eqs. ( 16)-( 18) yields The Zernike coefficients that pertain to a layer h in the stratified image region follow from the field in the exit pupil that includes the lens transmission and aberration function and from the functions t h (ρ) and r h (ρ).Although one could also include the lens mapping function T R (ρ) into the construction of the Zernike polynomials, we choose not to do so because the mapping is an entity that is basically detached from the lens quality and the object properties.Of course, in the case of lithographic imaging, there is virtually no choice left and the lens mapping is entirely defined by the Abbe sine condition.Regarding the functions t h,p (ρ), t h,s (ρ), r h,p (ρ) and r h,s (ρ), these can be obtained in a straightforward manner using, for instance, thin-film matrix theory.
We construct the required β-coefficients using T I t h,s , T I t h,p , T I r h,s and T I r h,p .At this point it is also straightforward to include birefringence of the optical system into the formalism by allowing non-identical aberration functions, T I,x and T I,y , acting on the xand ycomponents of the incident field, respectively.We then find the following sets of β-coefficients and similarly for the counter propagating field components For analysis of the field components in layer h, we only need these eight sets of coefficients.When going to another layer, new sets have to be constructed.Moreover, it will turn out that each sublayer also requires new diffraction integrals to be integrated and tabulated (V(r, f )-functions).The ENZ approach has to be compared with methods based on the direct summation of plane wave components in the sublayer h.The advantage of the ENZ-approach in the multilayer configuration is the availability of compact and accurate basic functions with the general shape V m n,j (r, f ) exp(imθ) that immediately yield the focused field in the entire layer by multiplication with the β-coefficients of Eqs. ( 22)-(29) for the forward and backward propagating fields.
With the Zernike expansions Eqs. ( 22)-( 25) and Eqs. ( 26)-( 29), expressions can be found for the forward and backward propagating Cartesian components of the electric and magnetic field.However, to continue the use of the V m n,j -type integrals that have been introduced in the earlier ENZ analysis with a homogeneous image space, we define the following sets of composite Zernike coefficients The Cartesian electric field components for each direction of incidence in layer h can now be evaluated, both for the forward and for the backward propagating field.Their directional or spectral dependence is determined by the normalized polar coordinates (ρ, θ) on the exit pupil sphere from where the propagating plane waves components in image space originate.The somewhat lengthy expressions for E t h,x (ρ, θ), E t h,y (ρ, θ) and E t h,z (ρ, θ), as well as the corresponding expressions for the back propagating components, are written in full in Appendix A.
Next, these Cartesian spectral components have to be inserted into the Debye diffraction integral.The propagation factor for the transmitted wave is exp{i[k h,x x + k h,y y + k h,z z]}; for the back propagating wave we have the exponential exp{−i[k h,x x + k h,y y + k h,z z]}.In both cases, the k z,hcomponent is obtained from the dispersion relation in the layer with index n h .Inserting the field components of forward and back propagating waves into the Debye integral, see [2,20], we find for the field vector E 2 in layer h the expression × exp {i2πrρ cos(θ−φ)} ρ dρdθ with and λ 0 the vacuum wavelength of the radiation.E t,r h , the socalled ray strength, represents the electric field components of a forward or backward propagating plane wave with wave vector ±(k x , k y , k z ), measured at a distance of 1 meter from the reference point in focus.The function E t h for the transmitted field in the integrand of the first integral has to be evaluated at the azimuth θ + π.This is because of the fact that a wave vector with azimuth θ originates in the exit pupil at an azimuthal position θ + π.In the second integral above for the reflected field, the azimuth θ of the integration variables and the azimuth of the incident wave vector that created the reflected wave are identical.
We may proceed now by applying ENZ theory to derive the fields E t,r 2 in the h-th layer of the focal region.Applying the formulae that were already derived in [20], the field in the focal region can now be written as (vector column notation between brackets) for the forward propagating contribution and in a similar fashion we get for the counter propagating contribution The integrals V n,j,t ± with j = −2, −1, 0, +1, +2 occurring in the expressions above, are given as and we also immediately find for absorption-free media that In Appendix B, a recipe is provided to efficiently compute the integrals given in Eq. (35).
To calculate the power flow and momentum flux in the focal region, we need the expression for the magnetic induction in a specific layer of the multilayer stack.In homogeneous space, for a plane wave, we have the relationship B = √ µ k × E, with k the unit propagation vector of the plane wave.In the optical domain, we are allowed to write µ = µ 0 with µ 0 the magnetic permeability of vacuum and 1/2 = n h ( 0 ) 1/2 .The unit propagation vector is measured in the layer of the stack with refractive index n h , for a certain value of the normalized pupil coordinates (ρ, θ) in the exit pupil, and is given by Using this quantity to obtain the magnetic induction on the exit pupil sphere and performing the Debye integral for the magnetic induction components we find the following quantities in the focal region (transmitted components in vector notation) with c the speed of light in vacuum, The backward propagating components of the magnetic induction are obtained in a similar way, with the unit vector s replaced by its negative counterpart for each reflected plane wave component, and we obtain

COMPUTATION OF THE POINT-SPREAD FUNCTION IN A MULTILAYER USING THE ENZ IMAGING METHOD
In Section 3 of this paper we have presented expressions, Eqs. ( 36)-(39), to calculate the forward and backward propagating electromagnetic field components in a given layer of a multilayer stack in the image region of an optical system.
In the present section, we will apply these new expressions to some optical systems in which layered configurations in image space play an important role and we will verify our results by comparing them with an existing method based on a numerical evaluation of the Richards and Wolf diffraction integral [14].

Multilayer effects in lithographic resist images
Our first example illustrates the importance of multilayer modelling in deep UV high numerical aperture optical lithography (λ = 193 nm, NA= 1.368 (water immersion)).In Figure 4, three different possible models for the wafer stack configuration are shown.In configuration a, it is assumed that the light has to cross a single material interface before image formation takes place in a half space of photoresist (n = 1.76).In configuration b, a more realistic wafer stack is considered in which the effect of a highly reflective silicon substrate (n = 0, 78 + 2, 46i) is taken into account.Configuration c shows a more advanced stack in which an anti-reflective coating (ARC) is included.This ARC acts as an optical gradient to minimize the light reflections coming from the substrate.For all three configurations we have computed the throughfocus point-spread function (PSF) that is formed in the resist layer.The geometric best-focus position was put one wavelength inside the layer of resist (z = 0 in Figure 5), meaning that focal shifts introduced by the layer transitions were not corrected for.The results are shown in Figure 5, where the first column contains data obtained with the ENZ method presented in this paper and, for comparison, the second column contains data for the same configurations, this time obtained by means of numerical integration.As can be seen from these figures, there is a very good agreement between the ENZ method and its numerical integration counterpart (a more detailed comparison of the accuracy of the two approaches is presented in Section 4.3).Comparing the results at different rows in Figure 5, we see that the simplified representation of the wafer stack given by configuration a accurately predicts the focal shift, relative to the geometrical best focus position, introduced by the medium transition.However, by comparing rows one and two we see that neglecting the influence of the substrate is not allowed.As a matter of fact, the standing wave pattern observed for configuration b is unacceptable from a lithographer's point of view and therefore in practice they often apply an anti-reflective coating (ARC).An ARC has been included in configuration c and it can be seen that it effectively reduces the standing wave non-uniformity down to an acceptable value.

Air-gap effects when imaging with a solid immersion lens
Our second example involves a proposed next-generation optical disk read-out system.In this system, a high refractive index 'solid' immersion lens is placed very closely to the disk in order to improve the read-out resolution by optical tunnelling of high spatial frequencies corresponding to evanescent plane waves.The thin layer of air between the lens and the spinning disk should be as small as possible for resolution purposes but large enough to avoid the risk of mechanical contact during play-back of a disk.The width of the air gap strongly influences the performance of the system.In Figure 6 we have plotted the through-focus intensity profile of the point-spread function that, in this example, serves as the read-out spot.Various effects introduced by the air-gap can be observed.A focal shift is present which should be accounted for in the actual design of the system.For large air-gaps, the average intensity of the light spot goes down and the full width at half maximum of the scanning light spot gradually increases.This The colour shading in each graph has been adapted to the maximum energy density occurring in the graph.
means that there is an optimum air-gap width for which the read-out beam still has a sufficiently small size whilst the spot intensity is adequate for read-out with a good signal-to-noise ratio.Figure 6 applies to a solid immersion lens with a numerical aperture of 1.45.It is seen that the full width at half maximum of the lateral intensity profile starts to increase substantially once the air gap width exceeds 40 nm.A practical compromise between disk storage capacity, optical transmission through the gap and mechanical robustness of the system was found to be a gap width of typically 25 nm [21].

Accuracy of the ENZ approach as compared with numerical integration
In previous work [17,22]- [24], an extensive study has been carried out on the correctness and the convergence of the ENZ semi-analytic expressions.In several cases, with exact analytic results available for comparison, it has been shown that the ENZ method can achieve an accuracy of 10 −6 in amplitude with a relatively modest number of terms included in the basically infinite summation series.In [22], it is shown that for a total defocus range of eight Rayleigh focal depths, the number of terms in each summation can be typically limited to 20.Much larger focal depths can be handled by applying a double Bessel expansion with respect to the lateral and the axial excursion from the centre of the point-spread function.The results for strongly varying exit pupil functions, both in amplitude and/or phase, are determined by the accuracy of the complex Zernike expansion of the pupil function.
In [24] and [20], it is shown how rather strongly oscillating pupil functions can be adequately matched by Zernike expansions with maximum radial and azimuthal indices of typically 20.The amplitude matching is generally correct down to a relative value 10 −6 ; only in the case of discontinuities in the pupil function itself or its derivatives, the residual error may amount to a value of, for instance, 10 −3 .The resulting ENZ calculations, chosen to be accurate themselves down to 10 −6 , will now be limited by the accuracy of the Zernike fit of the exit pupil field components.It was observed that the final accuracy of the field components in the image region was better than it would be expected on the basis of the accuracy of the Zernike fit of the field components in the exit pupil.This can be ascribed to an averaging effect when performing the Debye integral to obtain each field component in image space.
We have examined the differences between the results of the ENZ calculation scheme and of the numerical integration method.The geometry of the test problem was the one labelled c in Figure 4.The ENZ approach first requires a calculation-intensive tabulation of the V-integrals of Eq. ( 35).Once this preparatory work has been done, the Zernike coefficients of the field components in the exit pupil are calculated and the field distribution in the focal region is then obtained in a very fast way by using the tabulated values of the Vintegrals.In practice, this latter step is orders of magnitude faster than a numerical evaluation of the Debye diffraction integral.The summations in the ENZ calculations were truncated so as to achieve an accuracy of 10 −6 in the amplitude of the field components.The error tolerance of the numerical integration procedure was set to 10 −5 .This error tolerance for numerical integration was chosen because of memory and calculation time limitations.The differences between the results of the two methods for the dominating xand z-field components in the focal region were of the order of 1.10 −5 with respect to a normalized maximum amplitude of unity.This implies that the accuracy was not limited by the ENZ approach but determined by the error setting of the numerical integration procedure.

DISCUSSION AND CONCLUSIONS
We have shown that the Extended Nijboer-Zernike diffraction theory can be adapted to allow for the image formation in a stratified image space, as is encountered often in imaging applications like microscopy, optical lithography and optical data storage.With the aid of the angular plane wave spectrum of the transmitted and reflected wave field in a particular layer of the stratified medium, we are able to construct the effective transmitted and reflected fields in the exit pupil of the imaging system.We obtain the transmitted and reflected plane wave spectra in a specific layer by standard means, for instance, by using thin layer matrix theory.The transmitted and reflected fields are each used as input for the Debye diffraction integral that is then solved using the semi-analytic ENZ theory.With respect to standard numerical integration of the Debye integral, we have shown that the ENZ approach is highly accurate.Without loss of speed, the solution of the diffraction integral can easily be made accurate up to a level of 10 −6 in amplitude.Numerical methods often have a typical accuracy of 10 −3 up to 10 −4 in amplitude if reasonable calculation times are desired.The accuracy of the ENZ method is limited in practice by the accuracy of the Zernike fit of the transmitted and reflected fields in a sublayer, as projected backwards towards the exit pupil.This accuracy degrades if the effective pupil function shows sharp oscillations or discontinuities.Practical Zernike fit errors in our examples remained below a level of 10 −5 .It was shown that a numerical solution of the diffraction integral can be made accurate to 10 −5 .However, the calculation time becomes excessively large with respect to the ENZ approach.The use of numerical integration could be acceptable for a single field point evaluation.When repeated evaluations are needed like in object mask optimization for optical lithography, the ENZ approach will offer a superior performance, both in speed and accuracy.
Another advantage of the tabulation possibility of the ENZ method in a multilayer geometry is the following.For a given aperture angle in a selected sublayer in image space, the set of Zernike coefficients and V-integrals needs to be evaluated only once.The actual through-focus calculations in the layer are then obtained very quickly by using the prepared tabulated values of coefficients and integrals.This advantage is maintained for a specific sublayer when the aperture angle is left unchanged.Changes in other sublayers do not affect the tabulated data for the V-integrals in the layer under consideration.The changes only affect the plane wave spectra belonging to the forward and backward propagating waves in the sublayer and their corresponding Zernike coefficients.It is only when switching to another layer with a different aperture angle that a new tabulation is needed.The tabulation advantage is especially appreciated when many image calculations in the same medium are needed.Again, the optimization of object mask patterns to obtain a desired image pattern in the photoresist recording layer is a good example of such a repeated calculation with unchanged values of the tabulated V-integrals.
As we already pointed out in the introduction, an extra advantage a Zernike-based method is that accurate Zernike coefficients become available for the forward and backward propagating fields in a particular layer of the stratified medium in image space.The complex Zernike coefficients yield information about the imaging quality in the particular sublayer.Correction of the effective aberration state in a specific sublayer can then obtained by an appropriate design of the optical imaging system.For high-numerical-aperture imaging, the complex Zernike expansion define the polarizationdependent amplitude and phase deviations in a sublayer.These can then be corrected with the aid of birefringent means to obtain perfect imaging at a given position in the sublayer.
The stratified medium in image space may contain absorbing layers or layers where frustrated total reflection takes place via evanescent fields or plasmon creation.However, our implementation so far does not allow to calculate the field in a sublayer that itself shows appreciable absorption or where the plane wave spectrum comprises evanescent components.We actually work on an extension of the multilayer ENZ method so as to cover these cases as well.

A EXPRESSIONS FOR THE FORWARD AND BACKWARD PROPAGATING FIELD COMPONENTS IN LAYER h
In this Appendix we produce the expressions for the electric field components in layer h of the stratified medium in image space.These components follow from the substitution of the composite Zernike coefficients of Eq. (30) into the Eqs.( 19)- (21).We then split the Cartesian components into their forward and backward propagating parts, denoted by the upper indices t and r, respectively.Shorthand notation for goniometric quantities according to Eq. ( 9) is used to improve the readability of the expressions.
In this Appendix we present a method for obtaining a series expansion for the integral V m n,j,t ± (r, f ) given by We use the identity where σ = ±1 indicates both sign possibilities, to write Eq. ( 46) in the form Next, we follow a similar approach as in [20], Appendix A, to write the integral in Eq. (48) as a series expansion given by Furthermore, the coefficients g , f and B t , as well as C s in Eq. ( 49) can be obtained starting from Eq. ( 51) and following the recipe given in [20], Appendix A, Eqs.(A3-A22).The T k l given in Eq. ( 50) have already been computed in [18] for l − |k| ≥ 0, however in Eq. (49) values l − |k| < 0 can also occur.We therefore derive a computation scheme for T k l , that is valid for general integer values of l and k, below.We assume k ≥ 0 (note that J −k = (−1) k J k ).

1 FIG. 5
FIG. 5 Comparison between the ENZ method (left column) and a fully numerical integration method (right column).The figures show, from top to bottom, the intensity of the x-polarized point-spread distribution in the resist layer for the configurations a, b and c, presented in Figure 4.The geometric best focus in water would have been found at the axial position z = 0.

FIG. 6
FIG. 6 Point-spread formation in a layer of PC (n = 1.62) using a solid immersion lens (n =2.086) for various values of the air-gap thickness (N A = 1.45, λ = 405 nm).
The different layers have their respective refractive indexes, n h , (possibly complex) where the subscript refers to the h-th layer starting from the exit pupil.t h and r h are labels pertaining to the amplitudes of the forward and backward propagating plane wave components in the layer labeled h.The point F is the geometrical best focus position in case of a uniform image region and the numerical aperture NA is given by n 1 sin α with sin α = s 0 .
F FIG. 1 Schematic representation of imaging into a multi-layered image space by an optical system with magnification M.