High-NA aberration retrieval with the extended Nijboer-Zernike vector diffraction theory

The reconstruction of the exit pupil function of an optical system can basically be carried out by collecting intensity data in the focal region from a certain number of defocused image planes. In this paper we present the ﬁrst results of such a reconstruction operation for optical systems with a high numerical aperture using a point source in the object plane. The main feature of our approach is the use of the extended Nijboer-Zernike diffraction analysis that has been modiﬁed to incorporate vector diffraction effects. The quality of the optical system is expressed by means of a set of complex Zernike coefﬁcients that describe the phase and transmission variation in the exit pupil of the imaging system. The ’vector’ method will be compared to the more common scalar diffraction analysis. We also analyse the practical limits of the vector retrieval process regarding the maximum allowed aberration and the noise of the intensity data. The sensitivity of the method with respect to parameter settings (state of polarisation and value of numerical aperture) is also examined. [DOI: 10.2971/jeos.2006.06004]


I n t r o d u c t i o n
The quality assessment of an optical imaging system can be carried out in several ways. Impulse response measurement is a straightforward method where the intensity pattern in focus is detected in the presence of a point source in the object plane. Frequency-based methods rely on grating structures in the object plane and the measurement of the image contrast as a function of spatial frequency. In both cases, intensity measurements provide a quality measure of the system (impulse response, frequency transfer) but no direct access to the aberration function of the imaging system. For a characterisation or a 'repair' action of the imaging system, the aberration function is of chief importance because it produces a direct understanding of the nature of the defects in the imaging system. The big advantage of interferometric methods is that they allow a direct evaluation of the wavefront shape in the exit pupil [1]. Unfortunately, interferometry is also a rather elaborate method for lens quality measurement. Its implementation in practical situations can be cumbersome, among others because of the lack of adequate coherent sources at the wavelength of interest, extreme sensitivity to vibrations, etc. For that reason, methods have been developed to improve the reconstruction of the exit pupil aberration function from intensity measurement in the focal plane (inversion methods; [2,3,4] for an overview of these, see Ref. [5]). A first step has been to simultaneously study the intensity distribution in the focal plane and in the exit pupil. Another possibility is to measure the through-focus behaviour of the intensity of the image of a point source. The complex pupil function is then obtained by a numerical matching procedure between the measured intensity and the forward calculated intensity pattern obtained in a continuously improving cycle of updates. In practice, the numerical burden of these methods can be high and the unique final match is not always obtained. A twofold improvement has been presented in some recent publications of the present authors by introducing the following refinements: • representation of the exit pupil function of the optical system by means of a Zernike polynomial expansion with corresponding complex coefficients. Such a Zernike expansion is capable of representing rather complicated pupil functions with a moderate number of coefficients.
• representation of the complex amplitude in the focal region based on an analytic (truncated) Bessel series expansion related to the Zernike expansion of the exit pupil function (extended Nijboer-Zernike theory); the complex amplitude can be obtained for substantially large values of the defocus parameter f in the exponential exp{+i f ρ 2 } (or a comparable expression at high NA) that occurs in the diffraction integrals [6,7,8]. In this paper we will use the basic scheme for aberration reconstruction that has been developed in Ref. [9] for imaging systems with a relatively low numerical aperture. This scheme has been further developed to correctly include defocus effects at high numerical aperture [10]. The feasibility to include vector diffraction effects in the extended Nijboer-Zernike analysis has been demonstrated both for the forward calculation [11] and for the retrieval process [12]. With these extensions, a complete framework is available to retrieve aberrations (and birefringence effects) in high-NA systems. In Section 2 of this paper, we briefly recall the theoretical framework that was presented in Refs. [11,12] and we apply this material in Section 3 to some examples of aberrated point-spread functions. In this section we also develop insight into the size of the aberrations that can be handled by our method and the influence of noise on the retrieval process. The section also contains some examples that illustrate the sensitivity of the method for incorrect parameter settings such as the state of polarisation and the numerical aperture value of the lens. Section 4 concludes the paper with a description of the applicability of the extended Nijboer-Zernike retrieval method.

