Validity of the perturbation model for the propagation of MSF structures in 3D.

Mid-spatial frequency (MSF) structures on optical surfaces degrade system performance and a perturbation model is typically used to simplify the assessment of their effects. In this simple model, MSF phase structures are dragged along the nominal rays of a system to yield estimates of wavefronts in the exit pupil that may be used for further analysis. However, the validity of the perturbation model remains an open area of study. We extend our previous assessment of the validity of this model [K. Liang, Opt. Express 27, 3390-3408 (2019)] that was focused on the analysis of single-frequency MSF structures in two dimensions to now include error estimates for broad-spectra MSF structures in three dimensions.


Introduction
Mid-spatial frequency (MSF) structures are inevitable in most aspheric and freeform optical systems due to the subaperture tools that are used during the manufacturing process. Their characteristic frequencies lie between those of the common low-order aberrations and high-order scattering, and their detrimental effects on optical performance remain an active area of research. For example, there have been many efforts towards simplifying the tolerancing of optical parts afflicted with MSF [1][2][3][4][5][6]. To this end, the perturbation model is often used to cut down on the computation time needed to understand the propagation of MSF structures. This model, in which the MSF phase structure (which can vary significantly from part to part) is simply dragged along rays of the nominal system, is often used in order to avoid the need for new ray tracing for each MSF realization [7]. However, the validity of this perturbation model requires further treatment; its analysis in two dimensions was presented in Ref. 8 and those results are extended in Appendix A in a manner that we now generalize to three dimensions.
The mathematical framework used in this manuscript to estimate the error incurred by the perturbation model is based on an asymptotic analysis of the Helmholtz wave equation for the propagation of a monochromatic field in free space. A key step is the placement of the MSF structure at one asymptotic order beneath the nominal wavefront. This framework is similar to that in Ref. 8, but the solution now includes contributions from every asymptotic order (under appropriate approximations). This extension permits the analysis of MSF structures with broad spatial-frequency spectra.

