## Abstract

Assessment of the performance degradation caused by the mid-spatial frequency (MSF) structure on optical surfaces often relies on a perturbation method that dovetails with the familiar sequence of models based on geometrical and physical optics. In the case of imaging systems, the perturbative step yields estimates of wavefronts in the exit pupil which are, in turn, used to extract performance measures such as MTF, PSF, and Strehl ratio. To date, the validity of that perturbation appears to be poorly understood. We present methods to estimate the errors of this approach and thereby arrive at a rule of thumb for its accuracy: the error is approximately equal to the RMS of the MSF structure at its source multiplied by the square of the ratio between a particular Fresnel zone size and a characteristic length of the MSF structure.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Mid-spatial frequency (MSF) structure is inescapable in most aspheric and freeform optical systems due to the subaperture tools involved in manufacturing the optical components. Their characteristic frequencies lie between those of low-order aberrations and high-order scattering, and they affect optical performance in ways that continue to present challenges [1–4]. In particular, the tolerancing of parts is complicated by MSF structure and a perturbative model is often employed to simplify this process [5]. A similar approach was initially used to account for perturbations like misalignment [6–12]. This model is an approximation where perturbations, such as the MSF phase structure, are effectively dragged along the nominal rays of the system. In this way, it becomes straightforward to approximate the MSF phase impact in the exit pupil of the system for assessing system performance. An important advantage offered by the perturbation model is that, while the low-order aberrations are inherent to the system’s design, the MSF structure varies from part to part within a lot. The fact that perturbation employs the nominal rays removes the need to perform new ray tracing for each realization of the MSF structures. The validity of this perturbation model is not well understood, however.

In what follows, we derive simple estimates for the error associated with this perturbation. To do so, we first present an asymptotic approach for modeling the propagation of a monochromatic field in free space. This approach is based on using the rays determined by an initial nominal phase distribution, and modeling how these rays can carry the impact of an extra phase deviation that is taken to be at a lower asymptotic order (namely the MSF). Our goal is not to treat the details of how MSF surface structures become imprinted into the phase at the interfaces, but to focus on the propagation of that phase structure in homogeneous media. This asymptotic approach is then used to derive rules of thumb for understanding the limitation of the perturbation model when applied to imaging systems. For simplicity, the derivation is given for two dimensions; the results for three dimensions will be presented separately and augmented with practical design examples.

## 2. Asymptotic propagation estimate based on nominal rays

In two dimensions, a monochromatic scalar field Re[*U* (*x, z*) *e*^{−i}* ^{ωt}*] in a homogeneous transparent medium is governed by the Helmholtz equation

*k*=

*ω*/

*c*= 2

*π*/

*λ*, with

*λ*being the wavelength in the medium. This field is taken to be propagating towards larger values of

*z*and at some reference plane, say

*z*=

*z*

_{M}, it is nominally given by

*U*(

*x, z*

_{M}) =

*U*

_{0}

*A*(

*x*) exp[i

*kW*(

*x*)], where

*U*

_{0}is a constant with field units. Here,

*W*(

*x*) accounts for the nominal wavefront shape and

*A*(

*x*) is the nominal field amplitude, both assumed to vary slowly within the scale of the wavelength. Importantly, we now consider the case where, at this plane, the field also carries an additional MSF phase of the form exp[i

*ϕ*(

*x*)], where

*ϕ*(

*x*) is measured in radians and is typically smaller than

*π*. Without loss of generality, we take

*ϕ*to have a mean of zero. The goal in this section is to derive a simple estimate of how this MSF phase impacts the field under propagation in

*z*.

We begin by writing the field as *U*(*x, z*) = *U*_{0} exp[i*kϕ*(*x, z*)], where the complex function Φ(*x, z*) accounts for spatial variations in both the phase and amplitude. With this, Eq. (1) becomes

Equation (2) can be solved upon expressing Φ as an asymptotic series with parameter (i*k*)^{−1}:

By using Eq. (3) with Eq. (2) and separating terms of equal power of the asymptotic parameter *k*, we arrive at

The initial conditions discussed above can now be stated as Φ_{0}(*x, z*_{M}) = *W*(*x*), Φ_{1}(*x*, *z*_{M}) = ln[*A*(*x*)] + i*ϕ*(*x*), and Φ* _{N}*(

*x, z*

_{M})

*=*0 for

*N*≥ 2. It is now possible to work to progressively higher orders by integrating in

*z*at each order from these initial conditions.

This sequence starts with Eq. (4), the well-known Hamilton-Jacobi or Eikonal equation, which can be solved in terms of the nominal rays by using the following parametrization:

where*ξ*is the

*x*coordinate at the reference plane (

*z*=

*z*

_{M}) and

*s*is the arclength along the ray. The direction of each ray is given by the unit vector

*W*′(

*ξ*)

*, χ*(

*ξ*)], where $\chi \left(\xi \right)\triangleq \sqrt{1-{W}^{\prime}{\left(\xi \right)}^{2}}$, with ≜ denoting definition. It is shown in Appendix A that the solution to the Eikonal equation is simply given by its initial value at

*z*=

*z*

_{M}, namely

*W*(

*ξ*), plus the length of the ray: where an overline on any function

