Indrasen Bhattacharya

Introduction

Aberrations in the wavefront profile degrade the resolving power of imaging systems. It is highly desirable to measure and compensate any aberrations present in an imaging system. Such imaging systems may include microscopes, telescopes, lithography steppers, iPhone cameras, AR/VR goggles as well as the human eye. We shall consider the more traditional case of optical microscopes and lithography steppers in this work, but we do not rule out extensions to other imaging systems. The point spread function (PSF) at the image plane determines the resolving power of microscopes and the feature size in lithography systems. It is important to have a narrow PSF that is preferably diffraction limited by the numerical aperture (NA) of the imaging system. Wavefront aberrations reduce the sharpness of the PSF and lead to undesirable asymmetries and artifacts. The wavefront aberrations are typically expressed in an orthogonal radial basis on the unit circular pupil called the Zernike polynomials. The Zernike polynomials $R_{n}^{m}(\rho )$, together with azimuthal sinusoids, form a complete basis on the unit circle: any 2-dimensional phase function can be expressed as an appropriate linear combination of these polynomials. The weighting factors $\alpha _{n}^{m}$ of the Zernike polynomials are termed as Zernike coefficients. Since any wavefront aberration profile can be fitted to a series of Zernike coefficients: the pupil function and Zernike coefficients are equivalent representations. Conveniently, the Zernike polynomials can be interpreted as specific physical aberrations such as spherical aberration, coma, astigmatism, trefoil and others, which makes the Zernike coefficients an intuitive framework to work with.

More explicitly, if the phase due to aberrations is expressed in normalized Cartesian basis as: $e^{i\Phi (\kappa _{x},\kappa _{y})}=e^{i\Phi (\rho ,\theta )}$, the Zernike polynomial expansion of the aberration phase $\Phi$ is given by:

$\Phi (\rho ,\theta )=\Sigma _{n,m}\alpha _{n}^{m}R_{n}^{m}(\rho )cos(m\theta )$

where $0\leq \rho \leq 1,0\leq \theta \leq 2\pi$. For the purpose of this analysis, we consider only the cosine terms without loss of generality.

It is conceptually straightforward to calculate the point spread function from the wavefront profile: they are related by a Fourier transform where the aberrations occur in the phase of the pupil function. However, in this problem, we are interested in the reciprocal calculation of determining what is wrong with the system based on how the PSF appears. This is the inverse problem: determining the Zernike coefficients from the aberrated point spread function. It is important to explicitly state the inputs and outputs of the computational procedure we are hoping to undertake. If the wavefront aberration is given as a function of the NA-normalized pupil coordinates (in Cartesian basis) as $\Phi (\kappa _{x},\kappa _{y})$, the point spread function can be expressed as:

$U(x,y;f)={\frac {1}{\pi }}\int _{\kappa _{x}^{2}+\kappa _{y}^{2}\leq 1}\int e^{if(\kappa _{x}^{2}+\kappa _{y}^{2})}e^{i\Phi (\kappa _{x},\kappa _{y})}e^{[2\pi i\kappa _{x}x+2\pi i\kappa _{y}y]}d\kappa _{x}d\kappa _{y}$

Note that this formula neglects the Radiometric correction factor, which becomes relevant for high numerical aperture systems. It also approximates the defocus phase for the small NA case. However, these approximations typically suffice for NA smaller than 0.6, which captures many relevant systems.

We do not expect to be able to retrieve a complex phase from the intensity point spread function at a single imaging plane, unless we make additional prior assumptions regarding the phase. Instead, we use the approach proposed by Van Der Avoort in [2], which relies on using multiple point spread functions. In the approach proposed by these authors, we hope to acquire a focal progression of point spread functions: $U(x,y,-f_{max})...,U(x,y,-\delta f),U(x,y,0),U(x,y,\delta f),...U(x,y,f_{max})$ and perform computations that will let us closely recover the Zernike coefficients $\alpha _{n}^{m}$ for the phase function. The authors have shown that under a linearizing assumption, the aberration phase has a one to one correspondence with a focal progression of intensity point spread functions: it can be exactly inverted to obtain the Zernike coefficients $\alpha _{n}^{m}$. This approach is feasible for lithography and optical microscopy imaging systems, and indeed has been used in these applications. It will require some thought and experiment design to extend these ideas to the human eye - we leave this for future work.

