1. INTRODUCTION

I first met Emil Wolf in 1997 during a meeting on laser beam characterization at the University of Laval in Quebec City, Canada. Since then, my career, similarly to that of all scientists who started studying optics in the early nineties, has considerably been influenced by Wolf’s charisma as well as by the huge deal of work Wolf and his co-workers produced over four decades, especially in the field of classical coherence theory, of which he was universally considered the father. Wolf will also be remembered as one of the authors of the Optics Bible, namely, the monumental Principle of Optics (Principles henceforth), written together with Born [1]. A story about Principles Emil Wolf enjoyed telling his guests (I was not an exception during a lunch at Rochester’s Institute of Optics) was the following [2]:

In 1956, Born was already in retirement. I was on a visiting appointment at New York University, still working on our book [2]. One day, I received a letter from Born in which he asked me why the manuscript was not yet finished. I wrote back saying that the manuscript is almost completed, except for a chapter on partial coherence on which I was still working. Born replied at once saying, “Wolf, who apart from you, is interested in coherence? Leave the chapter out and send the manuscript to the printers.” I finished the chapter anyway, and our book was published in 1959, only a few months before the invention of the laser, and many of the reviews of our book, which were then appearing, stressed that Principles of  Optics contained an account of coherence theory, which had become of crucial importance to the understanding of some features of laser light.

I had the privilege of coauthoring a couple of works with Wolf, the first of which saw the light in the tea room of the cafeteria of the Department of Physics and Astronomy in Rochester, not far from a photographic portrait Olivia Newton-John, bearing the dedication:

To Emil, with love—Olivia

Not many know that Olivia Newton-John was Max Born’s granddaughter.

The present paper has been conceived as a tribute to Wolf’s memory and legacy. It deals with a classical topic of classical optics, namely, the Nijboer–Zernike diffraction theory of aberrations [3–7]. In 1951, Wolf published an important review paper on the subject [8], which later formed the core of Ch. IX in Principles.

Despite its formal elegance, Nijboer’s theory did not receive the attention it deserved. The main reasons have partly been indicated by Nijboer himself but were clarified 20 years ago by Janssen and his co-workers [9–12], who developed an important extension of the Nijober formalism aimed at overcoming three main drawbacks that, for nearly 80 years, have confined Nijboer–Zernike’s theory as a beautiful but nearly useless practical tool: (i) entity of the phase distortion; (ii) entity of defocusing; and (iii) amplitude uniformity of the wavefront. From a mathematical perspective, one of the main technical problems is the evaluation of several paraxial 1D diffraction integrals containing Bessel functions of the first kind. Such integrals are the computational bricks with which Nijboer’s perturbative aberration theory was built up. When the diffracted field is required across the focal plane, such integrals can exactly be evaluated in terms of suitable linear combinations of Bessel functions of the first kind [4–6,8]. In the presence of defocusing, a quadratic phase term appears inside the integrals and makes their computation a dramatically more demanding task. It should be stressed that, as far as the practical solution of several problems in lithography, microscopy, and astronomy, Janssen’s extension of the Nijboer–Zernike aberration theory constitutes an optimal answer. Interested readers are encouraged to consult the important review published on the subject in 2013 [13].

The scope of the present paper is different. Here, we want to show that the “old” Nijboer–Zernike theory has probably been abandoned too early, and that new computational strategies can be conceived in order to extend its practical feasibility at least to deal with the presence of arbitrary defocusing. To this end, a short mathematical tour will first be presented within the realm of Tchebychev’s polynomials. Unique among several families of orthogonal polynomials is their capability to provide excellent low-order polynomial approximants of continuous functions. In Section 2, the basics of Tchebychev’s polynomials approximation properties are briefly resumed and employed on the zeroth-order Bessel function of the first kind. Section 3 constitutes the core of the paper: Nijboer–Zernike’s theory is numerically reformulated in such a way that a quantitative study of primary optical aberrations, namely, spherical aberration, coma, and astigmatism, can be carried out as well as for arbitrarily high defocusing. To this end, all Nijboer’s diffraction integrals are converted into quickly converging series expansions whose evaluation can be effectively implemented through simple recursive algorithms. One of the most important practical consequences is that the numerical effort needed to evaluate the diffracted wavefield turns out to be practically independent of the defocusing amount, a fact that alone should render such an approach particularly attractive. In order to validate our algorithm, several numerical checks will be carried out using numerous previously published results by Nijboer himself [1,4–6,8] as well as by Janssen and co-workers [9,10,12,13].