*f*(

*x, z*) indicates that it is being expressed in terms of the ray parameters (

*ξ, s*) by using Eq. (6). That is, $\overline{f}\left(\xi ,s\right)\triangleq f[x(\xi ,s),z(\xi ,s)]$. Notice that, as a consequence of placing

*ϕ*one asymptotic order below

*W*, the rays used in this analysis are the nominal rays, not affected by the MSF structure. It is subsequently shown in Appendix A that Eq. (5) can also be solved in terms of this ray parametrization. For

*N*= 1, the parametrized solution is

*ξ, s*) is the Jacobian determinant of the coordinate transformation:

The first term of Eq. (8) accounts for the change in the amplitude due to the bunching or spreading of the rays under propagation; the second term indicates that, asymptotically, the effect of the MSF phase structure can be modeled by simply dragging this phase along the rays. This is precisely the perturbation model.

The leading contribution to the error incurred by using the perturbation model follows upon considering the next term in the asymptotic series. It is shown in Appendix A that this contribution is given by

*ϕ*, and

That is, all the dependence on *z* of this *ϕ*-dependent correction is as a quadratic in *ζ*. As is evident from Eqs. (10) and (11) and will become clearer in Sec. 4, the character of diffractive solutions means that the error of the perturbation model can remain finite in the limit of large propagation distances. It turns out that −*χ*^{3}/*W*″ is the displacement between *z*_{M} and the plane at which the ray in question crosses its neighbors, forming a caustic, i.e., the *z*-coordinate of this caustic point is *z*_{M} − *χ*^{3}(*ξ*)/*W*″(*ξ*). This follows from Eq. (9) or the observation that

*∂x*/

*∂ξ*being calculated for fixed

*z*.

In summary, the perturbation model can be expressed as follows:

*ϕ*-dependent component in $\overline{{U}_{\mathrm{P}}}$ is that contained in $\overline{{\Phi}_{1}}$. The leading correction to this approximation results from also considering the

*ϕ*-dependent terms in $\overline{{\Phi}_{2}}$ given in Eq. (10):

## 3. Simple field error estimates in a homogeneous medium

The relative root-mean-squared error (RMSE) of the perturbation model, *ϵ*, can be estimated as a function of propagation distance by integrating over the transverse plane the squared modulus of the difference between *U*_{P} and the corrected field estimate in Eq. (14). This is achieved by substituting *s* = (*z* − *z*_{M})/*χ*(*ξ*) in the parametrized estimates, and changing the variable of integration from *x* to *ξ* by using Eq. (12). Furthermore, we keep only the term linear in *ζ* of the exponential in Eq. (14) since it is initially dominant for small *z* − *z*_{M} and, as seen in what follows, it gives an adequate estimate in the domains of interest here. The squared RMSE is then

*z*, due to the fact that, when evanescent waves are not included, free propagation is a unitary operation. Upon defining the weighted average

*ϵ*

^{2}:

Further, both *A* and *χ* typically have a significantly slower dependence on *ξ* than *ϕ*. In this case, when their behavior is more or less homogeneous over the stop, an even simpler approximation is given by

The range of validity of the perturbation approximation can now be estimated by using Eq. (18). In particular, this approximation is useful when (*z, z*_{M}; *ϕ*) is small compared to the intrinsic impact of the MSF structure on the wavefield, referred to here as *φ* and defined according to

*φ i*s independent of z, and in the last step we assumed |

*ϕ*| ≪

*π*. Clearly, if

*ϵ i*s comparable to

*φ*, the error of the perturbation model is then comparable to that of just ignoring the MSF structure. We can therefore regard the perturbation model as valid whenever

*ϵ i*s sufficiently smaller than ϕ. For the sake of discussion, we can require

*ϵ < ϕ*/3.

As a first simple test to the RMSE estimate, consider a nominally uniform collimated beam [*W*(*ξ*) = 0, *A*(*ξ*) = 1, *χ*(*ξ*) = 1] incident on a surface at *z*_{M} that imparts on it a MSF phase *ϕ*. In this case, Eq. (11) becomes *ζ* = *z* − *z*_{M} and Eq. (18) reduces to

When the MSF phase structure is given by a simple sinusoid of the form

where*h*is the amplitude measured in radians and

*κ*is its spatial frequency, this RMSE estimate is then given by

*h*≪

*π*, which is justified by the fact that typical MSF phase structures have amplitudes that are much smaller than a wave,

It is convenient in what follows to normalize the RMSE by using *φ* of Eq. (19). Notice that in this case, $\phi ={\u3008|1-\mathrm{exp}(\mathrm{i}\varphi ){|}^{2}\u3009}_{1}=\sqrt{2[1-{J}_{0}(h)]}\approx h/\sqrt{2}$, where *J*_{0} is the zeroth order Bessel function of the first kind. Figure 1 shows plots of this normalized RMSE (NRMSE), defined as *ϵ*/*φ*, for *λ* = 632.8 nm and several values of *κ*, as a function of $|z-{z}_{\mathrm{M}}|/\mathcal{Z}$, where

with the assumption of small *h* being used in the last step. In these plots, the approximation of small *h* for both *φ* and $\mathcal{Z}$ is used for the normalization of the axes. To place these results in the same terms as those used later (corresponding to imaging systems), the results in Fig. 1 use *κ* = *C*/*L*, where *L* is the length of the aperture (set to 2 cm) and *C* is the number of cycles of the sinusoidal MSF phase across this aperture. Evidently, the NRMSE of the examples where *h* ≤ *π*/8 are well approximated by Eq. (22), for values of *z* for which *ϵ*(*z, z*_{M}; *ϕ*) is below *φ*, as the corresponding curves in Fig. 1 are nearly indistinguishable. In fact, the maximum deviation is less than 12% for *h <π*/4. That is, the error estimate in Eq. (22) is accurate for propagation by distances smaller than the characteristic distance $\mathcal{Z}$, which is the value of *z* at which the error estimate reaches an NRMSE of 1. Notice also that, in this weak MSF phase regime, this upper bound $\mathcal{Z}$ is independent of the magnitude of the MSF and, given the periodicity of the MSF, it happens to be the Talbot distance divided by 2*π*.