E X T E N D E D N I J B O E R -Z E R N I K E T H E O R Y O F V E C T O R D I F F R A C -T I O N A N D R E T R I E V A L
In this section we briefly recall the expressions for the Cartesian electric field vectors in the focal region of a high-NA imaging system in the case of a point source object, Figure 1. Assuming quasi-monochromatic radiation, the field in the entrance pupil is described by a coherent superposition of two orthogonally polarised linear states of polarisation according to E = (a, b) where a and b generally are complex numbers. The influence of the (non-perfect) high-NA optical system is the introduction of wavefront deformation and transmission changes, on top of an intrinsic amplitude distribution on the exit pupil sphere that is different for each Cartesian field component. Moreover, we observe the so-called radiometric effect. The various amplitude distributions for each field component are described in the basic paper by Richards and Wolf for an ideal optical system; the field components in the focal region are found by calculating three basic integrals, named I 0 , I 1 and I 2 [13]. The behaviour of an aberrated optical system is described by an expansion of the complex lens transmission function in terms of Zernike polynomials with complex coefficients that are supposed to be identical for each polarisation state. We suppose that each vector component of the electric field in the exit pupil has to be multiplied by the complex pupil transmission function where the β's are the Zernike coefficients of the expansion and the index m with |m| ≤ n takes on positive and negative values to describe the cosine-and sine-dependencies in P(ρ, θ).
In all practical applications we suppose that β 0 0 is the leading term; this will be the case for optical systems that are not too far away from the diffraction limit. Using the extended Nijboer-Zernike theory, the complex field vectors in the focal region can now be calculated and with reference to the aberration-free case, we now have modified integrals V m n,−2 , V m n,−1 , V m n,0 , V m n,+1 , V m n,+2 with the indices (n, m) pertaining to the Zernike polynomial expansion on the exit pupil function. The extended Nijboer-Zernike theory allows a semi-analytical evaluation of each of these integrals with sub-indices ranging through {−2, ..., 2} for each index combination (n, m), see Ref. [11]. The linear superposition of the two orthogonal polarisation components in the xand y-direction with complex weights a and b in the entrance pupil leads to a mixing in the high-NA focal region. The expression for the electric field components in the focal region for x-polarisation reads and a comparable expression for the field components for ypolarisation in the entrance pupil is given by The functions V m n,j that depend on the normalised radial coordinate r and the defocus parameter f are given by (j = −2, −1, 0, 1, 2) with series expansion expressions available to obtain quick and accurate values of the integral above [12]. The total field in the focal region is now given by Journal of the European Optical Society -Rapid Publications 1, 06004 (2006) S. van Haver,et. al.
in which the weighting factors a and b appear explicitly. In this expression the parameter f is the defocus parameter and s 0 stands for the numerical aperture in image space. The expressions for the constants γ and u 0 can be found in [12].

. 1 T h e e x p r e s s i o n f o r t h e e n e r g y d e n s i t y a n d i t s F o u r i e r c o mp o n e n t s
The energy density in the focal region is proportional to |E| 2 .
After some rearrangement and tedious manipulation [12], the expression for the energy density of an aberrated high-NA imaging system can be written where we have intentionally left out the influence of spatially varying birefringence of the imaging system. In the following we will use the property that, for our highquality imaging systems, β 0 0 is the dominant coefficient and we can linearise the general expression for G k,l according to (7) where νµ is unity for ν = µ = 0 and zero for all other values and where we also introduced the shorthand notation

. A b e r r a t i o n r e t r i e v a l s c h e m e
The retrieval scheme for obtaining the complex pupil function is based on a Fourier analysis of the measured and the analytically proposed intensity data. The Fourier decomposition is carried out with respect to the harmonics in the azimuthal dependence of the through-focus intensity distribution. To this goal we evaluate for the analytically calculated distribution. A comparable operation is performed on the measured intensity data, yielding functions Ψ m (r, f ).
With the aid of the expression given in Appendix A of [12] we arrive at the following expression for the Ψ m an -function derived from the analytically calculated intensity distribution of Eq.(6) Note that the expression above is not exact but applies to the linearised approximation for the G k,l -functions. The equations to be solved now read for each m-value and they can be merged into one large system of linearised equations. The practical solution procedure consists of taking inner products on both sides with the functions Ψ m n;k,l (r, f ) and to solve this new system of equations instead of finding a solution that provides an optimum match for each (r, f ) combination. The inner product above is defined by where the integration limits are determined by the axial and lateral range of the collected intensity data.