2. TCHEBYCHEV’S EXPANSION OF BESSEL FUNCTIONS

A. Approximation Problem Involving Tchebychev’s Polynomials

Tchebychev polynomials ${T_n}(x)$ constitute a well-known family of orthogonal functions across the real interval $[- 1,1]$, being ([14], Formula 8.949.3)

Property. Among all reduced polynomials of degree $n$, ${T_n}$ is the polynomial that deviates least from zero within the interval $[- 1,1]$.

By “reduced,” it should be mean a polynomial of the form ${x^n} + {a_{n - 1}}{x^{n - 1}} + \ldots + {a_0}$, i.e., such that the coefficient of the leading term be unitary. A simple proof of such property can be found, for instance, in [15] (Sec. 8.5.4). The “least deviation” property of Tchebychev’s polynomial is connected to its capability to build up excellent low-order polynomial approximations of continuous functions within the interval $[- 1,1]$. For the scope of the present paper, it is worth briefly illustrating such capabilities.

To this aim, consider the following approximation problem: given a continuous function $f(x)$, with $x \in [- 1,1]$, to find its best $n$ th-order polynomial approximation ${Q_n}(x)$ such that, given any other polynomial of degree $n$, say ${p_n}(x) \ne {Q_n}(x)$, the following inequality holds:

From what we’ve said here, the approximation provided by ${f_n}(x)$ will certainly be worse than that of the best one ${Q_n}(x)$, since ${\omega _n} \lt {{\cal S}_n}$, where ${{\cal S}_n}$ is defined by

Theorem. The ${L_\infty}$ norms ${\omega _n}$ and ${S_n}$ satisfy the following inequalities:

B. Tchebychev’s Expansion of ${J_0}$ Function

Bessel functions of the first kind play a dominant role within paraxial optics. In the Nijboer–Zernike aberration theory, the optical field generated by aberrated systems is represented via superposition of diffraction integrals whose integrands are products of Bessel functions, Zernike polynomials, and quadratic phase factors, the integration domain being (in suitable units) the interval [0, 1]. To overcome the difficulties of evaluating such highly oscillating (especially in the presence of nonnegligible defocusing) integrals, suitable representations of Bessel functions like ${J_\ell}(vx)$ (with $\ell \ge 0,v \ge 0$) for $x \in [0, 1]$ will be employed. To illustrate the main idea, consider first the zeroth-order Bessel function ${J_0}(vx)$, for which the following expansion holds {[16], Formula 9.3.6(3)}:

Equations (10) and (11) have the form in Eqs. (4) and (5), respectively, and, in this particular case, it is possible to analytically evaluate all quantities into Eq. (7) except, of course, ${\omega _n}$. This implies that it should be possible to quantitatively compare the performances of Tchebychev’s approximant into Eq. (11) and the (unknown) best approximant ${Q_n}(x)$. The upper limit ${\cal S}_n^ \star$, defined by Eq. (8), is first evaluated, yielding

In Fig. 1, an example of applications of the above inequalities is illustrated, also to give readers an idea about the powerfulness of Tchebychev’s approximation. In particular, in the left column, the behaviors of the quantities ${{\cal A}_n}$ and ${\cal S}_n^ \star$ are shown, as a function of $v$, as solid and dashed curves, respectively, for (a) $n = 4$, (c) $n = 8$, and (e) $n = 12$. The thin vertical line is located at the abscissa $v = 10$. It can be appreciated how, for a given $v$, on sufficiently increasing $n$ the curves representing ${A_n}$ and ${\cal S}_n^ \star$ becomes practically indistinguishable. According to Eq. (7), this implies that the $n$th-order polynomial approximant ${f_n}(v;x)$ given in Eq. (11) is expected to be close to the best $n$th-order polynomial approximant ${Q_n}(x)$ of the Bessel function ${J_0}(vx)$.

 

Fig. 1. Left column: behaviors of the quantities ${{\cal A}_n}$ (solid curve) and ${\cal S}_n^ \star$ (dashed curve), as a function of $v$, for (a) $n = 4$, (c) $n = 8$, and (e) ${n = 12}$. The thin vertical line is located at the abscissa $v = 10$. Right column: plots of polynomial approximations of the function ${J_0}(10x)$ (dashed curves) within the interval $x \in [0 ,1 ]$, (b) $n = 4$, (d) $n = 8$, and (f) $n = 12$. The solid curve denotes the exact profile of the Bessel function.