Asymptotic propagation estimate based on nominal rays
The propagation of a monochromatic scalar field, Re U(r)e −iωt , in a homogeneous medium is governed by the Helmholtz equation where k = ω/c = 2π/λ is the wavenumber in the medium. The field is taken to be propagating towards larger z and, at some reference plane z = z M , we take the initial value of the field to be given nominally by U(r ⊥ , z M ) = U 0 A(r ⊥ ) exp[ikW(r ⊥ )], where U 0 is a constant with field units and r ⊥ = (x, y) are the transverse coordinates. At z = z M , we also superpose an MSF phase factor of the form exp[iφ(r ⊥ )], where φ(r ⊥ ) is taken to have zero mean and the magnitude of its variation is less than π. Moreover, given its characterization as an MSF structure, φ(r ⊥ ) is assumed to vary more rapidly than either W(r ⊥ ) or A(r ⊥ ). The goal now is to derive an estimate, in a manner that is similar to that in Ref. 8, of how this MSF phase structure affects the field under propagation. We begin by writing U(r) = U 0 exp[ikΦ(r)], where Φ(r) is a complex quantity that accounts for spatial variations in both the phase and amplitude. With this, Eq. (1) becomes Eq.
(2) can be solved upon expressing Φ as an asymptotic series in the parameter (ik) −1 : By using Eq. (3) with Eq. (2) and separating terms of equal powers of k, we arrive at The initial conditions discussed above can now be stated as Φ 0 (r ⊥ , z M ) = W(r ⊥ ), Φ 1 (r ⊥ , z M ) = ln[A(r ⊥ )] + iφ(r ⊥ ), and Φ N (r ⊥ , 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. The image space of an imaging system, where the image plane (blue) is placed at z = 0. The exit pupil (red) and the image of the MSF phase structure (green) are located at z P and z M , respectively. Note that ξ ⊥ is the location of the intersection of the ray starting at r ⊥ in the plane z = z M . We begin with Eq. (4), the well-known Hamilton-Jacobi or Eikonal equation, which can be solved in terms of nominal rays by using the following parametrization involving ξ = (ξ, η, s): where ξ ⊥ = (ξ, η) are the transverse coordinates 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 denoting a definition and with ∇ ξ ⊥ (∂ ξ , ∂ η ). Figure 1 shows the relation between r ⊥ and ξ ⊥ for a nominally converging wavefront. It is shown in Appendix B that the solution to the Eikonal equation is simply where an overline on any function f (r) indicates that it is being expressed in terms of the ray parameters ξ by using Eq.
Notice that, as a consequence of placing φ one asymptotic order below W, the rays used in this analysis are the nominal rays, which are not affected by the MSF structure. It is furthermore shown in Appendix B that Eq. (5) can also be solved in terms of the parametrization in Eq. (6). For N = 1, the parametrized solution is where is the Jacobian determinant of the coordinate transformation given in Eq. (6), H W is the Hessian matrix of W, and θ GM is the Gouy-Maslov phase shift. This phase shift is a straightforward extension to three dimensions of that discussed in Appendix C of Ref. 8. As also discussed in Ref. 8, the first term of Eq. (8) accounts for the change in the amplitude due to the bunching or spreading of the nominal 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 (nominal) rays. This is precisely the perturbation model. In what follows, the method used for proceeding to larger values of N differs from that presented in Ref. 8. In order to appreciate the three-dimensional results, however, it is helpful to revisit the two-dimensional case and present the mathematical framework upon which the full three-dimensional treatment will follow by analogy, see Appendix A. As a reminder of the derivation in Ref. 8, recall that only the first correction to the perturbation model was analyzed. That is, the series in Eq. (3) is truncated at N = 2, hence the field is taken to be approximated by By using the field estimate in Eq. (10), the resulting rules of thumb for the error the perturbation model were ultimately found to be related to the fourth spectral moment of the MSF structure φ. This result inspired the development of a new family of rapidly decaying Fourier-like basis functions that yield finite spectral moments [9]. However, an alternative method can be used so that the field estimate contains contributions from every term on the right-hand side of Eq. (3). Although more approximations are necessary for this route, they are consistent with the assumptions regarding φ and are detailed in the full derivation shown in Appendix A. As a result, a better error estimate (not limited to N ≤ 2) for the perturbation model is obtained.
It is convenient, and necessary with regards to the derivations in Appendices A and B, to now work in image space in cases where, to a good approximation, a wavefront propagating from a point object source converges onto a point on the image plane. As in Ref. 8, we consider for simplicity only the on-axis object point, whose ideal image is located at the origin. It should be noted, however, that the analysis that follows can be used for off-axis object points as well since the choice of the origin was made out of convenience and similar methods can be applied to off-axis field points. Furthermore, we assume that the MSF content on each optical surface is adequately resolved in its corresponding conjugate plane in image space. Under these assumptions, the dominant error of the perturbation model is associated with the process of simply dragging these MSF structures along the nominally converging rays from their conjugate planes to the exit pupil. Henceforth, we will work with a single optical surface (with MSF structures) whose conjugate plane is located at z = z M . Furthermore, the locations of the exit pupil plane and the image plane are taken to be z = z P and z = 0, respectively, see Fig. 1. With this framework, the nominal (converging) wavefront and obliquity factor are given by To assess the error incurred by the perturbation model, one must go beyond the N = 1 term in Eqs. (3) and (5). It is shown in Appendix B, under the approximations that φ is small and that it varies more rapidly than the nominal quantities, that whereŴ is a differential operator, expressed in plane-polar coordinates (r, θ) where r = ξ 2 + η 2 and θ = arg(ξ + iη). This differential operator is discussed further in Appendix C. The perturbation model is given by We note that the only φ-dependent component in U P entered via Φ 1 . Furthermore, Ω represents contributions from N ≥ 2 that are independent of φ; aside from the stipulation that it is purely imaginary, its explicit form is unimportant here although it is discussed in Appendix B. The (approximate) correction to the perturbation model follows from including the first (φ-dependent) term in Eq. (12): Note that Eq. (15) represents a more complete expression for the field, beyond the perturbation model since it includes contributions from all N. This is in contrast with the equivalent twodimensional expression given in Eq. (14) of Ref. 8, which included only the first three terms (N ≤ 2) in Eq. (3). It should be noted that the approximations used in Appendix B for the purposes of obtaining Eq. (15) have the consequence of limiting our consideration to MSF structures with small amplitudes (|φ|<π), whereas the procedure employed in Ref. 8 explicitly produced a second-order correction term with respect to the amplitude of φ. However, our goal is to provide rule-of-thumb error estimates for residual MSF structures from typical processes in freeform manufacturing; these MSF phase structures generally have amplitudes that are a fraction of the wavelength.

Simple field error estimates in a homogeneous medium
The 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. (15). This is achieved by changing the variable of integration from (x, y) to (ξ, η) by using the differential area transformation where ∂r ⊥ /∂ξ ⊥ is the Jacobian matrix between (x, y) and (ξ, η) after substituting s = (z P − z M )/ χ(ξ ⊥ ) and holding (z P − z M ) constant. The Jacobian determinant is given by The squared RMSE is thus where a is the aperture in the initial reference plane and is the radius of the first Fresnel zone at z M , as illustrated in Fig. 2 . Note that the integral in the denominator of Eq. (19) is independent of z due to the fact that, when evanescent waves are not included, propagation through lossless media is a unitary operation. Upon defining the weighted average and expanding the (outer) exponential in Eq. (18), one can write the following simple approximation for 2 : As discussed in Sec. 4, there are situations when it is sufficient (and in some ways more insightful) to simply consider the first-order (of the argument within the sine function) approximation of Eq. (21). That is, Given thatŴ involves second derivatives, this type of truncation is what led to the aforementioned inspiration to develop the novel bases for the purpose of finding basis sets with finite fourth spectral moments [9]. Upon examination of the complete result in Eq. (21), however, it is clear that such considerations are not necessary.
In Section 4, we use both Eqs. (21) and (22) to obtain rules of thumb regarding the validity of the perturbation model. Although Eqs. (21) and (22) were formally derived for systems with arbitrary NA, the remainder of this work will focus on systems with low to moderate NA. This is because the appropriate analysis for high-NA systems should involve a vector field treatment and the scalar formalism described here is sometimes insufficient. Furthermore, the results for systems with low to moderate NA may be of more interest for manufacturers due to their ease of interpretation and utility. However, for those interested in the behavior of in the high-NA regime, a discussion of those results from this scalar treatment is included in Appendix C.

Rules of thumb for low to moderate NA
In this section, we provide rules of thumb for the error incurred by the perturbation model for an imaging system with low to moderate NA. In this case the expressions for in Eqs. (21) and (22) can be simplified by using χ ≈ 1 and A ≈ 1. Furthermore,Ŵ can be simplified to ∇ 2 ⊥ , the transverse Laplacian operator. With this, Eqs. (22) and (21) become respectively (where the subindices F and C stand for "first term" and "complete", respectively). In the following simulations, with λ = 632 nm, a low-NA system is used for demonstration. As explained in Ref. 8, the resulting rules of thumb are applicable regardless of whether the MSF is located before or after the aperture stop within the system.

Rules of thumb for sinusoidal MSF structures in the milled and turned geometries
For MSF structures whose spectra are well-localized, it turns out that ≈ F in Eq. (23) is sufficient. Note that Eq. (23) can be re-expressed as a normalized RMSE (NRMSE): where is a measure of the characteristic feature size of the MSF phase at z M and ϕ φ 2 1 is the RMS of the MSF structure. The expressions for milled and turned geometries are [in polar coordinates (r, θ)] given respectively by where 2h is the PV, κ is the spatial frequency of the MSF structure, β is a phase offset, and the milled groove pattern is perpendicular to the direction that makes an angle α to the x-axis. Both β and α will turn out to be irrelevant for determining error estimates.
where, as seen in Fig. 1, for a sufficiently large value of κ. That is, turned and milled MSF structures that are wellapproximated with a single frequency obey similar validity measures upon the use of the perturbation model. Note that, in the examples of Eqs. (27) and (28), For completeness and in anticipation of the discussion regarding MSF phases with broad spectra, Fig. 2(b) shows how NRMSE values from a turned MSF surface calculated numerically compare with both F and C of Eqs. (23) and (24), respectively, as functions of z M /z P . Figure 2(b) is an alternative way to view the same data as Fig. 2(a) without having to introduce the notion of r 1 , which may obscure the effects of what happens when, for example, the MSF is placed near the exit pupil or the focus. Near z M /z P = 1, F is accurate and, as is indicated by the previous discussion regarding Fig. 2(a), the perturbation model is valid there since /ϕ<1/3. Furthermore, we point out that it is possible for the perturbation model to be valid [for instance, the case of C = 5 in Fig. 2(b)] even for large values of z M /z P , such as those beyond the image plane. This is in keeping with the fact that /ϕ is a closed curve if one were to join the z M /z P = ±∞ edges of Fig. 2(b), as discussed in Ref. 8. For larger values of C, however, the NRMSE begins to oscillate, due to the Talbot effect, for values of z M /z P that are sufficiently far from unity (so that r 2 1 /R 2 is large) and this behavior is captured only by the complete error estimate. A notable feature of Fig. 2(b) is the region near z M /z P = 0, where the complete error estimate oscillates rapidly; it is evident that the perturbation model is not valid in this region and this fact is represented by the translucency of the plot. Recall that λ = 632 nm in these simulations and note that the effect of varying λ on the plots of Fig. 2(b) is to change the value of C to which they correspond, proportionally to λ −1/2 . For example, if λ were increased by a factor of 4, the plot for C = 20 would correspond instead to C = 10. Figure 2(a), on the other hand, is explicitly independent of λ.
The simple error estimate in Eq. (29) is the analogous rule of thumb, specifically for the milled and turned sinusoidal MSF groove geometries, to the one-dimensional version in Ref. 8. Although it works well for MSF structures in the form of Eqs. (27) and (28), it was observed in Ref. 8 that such an estimate appears to overestimate the error incurred by MSF structures that possess a broad spectrum; this was demonstrated with synthetic MSF structures with spectra that obeyed a power-decay law. For these specific examples, the simple analog of Eq. (24) proved to be sufficient since it accurately predicted the NRMSE behavior in the region 0< /ϕ<1/3, which is where the perturbation model was considered acceptable. However, it turns out that the consideration of MSF data requires an extension of the rule of thumb predicted by Eq. (25). For such MSF structures, it is necessary to consider the more complete NRMSE expression of Eq. (24).

Rules of thumb for MSF structures with broad spatial spectra
To begin, we present Fig. 3(a) so that it can be used as a reference for further discussion regarding the ineffectiveness of Eq. (25) when applied to MSF structures that are more complicated than those given by Eqs. (27) and (28). There, it is evident that such a simple estimate for /ϕ is useful only for small values of r 2 1 /R 2 . The behavior of the numerically calculated NRMSE departs from the simple estimate very quickly in some examples. Although it is fortunate that Eq. (25) overestimates the true NRMSE, it fails as a rule of thumb for MSF structures with broad spectra (such as those seen in Fig. 3). For instance, someone interested in using the perturbation model with an optical system with MSF structures similar to that color-coded purple in Fig. 3(a), operating at r 2 1 /R 2 ≈ 0.8 with an NRMSE threshold of 20%, would erroneously believe based on Eq. (25) that this model should not be used. Therefore, there is a need for a more complete estimate that is accurate for MSF structures with broad spectra; this is provided by Eq. (24), which can be rewritten in the Fourier domain, after normalization by ϕ 2 , as where the Fourier transform of φ(ξ ⊥ ) is taken to be defined bỹ Eq. (30) shows that the validity of the perturbation model is best understood in the Fourier domain and matches the behavior of the numerically calculated NRMSE values in Fig. 3(a) very well. As done in Fig. 2, an alternative view of Fig. 3(a) is given by Fig. 3(b), which shows the NRMSE as a function of z M /z P . For the numerical simulations performed with the MSF structures shown in Fig. 3, the MSF data were pre-processed to have zero mean and normalized to have RMS of π/8. As was done with the sinusoidal examples in Sec. 4.1, field propagation was modeled by using the angular spectrum. Specifically, the nominal field is initially generated over an array of points at z P and propagated to the location z M , where it is then multiplied by an exponential containing φ, the MSF structure. Depending on the value of the ratio |z M /z P |, the dimensions of the MSF array must be scaled in order to ensure that the converging beam sees the same MSF phase over its diameter. Such a consideration was technically less cumbersome for the analysis in Sec. 4.1 because the MSF structures considered there are periodic and new arrays for φ at any z M can be sampled directly from an analytic sinusoidal function with the appropriate number of cycles. Although the calculated NRMSE does not oscillate quickly near z M /z P = 0 for MSF phases with broad spectra [compare the complete estimates in Fig. 2(b) and Fig. 3(b)], this region is similarly marked with translucency (to indicate the invalidity of the perturbation model regardless of numerically calculated results) because the simulation procedure described earlier in this process becomes invalid as the MSF phase is placed near the vicinity of the caustic at z = 0. Figure 3 shows that, in the paraxial regime, the more complete estimate in Eq. (30) accurately predicts the NRMSE due to the perturbation model. As mentioned earlier, the simple estimate in Eq. (25) always overestimates the NRMSE predicted by Eq. (30). That is, the perturbation model appears to be more valid for MSF structures with broad spectra, as compared to those that are approximated well with a single frequency. One can understand this by noting that the real MSF data used in the analysis is dominated by low-frequency contributions. As a result, the plots shown in Fig. 3(b) are more akin to the plots in Fig. 2(b) with small values of C (such as C = 5 and C = 10) that display a less oscillatory behavior. In other words, the calculated NRMSE for MSF data with broad spectra do not display the oscillatory Talbot re-imaging behavior seen in Fig. 2(b) for the larger C values because the MSF data contain multiple frequencies that wash out the Talbot effect (in particular, validity is mainly governed by the low spatial frequencies).  Fig. 2(b) within a region that shrinks with larger C (or κ). However, despite this shrinking region of accuracy, Eq. (25) is sufficiently accurate for 0 ≤ F /φ<1, which is the region that is relevant for assessing the validity of the perturbation model. The same cannot be said for the MSF examples in Fig. 3(b), whose spectra include both small and large values of κ. Although the presence of large spatial frequencies consistently explains the small region of accuracy of Eq. (25), it is evident that this simple error estimate fails far before /ϕ ≈ 1.

Concluding remarks
We investigated the validity of the perturbation model in three dimensions by using an asymptotic framework like that in Ref. 8. However, we expanded upon the findings in Ref. 8 not only in the consideration of one extra spatial dimension, but also in the completeness of estimates for the error incurred by the perturbation model. That is, through further consistent approximations within the asymptotic framework, it is possible to solve for the complete correctional term [not just the first correction seen in Eq. (22)], which proves to be necessary when considering realistic MSF structures. This more complete approach gives accurate rules of thumb for the validity of the perturbation model. In particular, for the case of imaging systems with low to moderate NA, a general rule of thumb for the validity of the perturbation model for small-amplitude MSF structures is provided in Eq. (30). This more complete error estimate replaces the one in Ref. 8, which involved a troublesome fourth-order spectral moment of the MSF phase structure. These possibly-divergent moments are evidently replaced by the well-behaved integral of Eq. (30). As observed in Ref. 8, when an optical system has more than a single instance of MSF, the total (mean squared) error incurred upon using the perturbation model is simply given by the sum of the (mean squared) error for each individual MSF structure, provided that the MSF structures are statistically uncorrelated with each other. Furthermore, we mention that our results, which are derived by imaging the MSF structure to its conjugate location in image space and then propagating to the exit pupil, would be the same if we had instead imaged the MSF to the neighborhood of the aperture stop and performed the propagation steps there before imaging to the pupil. This invariance is discussed in Appendix B of Ref. 8.
The general framework in Sec. 2 was specialized to that of imaging systems (in fact, the complete error estimates are accurate only for a spherical nominal wavefront in image space), where only three locations in image space are relevant: the image plane, the exit pupil plane, and the plane that is conjugate to the MSF interface itself. The simplest error estimate, F , is suitable for imaging systems with low to moderate numerical apertures and MSF structures that are well described with a single spatial frequency with /ϕ<1. This error estimate is subsumed by the more general one described in Eq. (30), which is nicely represented in the Fourier domain; this estimate was shown to be sufficient in estimating the validity of the perturbation model for MSF structures with broad spectra. The real MSF data used in the analysis contained low spatial frequencies, which dominated the general behavior of the validity of the perturbation model. This, along with the washing-out of the Talbot effect, accounts for why the NRMSE, for many of the real MSF data used in this analysis, is low and not oscillatory when compared with the plots in Fig. 2. Moreover, the simple error estimate of Eq. (25) fails for both the MSF examples in Figs. 2 and 3 outside a domain that shrinks with the presence of higher spatial frequencies. However, the estimate is always sufficiently accurate within the range F /ϕ<1 for MSF structures whose spectra are well-localized. For MSF structures with broad spatial spectra, one must instead use the complete error estimate, C , given by Eq. (30). Even though it was derived in the context of low to moderate NA, we remark that the results in Appendix C lead us to expect that Eq. (30) is, to within a factor of 2 or so, useful for numerical apertures up to 0.8. Therefore, even in the analysis of systems with appreciable NA, it is possible to use Eq. (30) to obtain a rough error estimate of the perturbation model rather than involving the more complicated, and incomplete, discussion regarding systems with high NA in Appendix C.

A. Extended derivation for two spatial dimensions
In this appendix, an alternative approximation for the field, U(x, z), in two spatial dimensions is presented; the definitions of symbols correspond to those given in Appendix A of Ref. 8. One begins by considering the general differential equation for Φ N : where the quantities with an overbar are functions of the ray coordinates (ξ, s). In Ref. 8, the equation for N = 2, which represents the first correction of using the perturbation model, was solved exactly. It turns out, however, that progress can be made by dropping the last (non-linear) term on the right-hand side of Eq. (32), leaving This approximation is justified by the fact that φ, which appears linearly in Φ 1 , is small. Equation (33) can be expressed explicitly in terms of derivatives with respect to the ray coordinates as where J −1 il is the (i, l) element of the J −1 matrix defined in Eq. (37) of Ref. 8. At this point, W is taken to be a converging nominal wavefront, centered at z M : This early specification of W is another discrepancy between the method presented here and that in Ref. 8; however, since imaging systems are of interest, where the wavefront converges to a nominal point in image space, this loss of generality is insignificant for our purposes. An approximation is now made with regards to Eq. (34): the only term retained on the right-hand side is the one that involves ∂ 2 ξ . The reason for making this simplification, which is discussed more in Appendix B, is ultimately due to the fact that the derivatives of φ (which varies quickly under its characterization as MSF) are large when compared to the other terms. These other terms contain nominal quantities and their transverse derivatives. By doing this, and explicitly using Eq. (35) in the definition of J −1 , Eq. (34) is approximated as We propose a change of variables from s tos s/( χ∆), where Recall that which provides the initial right-hand side for Eq. (38) (for the case N = 2). Upon inserting Eq. (39) into Eq. (38) for N = 2, it is possible to drop the contribution due to the second term of Eq. (39) under the assumption that the nominal field quantities vary significantly more slowly than φ. In actuality, this φ-independent term is ultimately included in the perturbation model and should be carried along in the following derivation. However, its explicit form (and subsequently derived expressions) is not important and can henceforth be dropped. Equation (38) now leads to This process can be iterated to give Note that, Eq. (41) is a further approximation sinces depends on ξ and strictly cannot be pulled out of the derivatives in ξ. However all the ξ-dependence ofs is through nominal field quantities and, as was mentioned after Eq. (39), the derivatives of φ vary much more quickly. That is, when integrated over the domain of interest, we assume As a result, in the iterative process of approximately solving for Γ N , it is possible to pull the factors ofs out of the derivatives in ξ.
With Eq. (41), it is now possible to calculate The method presented here differs from that shown in Appendix A of Ref. 8 mainly in two ways. First, the new derivation includes every term in the summation seen in Eq. (43); in Ref. 8, this summation was truncated at N = 2. Second, the specification of W to be a converging spherical wavefront was used. These two differences, along with approximations that φ is small and varies more quickly than the nominal quantities, W and A, allows for a field estimate that gives rise to a complete error estimate upon using the perturbation model that includes contributions from all N. Although the method in Ref. 8 included higher-order corrections in Φ 2 that accounted for larger φ, these terms were ultimately discarded in the development of a simple rule of thumb.

B. MSF-independent rays derivation in three spatial dimensions
Equation (4) can be solved with the method of characteristics, which leads to solutions given in the parametrization of Eq. (6). When making this change of variables, it is convenient to use the transpose of the Jacobian matrix With this, the derivatives in Cartesian coordinates r = (x, y, z) can be written in terms of derivatives in the ray parameters ξ = (ξ, η, s) according to the chain rule, where ∇ ξ (∂ ξ , ∂ η , ∂ s ).
To begin, we show that Eq. (7) is the solution to the Eikonal equation in Eq. (4). By using Eq. (45), and with some simplification, we find that It is simple to see that the dot product of Eq. (46) with itself gives unity, therefore satisfying Eq. (4). Furthermore, by using both Eqs. (45) and (46), we observe that Equation (47) allows us to rewrite the left-hand side of Eq. (5) as a derivative in s. The Eikonal can then be expressed as For N = 1, Eq. (5) reduces to By parameterizing in terms of ξ and using Eq. (47), the left-hand side of Eq. (49) simply becomes ∂ s Φ 1 . For the right-hand side, the parametrization leads to where ∆ = det(J), as given by Eq. (9). Equation (49) therefore becomes which has a simple solution that satisfies the initial conditions of Φ 1 (ξ, η, 0) = ln[A(ξ, η)]+iφ(ξ, η): where the s-dependence is fully encapsulated in ∆. Note that Eq. (52) has a form that is similar to that of the corresponding quantity in the two-dimensional analysis and the logarithmic portion is an amplitude factor that accounts for the bunching of the rays under propagation. This factor diverges at the caustics of these nominal rays.
It is useful to restate the remaining equations for N ≥ 2. Once again, the left-hand side of Eq. (5) can be simplified to ∂ s Φ N . That is, Although progress can be made with Eq. (53), particularly with the case of N = 2 (as discussed separately later), as it is presented, we seek an expression for general N. In order to proceed in this direction, several approximations are made; to begin, we neglect the second term on the right-hand side of Eq. (53), which is a nonlinear term for preceding values of Φ n . This can be justified by the smallness of φ. What remains is where we use the Einstein implicit summation convention. We now consider only MSF structures φ for which there is an appreciable number of cycles across the aperture. This means the quantities within the derivatives and so only the upper-left 2 × 2 submatrix of J T J is relevant. With this, we arrive at which can now be directly integrated to give The processes between Eqs. (61) and (63) can be iterated to give [using an approximation akin to that leading to Eq. (41) in Appendix A] where Tr represents a matrix trace. Note that whereŴ is a differential operator defined aŝ for the plane polar coordinates r = ξ 2 + η 2 and θ = arg(ξ + iη). With Eq. (65), we can rewrite Eq. (64) as With Eq. (67), it is now possible to find Φ through an infinite sum.
C. Discussion of the high NA regime In this section, we discuss the NRMSE expressions given by Eqs. (21) and (22) for systems with high NA. In particular, we examine how going into the high NA regime complicates the simple rules of thumb seen in Sec. 4. The discussion in this section highlights a non-trivial aspect of the generalization to three dimensions; the rule-of-thumb expressions in Eqs. (21) and (22). We restate here that the following results may be inappropriate for a rigorous treatment of high NA systems, where polarization effects can no longer be ignored for field calculations. As observed earlier, the main distinction between the low/moderate NA and high NA error estimates is comprised of the appearance of the obliquity factor χ and the differential operatorŴ (in place of ∇ 2 ⊥ ) for the latter. For convenience, we re-state the approximate NRMSE formulas, taken from Eqs. (22) and (23), as for the low/moderate NA and high NA regimes, respectively. In order to quantify their differences, we chose the nominal amplitude to be unity (A = 1) and consider the following ratio From Sec 4, it was evident that the first-order NRMSE expressions in Eq. (69) are accurate for MSF structures that are well represented by a single frequency (such as those with milled and turned geometries; it turns out to be satisfactory for spoked geometry as well, as shown in Fig. 4). Figure 4 shows Q for the case of a turned, milled, and spoked surface with 10 cycles across the aperture; the behavior of Q is fairly insensitive to the value of C. Note that, although we do not show it explicitly in Sec. 4 and Fig. 2, the first-order NRMSE expression in Eq. (25) works well in the domain of /ϕ<0.6 for the spoked structure (shown in green) in Fig. 4. Furthermore, Fig. 4   To make a connection with the rules of thumb found in Ref. 8 regarding the high NA regime, it is prudent to ask which part of h /ϕ is responsible for most the behavior of Q as the numerical aperture of the system is increased. From Eq. (69), it is clear that the two possible sources are the factor of χ −6 (in the angled brackets of the numerator) and the differential operatorŴ. For MSF structures φ with more than a few cycles across the aperture, we can further approximate h /ϕ by writing the average of a product [ χ −6 and (∇ 2 ⊥ φ) 2 ] as a product of averages. This approximation is justified because χ −6 varies much more slowly across the aperture, even for large numerical apertures, when compared with φ. Furthermore, as can be seen from Eq. (13), the differential operatorŴ can be approximated by ∇ 2 ⊥ so long as the dominant variation of φ is in r rather than θ. However, if the variation in θ dominates, thenŴ ≈ χ 2 ∇ 2 ⊥ . The NRMSE of the perturbation model, for φ that vary dominantly in r or θ, is given by respectively. The corresponding approximate forms of Q are then given by which are both independent of the choice of φ (aside from the inherent assumption regarding its geometry). Figure 4 shows that Q I matches well with the true value predicted by Eq. (70) for the case where φ is a turned MSF structure; for the milled case, there is a notable discrepancy. Once again, this is because a milled MSF structure has non-negligible θ-dependent variations. An example with an even larger discrepancy is that of the spoked geometry, where the variation in θ is much more significant that that in r; in this case, Q II is a much more accurate prediction. For (most) other types of MSF structures, the value of Q (which gives the first-order approximation of the RMS difference between the low/moderate-NA and high-NA NRMSE rules of thumb) will lie in between the values predicted by the two expressions in Eq. (72) (between the black and gray curves in Fig. 4). There are unique exceptions; for instance, one may consider MSF structures of the form φ(r, θ) = h(r/r 0 ) C cos(Cθ), where C is an integer and r 0 is some constant value. For this particular example, ∇ 2 ⊥ φ = 0 butŴ 0. Therefore, Q → ∞ for any numerical aperture. One should keep in mind, though, that this example is very specific and, althoughŴφ 0, it is very near zero [as can be reasoned from the approximations discussed before Eq. (71)].
As a final comment on the comparison of the expressions in Eq. (69), it should be noted that, although Fig. 4 illustrates an intriguing geometrical dependence of Q, its actual value is very close to unity for systems with a moderately large numerical aperature. For an imaging system with a numerical aperture of 0.6, for example, it can be seen that Q 1.2. Therefore, for the purposes of obtaining a rule-of-thumb error estimate for using the perturbation model in systems with moderate numerical apertures, it may be sufficient to use lm /ϕ in Eq. (69).

Disclosures
The authors declare that there are no conflicts of interest related to this article.