Keep in mind that the MSF impact is roughly proportional to *h*. However, our analysis is aimed at analyzing the validity of the perturbation model and, for small *h*, validity is found to be asymptotically independent of *h*.

## 4. Application to imaging systems

The ideas presented above are now used to derive rules of thumb applicable to imaging systems for which it is convenient to work solely in image space. In such cases, to a good approximation, the wave field generated by an object point now converges onto a point on the image plane. Akin to the standard practice for the stop and pupil, it is effective to assume that the frequency band of interest in the MSF phases generated by each optical surface in the system is faithfully resolved at its associated conjugate plane. Those conjugate planes can sit before or after either the image plane and/or the exit pupil. Given this assumed fidelity in those images, the dominant error of the perturbation model is then generated by dragging these MSF phases along the nominal rays from their respective conjugate planes to the exit pupil. (It is shown in Appendix B that working in image space is equivalent to working in any other conjugate space, e.g. at the stop.) By using the results derived above, we are now in a position to estimate the accuracy of that step. In this way, the effects of the MSF phase as well as of the nominal aberrations and diffraction from the aperture stop can all ultimately be estimated by using a single wave propagation integral from the exit pupil to the image plane.

Consider first the MSF introduced by a single optical surface within the system whose conjugate is at *z* = *z*_{M} and where the exit pupil plane is at *z* = *z*_{P} with coordinates chosen such that the system’s image plane is at *z* = 0, as shown in Fig. 2. To simplify the treatment, we consider the on-axis object point, whose ideal image is at the origin, and assume that the nominal wavefronts are perfectly circular and the field amplitude across each wavefront is uniform. (In reality, for precision optics where MSF is critical, the nominal aberrations are sufficiently small that the error estimates found in what follows are expected to remain useful.) The equation for the nominal rays is then simply *x* = *ξ z*/*z*_{M}. The initial conditions for the nominal phase and amplitude are

_{0}and Φ

_{1}follow from Eqs. (7)-(9):

*z*= 0 (although they break in the vicinity of this plane).

The perturbation model amounts to multiplying the nominal field estimate at the exit pupil plane by the exponential of the term at the end of Eq. (25b). The estimate for the error involved in this approximation, given in Eq. (18), simplifies in this case because the ray caustic reduces to a point and it follows from Eq. (11) that *ζ* is independent of *ξ* for a nominally perfect converging wave:

*r*1 of the first Fresnel zone [see illustration in Fig. 2(b)], defined as

Equation (18) then becomes

*ϕ*

^{″2}»

*ϕ*

^{′4}for small MSF, as discussed after Eq. (22).

For moderate numerical apertures, we can drop the obliquity factor of $\sqrt{{\u3008{\chi}^{-6}\u3009}_{A}}$ in Eq. (28) since it is of the order of unity (e.g. this factor makes a change of the order of 3% and 10% when the NA reaches 0.25 and 0.4, respectively). With this, is useful to rewrite Eq. (28) in the form

*z*

_{M}, defined as

The rule of thumb in Eq. (29) is the main result of this work. It states that the error of the perturbation model scales with (*r*_{1}/ℛ)^{2}, and that, as discussed after Eq. (19), once this factor is comparable to unity (say, greater than 1/3), the perturbation model loses its usefulness. While *r*1 is a familiar and intuitive entrant here, the explicit definition and role of the characteristic length of the MSF, namely ℛ, is our key result. Note that for exit-telecentric systems where propagation across large distances may be expected to be problematic we simply have *z*_{P} →∞ and therefore ${r}_{1}^{2}\to \lambda \left|{z}_{\mathrm{M}}\right|$. Similarly, the limit of large *z*_{M} can be handled by dividing both ℛ and *r*_{1} by *z*_{M} to convert those entities to the tangents of angles subtended at the origin (see Fig. 2) while effectively leaving Eq. (29) intact.

To illustrate these ideas, consider the simple sinusoidal MSF phase structure given in Eq. (21) with *κ* = *C*/*L*_{M}, where *L*_{M} ≜ *L*|*z*_{M}/*z*_{P}| is the size of the beam footprint in the part’s conjugate plane and *C* is the number of MSF cycles across this footprint. Because *A* varies significantly more slowly than *ϕ* and its derivatives, its value within the integral of Eq. (16) can be approximated by a constant, and it follows that ℛ is approximately equal to the MSF’s period divided by $\sqrt{\pi}$, i.e. $\mathcal{R}\approx 1/(\sqrt{\pi}\kappa )$. Figure 3(a) shows a contour plot that illustrates the validity of Eq. (28) for sinusoidal MSF structures, which includes the average of the obliquity factor, namely,

*η*is the numerical aperture in image space. The contours are the values of

*C*such that

*ϵ*=

*φ*/3, and the square of this value is written

