On the computation of the Nijboer-Zernike
aberration integrals at arbitrary defocus
A.J.E.M. Janssen
1)
, J.J.M. Braat
2)
and P. Dirksen
1)
1)
Philips Research Laboratories, 5656 AA Eindhoven, The Netherlands,
e-mail: {a.j.e.m.janssen,peter.dirksen}@philips.com
2)
Optics Research Group, Department of Applied Sciences, Delft Univer-
sity of Technology, 2628 CJ Delft, The Netherlands,
e-mail: j.j.m.braat@tnw.tudelft.nl
Short title: Nijboer-Zernike aberration integrals
PACS codes: 42.15.Fr, 42.25.Fx, 42.30.Va
Abstract
We present a new computation scheme for the integral expressions describing
the contributions of single aberrations to the diffraction integral in the con-
text of an extended Nijboer-Zernike approach. Such a scheme, in the form
of a power series involving the defocus parameter with coefficients given ex-
plicitly in terms of Bessel functions and binomial coefficients, was presented
recently by the authors with satisfactory results for small-to-medium-large
defocus values. The new scheme amounts to systemizing the procedure pro-
posed by Nijboer in which the appropriate linearization of products of Zernike
polynomials is achieved by using certain results of the modern theory of or-
thogonal polynomials. It can be used to compute point-spread functions of
general optical systems in the presence of arbitrary lens transmission and
lens aberration functions and the scheme provides accurate data for any,
small or large, defocus value and at any spatial point in one and the same
format. The cases with high numerical aperture, requiring a vectorial ap-
proach, are equally well handled. The resulting infinite series expressions for
these point-spread functions, involving products of Bessel functions, can be
shown to be practically immune to loss of digits. In this respect, because
of its virtually unlimited defocus range, the scheme is particularly valuable
in replacing numerical Fourier transform methods when the defocused pupil
functions require intolerably high sampling densities.
1
1 Introduction
The computation of strongly defocused amplitude distributions has been con-
sidered by several authors. Their effort has been directed towards the stable
evaluation of the diffraction integral in the presence of a strongly oscillating
defocusing phase factor in the integrand of the basic diffraction integral. We
refer to the end of this section for a short survey of numerical and analytic
approaches to solve the defocus problem that can be found in the literature
[1]-[5]. In this paper we treat the strong defocus problem in the framework
of the recently developed extension of the Nijboer-Zernike approach to the
computation of optical point-spread functions of general aberrated optical
systems [6]-[10]. In the extension to the Nijboer-Zernike approach, power
series expressions involving the defocus parameter f , with coefficients explic-
itly given in terms of Bessel functions and binomial coefficients, were given
for the contribution to the diffraction integral of a single aberration term
β
nm
R
m
n
(ρ)cosmϑ with R
m
n
(ρ) a Zernike polynomial, see [11], Sec. 9.2. We
follow the developments as given in [7], Secs. 1–2. Thus, given the Zernike
expansion
n,m
β
nm
R
m
n
(ρ)cosmϑ of the pupil function A exp[iΦ], with the
transmission function A and the aberration phase Φ assumed to be symmet-
ric in the angular coordinate ϑ, the point-spread function U can be written
as
U(x, y)=2
n,m
β
nm
i
m
V
nm
(r, f )cosmϕ . (1)
Here we have used Cartesian coordinates x, y in the image plane that can be
transformed to polar coordinates r, ϕ according to x + iy = r exp(iϕ). The
V
nm
in (1) are the basic integrals
V
nm
(r, f )=
1
0
e
ifρ
2
R
m
n
(ρ) J
m
(2πρr) ρdρ , (2)
which should be considered for all integer n, m ≥ 0withn −m ≥ 0 and even.
For these V
nm
there holds the power series expansion, see [6]–[7],
V
nm
(r, f )=e
if
∞
l=1
(−2if)
l−1
p
j=0
v
lj
J
m+l+2j
(2πr)
l(2πr)
l
, (3)
with
v
lj
=(−1)
p
(m + l +2j)
m+j+l−1
l−1
j+l−1
l−1
l−1
p−j
q+l+j
l
(4)
2
for l =1, 2, ... , j =0, 1, ..., p in which
p =
1
2
(n − m) ,q=
1
2
(n + m) . (5)
This approach of computing point-spread functions has been assessed in [7]
from an optical and numerical point of view. As a rule of thumb one should
include in the infinite series over l some 3f terms to get sufficiently accu-
rate results. The approach has been extended further in [8] so as to cover
the cases of high numerical aperture which requires computation of field
components as well as inclusion of the radiometric effect and the state of po-
larization. Application of the method in a lithographic context is considered
in [9]–[10], where the inverse problem of estimating the coefficients β
nm
from
measurements of the (intensity) point-spread function in the focal region is
solved. In this kind of applications the focus variable f is often taken to
be complex-valued so that illuminated objects of small but finite size can be
accommodated
As said, the series expression in (3) yields accurate results when some 3f
terms are included in the series over l. The evaluation of the required Bessel
functions is normally no problem since one can exploit recursion formulas as
in [12], expression 9.1.27 on p. 361, when the efficient computation of the
Bessel functions is not already available in the software environment of the
user. A much more serious problem is posed by the fact that large values of
f lead to a considerable loss of digits in the series in (3). Typically, one has
terms of the order of magnitude |f|
l
/l! in (3) while the V
nm
’s themselves are
of the order of unity. Practically, the use of the series in (3) is limited to a
range like |f|≤5π, so that an axial range of the order of typically ten focal
depths can be handled.
A similar problem was noted in [13] where the attention is limited to
radially symmetric aberrations. Here for the computation of the integrals
T
2p,0
=
1
0
e
ifρ
2
ρ
2p
J
0
(2πρr) ρdρ , (6)
the exp(ifρ
2
) is expanded in powers of f so that the so-called generalized
Jinc functions,
Jinc
n
(u)=
1
u
2n+2
u
0
v
2n+1
J
0
(v) dv , n =0, 1, ... , (7)
appear for which explicit, finite Bessel series are given in [13]. Also in this
case one should limit f to a range like |f |≤5π.
3
In optical problems, axial excursions beyond the described limits are fre-
quently encountered, and we list here some examples:
– the self-imaging by periodic structures manifests itself far outside the focal
region,
– amplitude oscillations close to the geometrical boundary are encountered
in the Fresnel diffraction regime,
– certain optical problems use the fractional Fourier transform in which light
distributions occur that are intrinsically quite remote from the standard
Fraunhofer pattern,
– in optical recording, the extension to volumetric storage using several
recording layers in depth requires propagation of a focused beam through
several strongly defocused information layers.
In all these cases a reliable analytical method is required that allows the
calculation of strongly defocused, aberrated optical fields. In this paper it is
shown how the analysis developed in the framework of our extended Nijboer-
Zernike theory can be used to achieve this goal. Below we present the basic
features of our approach.
In [14] a method has been analyzed to compute Lommel’s functions of
two variables without loss of digits (the first two Lommel functions can be
expressed in terms of T
00
= V
00
, see (2) and (6)). This method was developed
by Nijboer and Zernike [15]–[16] and uses Bauer’s identity
e
ifρ
2
= e
1
2
if
∞
k=0
(2k +1)i
k
j
k
(
1
2
f) R
0
2k
(ρ) , (8)
see [12], formula 10.1.47 on p. 440 (observe that R
0
2k
(ρ)=P
k
(2ρ
2
− 1)) or
[11], formula (10) on p. 534, with
j
k
(z)=
π
2z
J
k+
1
2
(z) ,k=0, 1, ... , (9)
the spherical Bessel functions of the first kind, see [12], Ch. 10. Applying (8)
for the computation of T
00
in (6) we get
T
00
= e
1
2
if
∞
k=0
(2k +1)i
k
j
k
(
1
2
f)
1
0
R
0
2k
(ρ) J
0
(2πρr) ρdρ . (10)
4
The remaining integrals are computed using the basic result
1
0
R
m
n
(ρ) J
m
(2πρr) ρdρ=(−1)
1
2
(n−m)
J
n+1
(2πr)
2πr
(11)
from the “classical” Nijboer-Zernike theory, see [11], Sec. 9.2, formula (9) on
p. 525. Due to certain bounds on the Bessel functions j
k
(
1
2
f)andJ
n+1
(2πr),
it can be shown, see [14], that the resulting method to compute T
00
does not
suffer from loss of digits.
In this paper we extend the method of Nijboer and Zernike to the compu-
tation of the V
nm
in (2). Thus we insert Bauer’s formula (8) into the integral
at the right-hand side of (2) and interchange integral and summation. Then
one is faced with the problem of writing products R
0
2k
R
m
n
as a linear combi-
nation of Zernike polynomials with upper index m so as to be able to apply
(11). This problem was attacked by Nijboer in [15] for modest values of k,
m, n by employing recursion formulas for the Zernike polynomials. A sys-
tematic procedure for this does not seem to have been devised by Nijboer or
thereafter (also see [11], top of p. 535, for this point). By using results of the
relatively recent modern theory of orthogonal polynomials, as can be found
in [17], Secs. 6.8 and 7.1, we are able to find finite series expressions, with
favourable properties from a computational point of view, for the coefficients
needed in the “linearization” of the products R
0
2k
R
m
n
.
This paper has been organised as follows. In Sec. 2 we present the main
result and its derivation, and we comment on the nature and magnitude of
the linearization coefficients that yield the main result. In Sec. 3 we show
how the main result can be used to evaluate integrals T
nm
of the type that
occurs in (6). In Sec. 4 we briefly comment on how to extend the method to
the computation of integrals that occur in the case of high numerical aperture
for which a vector formalism as well as inclusion of the radiometric effect and
the state of polarization is required. In Sec. 5 we present some examples and
results of computations to compare the new scheme with the one based on
(3) and the one as follows from the results in [13] for T
2p,0
.
We conclude this section with some comments on different approaches
that can be found in the literature. Formally, the phase factor in the diffrac-
tion integral is proportional to the projection of the defocus distance onto
the direction of the plane wave contribution in the integral. In [1]-[2], the
evaluation of the diffraction integral is carried out using numerical Fourier
transform techniques, and the effect of the phase factor is mitigated by in-
tentionally introducing a compensating quadratic phase factor that can be
incorporated in the FT-scheme. An analytic approach to obtain stable ex-
pressions for strongly defocused fields is found in [3]-[4]. The radiometric
5