As a concrete example, a focal progression is displayed in the figure below. The focus increases in steps of 1 normalized unit (see 'Notations and Conventions') starting at -6 on the top left corner, and increasing rightwards. This was generated using a numerical implementation of the fourier transform formula above, using the DFT method. The units for the X and Y axes are normalized according to the conventions specified in the section on 'Notations and Conventions'. The system was simulated to be aberrated by vertical astigmatism and horizontal coma. It can be qualitatively observed that the impacts of the aberration become much more apparent at a higher defocus, particularly the coma.

Fig. 1: Focal progression of intensity for an optical system with two higher order aberrations: vertical astigmatism $(\alpha _{2}^{2}=0.5)$ and horizontal coma $(\alpha _{3}^{1}=0.3)$

Notation and Conventions

In this project, we have focused on one particular method: the one proposed in [2]. We have used python for a computational implementation of the technique, following the numerical conventions used in the paper. In particular, all spatial units $(X,Y,Z)$ have been normalized to the numerical aperture of the system as follows:

$x=X\times {\frac {2\pi NA}{\lambda }},y=Y\times {\frac {2\pi NA}{\lambda }},f=Z\times {\frac {\pi NA^{2}}{\lambda }}$

These are the spatial quantities used in the Fourier transform expression above.

We are also using the convention for Zernike polynomials used in reference [1], appendix VII and section 9.2. Please note that this convention is different from the normalized convention used on the Wikipedia page. As is the case for most conventions, we use m to denote the azimuthal degree of the Zernike polynomial and n the radial degree. n takes values m, m+2, m+4, ... Essentially, a certain azimuthal order can only occur for a sufficiently high radial order. Also, the difference between the radial and azimuthal order of the polynomials is always even. These are certain mathematical aspects of the Zernike polynomial expansion of the pupil phase that are useful to have at reference.

Certain choices regarding discretization and infinite series truncation have also been made in this implementation. Where these choices are different from the implementation in the original paper, I have clarified in the text in the relevant sections below.

Analytical Forward Model

An important aspect of the inverse problem solution is the development of a linearized model for the field point spread function. As a result, it will turn out that the field PSF is linear in the Zernike coefficients and the intensity PSF is quadratic. The error in the field PSF is of third order in the Zernike coefficients. This linearization will later turn out to be greatly beneficial for solving the inverse problem. We omit the derivation here but specify the key assumptions and results.

If the aberration phase is small, we can linearize it as follows: $e^{[i\Phi (\rho ,\theta )]}\sim 1+i\Phi (\rho ,\theta )=1+i\Sigma _{n,m}\alpha _{n}^{m}R_{n}^{m}(\rho )cos(m\theta )$

Substituting this simplification into the Fourier transform formula and using certain mathematical identities, we obtain a closed form solution for the field point spread function in normalized radial coordinates:

$U(r,\phi ;f)=2V_{0}^{0}(r,f)+2i\Sigma _{n,m}i^{m}\alpha _{n}^{m}V_{n}^{m}(r,f)cos(m\phi )$

We note that this formula is linear in the Zernike coefficients $\alpha _{n}^{m}$. The special functions $V_{n}^{m}$ used here are given by:

$V_{n}^{m}(r,f)=\int _{0}^{1}\rho e^{(if\rho ^{2})}R_{n}^{m}(\rho )J_{m}(2\pi r\rho )d\rho$, where $J_{m}(x)$ is the Bessel function corresponding to the azimuthal order.

The authors further express the special functions as an infinite summation using Bessel functions. We omit that expression here, but it can be found in [2], equation (10). This closed form analytical expression has been implemented in python. Please note that the infinite summation in equation (10), using the indexing variable l requires termination for a computational implementation. We will discuss this in the next section.

Note that the intensity PSF is simply given as: $I=|U|^{2}$ and can be written based on the above special functions as well. The intensity PSF will turn out to be quadratic in the Zernike coefficients, under this linearization assumption.

Verification of ENZ Forward Model

The authors of [2] have termed the analytical closed form solution for the field as the extended Nieboer-Zernike approach (ENZ approach). Since the inverse problem solution relies on these special functions and the accuracy of the forward model, we decided to independently test the forward model first. This will also provide some insight on how to terminate the infinite summation used in equation (10) of [2] for expressing the special functions $V_{n}^{m}(r,f)$.