Download Full Size | PDF

This is confirmed, at least at a visual level, by the right column, where plots of different polynomial approximations of the function ${J_0}(10x)$ are shown (dashed curves) for the same values of the truncation $n$, namely, (b) $n = 4$, (d) $n = 8$, and (f) $n = 12$. In the next section, such a capability of Tchebychev’s polynomials to generate excellent low-order approximations will be employed to numerically revisit and improving the Nijboer–Zernike aberration perturbative theory.

3. APPLICATION TO NIJBOER–ZERNIKE’S DIFFRACTION THEORY OF ABERRATION

A. Brief Résumé of Nijboer’s Formalism

In a perfect optical system, the light waves, which proceed from each element of the object, emerge in the image space as convergent waves spherical in form. Such an instrument represents an idealization that cannot be realized in practice; in the image space of an actual instrument, the wave surfaces will, as a rule, be of a more complicated form.

These are the opening words of an important historical review of the subject of aberration theory, written in a 1951 paper by Wolf [8]. The paper formed the core of Ch. IX in Principles, one of the most beautiful of the whole opera. In the rest of paper, the Nijboer–Zernike diffraction aberration theory [4–6] will be quantitatively revisited in the light of that reported in the previous section. To this end, readers are encouraged to first go through ([1], Ch. IX) or even through [4–6,8] in order to familiarize with the notations as well as the meaning of the symbols, which will be employed in the rest of the paper.

Given a centered imaging optical system, the deformation from the perfectly spherical convergent wavefront in the region of the exit pupil, measured in reduced wavelength units $\rlap{--} \lambda = \lambda /2\pi$, will be denoted $\Phi (\tau ,\vartheta)$, where $(\tau ,\vartheta)$ represent polar coordinates upon the exit pupil plane, normalized so that $\tau$ ranges from 0 to 1.

Further, the image space will be parametrized by dimensionless cylindrical coordinates $(u,v,\varphi)$: the variable $u$, related to the Fresnel number, represents a dimensionless measure of the position of the transverse observation plane, with respect to the focal one (corresponding to $u = 0$). The couple $(v,\varphi)$ represents (dimensionless) polar coordinates across the typical plane $u = {\rm const}$. The diffracted wavefield, say $\Psi (u,v,\varphi)$, will be expressed, within paraxial approximation and apart from unessential overall phase and amplitude factors, by the following integral [8]:

For an optical system free of any aberration, i.e., when $\Phi \equiv 0$, the integral in Eq. (14) turns out to be equal to $2{\cal L}_0^0(u,v)$, where [1]

The perturbative Nijboer approach to the diffraction theory of aberrations was conceived to deal with “small” values of $|\Phi |$, typically of the order of unity. To this end, the aberration function $\Phi$ is first written in the following form:

B. Numerical Computation of Nijboer’s Diffraction Integrals

The numerical computation of the contributions into the series expansion in Eq. (18) is the central practical problem of such perturbative approach to aberration theory. Nijboer developed a clever analytical treatment based on basic properties of Zernike’s polynomials in order to represent each integral in Eq. (B1) through a suitable series of Bessel functions. In particular, such series terminate when $u = 0$, i.e., across the focal plane, thus allowing the sequence $\{{{\cal C}_k}(0,v;n,m)\} _{k = 0}^4$ to be exactly found. The same is no longer true when the computation is required out of the focal plane, i.e., for $u \ne 0$.

In order for the practical feasibility of Nijboer’s theory to be extended for dealing with any couple $(u,v)$, the following idea will be pursued: given a primary aberration described by Eq. (17), each power of the Zernike polynomial $R_m^n(\tau)$ into the diffraction integrals of Eq. (B1) will be expanded as powers of $\tau$. In this way, each term of the sequence $\{{{\cal C}_k}(u,v;n,m)\} _{k = 0}^4$ will be recast as a linear combination of generalized Lommel integrals, as shown in Eq. (16). The latter can be converted into rapidly converging series by using the generalization of Eq. (10) to expand Bessel functions ${J_\ell}(v\tau)$, with $\ell \ne 0$, in terms of Tchebychev polynomials. In particular, the following expansions will be employed ([16], Sec. 9.3.6):

For $\ell \ne 0$, the lower bound ${{\cal A}_n}$ of the maximum deviation ${{\cal S}_n}$, defined in Eq. (9), can still be exactly computed for any $\ell$, as detailed in Appendix A, which yields

