Validity of the perturbation model for the propagation of MSF structure in 2D

: 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.


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 where 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[ikW(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[ikΦ(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 (ik) −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: x(ξ, s) = ξ + sW (ξ), z(ξ, s) = z M + s χ(ξ), 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 χ(ξ) 1 − W (ξ) 2 , with denoting a 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: 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 where ∆(ξ, 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 where Ω 2 is independent of φ, 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 with ∂ x/∂ξ being calculated for fixed z.
In summary, the perturbation model can be expressed as follows: where the first factor on the right-hand side of the first line is simply the underlying field in the absence of MSF phase structure. That is, the only φ-dependent component in U P is that contained in Φ 1 . The leading correction to this approximation results from also considering the φ-dependent terms in Φ 2 given in Eq. (10):

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 where a is the aperture in the initial reference plane. In the second line of Eq. (15) we used Eqs. (6), (9), and (13), along with the fact that Ω 2 is real-valued, to see that Note that the integral in the denominator is independent of z, due to the fact that, when evanescent waves are not included, free propagation is a unitary operation. Upon defining the weighted average and expanding the exponential in Eq. (15), one can write a simple approximation for 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 where we used the field's initial conditions and the linearity of the Helmholtz equation as well as the unitarity of free propagation to see that ϕ is independent of z, and in the last step we assumed |φ| π. Clearly, if is 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 is 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 where in the last step we assumed 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, 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 M |/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 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 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 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.

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 from which the solutions for Φ 0 and Φ 1 follow from Eqs. (7)- (9): where sgn(·) is the signum function. Note that, following the considerations in Appendix C which dictate the correct choice for the Gouy-Maslov phase shift, these expressions are asymptotically valid at either side of the image plane 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 The image space of an imaging system where the origin is placed at the image plane (blue). The image of the MSF phase (green) and the exit pupil (red) are located at z M and z P , respectively. Here, z M < z P < 0 but either or both could be positive. (a) We assume in this work that the lion's share of the error in the perturbation model can be accounted for by just the step of carrying the MSF phase structure along the nominal rays (gray lines) from z M to z P . Note that ξ is the value of the x intercept at z M . The size of the exit pupil is L, and that of the beam footprint at z M is L M . (b) Definition of r 1 as the half-width of the first Fresnel zone at z M , that is, as the height at which two circles, centered at the center of the exit pupil and the image point, respectively, and touching the axis at z M , have a separation of λ/2 in the Fresnel approximation.
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: where this quantity is expressed in terms of the radius r 1 of the first Fresnel zone [see illustration in Fig. 2(b)], defined as Equation (18) then becomes where in the last step we used the fact that φ 2 φ 4 for small MSF, as discussed after Eq. (22). For moderate numerical apertures, we can drop the obliquity factor of χ −6 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 where R is a measure of the characteristic feature size of the MSF structure at 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 /R) 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 R, 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 2 1 → λ|z M |. Similarly, the limit of large z M can be handled by dividing both R 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 R is approximately equal to the MSF's period divided by √ π, i.e. R ≈ 1/( √ πκ). 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, where η 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 R 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 2 1 /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 2 1 /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 before or after 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 R = 1 √ π κ 4 −1/4 , where κ 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 2 1 /R 2 is less than about a third. (Appendix D provides a description of this error for larger values of r 2 1 /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 m . The perturbation approximation then involves multiplying the nominal field at the exit pupil by exp i M m=1 φ m (xz m /z P ) . The generalization of the dominant part of Eq. (14) results from simply introducing a summation over these contributions: where, in keeping with Eq. (26), ζ 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 Fig. 5. The normalized RMSE plot, as a function of z M /z P , is shown in (a) for randomly generated MSF structures that possess a power-decay spectra, whose extent over the aperture is L = 2 cm. The solid curves are generated from Eq. (29), and the numerically calculated values are shown as dots whose color corresponds to the corresponding structure, shown in the subplot. Part (b) shows the same information but plotted against r 2 1 /R 2 , where the single black line is the estimate. The green region indicates ≤ ϕ/3.
Note that (z P , z m ; φ m ) < ϕ m /3 should be checked individually, where ϕ m is given by Eq. (19) for the corresponding φ m .

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, where ∇ ξ = (∂/∂ξ, ∂/∂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 ∂Φ 1 /∂s. The corresponding parametrization of the right-hand side gives with ∆ = Det(J) = ( χ 2 + sW )/ χ. With these results, Eq. (41) now becomes which has a simple solution that satisfies Φ 1 (ξ, 0) = ln[A(ξ)] + iφ(ξ): 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 where in the last part we use the convention of implicit summation over repeated subindices. By using the product rule and some simplification, this expression can be separated into three types of contribution, depending on how many of the two derivatives act on the exponential: where 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 Ω 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 where 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 where 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 n f 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 δ i , δ m , and δ p 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 L m , and the radius of the first Fresnel zone when propagating a nominally spherical wave centered at a point δ i from δ m to δ p . 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 n λ where the inverse of the magnification between the MSF-bearing planes, namely n f /δ m = δ m /(n f ), was used to convert the characteristic length. The fact that the final expression of Eq. (53) matches the original expression on the left-hand side, except it now has no primes, means that our ratio is invariant under propagation through any part, or all of, an optical system. That is, for example, we'll get the same estimate of the error in accounting for the MSF on a given surface whether we first propagate it to the stop and image that result to the exit pupil or instead begin by imaging it to its conjugate plane in image space and then propagate from there to the exit pupil.

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,Ũ, by whereŨ is related to the initial field U(ξ, z M ) bỹ For the specific case of a converging wave that is perturbed by a MSF structure, , 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 with β = γ (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[ikΦ 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:K z , which denotes free propagation by a distance z; M M , which denotes the MSF phase factor at z M , the image plane for the surface in question;M P , which denotes the MSF written at the pupil plane z P according to the perturbation approximation; andâ 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 ) = x|K z |x . This field is a converging wave if z < 0 and a diverging one if z > 0. Note that x|K z |x 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 The effect of the MSF is accounted for by multiplying this field by the MSF phase structure, , where we use the approximation exp[iφ(ξ)] ≈ 1 + iφ(ξ) = 1 + ih sin(2πκξ). We can now back-propagate this field to the pupil plane through an (inverse) Fresnel propagation integral: where 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 )], and the back-propagated field U ii (x, z P ) must be multiplied by 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 some simplification, this normalized integral can be written as for 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 I ϕ 2 ≈ 4 [1 − sinc (2πC)] sin 2 πκ 2 r 2 where we used ϕ 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 κ 2 r 2 1 , 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 2 I and 2 II give the simpler (and more representative) results 2 I /ϕ 2 ≈ 4 sin 2 πκ 2 r 2 1 /2 , 2 II /ϕ 2 ≈ 4 sin 2 πκ 2 r 2 1 /2 + 2 κ 2 r 2 1 C cos πκ 2 r 2 1 .
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.
where we used ϕ 2 ≈ m |a m | 2 . If κ 2 m r 2 1 1 for all meaningful MSF spatial frequencies, the sine function can be approximated by its argument, leading to the result in Eq. (29). The departure of the plots in Fig. 5 from this linear approximation is explained by the nonlinearity of Eq. (75).

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.