It is worth noting that these contours are symmetric about both *z*_{M}/*z*_{P} = 0 and *z*_{M}/*z*_{P} = 1, and that if this plot is wrapped onto a cylinder, the contours are smoothly continuous at the join corresponding to *z*_{M}/*z*_{P} = ±∞. Also note that *ϵ* of Eq. (29) is proportional to *C*^{2} (because ℛ is inversely proportional to *C*) so the perturbation model is useful whenever *C < C*_{max}. It follows that Fig. 3(a) can be interpreted as defining the frequency passband in which the perturbation model is valid as well as providing an accuracy estimate. Note too that modifying the value of *L* or *λ* simply rescales the numbers of the rainbow-like contour legend; the plot itself is unchanged.

As an example, for a system with *η* = 0.4, it follows upon inspection of the associated horizontal dotted gray line in Fig. 3(a) that the perturbation model is valid for frequencies up to 20 cycles except within the interval *z*_{M}/*z*_{P} ∈ [−0.17*,* 0.13] (indicated by the red line segment). That is, the perturbation model is unhelpful when the conjugate plane approaches the image plane. Unsurprisingly, higher frequencies can be resolved when the conjugate plane is nearer the pupil plane. Again with *η* = 0.4 as an example, it is possible to handle more than 160 cycles when *z*_{M}/*z*_{P} ∈[0.9*,* 1.1] (indicated by the blue line segment). Also note that smaller *η* values have larger exclusions. For example, for *η* = 0.05 and frequencies beyond 20 cycles, the perturbation model is invalid for all *z*_{M}/*z*_{P} *<* 0.6 (indicated by the orange line segment). This includes all cases where the conjugate plane is closer to the image than to the exit pupil (hence, for the perturbation model to be of value for this number of cycles, the conjugate plane cannot sit on the far side of the image plane). Figure 3(b) shows the NRMSE as a function of ${r}_{1}^{2}/{\mathcal{R}}^{2}$ for *η* = 0.05, and this is generated by varying *z*_{M}. The estimate in Eq. (29) is seen to be accurate within the region of interest provided the amplitude *h* is sufficiently small. This plot is strikingly similar to that in Fig. 1. A different perspective is presented in Fig. 4 which, for *η* = 0.2, shows NRMSE plots for various values of *h* and *C*. Figures 4(a) and 4(c) are essentially squared cross-sections of Fig. 3(a) at *η* = 0.2, and the regions of validity of Eq. (29) can be seen to correspond to those indicated by Fig. 3(a). The NRMSE plots as a function of ${r}_{1}^{2}/{\mathcal{R}}^{2}$ are shown in Figs. 4(b) and 4(d). Again it is clear that the small *h* approximation gives accurate estimates for *h* ≤ *π*/8. In fact, the maximum deviation between the numerically calculated NRMSE and the simple estimate from Eq. (29) is less than 12% provided that *h* ≤ *π*/4.

At this point it is important to make a distinction between two situations. As mentioned earlier, it is possible for *z*_{M}, the axial position of the image of the optical surface introducing the MSF, to be to the left or right of the exit pupil, *z*_{P}, and for either of these to be to the left or right of the system’s image plane at *z* = 0. However, for an optical system in which there are internal images, this order does not necessarily reflect the order in which the corresponding actual optical elements are encountered in the system. The order in the system between the surface introducing the MSF phase and the aperture stop is important when more rigorously simulating the field at the image plane. Even when working in image space, this order can be used to determine whether diffraction at the exit pupil takes place befor*e* or afte*r* the MSF phase is applied, regardless of the order of z_{P} and *z*_{M}. A theoretical argument of why the RMSE estimate is valid for these two situations is provided in Appendix D. The validity of the RMSE for both orderings was also verified numerically.

It is possible to gain added intuition about Eq. (28) upon expressing the MSF phase in terms of some type of frequency spectrum. Although orthogonal polynomials also represent an attractive option [4], the essential result is perhaps simplest to appreciate by contemplating a band-limited Fourier series containing only MSF frequencies of interest, say

where*κ*=

_{m}*m*/

*L*

_{M}. If the field amplitude is roughly constant, it now follows that $\mathcal{R}=\frac{1}{\sqrt{\pi}}{\left(\overline{{\kappa}^{4}}\right)}^{-1/4}$, where $\overline{{\kappa}^{4}}$ is a normalized fourth moment of the spatial frequency spectrum, i.e.

Note that, while the spectrum discussed here uses the exit pupil as its domain, i.e. the projected beam footprint rather than the entire aperture of each surface of interest, the entire apertures are included within the analysis upon averaging over the field to extract an overall measure of the MSF’s impact.

Figure 5 presents plots of the NRMSE for five nonsinusoidal MSF phase structures, shown in the inset in Fig. 5(a), with equal RMS value and whose power spectral densities decay according to a power law. Again, these error plots are generated by varying *z*_{M} while keeping constant the projection of the MSF phase structure (in keeping with the perturbation model) across the exit pupil. The wavelength, exit pupil size and NA are the same as in the previous example. Notably, the error estimate in Eq. (29) approximates well all these errors within the regime in which the perturbation model is useful, namely when ${r}_{1}^{2}/{\mathcal{R}}^{2}$ is less than about a third. (Appendix D provides a description of this error for larger values of ${r}_{1}^{2}/{\mathcal{R}}^{2}$).