In the following, each generalized Lommel integral in Eq. (16) will be replaced by its $N$th-order approximant, defined by

C. Preliminary Numerical Checks

Since the central idea is to replace, into Eq. (16), the Bessel function ${J_\ell}(v\tau)$ via its $N$ th-order approximant ${f_N}(v,\tau)$, we introduce the absolute error, say ${\cal E}$, i.e., the difference between the exact value in Eq. (16) and the approximated one in Eq. (24), which can be recast as follows:

 

Fig. 2. Behavior of the absolute error in Eq. (28) obtained by evaluating Lommel’s function ${\cal L}_0^0(0,v)$ via Eq. (24) as a function of the truncation order $N$ (dots), for (a) $v = 10$, (b) $v = 20$, (c) $v = 30$, and (d) $v = 40$. The open circles represent the values of the upper bound given by Eq. (29).

Download Full Size | PDF

For instance, when $u = 0$, each of the integrals into Eq. (B1) have been evaluated in closed form by Nijboer in the case of primary aberrations, e.g., spherical aberration, coma, astigmatism. It is then rather natural to use Nijboer’s own results to carry out preliminary checks about the numerical efforts required by the present method.

The simplest test consists in evaluating Lommel’s integral ${\cal L}_0^0(u,v)$, which, across the focal plane, $u = 0$ reduces to

Another example of exact evaluation of Lommel’s function ${\cal L}_0^0$ is when $u = v$, which yields [17]

Similarly, as done in Fig. 2, in Fig. 3 the behavior of the absolute error in Eq. (28) obtained by evaluating Lommel’s function ${\cal L}_0^0(v,v)$ via Eq. (24) is plotted against the truncation order $N$ (dots) for (a) $v = 10$, (b) $v = 20$, (c) $v = 30$, and (d) $v = 40$. The open circles represent again the values of the upper bound given by Eq. (29).

 

Fig. 3. Behavior of the absolute error in Eq. (28) obtained by evaluating Lommel’s function ${\cal L}_0^0(v,v)$ via Eq. (24) as a function of the truncation order $N$ (dots), for (a) $v = 10$, (b) $v = 20$, (c) $v = 30$, and (d) $v = 40$. The open circles represent the values of the upper bound given by Eq. (29).

Download Full Size | PDF

Per Figs. 2 and 3 (and for several others’ error behaviors not show here), it is possible to understand why the upper bound in Eq. (29) is not numerically significative to guarantee an efficient choice of the truncation order $N$.

Nevertheless, it can be appreciated, within the behavior of the absolute error, a sort of plateau, which extends for values of the truncation orders up to, roughly, $v/2$. When $N \gt v/2$, the error reduction considerably increases, according to the exponential law represented in figures by thick black lines. The following rule for the truncation error choice will then be adopted:

D. Recurrence Relation for von Lommel’s Integrals

The numerical evaluation of the generalized von Lommel integrals ${\cal L}_\ell ^m(u,v)$ can further be improved by deriving a simple but computationally important recurrence relation, based on a well-known property of Bessel functions. To this end, consider Eq. (16) written for $\ell + 1$, i.e.,

 

Fig. 4. Recurrence relation in Eq. (35).

Download Full Size | PDF

It is worth giving readers at least an idea about the computational times of the present method, compared to those pertinent to standard integration packages of Mathematica. To this end, the evaluation of the quantity ${\cal L}_{25}^{- 1}(60,2\pi)$ has been carried out on using both Eq. (24) and the Mathematica command NIntegrate. The result of the computational experiment is reported in Table 1. The first column contains the value of the parameter $K$ to be inserted into Eq. (32) in such a way a relative error of the order of magnitude given in the fifth column can be reached. The third column contains the value of the corresponding WorkingPrecision parameter to be set into the NIntegrate command in order for the same relative error order to be achieved. The second and the fourth columns contain the computational times (expressed in seconds). The relative errors have been computed with respect to an “exact” value of ${\cal L}_{25}^{- 1}(60,2\pi)$ evaluated again via Mathematica NIntegrate command, but on setting WorkingPrecision to 200. From the table it is possible to realize that the computational speed of the algorithm of Eq. (24) turns out to be, for a given relative error order, about one order of magnitude greater than the speed of the standard Mathematica integration package, regardless of the level of accuracy required. The results of several other numerical experiments, not reported here, fully confirmed such a speed difference.

Table 1. Evaluation of the Quantity ${\cal L}_{25}^{- 1}(60,2\pi)$^a