In order to obtain a ground truth comparison, we have implemented the Fourier transform expression for the PSF using in-built DFT packages in python:

$U(x,y;f)={\frac {1}{\pi }}\int _{\kappa _{x}^{2}+\kappa _{y}^{2}\leq 1}\int e^{if(\kappa _{x}^{2}+\kappa _{y}^{2})}e^{i\Phi (\kappa _{x},\kappa _{y})}e^{[2\pi i\kappa _{x}x+2\pi i\kappa _{y}y]}d\kappa _{x}d\kappa _{y}$

The phase and resulting PSF are shown in Fig. 2 below. This was for a case with both astigmatism, coma and defocus. The amount of defocus is 3 normalized units; astigmatism and coma are same as in Fig. 1.

Fig. 2: Ground truth PSF calculation using the DFT approach. Top panel: pupil phase and 1D cuts, Bottom panel: resulting intensity PSF in real space and 1D cuts

In order to verify the analytical formula based on the special functions, we have implemented this and compared it with the ground truth in Fig. 2. For the same defocus and aberration parameters, the resulting point spread function image is shown below in Fig. 3. The difference with the ground truth, normalized to the maximum of the ground truth PSF is shown in Fig. 4. As observed, the intensity PSF resulting from the analytical approach is a good approximation for the DFT based calculation. The error is on the order of 5% close to the peak of the PSF. This is a good sanity check to show that the analytical approximation is accurate enough for us to continue with this approach.

Fig. 3: PSF resulting from the closed form analytical solution, using the same aberration and defocus parameters as Fig. 2.

Fig. 4: Difference in PSF, normalized by maximum of the PSF. As can be seen, errors on the order of 5% can be observed.

We note that this small error was able to be achieved by terminating the infinite summation in equation (10) of [2] with 10 terms. Therefore, we continue with this choice for the inverse problem solution. Please note that the aberrations specified in Zernike coefficients are in radians (vertical astigmatism $(\alpha _{2}^{2}=0.5)$ and horizontal coma $(\alpha _{3}^{1}=0.3)$). These can be converted to wavelength units, resulting in astigmatism of 79.6 milliwaves, and horizontal coma of 47.7 milliwaves. These are fairly realistic numbers and this case (among other sanity checks) will be used for the rest of this study.

Inverse Problem Solution

Here we discuss how to retrieve the pure phase aberrations from a defocus series of PSFs $I(r,\phi ;f)$. In a real world situation, these PSFs would be obtained empirically. However, for the sake of our computational experiments, I have used the Fourier transform method to generate a ground truth PSF progression. Next, we summarize the key steps in retrieving the Zernike coefficients. For a flavor of the derivation, please see section 3 of reference [2]. The key details involve:

1. Expressing the intensity PSF in the analytical form $I=|U|^{2}$

2. Inner product of the measured PSFs $I_{meas}(r,\phi ;f)$ with the azimuthal sinusoid in order to obtain azimuth-free special functions $\Psi _{meas}^{m}(r,f)$

3. Calculation of intensity basis functions $\Psi _{n}^{m}(r,f)$ from the original special functions $V_{n}^{m}(r,f)$

4. For each azimuthal order $m$, the formulation of a system of equations for obtaining the Zernike coefficients for all valid radial orders corresponding to that azimuthal order: $\alpha _{n}^{m},n=m,m+2,m+4,....m+2M$; where $M$ is used to express an appropriately chosen termination of the Zernike polynomial radial order

5. The inversion of the system of equations (by inverting the corresponding Gram matrix) in order to obtain the Zernike coefficients

The procedure appears to be rather involved, but it can be summarized as follows. The measured intensity PSF is written in terms of unknown Zernike coefficients. Then, the formula is projected into the basis of the azimuthal sinusoids and the intensity basis functions. Once this is done, the quadratic term is rigorously found to be zeroed out, leaving only terms linear in the Zernike coefficients. The linear expression can be solved for the Zernike coefficients just like any other system of linear equations.

Here we state the formula for the measurement projection (equation 23 in [2]):

$\Psi _{meas}^{m}(r,f)={\frac {1}{2\pi }}\int _{0}^{2\pi }I_{meas}(r,\phi ;f)cos(m\phi )d\phi$

The intensity basis functions can be pre-calculated and are given by (equation 24 in [2]):