Finally, we discuss the generalization of these results to the case of *M* optical surfaces, each introducing an MSF phase structure *ϕ _{m}*, for

*m*= 1

*,*2

*,…, M*. Let the conjugates of these surfaces in image space sit at

*z*=

*z*. The perturbation approximation then involves multiplying the nominal field at the exit pupil by $\mathrm{exp}\left[\mathrm{i}{\displaystyle {\sum}_{m=1}^{M}{\varphi}_{m}(x{z}_{m}/{z}_{\mathrm{P}})}\right]$. The generalization of the dominant part of Eq. (14) results from simply introducing a summation over these contributions:

_{m}*ζ*= (

_{m}*z*

_{P}−

*z*)

_{m}*z*/

_{m}*z*

_{P}. By assuming that the MSF phases of the different surfaces are statistically uncorrelated, the total relative RMSE,

*ϵ*

_{T}, can be estimated as the sum of the squares of the RMSEs for each surface as

Note that *ϵ*(*z*_{P}*, z _{m}*;

*ϕ*)

_{m}*< φ*/3 should be checked individually, where

_{m}*φ*is given by Eq. (19) for the corresponding

_{m}*ϕ*.

_{m}## 5. Concluding remarks

By investigating the propagation in a homogeneous medium of a general nominal wavefield that is perturbed by an MSF phase structure, we have constructed an asymptotic framework that can be used to assess the validity and accuracy of the perturbation model. For this framework to be consistent with the perturbation model, the effects of the MSF structure are placed one asymptotic order lower than the nominal wavefront. This leads to a solution that employs the nominal rays but modifies the phase and amplitude of their contributions. This framework allows not only the derivation of the perturbation model itself but also of a series of corrections, the leading term of which is used to find a simple RMS error estimate. Notice that, while not the goal of the current work, these corrections can be used not only to estimate the error of the perturbation model but to actually provide a better approximation of the field if desired. Importantly, our results give explicit estimates of both when the perturbation model retains validity (i.e. leads to errors that are less than some agreed fraction, say a third, of the impact caused by the MSF structure) as well as of the model’s accuracy within its domain of validity.

As a specific demonstration of these ideas we showed that, under an additional assumption, significant progress can be made for optical imaging systems. Our process for considering the impact of MSF on a particular interface does not require knowledge of all the system details, but involves only three locations of interest in image space: the system’s image plane and exit pupil plane as well as the plane that is conjugate to the interface itself. Just as the imagery from stop to pupil is generally taken for granted, we assumed that the frequency band of interest in the interface’s MSF is resolved at its associated conjugate plane with higher fidelity than in the free-space propagation to the exit pupil. While it would be straightforward to investigate this assumption for a variety of standard system types (perhaps using nothing more than resolution estimates based on geometrical optics), that falls outside the scope of this more general contribution. Our asymptotic treatment led to the intuitive rule of thumb in Eq. (29) for the RMS error caused by the perturbation model in terms of the size of the first Fresnel zone when propagating to the exit pupil from the image of the surface introducing the MSF phase structure: the RMS error is simply proportional to the inherent impact of the MSF phase, times a scaling factor given by the square of the ratio between the first Fresnel zone’s half-width and a specific characteristic length of the MSF phase. (This characteristic length can be expressed either in terms of an average second derivative or a fourth spectral moment.) That is, whenever the central Fresnel zone is sufficiently smaller than this characteristic length, the perturbation model remains valid. Finally, these relations are extended to multiple surfaces and are shown to be equally valid regardless of the ordering between the aperture stop and the surface introducing the MSF structure, as well as of the order of their images in image space.

## A. MSF-independent rays derivation

The method of characteristics leads to the solution to Eq. (4) in terms of the parametrization in Eq. (6). When changing variables according to Eq. (6), it is convenient to introduce the transpose of the Jacobian matrix:

With this, derivatives in the Cartesian coordinates **r** = (*x, z*) can be written in terms of derivatives in the ray parameters ** ξ** = (

*ξ, s*) according to the chain rule,

*= (*

_{ξ}*∂*/

*∂ξ, ∂*/

*∂s*).

We start by showing that Eq. (7) is a solution to the Eikonal equation in Eq. (4). By using Eq. (38), the parametrized gradient of Eq. (7) gives, after some simplification,

It is then straightforward to show that the dot product of Eq. (39) with itself gives unity, hence satisfying Eq. (4). Notice also that Eq. (39) combined with Eq. (38) allow turning the left-hand side of Eq. (5) into a simple derivative in *s*:

For *N* = 1, Eq. (5) reduces to

By changing variables to ** ξ** and using Eq. (40), the left-hand side of Eq. (41) transforms into $\partial \overline{{\Phi}_{1}}/\partial s$. The corresponding parametrization of the right-hand side gives

The logarithmic portion of Eq. (44) leads to an overall amplitude factor that accounts for the bunching of the rays under propagation. This factor diverges at the caustics of these nominal rays.