R E T R I E V A L E X A M P L E S U S -I N G N I J B O E R -Z E R N I K E V E C -T O R D I F F R A C T I O N T H E O R Y
In this section we will first present an example that shows the inadequacy of a retrieval method based on scalar diffrac-tion theory when dealing with high-NA imaging systems. Secondly, we will show the ranges of aberration and transmission defects that can be handled by the linearised system of equations given above. A wider range of retrieved values can be realised once we use the so-called predictor-corrector extension when solving the linearised equations. Subsequently, the error is evaluated that is introduced in the retrieval process when system parameters, such as the numerical aperture and incident polarisation, are not exactly known. Finally, we present a simulated retrieval example performed in such a way that it closely resembles the treatment of experimental data. In order to show the necessity of treating the full vectorial case when performing retrieval upon systems of high NA, we present the following simulation. A through-focus intensity distribution (containing five axially displaced through-focus images, f = −2, −1, 0, 1, 2 units in dimensionless axial coordinates) is calculated for an aberration-free optical system having an NA of 0.95 that is illuminated by purely x-polarised light. We use the forward-calculation scheme embodied by Eqs. (2)(3)(4)(5). Thereupon, the high-NA data set acquired through this operation, is analysed using a scalar version of the ENZretrieval scheme. This gives rise to the set of retrieved βcoefficients that can be found in the right column of Table 1. The fit imposed by this limited set of β-coefficients is remarkably good as can be derived from Figure 2. Nevertheless, on comparing the β-coefficients obtained through the scalar retrieval operation with the actual β's (aberration-free case) used for the simulation, we see poor correspondence. We observe that, scalar ENZ-theory can fit high-NA intensity distributions, but the β-coefficients resulting from this approach no longer have a direct physical relevance. Scalar ENZ-retrieval needs to introduce a strongly amplitude-deformed pupil function in order to describe the high-NA and vectorial effect that are not included in the scalar model. Note that the strong amplitude deformation even includes a region in the pupil where the amplitude is negative (see Figure 3). This effectively means that this region should be subject to a phase shift of π.

Input Retrieved
The above example shows that the β-coefficients that are found when applying scalar theory to intensity distributions governed by the vectorial regime, do not have direct physical relevance and no longer give a correct description of the system under consideration. If we want to retrieve the physically relevant β-coefficients for a high-NA optical system it is therefore mandatory to treat the full vectorial case. The high-NA ENZ-retrieval formalism was introduced in Section 2 and we will examine some of its characteristics in the subsequent sections.

. 2 A c c u r a c y o f t h e h i g h -N A E N Zr e t r i e v a l s c h e m e
When applying the high-NA ENZ-retrieval scheme, the retrieved β-coefficients are generally not exact. Only when retrieval is performed upon a system that is aberration-free (which effectively means we simulate a system described by only one β-coefficient, β 0 0 = 1), the β-coefficients obtained will be exact. In all other cases, when one or more aberrations are present in the system (described by additional β-coefficients), an error in the retrieved β-coefficients is present. This error originates from the linearision applied when deriving Eq. (10).
In order to obtain an indication of the magnitude of the error present in the retrieved β-coefficients the following simulations are performed. Starting from a perfect (aberration-free) system, represented by a single β-coefficient, β 0 0 = 1, we introduce one additional non-zero β-coefficient (β 2 2 = 0). Next, this pair of β-coefficients is used to simulate a through-focus intensity distribution that serves as input for the ENZ-retrieval operation. The ENZ-retrieval process will now generate an estimate for this pair of β-coefficients describing the system. These estimates, which we shall denote as β', are not exact and include a certain error. The above process is repeated for ever increasing sizes of β 2 2 , which leads to ever increasing sizes of the error present in β 2 2 '. The results can be found in Figure 4, where we have plotted the maximum error in the retrieved estimates for the β-coefficients versus the magnitude of β 2 2 used in the simulation.    2 2 . This is as expected, as they are exactly those terms, dependent on (β) 2 , that we omit when applying the linearised expression for G k,l (Eq. (7)). Because of this relation, ENZretrieval is very accurate for well-corrected optical systems (β ≤ 10 −1 ). On the other hand, if we have a system influenced by larger aberrations (β of the order of 1), the errors present in the retrieved β-coefficients will be of the same order, which means that the quality of the ENZ-retrieval will be very poor. Fortunately, the ENZ-formalism provides an elegant correction scheme that improves retrieval results and makes the retrieved β-coefficients converge to their exact values. This correction scheme, known as the "predictor-corrector" method, will be treated in somewhat more depth below.

. P r e d i c t o r -c o r r e c t o r :
i m p r o v i n g r e t r i e v a l q u a l i t y As already stated above, basic ENZ-retrieval is not exact and this becomes problematic for systems containing medium-tolarge aberration. This was also encountered when scalar ENZretrieval was treated. In order to improve scalar retrieval results, a so-called predictor-corrector iteration scheme was proposed and thoroughly tested in Ref. [14]. We now propose an equivalent iterative procedure for the high-NA case, which is formulated in Appendix A. This predictor-corrector procedure is based on predicting the error introduced by the linearision and correcting the intensity distribution for this error. This operation will improve the retrieval quality and, when applied in an iterative manner, will make that retrieval results, obtained from simulated intensity data, converge to their exact values. to various aberration terms that were either initially present or that were erroneously detected at the start of the iterative retrieval process. The end value is determined by machine precision.
In Figure 5 results for high-NA ENZ-retrieval with the predictor-corrector procedure are shown for an optical system subject to astigmatism in the x-direction (β 0 0 = 1, β 2 2 = β −2 2 = 0.5i, NA= 0.95, x-polarisation). One can clearly observe the steady decrease of the error present in the retrieval result for an increased number of iterations. For synthetic data, as used in this example, the error eventually goes down to machine precision of the calculation software used (MatLab). Note that accuracies, customary for practical applications, are reached within less than 10 cycles.
For real experimental data, when numerous inaccuracies (such as noise) are inevitable, the attainable precision will be limited. Still in that case, the residual errors in the retrieved β-values, obtained through the predictor-corrector procedure, will be small and of the same (or lower) order of magnitude than the noise present in the data. We have observed such a performance for aberration values that can be as large as twice the diffraction-limit, e.g. up to 0.15 λ rms wavefront deviation. In Subsection 3.5 a retrieval example is treated for data obtained from a simulated general optical system in which noise is present. But first we will investigate the effect on the retrieval quality of system parameters that are not exactly known.

. 4 U n c e r t a i n t y i n s y s t e m p a r a m et e r s
One of the possible complications encountered when going from simulated data to experimental data obtained from a real optical system, is that certain system parameters are not exactly known. In this subsection we will investigate the effect on the retrieval quality when incorrect values for the NA and the orientation of the linear state of polarisation are assumed.
To this end, we simulate a through-focus intensity distribution according to the β-values found in Table 2 system defined in Table 2 is subjected to ENZ-retrieval while assuming a range of different values for the numerical aperture. During this process the error in several β-coefficients is monitored and the results can be found in Figure 6. A striking observation is that the residual error in all β-coefficients is minimal for the correct value of the numerical aperture. This not only makes it possible to tune the retrieval process to the exact value for the numerical aperture, but also suggests a procedure that can accurately determine the numerical aperture of an unknown system. Secondly, an equivalent simulation was performed for several supposed polarisation states while the numerical aperture is known. Starting from an exclusively x-polarised illumination state the parameters a and b are varied while making sure that (|a| 2 + |b| 2 ) = 1 (a and b real) because of normalisation purposes. This leads to a rotation of the polarisation state to a certain angle relative to the x-axis. This presumed polarisation state, different from the actual state of the system leads to an additional error present in the retrieved set of β-coefficients. The results of these simulations can be found in Figure 7, where the error present in the set of retrieved βcoefficients is plotted versus the angle between the supposed and actual polarisation direction. From Figure 7 one observes that it is also important to have accurate knowledge about the polarisation state of the system under consideration in order to obtain good retrieval quality.
The example above illustrates an interesting relation between the possible inaccuracy in the polarisation state and the error in the retrieved β-coefficient. For an optical system of which the polarisation state is approximately known an equivalent operation, as was used for generation of Figure 7, can be used to determine the polarisation state with great accuracy. For the special case that one has no knowledge whatsoever about the state of polarisation, the above procedure is no longer applicable. This is caused by the fact that the predictor-corrector procedure is not valid for very large deviations from the real state of polarisation. In this case, basic ENZ-retrieval can be used to obtain an approximate polarisation state after which the predictor-corrector method can be applied to determine the state of polarisation with great accuracy. , is plotted (dotted line).