View Table

 

Fig. 5. Behavior of the relative error obtained by evaluating the diffraction integral $V_n^m(u,v)$ defined through Eq. (36) via the expansions given in Eq. (37) (solid curve) and in Eq. (40) (dots). The chosen values of the parameters are $n = 25$, $u = 60$, and $v = 2\pi$, while (a) $m = 1$, (b) $m = 5$, (c) $m = 15$, and (d) $m = 25$. The “exact” value has been evaluated through the standard integration package of Mathematica.

Download Full Size | PDF

E. Numerical Comparison with Janssen’s Extended Nijboer–Zernike Theory

The present section will be devoted to present numerical comparisons between our algorithm and Janssen’s extended Nijboer–Zernike [9,10,13]. Although the evaluation of generalized von Lommel’s integrals in Eq. (16) was also considered in [10], here we prefer to employ our algorithm to evaluate directly the following integral:

In [13], Haver and Janssen pointed out how the convergence of the series into Eq. (37) becomes numerically problematic when $u$ becomes greater than about 50. This is due in part to the fact that the truncation order of Janssen’ series grows with the modulus of $u$. To overcome such a problem, Janssen and co-authors gave the series in Eq. (37) new, computationally more efficient forms, which will not be considered here [12,13]. We shall limit ourselves to work on a single example, namely, the computation of $V_n^m(u,v)$, with $n = 25$, $u = 60$, and $v = 2\pi$, and for different (odd) values of $m$.

To this end, it is worth recalling the explicit definition of Zernike polynomial’s $R_n^m(\tau)$, namely,

In Fig. 5, the behavior of the relative error, against the truncation orders, obtained by evaluating the diffraction integral $V_{25}^m(60,2\pi)$ via the expansions given in Eq. (37) (solid curve) and in Eq. (40) (dots), is plotted for (a) $m = 1$, (b) $m = 5$, (c) $m = 15$, and (d) $m = 25$. The “exact” value has again been evaluated through the standard integration package of Mathematica. It should be appreciated how the numerical effort required to Eq. (40) to evaluate the integral is substantially independent of $m$, being the truncation order given by Eq. (32), with $K$ not exceeding 23. Moreover, the error reduction is monotonically decreasing on growing the truncation order, i.e., the parameter $K$. The above are the main differences with respect Janssen’ series, whose truncation order displays a clear dependance on $m$, and the partial sums into Eq. (37) do not monotonically decrease on increasing their order.

 

Fig. 6. Isophotes in the presence of third-order spherical aberration of the amount $\beta = 1/2$. The figure should be visually compared with Fig. 3 in [8].

Download Full Size | PDF

 

Fig. 7. Same as in Fig. 6, but for a fifth-order spherical aberration. The figure should be visually compared with Fig. 4 in [8].

Download Full Size | PDF

F. Third- and Fifth-Spherical Aberration

The first example of application of our algorithm to Nijboer–Zernike’s theory concerns spherical aberration, characterized by the following phase function inside the integral in Eq. (14):

Consider now third- and fifth-order spherical aberrations, corresponding to set into Eq. (41) $n = 4$ and $n = 6$, respectively. The parameter $\beta$ will be set to 1/2, in order to reproduce Figs. 3 and 4 in [8]. To this end, it is sufficient to expand the factors ${[R_n^0(\tau)]^k}$ into the integrals in Eq. (42) and then to apply the definition in Eq. (16) to recast them as a linear combination of generalized Lommel integrals ${\cal L}_0^m(u,v)$. For space reasons, the explicit expressions of ${{\cal C}_k}$’s will be not reported here. Each function ${\cal L}_0^m(u,v)$ is then numerically evaluated by suitably choosing the truncation order $N$ via Eq. (32). Results are shown in Figs. 6 and 7, which are aimed at reproducing the isophotes of Figs. 3 and 4 in [8], respectively. The visual agreement is perfect, and several numerical tests (not shown here) were performed directly on the exact analytical expressions of the focal field, for further checking purposes.

An important comment, of a numerical nature, in regard to the above figure: once the truncation order $N$ has been fixed, our algorithmic efforts are independent of the variabile $u$, i.e., of the position of the observation plane. This represents, in our opinion, one of the most important features of the present approach, which is basically due to the marriage between the simplicity and elegance of Nijboer’s analytical treatment and the extraordinary efficiency of Bessel function Tchebychev’s expansions, which allows detailed maps of the diffracted field to be produced in a few minutes. To give an idea, the time needed to generate Figs. 6 and 7, both containing about 20,000 sampled points, was about 7 min, working with an unoptimized Mathematica code running onto a commercial laptop equipped with a 3 GHz Intel Core i7 processor and 16 GB RAM.

 