The *N* = 2 case of Eq. (5) can be written as

By parametrizing in terms of ** ξ** and using Eq. (40) to simplify the left-hand side we get

*a*is the contribution that is independent of the MSF perturbation

*ϕ*, and is given by

Since this contribution is unrelated to the MSF structure, we will not write explicitly its integral in *s* and will just refer to the resulting expression as $\overline{{\Omega}_{2}}$. To integrate the rest, we use the simple results

Because derivatives in *ξ* and integration in *s* commute, we can integrate Eq. (47) to find

*ζ*=

*χ*

^{2}

*s*/Δ = (

*z*−

*z*

_{M})/[1 + (

*z*−

*z*

_{M})

*W*

^{″}(

*ξ*)/

*χ*

^{3}(

*ξ*)]. Notice that a factor of 1 = Δ/

*χ*−

*sW*

^{″}/

*χ*

^{2}was introduced in the third line, which allowed simplifying the expression.

## B. Invariance of the error estimate in conjugate spaces

For any optical system (or subsystem) with non-zero focal power, the angle characteristic [13] can be written to second order as

*f*is an effective focal length, primes distinguish entities in exit space,

*n*is the refractive index,

*p*is an optical direction cosine (i.e. refractive index times transverse direction cosine), and

*δ*is a displacement from the focal point. The front and rear focal lengths are

*nf*and

*n*

^{′}

*f*, respectively. From the basic equations of Hamiltonian optics [13], namely

*y*= −

*∂T*/

*∂p*and

*y*

^{′}=

*∂T*/

*∂p*

^{′}, it follows that the displaced reference planes are conjugate when

*δδ′*= −

*nn′*

*f*

^{2}, and the magnification is then

*δ′*/(

*n′*

*f*). Now, suppose we know three locations in exit space (or conjugates to these locations), namely ${\delta}_{\mathrm{i}}^{\prime},{\delta}_{\mathrm{m}}^{\prime}$

*,*and ${\delta}_{\mathrm{p}}^{\prime}$ for the image, MSF-bearing, and pupil planes respectively. The unprimed entities are their respective conjugates in entrance space with

*δ*= −

*nn*

^{′}

*f*

^{2}/

*δ*

^{′}. Our estimate of error in the perturbation model involves the ratio of a characteristic length in the MSF-bearing plane, say ${\mathcal{L}}_{\mathrm{m}}^{\prime}$, and the radius of the first Fresnel zone when propagating a nominally spherical wave centered at a point ${\delta}_{\mathrm{i}}^{\prime}$ from ${\delta}_{\mathrm{m}}^{\prime}$ to ${\delta}_{\mathrm{p}}^{\prime}$. The square of this ratio is given by

