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 $Insertformulahere$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.

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 for this project can be found on this iPython notebook: Code_okkeun.zip