$\Psi _{n}^{m}(r,f)=-8\epsilon _{m}^{-1}Im[i^{m}V_{n}^{m}(r,f)V_{0}^{0*}(r,f)]$

And each linear equation for solving for the Zernike coefficients is given simply by:

$\Sigma _{n=m,m+2,...m+2M}\alpha _{n}^{m}(\Psi _{n}^{m},\Psi _{n'}^{m})=(\Psi _{meas}^{m},\Psi _{n'}^{m})$

The full system of equations can be expressed by vectorizing over the radial order $n'$, leading to a matrix formulation. There is one matrix for each azimuthal order, and the size of the matrix is [M, M], given by how we originally chose to truncate the Zernike polynomial radial order. By inverting this matrix equation we obtain a vector of size [M,1] containing all the relevant Zernike coefficients for that azimuthal order. This needs to be done for each azimuthal order that we choose to model. I have explained the discretization and truncation choices in the next section.

Results

We have implemented the ENZ inverse calculation on a small scale problem and compared with ground truth. The following three cases were tested:

Case 1: Sanity check with no aberrations, all Zernike coefficients expected to be zero

Case 2: Only vertical astigmatism ($\alpha _{2}^{2}$) is present

Case 3: Vertical astigmatism and horizontal coma ($\alpha _{3}^{1}$) are both present

The results below (with retrieved aberration in shown as red) demonstrate that there is a reasonable qualitative correspondence in the retrieved aberrations and ground truth. However, the quantitative correspondence was found to be lacking. Some of the reasons for this are explained at the end of this section. Please note that the aberrations are specified by both the descriptive name and the (m, n) order of the Zernike coefficient in the comparison graphs below.

Fig. 5: Results for Case 1 (sanity check with no aberrations)

Fig. 6: Results for Case 2 (only one form of aberration)

Fig. 7: Results for Case 3 (two dominant aberrations: same as the case shown in the PSF progression)

Due to time and resource constraints, the simulation settings were rather coarse. For instance, the choice of the quantity M, the number of radial orders per azimuthal order could only be 2. Ideally, this number would be on the order of 10 in order to consider a sufficiently large number of radial orders. The Gram matrix method attempts a best fit of the data to the basis functions included in the projection, up to M. Hence, the results may be poor if M is rather low. Another constraint was the spatial discretization, which could not be too fine due to memory and run time constraints on my personal computer. It would be very interesting to perform a more detailed study exploring the impact of focus and spatial discretization, as well as M on the accuracy of the results.

Limitations of Method and Future Work

There are some well known limitations of the approach as presented in this write-up (which have been addressed in later publications [3]-[5]):

1. The approach is valid only for pure phase aberrations. There is an extension to phase and amplitude aberrations using complex Zernike coefficients that may be worth investigating.

2. The radiometric correction factor has been neglected in the current formulation.

3. A vectorial formulation is required for higher numerical aperture systems.

4. The effect of finite source extent has been neglected, although it is fairly straightforward to incorporate.

Ideally, a little more investigation into the accuracy, numerical stability and a few more test cases would have been great - but hopefully this is a good starting point.

References

[1] M. Born and E. Wolf, Principles of Optics, (Cambridge University Press, Cambridge, 2002)

[2] C. Van Der Avoort, J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, "Aberration retrieval from the intensity point-spread function in the focal region using the extended Nijboer-Zernike approach", J. Mod. Opt., 52, 12, 2005, pp. 1695-1728

[3] A. J. E. M. Janssen, "Extended Nijboer-Zernike approach for the computation of optical point-spread functions", J. Opt. Soc. Am. A, 19, 5, 2002

[4] J. J. M. Braat, P. Dirksen, A. J. E. M. Janssen, A. S. van des 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, 12, 2003

[5] A. J. E. M. Janssen, J. J. M. Braat, P. Dirksen, "On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus", J. Mod. Opt., 51, 5, 2004, pp. 687-703

[6] P. Dirksen, J. Braat, A. J. E. M. Janssen, "Estimating resist parameters in optical lithography using the extended Nijboer-Zernike theory", J. Microlith., Microfab., Microsyst. 5(1), 013005, 2006

[7] Extended Nijboer-Zernike (ENZ) Analysis and Aberration Retrieval research site (Netherlands group): https://www.nijboerzernike.nl/

Appendix

The computational implementation (jupyter notebook) and results for this project can be found in this zipped folder: File:Indrasen final.zip