. 5 H i g h -N A E N Z -r e t r i e v a l e x a m p l e
At the time this publication was written, no experimental data was available that could serve as input for the high-NA retrieval operation described in this paper. Therefore the ultimate test, retrieval of the aberrations for a real-world high-NA optical system, was not possible. As the best available alternative we present the following numerical experiment.
A through-focus intensity distribution, corresponding to the exit-pupil field defined in Table 3, is simulated (top row Figure  8). Next, in order to simulate an experiment, we add noise to this distribution. The result can be found in the second row of Figure 8. Now this distribution serves as the input of the ENZretrieval process while assuming that the NA and polarisation state of the system are known with great accuracy, meaning we will not have to bother about the problems discussed in the preceding subsection. Results for basic ENZ-retrieval (a single retrieval cycle) can be found in the third row and the result obtained through the predictor-corrector method, which uses the full power of the ENZ-formalism, is shown in the bottom row of Figure 8.

06004-6
Journal of the European Optical Society -Rapid Publications   Table 3. The upper row is the actual distribution, the second row is what results after adding noise with a SNR of 10, the third row is the distribution defined by the first β-estimates and the last row is the distribution resulting after the predictor-corrector procedure. Note that all images have been scaled according to their maximum value in order to show maximum detail. The SNR-value of 10 applies to the in-focus distribution; from the pictures in the second row it is clear that the SNR-value is much lower for the out-of-focus intensity distributions.

C O N C L U D I N G R E M A R K S
In this paper we have shown the feasibility of doing aberration retrieval according to the full vector diffraction version of the Extended Nijboer-Zernike theory. It is shown how to construct a linear system of equations that, upon solving, results in an unambiguous set of Zernike coefficients describing the optical system under consideration. At the time this paper was written, no experimental high-NA through-focus intensity data was available, therefore numerous test on simulated data were performed and presented in Section 3. It became obvious that it is mandatory to treat the complex vectorial case, rather than the far more simple scalar case, in order to acquire meaningful retrieval results. In addition, we discussed the fact that basic ENZ-retrieval is generally not exact and that this becomes problematic for medium-to-large aberrations. Recognising this inadequacy of basic ENZ-retrieval, we inserted the so-called predictor-corrector procedure which significantly improves the possible retrieval range and quality. Looking forward to experimental high-NA data sets to become available we investigated some possible complications that are inextricably connected to dealing with experimental data. There we found that it is very important to have accurate knowledge about the numerical aperture and state of polarisation of the system under consideration in order to assure good retrieval quality and stability of the predictor-corrector procedure. Finally, we concluded this paper with a simulation that illustrates the full power and versatility of the ENZ-retrieval scheme.

A T H E P R E D I C T O R -C O R R E C T O R P R O C E D U R E
The predictor-corrector procedure has been described and tested in simulation in Ref. [14], Section 4, for the case of relatively low-NA optical systems that allow a scalar treatment of the image formation. The extension to the high-NA vectorial case is rather straightforward, the basic principles being identical, and so we give here only a brief outline.
We assume that a measured through-focus intensity distribution I m is available that is to be represented in the form Here χ 0,0 0,0 and χ m,0 n,0 pertain to the dominant aberration-free auto-term and the dominant cross-terms that arise in accordance with Eqs. (6)(7)(8), and χ m,m n,n is an elaborate term, that involves products V m n;j V m * n ;j , pertaining to small cross terms. The -signs in Eq.(A.1) indicate that the terms with n = m = 0 and n = m = 0 should be deleted. In the basic retrieval scheme, we choose the β's in the small cross-term deleted version 2) and I m is maximal; this is done in accordance with Eqs. (10-12). The re-