This can be expressed in terms of the locations in entrance space by using *δ*^{′} = −*nn*^{′}*f* ^{2}/*δ* to find

## C. Perturbation model in angular spectrum domain

The asymptotic series in Eq. (3) provides a relatively simple method to approximate solutions to the Helmholtz equation. By using these solutions, we obtained RMSE results that give insight into the validity of the perturbation model. However, these results are limited because of their failure at nominal caustics. Here we present an alternative approach, based on the angular spectrum and two applications of the stationary phase method. This approach leads to the same algebraic results for the field approximations, but includes an unambiguous Maslov-Gouy phase shift factor associated with the passage through caustics and therefore gives a result of extended validity.

A monochromatic 2D field, *U*, is related to its angular spectrum, $\tilde{U}$, by

*U*(

*ξ, z*

_{M}) by

For the specific case of a converging wave that is perturbed by a MSF structure, *U*(*ξ, z*_{M}) = *U*_{0} *A*(*ξ*) exp[i*ϕ*(*ξ*)] exp[i*kW*(*ξ*)], where *A* and *W* are given in Eq. (24). Like the process in Sec. 2, the stationary phase approximation relies on the assumption that *k* is a “large” asymptotic parameter, i.e. that *λ* is much smaller than any characteristic length of the field at the plane in question. To perform this approximation, Eq. (55) is written as

*a*(

*ξ*) =

*U*

_{0}

*A*(

*ξ*) exp[i

*ϕ*(

*ξ*)] is a slowly-varying complex amplitude,

*kγ*(

*ξ*) =

*k*[

*W*(

*ξ*)−

*αξ*]is the phase associated with the rapid oscillations, and in the second line we changed variables to

*ξ*=

*ξ*

_{0}+

*τ*, where

*ξ*

_{0}is the stationary point of

*γ*, i.e.

*γ*

^{′}(

*ξ*

_{0}) = 0. (For the focused wave we are considering, there is only one such stationary point.) The factor in large parentheses within the second line of Eq. (56) is expanded about

*τ*= 0. Each of the resulting terms can be integrated analytically by using

*β*=

*γ″*(

*x*0). Hence, from Eq. (24), sgn(

*β*) = sgn(

*z*

_{M}). For the current purposes it is sufficient to keep only the two leading terms, proportional to

*k*

^{−1/2}and

*k*

^{−3/4}, respectively. The resulting approximate expression is then substituted into Eq. (54), and the integral is again evaluated by using the method of stationary phase, keeping the two leading terms according to the powers of

*k*they contain. After a long but straightforward calculation, it is found that

*U*(

*x, z*) ≈

_{P}*U*

_{0}exp[i

*k*Φ

_{0}(

*x, z*

_{P}) + Φ

_{1}(

*x, z*

_{P})], with Φ

_{0}and Φ

_{1}given in Eqs. (25a) and (25b), which includes the correct Maslov-Gouy phase terms.

## D. Equivalence of error estimate regardless of order between stop and MSF

We now give a justification of why the RMSE derived in the main body is valid whether the surface introducing the MSF is before or after the aperture stop. For brevity we use operator/Dirac notation. The operators used in this proof are: ${\widehat{K}}_{z}$, which denotes free propagation by a distance *z*; ${\widehat{M}}_{\mathrm{M}}$, which denotes the MSF phase factor at *z*_{M}, the image plane for the surface in question; ${\widehat{M}}_{\mathrm{P}}$, which denotes the MSF written at the pupil plane *z*_{P} according to the perturbation approximation; and ${\widehat{a}}_{\mathrm{P}}$, which denotes the transmission function of the aperture image at the exit pupil, normalized by the extent of the pupil. The ideal field at any plane *z* focused towards an image point *x*′ (ignoring the aperture or MSF) can be approximated, to within an unimportant constant amplitude, as $U(x,z;{x}^{\prime})=\u3008x\left|{\widehat{K}}_{z}\right|{x}^{\prime}\u3009$. This field is a converging wave if *z <* 0 and a diverging one if *z >* 0. Note that $\u3008x\left|{\widehat{K}}_{z}\right|{x}^{\prime}\u3009$ also corresponds to the Rayleigh-Sommerfeld propagation kernel.

Let us study first the simpler case of a field focused at the image point *x*′ that aquires a MSF phase at *z*_{M} and then travels to the exit pupil, where it is diffracted. After propagating to the image plane *z* = 0 this field given by

Notice that since both ${\widehat{M}}_{\mathrm{P}}$ and ${\widehat{a}}_{\mathrm{P}}$ are multiplicative, they commute. The square of the RMSE for this first case is then

Now consider the case in which the field is first diffracted by the stop and then acquires a MSF phase. After propagating to the image plane this field is given by

It is worth noting that Eqs. (60) and (63) have very different forms, the latter having the two aperture operators separated. The symmetry between these two expressions can be restored by averaging over all object points, that is, by integrating in *x*^{′} over a sufficiently large range, replacing the *x*^{′}-dependent normalization factor $\mathcal{N}$ by $\mathcal{M}={\displaystyle \iint |U(x,0;{x}^{\prime}){|}^{2}\mathrm{d}x\mathrm{d}{x}^{\prime}}$. Assuming that the range of integration is sufficiently large to approximate $\int |{x}^{\prime}\u3009\u3008{x}^{\prime}|\mathrm{d}{x}^{\prime}}\approx \widehat{1$, we arrive at

*x′*, we then can expect that the asymptotic behavior of the error is the same independently of the order within the system between the aperture stop and the surface introducing the MSF.

To appreciate why the order of diffraction by the stop and MSF does not affect significantly the error estimate, even for a single object/image point, consider a simple case that allows a closed-form approximate expression: a sinusoidal MSF phase of the form in Eq. (21) for *h* ≪ *π*/2, and under the paraxial approximation. It is convenient to evaluate the errors at the exit pupil plane, *z*_{P}. Let us begin by considering some nominal initial field *U*_{i}(*x, z*_{P}), and propagating it to the plane conjugate to the MSF structure according to the Fresnel propagation formula:

The effect of the MSF is accounted for by multiplying this field by the MSF phase structure, *U*_{ii}(*ξ, z*_{M}) = *U*_{i} (*ξ, z*_{M}) exp[i*ϕ*(*ξ*)], where we use the approximation exp[i*ϕ*(*ξ*)] ≈ 1 + i*ϕ*(*ξ*) = 1 + i*h* sin(2*πκξ*). We can now back-propagate this field to the pupil plane through an (inverse). Fresnel propagation integral:

*x*≜

_{κ}*κλ*(

*z*

_{M}−

*z*

_{P}). For the case in which the MSF phase is acquired before diffraction at the stop, the initial field is a perfect paraxial converging wave,

*U*

_{i}(

*x, z*

_{P}) =

*U*

_{0}exp[i

*πx*

^{2}/(

*λz*

_{P})], prior by rect and the back-propagated field

*U*

_{ii}(

*x, z*

_{P}) must be multiplied a pupil function rect(

*x*/

*L*) prior to comparison with the perturbation model. On the other hand, for the case in which the MSF phase follows diffraction by the stop, the initial field already includes the aperture function, i.e.,

*U*

_{i}(

*x, z*

_{P}) =

*U*

_{0}exp[i

*πx*

^{2}/(

*λz*

_{P})] rect(

*x L*). In both cases, the resulting back-propagated fields at the pupil must be compared with the perturbation model given by

The subtraction of *U*_{P} from either back-propagated field leads to the cancellation of the leading terms, such that the resulting difference is proportional to *h*. The integral of the modulus squared of this difference, normalized by *L U*_{0}|^{2}, gives the square of the RMSE. After simplification, this normalized integral can be written some

*j*= I,II (representing the two cases, namely the MSF phase being acquired before or after diffraction at the stop), with r

_{±}

_{,}_{I}(

*x*) ≜ rect(

*x*/

*L*) and r