Fig. 8. Isophotes at the plane $u = 0$ in presence of a coma in Eq. (43) of amount (a) $\beta = 1$ and (b) $\beta = 3$. The figures should be compared with Figs. 5 and 6 in [8], respectively.

Download Full Size | PDF

G. Primary Coma

The coma aberration phase is [1,5]

In Fig. 9, a slightly modified version of the same code has been used to reproduce the numerical results presented in [9]. In particular, to adhere as much as possible with Janssen’s simulations, the aberration phase in Eq. (43) has been changed into ${\Phi _{{\rm coma}}} = \beta {\tau ^3}\cos \vartheta$, with $\beta = 1$. Isophotes of the diffracted field have been evaluated across the transverse planes (a) $u = 0$, (b) $u = \pi /2$, (c) $u = 2\pi$, and (d) $u = 4\pi$. Each figure contains a matrix of about $300 \times 200$ sampled points, with computational times of about 30 s per figure. Also in that case, the visual agreement with the results obtained by Janssen is excellent, once the considerably different densities of sampled points of the figures have been taken into account.

H. Primary Astigmatism

The last example of application deals with astigmatism, whose phase model ${\Phi _{{\rm ast}}}$ is

 

Fig. 9. Isophotes in presence of coma ${\Phi _{{\rm coma}}} = 3\beta {\tau ^3}\cos \vartheta$, for $\beta = 1/3$ at different transverse planes: (a) $u = 0$, (b) $u = \pi /2$, (c) $u = 2\pi$, and (d) $u = 4\pi$. The figure should be compared with Fig. 2 in [9].

Download Full Size | PDF

 

Fig. 10. Isophotes at the plane $u = 0$ in the presence of an astigmatism in Eq. (44) of amount (a) $\beta = 1$ and (b) $\beta = 3$. The figures should be compared with Figs. 7 and 8 in [8], respectively.

Download Full Size | PDF

Finally, Fig. 11 shows examples of field computation at different transverse planes as shown for $\beta = 1$. Isophotes of the diffracted field have been produced for (a) $u = 0$, (b) $u = \pi /2$, (c) $u = 2\pi$, and (d) $u = 4\pi$. Each figure contains a matrix of about $300 \times 200$ sampled points, with computational times of about 30 s per figure.

 

Fig. 11. Isophotes in the presence of astigmatism in Eq. (44) of amount $\beta = 1$ at different transverse planes: (a) $u = 0$, (b) $u = \pi /2$, (c) $u = 2\pi$, and (d) $u = 4\pi$.

Download Full Size | PDF

4. CONCLUSION

The present paper has been conceived as a tribute to Emil Wolf’s memory and legacy in respect to his authorship of Principle of Optics. The paper’s content was also intended to honor Wolf as a mathematician (he earned his PhD in math from Bristol University in 1948).

What has been celebrated is the marriage of two mathematical jewels: the Nijboer–Zernike aberration theory and the Tchebychev polynomial approximation theory. The result is a powerful computational tool able to evaluate, up to the desired degree of accuracy, all terms of classical perturbative treatment of primary aberrations for, in principle, any values of the physical parameters in the observation space. In particular, the computational effort required to achieve such a task turns out to be independent of the defocusing, a fact that could greatly improve the practical applicability of Nijboer–Zernike’s theory to a variety of scenarios, from vectorial approaches to aberration theory [11], to the numerical evaluation of Fresnel transforms [18,19]. We believe that the present paper could also open mathematical questions worthy of consideration, e.g., the use of different families of polynomials than Tchebychev to approximate continuous functions into the interval $[- 1,1]$, for instance, Legendre polynomials [20]. Another concerns the summability of the complete Nijboer’s perturbative series, in order to also deal with nonsmall values of $\beta$. To this end, an (symbolic language) algorithm for generating functions ${{\cal C}_k}(u,v;n,m)$ of arbitrary orders $k$ should first be conceived to also grasp information about their asymptotic behavior for $k \gg 1$. Then, suitable convergence acceleration algorithms, for instance, based on Padé approximants or on Levin-type nonlinear sequence transformations, could be used to practically sum Nijboer’s series.

In conclusion, no words would be better: