Computing Zernike polynomials of arbitrary degree using the discrete Fourier transform

The conventional representation of Zernike polynomials Rn (ρ) gives unacceptable numerical results for large values of the degree n. We present an algorithm for the computation of Zernike polynomials of arbitrary degree n. The algorithm has the form of a discrete Fourier (cosine) transform which comes with advantages over other methods in terms of computation time, accuracy and ease of implementation. As an application we consider the effect of NA-scaling on the lower-order aberrations of an optical system in the presence of a very high order aberration. [DOI: 10.2971/jeos.2007.07012]


INTRODUCTION
The (radial part of the) Zernike polynomials R m n (ρ) are widely used in the representation of the aberrations of optical systems and in the computation of the diffraction integral defining the point-spread function of these systems [1]- [5].When we are dealing with smooth exit pupil functions, it is, in general, sufficient to consider the R m n for modest values of the degree n and azimuthal order m, say, n + m ≤ 12.For such pupil functions, the conventional polynomial representation [ can be used to calculate the Zernike polynomials.Some low order Zernike polynomials are shown in Table 1.In the case that the exit pupil function contains discontinuities, or is roughly behaved in a more general sense, it is necessary to consider Zernike polynomials of much higher degree and order.For instance, when the pupil function has a central obstruction (0 and 1 on two concentric sets 0 ≤ ρ < a and a ≤ ρ ≤ 1), the coefficient of R 0 n (ρ) in the Zernike expansion of the pupil function can be shown to decay only like n −1/2 .Then Eq. (1) becomes cumbersome because of the high-order factorials that are required.Also, for m = 0, it can be shown from the Stirling's formula that the largest coefficient of ρ n−2s occurring in the series in Eq. (1) behaves like (1 + √ 2) n .Accordingly, when computing with d digits, Eq. (1) produces errors of the order of unity or larger from n = d/log(1 + √ 2) onwards.Hence, for the commonly used 15 digits precision, one has serious problems from n = 40 onwards, as shown in Figure 1, top, right.An alternative to compute Zernike polynomials is to use recursions for them such as those found in Ref. [6].These recursion schemes are, however, computationally more expensive and less direct than a formula like Eq. ( 1) and their accuracy due to error propagation is also an issue.
In this paper, we present a new computation scheme in which one can allow degrees as large as 10 5 without problems.This new algorithm is of the discrete-cosine transform (DCT) type, and is direct and transparent.Furthermore, the computation can be done using the FFT-algorithm which comes with the following advantages [7,8]: • Very favorable and well-established accuracy • Simultaneous computation of all Zernike polynomials of the same degree n in as few as O(n logn) operations As an application we consider the effect of NA-scaling on the lower-order aberrations of an optical system in the presence of a very high order aberration.For this we use a recently found formula [9], entirely in terms of Zernike polynomials, for the Zernike coefficients of scaled pupils.

DCT FORMULA FOR ZERNIKE POLYNOMIALS
In Appendix A we show that where N is any integer > n + m.In Eq. ( 2) we have integer n, m ≥ 0 with n − m even and ≥ 0 (as usual), and U n is the Chebyshev polynomial of the second kind and of degree n.No matter how large n is, the evaluation of U n (x) is no problem since we have Eq. ( 2) is a consequence of the formula (A.10) that represents R m n (ρ) as an integral of a trigonometric polynomial of degree n + m over a periodicity interval.Such an integral can be computed error-free as a series when the number of equidistant sample points N exceeds the degree n + m.It also follows from this that the right-hand side of Eq. ( 2) (and that of Eq. (A.10)) vanish when m > n or when n and m have different parity (with again N > n + m).
Eq. ( 2) gives R m n (ρ) as the m th component of the DCT of the sequence (U n (ρ cos 2πk/N)) k=0,1,...,N−1 , hence we get all R m n (ρ), with m ≥ 0 and m = n, n − 2, ... , using O(N log N) operations.Since m ≤ n, it is sufficient to take N any integer > 2n.
In Figure 1, bottom, we show R m n (ρ), computed according to Eq. ( 2), with m = 0 and n = 10000 and ρ very close to 1.We see that R m n (ρ = 1) = 1 which is in agreement with the theory [1].

AN APPLICATION: HIGH-ORDER ABERRATIONS AND SCALING
In lithographic imaging systems, the numerical aperture (NA) is varied intentionally below its maximum value so as to optimize the performance for the particular object to be imaged.In Ref. [9] the effect of NA-scaling on the Zernike coefficients describing the optical system has been concisely expressed in terms of Zernike polynomials.Thus we consider a pupil function in polar coordinates with real phase Φ, and we assume that Φ is expanded as a Zernike series according to where the summation is over n = n, n + 2, ... (R n+2 n ≡ 0).In case of a non-smooth phase function Φ, one should expect significant values of α m n for very high degrees n .Also, scaling is normally done using values of ε close to its maximum 1, where Eq. (1) produces the largest numerical error.Thus, formula ( 6) is not practicable in these cases when Eq. ( 1) is used to evaluate R n n (ε) − R n+2 n (ε), but becomes so when Eq. ( 2) is used instead.As an example, we consider the effect of a single high order aberration term α m n on the totality of α m n with n = m, m + 2, ..., n while scaling to relative size ε.We take α m n = 0 for n = n and α m n = 1, and get from Eq. ( 6)

FIG. 2
FIG.2The disturbance α 0 n (ε) of the aberration of order n = 0, 2, • • • , 100, due to the presence of an aberration of amplitude 1 and of the order n = 100 when the system is scaled to relative size ε = 0.50 (left) and ε = 0.98 (right).

TABLE 1
Low order Zernike polynomials