_{±}

_{,}_{II}(

*x*) ≜ rect[(

*x*±

*x*)/

_{κ}*L*] and where we used the fact that

*κ |z*

_{M}/

*z*

_{P}| =

*C*/

*L.*That is, the only difference between both cases is a slight shift by ±

*x*of the regions in which some of the terms contribute. For

_{κ}*L |x*, this shift has little effect on the result. In either case, the integral can be carried out analytically, leading to

_{κ}|*ϕ*

^{2}≈

*h*

^{2}/2. Notice that these expressions depend on only two dimensionless parameters: the number of cycles

*C*, which is significantly greater than unity, and the product ${\kappa}^{2}{r}_{1}^{2}$, which is smaller than unity for cases in which

*ϵ*

_{I}

_{,}_{II}

*< φ*according to Eq. (29). It can be easily shown that the relative difference between these two quantities is small as long as the following condition is satisfied:

Note that in deriving Eqs. (70) and (71) we assumed that the sinusoidal MSF in Eq. (21) is exactly antisymmetric around the axis. If this sinusoidal were shifted by some distance *ξ*_{0}, Eqs. (70) and (71) would change. This change is particularly simple for Eq. (71), where the only modification is that the contribution proportional to sinc (2*πC*) woud be multiplied by a factor of cos(4*πκξ*_{0}). When considering the level of error of a generic sinusoidal MSF structure, perhaps what is more meaningful is the average of the error over all possible shifts *ξ*_{0} of the sinusoidals, especially because this has a similar effect to the averaging over different object points mentioned earlier. These averages for both ${\u03f5}_{\mathrm{I}}^{2}$ and ${\u03f5}_{\text{II}}^{2}$ give the simpler (and more representative) results

Note that, again, the relative difference between these two errors is negligible if Eq. (72) is met.

Finally, consider MSF structure containing many spatial frequencies. For simplicity, we consider explicitly only the case where diffraction at the stop follows the acquisition of the MSF phase. We again use the approximation exp(i*ϕ*) ≈ 1 + i*ϕ*, with *ϕ* now being the multi-frequency phase in Eq. (33) instead of the sinusoidal in Eq. (21). Due to the linear independence of the different frequencies, the result is just a weighted sum of contributions:

## Funding

National Science Foundation (NSF) I/UCRC Center for Freeform Optics (IIP-1338877).

MAA received funding from the Excellence Initiative of Aix-Marseille University - A^{∗}MIDEX, a French “Investissements d’Avenir” programme.

## Acknowledgments

We thank the members of CeFO for their support and insights, and acknowledge in particular useful discussions with James Fienup.

## References

**1. **R. J. Noll, “Effect of mid- and high-spatial frequencies on optical performance,” Opt. Eng. **18**(2), 182137 (1979). [CrossRef]

**2. **D. Aikens, J. E. DeGroote, and R. N. Youngworth, “Specification and control of mid-spatial frequency wavefront errors in optical systems,” in *Frontiers in Optics 2008/Laser Science XXIV/Plasmonics and Metamaterials/Optical Fabrication and Testing, OSA Technical Digest (CD)* (Optical Society of America, 2008), paper OTuA1. [CrossRef]

**3. **J. M. Tamkin, T. D. Milster, and W. Dallas, “Theory of modulation transfer function artifacts due to mid-spatial-frequency errors and its application to optical tolerancing,” Appl. Opt. **49**, 4825–4835 (2010).

**4. **G. W. Forbes, “Never-ending struggles with mid-spatial frequencies,” Optical Measurement Systems for Industrial Inspection IX, Proc. of SPIE **9525**, 95251B (2015).

**5. **R. N. Youngsworth and B. D. Stone, “Simple estimates for the effects of mid-spatial-frequency surface errors on image quality,” Appl. Opt. **39**, 2198–2209 (2000). [CrossRef]

**6. **P. A. Sturrock, “Perturbation characteristic functions and their application to electron optics,” The Royal Society. **210–1101**, 269–289 (1951).

**7. **H. A. Buchdahl, “Perturbed characteristic functions, II second-order perturbation,” International Journal of Theoretical Physics. **24**, 457–465 (1985). [CrossRef]

**8. **E. Garbusi and W. Osten, “Perturbation method in optics: application to the interferometric measurement of surfaces,” J. Opt. Soc. Am. A. **26**2538–2549 (2009). [CrossRef]

**9. **H. H. Hopkins and H. J. Tiziani, “A theoretical and experimental study of lens centring errors and their influence on optical image quality,” Br. J. Appl. Phys. **17**33–55 (1966). [CrossRef]

**10. **H. A. Buchdahl, “Perturbed characteristic functions: application to the linearized gravitational field,” Aust. J. Phys. **32**405–410 (1979). [CrossRef]

**11. **B. D. Stone, “Perturbations of optical systems,” J. Opt. Soc. Am. A. **14**2837–2849 (1997). [CrossRef]

**12. **M. Rimmer, “Analysis of perturbed lens systems,” Appl. Opt. **9**533–537 (1970). [CrossRef] [PubMed]

**13. **G.W. Forbes and B.D. Stone, “Restricted characteristic functions for general optical configurations,” J. Opt. Soc. Am. A **10**, 1263–1269, (1993). [CrossRef]