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

Kevin Liang; G. W. Forbes; Miguel A. Alonso

doi:10.1364/OE.27.003390

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⁻ⁱ^ωt] in a homogeneous transparent medium is governed by the Helmholtz equation

\nabla^{2} U (x, z) + k^{2} U (x, z) = 0,

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₀A(x) exp[ikW(x)], where U₀ 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₀ 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

{| \nabla Φ |}^{2} - 1 = - \frac{1}{i k} \nabla^{2} Φ .

Equation (2) can be solved upon expressing Φ as an asymptotic series with parameter (ik)⁻¹:

Φ (x, z) = \sum_{n = 0}^{\infty} \frac{Φ_{n} (x, z)}{{(i k)}^{n}} .

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

{| \nabla Φ_{0} |}^{2} - 1 = 0,

\nabla Φ_{0} \cdot \nabla Φ_{N} = - \frac{1}{2} (\nabla^{2} Φ_{N - 1} + \sum_{n = 1}^{N - 1} \nabla Φ_{n} \cdot \nabla Φ_{N - n}), N \geq 1 .

The initial conditions discussed above can now be stated as Φ₀(x, z_M) = W(x), Φ₁(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) = ξ + s W^{'} (ξ), 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

χ (ξ) ≜ \sqrt{1 - W^{'} {(ξ)}^{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:

{\bar{Φ}}_{0} (ξ, s) = W (ξ) + s,

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,

\bar{f} (ξ, s) ≜ f [x (ξ, s), z (ξ, 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

{\bar{Φ}}_{1} (ξ, s) = \ln [A (ξ) \sqrt{\frac{χ (ξ)}{Δ (ξ, s)}}] + i ϕ (ξ),

where ∆(ξ, s) is the Jacobian determinant of the coordinate transformation:

Δ (ξ, s) ≜ \frac{\partial (x, z)}{\partial (ξ, s)} = χ (ξ) + \frac{s W^{″} (ξ)}{χ (ξ)} .

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

\bar{Φ_{2}} (ξ, s) = \bar{Ω_{2}} (ξ, s) + \frac{1}{2} [\frac{{ϕ^{'}}^{2}}{χ^{2}} - \frac{i}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ^{3}})}^{'}] ζ + \frac{i ϕ^{'}}{2 χ^{3}} {(\frac{W^{″}}{χ^{2}})}^{'} ζ^{2},

where

\bar{Ω_{2}}

is independent of ϕ, and

ζ ≜ \frac{χ^{2} (ξ) s}{Δ (ξ, s)} = \frac{z - z_{M}}{1 + (z - z_{M}) W^{″} (ξ) / χ^{3} (ξ)} .

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 −χ³/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 − χ³(ξ)/W″(ξ). This follows from Eq. (9) or the observation that

d x = \frac{\partial x}{\partial ξ} d ξ = [1 + (z - z_{M}) W^{″} (ξ) / χ^{3} (ξ)] d ξ,

with ∂x/∂ξ being calculated for fixed z.

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

\begin{array}{l} {\bar{U_{P}} (ξ, s) = \bar{U} (ξ, s) |}_{ϕ (ξ) = 0} \exp [i ϕ (ξ)] \\ = A (ξ) \sqrt{\frac{χ (ξ)}{Δ (ξ, s)}} \exp {i k [W (ξ) + s] - i \frac{\bar{Ω_{2}} (ξ, s)}{k}} \exp [i ϕ (ξ)], \end{array}

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

\bar{U_{P}}

is that contained in

\bar{Φ_{1}}

. The leading correction to this approximation results from also considering the ϕ-dependent terms in

\bar{Φ_{2}}

given in Eq. (10):

\bar{U} (ξ, s) \approx \bar{U_{P}} (ξ, s) \exp {- \frac{1}{2 k} [i \frac{{ϕ^{'}}^{2}}{χ^{2}} + \frac{1}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ^{3}})}^{'}] ζ + \frac{ϕ^{'}}{2 k χ^{3}} {(\frac{W^{″}}{χ^{3}})}^{'} ζ^{2}} .

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

\begin{array}{l} ϵ^{2} (z, z_{M}; ϕ) ≜ \frac{\int_{a} {| \bar{U_{P}} [ξ, (z - z_{M}) / χ (ξ)] - \bar{U} [ξ, (z - z_{M}) / χ (ξ)] |}^{2} [1 + (z - z_{M}) W^{″} (ξ) / χ^{3} (ξ)] d ξ}{\int_{a} {| \bar{U} [ξ, (z - z_{M}) / χ (ξ)] |}^{2} [1 + (z - z_{M}) W^{″} (ξ) / χ^{3} (ξ)] d ξ} \\ = \frac{\int_{a} A^{2} {| 1 - \exp {- \frac{1}{2 k} [i \frac{ϕ^{'}}{χ^{3}} + \frac{1}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ^{3}})}^{'}] ζ (z, ξ) + O (ξ^{2})} |}^{2}}{\int_{a} A^{2} d ξ}, \end{array}

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

\bar{Ω_{2}}

is real-valued, to see that

[1 + (z - z_{M}) W^{″} / χ^{3}] {| \bar{U_{P}} [x, (z - z_{M}) / χ] |}^{2} = A^{2}

. 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

{〈 Q 〉}_{A} ≜ \frac{\int_{a} Q A^{2} d ξ}{\int_{a} A^{2} d ξ},

and expanding the exponential in Eq. (15), one can write a simple approximation for ϵ²:

ϵ^{2} (z, z_{M}; ϕ) \approx \frac{1}{4 k^{2}} {〈 {{[\frac{1}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ^{3}})}^{'}]}^{2} + {(\frac{{ϕ^{'}}^{2}}{χ^{3}})}^{2}} ζ^{2} 〉}_{A} .

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

ϵ^{2} (z, z_{M}; ϕ) \approx \frac{1}{4 k^{2}} {〈 χ^{- 6} ζ^{2} 〉}_{A} {〈 {ϕ^{″}}^{2} + {ϕ^{″}}^{4} 〉}_{A} .

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

\begin{array}{l} φ^{2} ≜ \frac{\int {| U (x, z) |_{ϕ = 0} - U (x, z) |}^{2} d x}{\int {| U (x, z) |_{ϕ = 0} |}^{2} d x} = \frac{\int {| U (x, z_{M}) |_{ϕ = 0} - U (x, z_{M}) |}^{2} d x}{\int {| U (x, z_{M}) |_{ϕ = 0} |}^{2} d x} \\ = {〈 | 1 - \exp (i ϕ) |^{2} 〉}_{A} = 4 {〈 \sin^{2} (ϕ / 2) 〉}_{A} \approx {〈 ϕ^{2} 〉}_{A}, \end{array}

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

ϵ (z, z_{M}; ϕ) \approx \frac{| z - z_{M} |}{2 k} \sqrt{{〈 {ϕ^{″}}^{2} + {ϕ^{'}}^{4} 〉}_{1}} .

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

ϕ (ξ) = h \sin (2 π κ ξ),

where h is the amplitude measured in radians and κ is its spatial frequency, this RMSE estimate is then given by

ϵ (z, z_{M}; ϕ) \approx | z - z_{M} | π λ h κ^{2} \sqrt{\frac{1}{2} + \frac{3 h^{2}}{8}} \approx \frac{| z - z_{M} | π λ h κ^{2}}{\sqrt{2}},

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, $φ = {〈 | 1 - \exp (i ϕ) |^{2} 〉}_{1} = \sqrt{2 [1 - J_{0} (h)]} \approx h / \sqrt{2}$ , where J₀ 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

Z ≜ \frac{1}{π λ κ^{2}} \frac{φ}{h \sqrt{\frac{1}{2} + \frac{3 h^{2}}{8}}} \approx \frac{1}{π λ κ^{2}},

Fig. 1 NRMSE plots for a sinusoidal MSF phase structure on a collimated beam, as a function of $| z - z_{M} | / Z$ , for various values of h and C, with L = 2 cm. The lines indicated in the legend are the NRMSE between the perturbation model and numerically calculated exact fields. These NRMSE are compared with Eq. (22), which is shown as a thick black line.

Download Full Size | PDF

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.

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

W (ξ) = z_{M} \sqrt{1 + \frac{ξ^{2}}{z_{M}^{2}}}, A (ξ) = \frac{| z_{P} |^{1 / 2}}{{(z_{M}^{2} + ξ^{2})}^{1 / 4}},

from which the solutions for Φ₀ and Φ₁ follow from Eqs. (7)-(9):

Φ_{0} (x, z_{P}) = z_{P} \sqrt{1 + \frac{x^{2}}{z_{P}^{2}}},

Φ_{1} (x, z_{P}) = \ln [\sqrt{| \frac{z_{P}}{Φ_{0} (x, z_{P})} |}] + i \frac{π}{4} [sgn (z_{P}) - sgn (z_{M})] + i ϕ (x \frac{z_{M}}{z_{P}}),

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

Fig. 2 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₁ 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.

Download Full Size | PDF

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:

| ζ |_{z = z_{P}} | = | \frac{z_{P} - z_{M}}{z_{P} / z_{M}} | = \frac{r_{1}^{2} (z_{P}, z_{M})}{λ},

where this quantity is expressed in terms of the radius r1 of the first Fresnel zone [see illustration in Fig. 2(b)], defined as

r_{1} (z_{P}, z_{M}) ≜ \sqrt{λ | \frac{(z_{M} - z_{P}) z_{M}}{z_{P}} |} .

Equation (18) then becomes

ϵ (z_{P}, z_{M}; ϕ) \approx \frac{r_{1}^{2} (z_{P}, z_{M})}{4 π} \sqrt{{〈 χ^{- 6} 〉}_{A}} \sqrt{{〈 {ϕ^{″}}^{2} + {ϕ^{'}}^{4} 〉}_{A}} \approx \frac{r_{1}^{2} (z_{P}, z_{M})}{4 π} \sqrt{{〈 χ^{- 6} 〉}_{A}} \sqrt{{〈 {ϕ^{″}}^{2} 〉}_{A}},

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 $\sqrt{{〈 χ^{- 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

ϵ (z_{P}, z_{M}; ϕ) \approx \frac{r_{1}^{2} (z_{P}, z_{M})}{ℛ^{2}} φ,

where

ℛ

is a measure of the characteristic feature size of the MSF structure at z_M, defined as

ℛ ≜ 2 \sqrt{π} {(\frac{{〈 ϕ^{2} 〉}_{A}}{{〈 {ϕ^{″}}^{2} 〉}_{A}})}^{1 / 4} .

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₁/ℛ)², 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 r1 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 λ | z_{M} |$ . Similarly, the limit of large z_M can be handled by dividing both ℛ and r₁ 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{π}$ , i.e. $ℛ \approx 1 / (\sqrt{π} κ)$ . 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,

\sqrt{{〈 χ^{- 6} 〉}_{A}} = \sqrt{\frac{35 - 70 η^{2} + 56 η^{4} - 16 η^{6}}{35 {(1 - η^{2})}^{3}}},

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

C_{\max}^{2} ≜ \frac{2 η L (1 - η^{2})}{3 π λ} | \frac{z_{M} / z_{P}}{z_{M} / z_{P} - 1} | \sqrt{\frac{35}{35 - 70 η^{2} + 56 η^{4} - 16 η^{6}}} .

Fig. 3 (a) C_max as a function of η and z_M/z_P given by Eq. (32) for L = 20 mm and λ = 632.8 nm. Note that the horizontal axis is scaled such that it is proportional to tan⁻¹(2z_M/z_P − 1). (b) NRMSE for η = 0.05 as a function of $r_{1}^{2} / ℛ^{2}$ , for various values of C and h. The black line is the estimate from Eq. (29), and is seen to agree well with the numerically calculated values, shown as colored dots, particularly for small values of h.

Download Full Size | PDF

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² (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} / ℛ^{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} / ℛ^{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.

Fig. 4 Plots for an imaging system with η = 0.2, L = 20 mm, and a sinusoidal MSF structure. (a) NRMSE as a function of z_M/z_P for h = π/8 and various values of C given by the legend. The solid curves are the estimate from Eq. (29) and the dots are numerically calculated actual values. The green region indicates ϵ ≤ φ/3. (b) The same information, plotted against $r_{1}^{2} / ℛ^{2}$ , where the single black line is the estimate. (c/d) show similar information to (a/b), but with a fixed value of C = 20 cycles and various values of h given by the legend. Note that in (c), unlike (a), the estimate is represented by a single black curve for all examples since they all have the same value of C.

Download Full Size | PDF

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

ϕ (ξ) = \sum_{m} a_{m} \exp (2 π i κ_{m} ξ),

where κ_m = m/L_M. If the field amplitude is roughly constant, it now follows that

ℛ = \frac{1}{\sqrt{π}} {(\bar{κ^{4}})}^{- 1 / 4}

, where

\bar{κ^{4}}

is a normalized fourth moment of the spatial frequency spectrum, i.e.

\bar{κ^{4}} ≜ \frac{\sum_{m} | a_{m} |^{2} κ_{m}^{4}}{\sum_{m} | a_{m} |^{2}} .

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} / ℛ^{2}$ is less than about a third. (Appendix D provides a description of this error for larger values of $r_{1}^{2} / ℛ^{2}$ ).

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_{1}^{2} / ℛ^{2}$ , where the single black line is the estimate. The green region indicates ϵ ≤ φ/3.

Download Full Size | PDF

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 \sum_{m = 1}^{M} ϕ_{m} (x z_{m} / z_{P})]$ . The generalization of the dominant part of Eq. (14) results from simply introducing a summation over these contributions:

U (x, z_{P}) \approx U_{P} (x, z_{P}) \exp [- \frac{1}{2 k} \sum_{m = 1}^{M} \frac{ϕ_{m}^{″} (x z_{m} / z_{P})}{χ^{3} (x z_{m} / z_{P})} ζ_{m}],

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

ϵ_{T}^{2} \approx \sum_{m = 1}^{M} ϵ^{2} (z_{P}, z_{m}; ϕ_{m}) .

Note that ϵ(z_P, z_m; ϕ_m) < φ_m/3 should be checked individually, where φ_m is given by Eq. (19) for the corresponding ϕ_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:

J ≜ \frac{\partial r}{\partial ξ} = [\begin{matrix} \frac{\partial x}{\partial ξ} & \frac{\partial z}{\partial ξ} \\ \frac{\partial x}{\partial s} & \frac{\partial z}{\partial s} \end{matrix}] = [\begin{matrix} 1 + s W^{″} (ξ) & - s \frac{W^{'} (ξ) W^{″} (ξ)}{χ (ξ)} \\ W^{'} (ξ) & χ (ξ) \end{matrix}] = [\begin{matrix} 1 + s W^{″} (ξ) & s χ^{'} (ξ) \\ W^{'} (ξ) & χ (ξ) \end{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,

\nabla F [x (ξ), z (ξ)] = \bar{\nabla F} (ξ, s) = J^{- 1} \cdot \nabla_{ξ} \bar{F} (ξ, s),

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,

\bar{\nabla Φ_{0}} = (0 1) \cdot J = (W^{'}, χ) .

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:

\bar{\nabla Φ_{0}} \cdot \bar{\nabla Φ_{N}} = (0 1) \cdot \nabla_{ξ} \bar{Φ_{N}} = \frac{\partial}{\partial s} \bar{Φ_{N}} .

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

\nabla Φ_{0} \cdot \nabla Φ_{1} = - \frac{1}{2} \nabla^{2} Φ_{0} .

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

\bar{\nabla_{r}^{2} Φ_{0}} = (J^{- 1} \cdot \nabla_{ξ}) \cdot \bar{\nabla Φ_{0}} = (J^{- 1} \cdot \nabla_{ξ}) \cdot [(0 1) \cdot J] = Tr (J^{- 1} \cdot \frac{\partial}{\partial s} J) = \frac{\partial}{\partial s} \ln (Δ) .

with

Δ = Det (J) = (χ^{2} + s W^{″}) / χ

. With these results, Eq. (41) now becomes

\frac{\partial}{\partial s} \bar{Φ_{1}} = - \frac{1}{2} \frac{\partial}{\partial s} \ln (Δ),

which has a simple solution that satisfies

\bar{Φ_{1}} (ξ, 0) = \ln [A (ξ)] + i ϕ (ξ)

:

\bar{Φ_{1}} (ξ, s) = \ln (A \sqrt{\frac{χ}{Δ}}) + 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

\nabla Φ_{0} \cdot \nabla Φ_{2} = - \frac{1}{2} \frac{\nabla^{2} \exp (Φ_{1})}{\exp (Φ_{1})} .

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

\begin{matrix} \frac{\partial \bar{Φ_{2}}}{\partial s} = - \frac{1}{2} \exp (- \bar{Φ_{1}}) \bar{\nabla^{2} \exp (Φ_{1})} \\ = - \frac{1}{2} \frac{1}{A \exp (i ϕ)} \sqrt{\frac{Δ}{χ}} J_{i l}^{- 1} \frac{\partial}{\partial ξ_{l}} {J_{i j}^{- 1} \frac{\partial}{\partial ξ_{j}} [A \sqrt{\frac{χ}{Δ}} \exp (i ϕ)]}, \end{matrix}

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:

\begin{matrix} \frac{\partial \bar{Φ_{2}}}{\partial s} = a - \frac{i ϕ^{'}}{2 A^{2} χ} (1, - W^{'}) \cdot \nabla_{ξ} (\frac{A^{2} χ}{Δ^{2}}) + \frac{{ϕ^{'}}^{2} - i ϕ^{″}}{2 Δ^{2}} \\ = a - \frac{i}{2 A^{2} χ} {(\frac{A^{2} χ ϕ^{'}}{Δ^{2}})}^{'} + \frac{i χ^{'} ϕ^{'}}{Δ^{3}} + \frac{{ϕ^{'}}^{2}}{2 Δ^{2}}, \end{matrix}

where a is the contribution that is independent of the MSF perturbation ϕ, and is given by

a = - \frac{1}{2 A} \sqrt{\frac{Δ}{χ}} (J^{- 1} \cdot \nabla_{ξ}) \cdot [J^{- 1} \cdot \nabla_{ξ} (A \sqrt{\frac{χ}{Δ}})] .

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 $\bar{Ω_{2}}$ . To integrate the rest, we use the simple results

\int_{0}^{s} \frac{d s}{Δ^{2}} = \frac{s}{χ Δ}, \int_{0}^{s} \frac{d s}{Δ^{3}} = \frac{(χ + Δ) s}{2 χ^{2} Δ^{2}} .

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

\begin{matrix} \bar{Φ_{2}} (ξ, s) = \bar{Φ_{2}} (ξ, s) - \frac{i}{2 A^{2} χ} {(\frac{A^{2} ϕ^{'} s}{Δ})}^{'} + \frac{i χ^{'} ϕ^{'} (χ + Δ) s}{2 χ^{2} Δ^{2}} + \frac{{ϕ^{'}}^{2} s}{2 χ Δ} \\ = \bar{Φ_{2}} (ξ, s) - \frac{i}{2 A^{2}} {(\frac{A^{2} ϕ^{'}}{Δ})}^{'} \frac{s}{Δ} + \frac{i ϕ^{'} (χ^{'} + Δ^{'}) s}{2 χ Δ^{2}} + \frac{{ϕ^{'}}^{2}}{2 χ} \frac{s}{Δ} \\ = \bar{Φ_{2}} (ξ, s) + \frac{1}{2} [\frac{{ϕ^{'}}^{2}}{χ} - \frac{i}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ})}^{'}] \frac{s}{Δ} + \frac{i χ^{'} ϕ^{'}}{χ} \frac{s}{Δ^{2}} (\frac{Δ}{χ} - \frac{W^{″} s}{χ^{2}}) + \frac{i ϕ^{'}}{2 χ} {(\frac{W^{″}}{χ})}^{'} \frac{s^{2}}{Δ^{2}} \\ = \bar{Φ_{2}} (ξ, s) + \frac{1}{2} [\frac{{ϕ^{'}}^{2}}{χ} - \frac{i}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ})}^{'} + 2 \frac{i χ^{'} ϕ^{'}}{χ^{2}}] \frac{s}{Δ} + \frac{i ϕ^{'}}{2 χ} [{(\frac{W^{″}}{χ})}^{'} - 2 \frac{χ^{'} W^{″}}{χ^{2}}] \frac{s^{2}}{Δ^{2}} \\ = \bar{Ω_{2}} (ξ, s) + \frac{1}{2} [\frac{{ϕ^{'}}^{2}}{χ^{3}} - \frac{i}{A^{2}} {(\frac{A^{2} ϕ^{'}}{χ^{3}})}^{'}] ζ + \frac{i ϕ^{'}}{2 χ^{2}} {(\frac{W^{″}}{χ^{3}})}^{'} ζ^{2}, \end{matrix}

where ζ = χ²s/Δ = (z − z_M)/[1 + (z − z_M)W^″ (ξ)/χ³(ξ)]. Notice that a factor of 1 = Δ/χ − sW^″ /χ² 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

T (p, p^{'}) = f p p^{'} + \frac{δ}{2 n} p^{2} - \frac{δ^{'}}{2 n} p^{2} - \frac{δ^{'}}{2 n^{'}} {(p^{'})}^{2},

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 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 ², 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 ²/δ^′. Our estimate of error in the perturbation model involves the ratio of a characteristic length in the MSF-bearing plane, say

ℒ_{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

\frac{n^{'}}{λ} (\frac{1}{δ_{m}^{'} - δ_{i}^{'}} + \frac{1}{δ_{p}^{'} - δ_{m}^{'}}) {(ℒ_{m}^{'})}^{2} = \frac{n^{'}}{λ} \frac{δ_{p}^{'} - δ_{i}^{'}}{(δ_{m}^{'} - δ_{i}^{'}) (δ_{p}^{'} - δ_{m}^{'})} {(ℒ_{m}^{'})}^{2}

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

\begin{matrix} \frac{n^{'}}{λ} \frac{δ_{p}^{'} - δ_{i}^{'}}{(δ_{m}^{'} - δ_{i}^{'}) (δ_{p}^{'} - δ_{m}^{'})} {(ℒ_{m}^{'})}^{2} = - \frac{1}{λ n f^{2}} \frac{1 / δ_{p} - 1 / δ_{i}}{(1 / δ_{m} - 1 / δ_{i}) (1 / δ_{p} - 1 / δ_{m})} {(ℒ_{m}^{'})}^{2} \\ = \frac{n}{λ} \frac{δ_{p} - δ_{i}}{(δ_{m} - δ_{i}) (δ_{p} - δ_{m})} {(\frac{δ_{m}}{n f} ℒ_{m}^{'})}^{2} = \frac{n}{λ} \frac{δ_{p} - δ_{i}}{(δ_{m} - δ_{i}) (δ_{p} - δ_{m})} ℒ_{m}^{2}, \end{matrix}

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, $\tilde{U}$ , by

U (x, z_{P}) = \frac{1}{λ} \int \tilde{U} (α) \exp {i k [α x + \sqrt{1 - α^{2}} (z_{P} - z_{M})]} d α,

where

\tilde{U}

is related to the initial field U(ξ, z_M) by

\tilde{U} (α) = \int U (ξ, z_{M}) \exp (- i k α ξ) d ξ .

For the specific case of a converging wave that is perturbed by a MSF structure, U(ξ, z_M) = U₀ A(ξ) exp[iϕ(ξ)] exp[ikW(ξ)], 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

\begin{matrix} \tilde{U} (α) = \int a (ξ) \exp [i k γ (ξ)] d ξ \\ = \exp [i k γ (ξ_{0})] \int (a (ξ_{0} + τ) \exp {i k [γ (ξ + τ) - γ (ξ_{0}) - γ^{″} (ξ_{0}) \frac{τ^{2}}{2}]}) \\ \times \exp [i k γ^{″} (ξ_{0}) \frac{τ^{2}}{2}] d τ, \end{matrix}

where a(ξ) = U₀ 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 ξ = ξ₀ + τ, where ξ₀ is the stationary point of γ, i.e. γ^′(ξ₀) = 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

\int τ^{2 n} \exp (i \frac{k}{2} τ^{2} β) d τ = k^{- (n + 1 / 2)} \sqrt{\frac{2 π}{| β |}} \frac{(2 n - 1)!!}{| β |^{n}} \exp [i \frac{π}{4} sgn (β)],

with β = γ″ (x0). 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₀ exp[ikΦ₀(x, z_P) + Φ₁(x, z_P)], with Φ₀ and Φ₁ 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: ${\hat{K}}_{z}$ , which denotes free propagation by a distance z; ${\hat{M}}_{M}$ , which denotes the MSF phase factor at z_M, the image plane for the surface in question; ${\hat{M}}_{P}$ , which denotes the MSF written at the pupil plane z_P according to the perturbation approximation; and ${\hat{a}}_{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 | {\hat{K}}_{z} | x^{'} 〉$ . This field is a converging wave if z < 0 and a diverging one if z > 0. Note that $〈 x | {\hat{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

U (x, 0; x^{'}) = 〈 x | {\hat{K}}_{- z_{P}} {\hat{a}}_{P} {\hat{K}}_{z_{P}} - z_{M} {\hat{M}}_{M} {\hat{K}}_{z_{M}} | x^{'} 〉,

where the argument following the semicolon indicates the point at which the wave is focused. The perturbation approximation to this result, on the other hand, is given by

U_{P} (x, 0; x^{'}) = 〈 x | {\hat{K}}_{- z_{P}} {\hat{a}}_{P} {\hat{M}}_{P} {\hat{K}}_{z_{P}} | x^{'} 〉 = 〈 x | {\hat{K}}_{- z_{P}} {\hat{a}}_{P} {\hat{M}}_{P} {\hat{K}}_{z_{P} - z_{M}} | x^{'} 〉 .

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

ε_{I}^{2} (x^{'}) = N^{- 1} \int | U_{P} (x, 0; x^{'}) - U (x, 0, x^{'}) |^{2} d x = N^{- 1} 〈 x^{'} | {\hat{K}}_{- z_{M}} {\hat{T}}_{I}^{†} {\hat{a}}_{P}^{2} {\hat{T}}_{I} {\hat{K}}_{z_{M}} | x^{'} 〉,

with

N (x^{'}) = \int | U (x, 0; x^{'}) |^{2} d x

,

{\hat{T}}_{I} = {\hat{K}}_{z_{P} - z_{M}} {\hat{M}}_{M} - {\hat{M}}_{P} {\hat{K}}_{z_{P} - z_{M}}

, and where we used the fact that

{\hat{a}}_{P}^{†} = {\hat{a}}_{P}

,

{\hat{K}}_{z}^{†} = {\hat{K}}_{- z}

, and

\int | x 〉 〈 x | d x = \hat{1}

. Note also that for a binary aperture,

{\hat{a}}_{P}^{2} \propto {\hat{a}}_{P}

, but we do not make this simplification so that the result applies also to apodized apertures.

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

U (x, 0; x^{'}) = 〈 x | {\hat{K}}_{- z_{M}} {\hat{M}}_{M} {\hat{K}}_{z_{M} - z_{P}} {\hat{a}}_{P} {\hat{K}}_{z_{P}} | x^{'} 〉,

and its corresponding perturbation approximation is given by

U_{P} (x, 0; x^{'}) = 〈 x | {\hat{K}}_{- z p} {\hat{M}}_{P} {\hat{a}}_{P} {\hat{K}}_{z_{P}} | x^{'} 〉 = 〈 x | {\hat{K}}_{- z_{M}} {\hat{K}}_{z_{M} - z_{P}} {\hat{M}}_{P} {\hat{a}}_{P} {\hat{K}}_{z_{P}} | x^{'} 〉,

so that, by using

{\hat{T}}_{II} = {\hat{K}}_{z_{P} - z_{M}} {\hat{M}}_{M}^{†} - {\hat{M}}_{P}^{†} {\hat{K}}_{z_{P} - z_{M}}

, the square of the RMSE can be written as

ϵ_{II}^{2} (x^{'}) = N^{- 1} 〈 x^{'} | {\hat{K}}_{- z_{P}} {\hat{a}}_{P} {\hat{T}}_{II} {\hat{T}}_{II}^{†} {\hat{a}}_{P} {\hat{K}}_{z_{P}} | x^{'} 〉,

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 $N$ by $ℳ = \iint | U (x, 0; x^{'}) |^{2} d x d x^{'}$ . Assuming that the range of integration is sufficiently large to approximate $\int | x^{'} 〉〈 x^{'} | d x^{'} \approx \hat{1}$ , we arrive at

\bar{ϵ_{I}^{2}} = ℳ^{- 1} \int N ϵ_{I}^{2} d x^{'} \approx ℳ^{- 1} Tr ({\hat{K}}_{- z_{M}} {\hat{T}}_{I}^{†} {\hat{a}}_{P}^{2} {\hat{T}}_{I} {\hat{K}}_{z_{M}}) = ℳ^{- 1} Tr ({\hat{a}}_{P}^{2} {\hat{T}}_{I} {\hat{T}}_{I}^{†}),

\bar{ϵ_{II}^{2}} = ℳ^{- 1} \int N ϵ_{II}^{2} d x^{'} \approx ℳ^{- 1} Tr ({\hat{K}}_{- z_{M}} {\hat{a}}_{P} {\hat{T}}_{II} {\hat{T}}_{II}^{†} {\hat{a}}_{P} {\hat{K}}_{z_{P}}) = ℳ^{- 1} Tr ({\hat{a}}_{P}^{2} {\hat{T}}_{II} {\hat{T}}_{II}^{†}),

where we used the property

Tr (\hat{A} \hat{B}) = Tr (\hat{B} \hat{A})

. The error estimates in Eqs. (64) and (65) coincide if the MSF phase structure of the first equals the negative of that of the second, since the only difference between

{\hat{T}}_{I}

and

{\hat{T}}_{II}

is that

{\hat{M}}_{M}

in the former is replaced by

{\hat{M}}_{M}^{†}

in the latter. Since in each situation the error measures

ϵ_{I, II}^{2} (\geq 0)

are probably not strongly dependent on the image points 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:

U_{i} (ξ, z_{M}) = \frac{1}{\sqrt{i λ (z_{M} - z_{P})}} \int_{- \infty}^{\infty} U_{i} (x, z_{P}) \exp [\frac{i π {(ξ - x)}^{2}}{λ (z_{M} - z_{P})}] d x .

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 + ih sin(2πκξ). We can now back-propagate this field to the pupil plane through an (inverse). Fresnel propagation integral:

\begin{array}{l} U_{ii} (x, z_{P}) \approx \frac{1}{\sqrt{- i λ (z_{M} - z_{P})}} \int_{- \infty}^{\infty} U_{i} (ξ, z_{M}) [1 + i h \sin (2 π κ ξ)] \exp [- \frac{i π {(x - ξ)}^{2}}{λ (z_{M} - z_{P})}] d ξ \\ = U_{i} (z, z_{P}) + \frac{h}{2} \exp (i π κ x_{κ}) [U_{i} (x + x_{κ}, z_{P}) \exp (i 2 π κ x) - U_{i} (x - x_{κ}, z_{P}) \exp (- i 2 π κ x)], \end{array}

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₀ exp[iπx²/(λ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₀ exp[iπx²/(λ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

U_{P} (x, z_{P}) \approx U_{0} \exp (\frac{i π x^{2}}{λ z_{P}}) rect (\frac{x}{L}) [1 + i h \sin (2 π κ \frac{z_{M}}{z_{P}} x)] .

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₀|², gives the square of the RMSE. After simplification, this normalized integral can be written some

ϵ_{j}^{2} \approx \frac{h^{2}}{4 L} \int {| \sum_{σ = \pm} σ \exp (σ i 2 π \frac{C}{L} x) [r_{σ, j} (x) \exp (i π \frac{C}{L} x_{κ}) - rect (\frac{x}{L})] |}^{2} d x,

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

\frac{ϵ_{I}^{2}}{φ^{2}} \approx 4 [1 - sinc (2 π C)] \sin^{2} (\frac{π κ^{2} r_{1}^{2}}{2}),

\frac{ϵ_{II}^{2}}{φ^{2}} \approx 2 {1 + \frac{κ^{2} r_{1}^{2} - C}{C} \cos (π κ^{2} r_{1}^{2}) + \frac{1 + 2 \cos (π κ^{2} r_{1}^{2})}{π C} \sin^{2} (\frac{π κ^{2} r_{1}^{2}}{2}) \sin [2 π (C - κ^{2} r_{1}^{2})]},

where we used ϕ² ≈ h²/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_{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:

C ≫ Max [1, \frac{π κ^{2} r_{1}^{2}}{2} \sin^{- 2} (\frac{π κ^{2} r_{1}^{2}}{2})] .

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 ξ₀, 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πκξ₀). 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 ξ₀ of the sinusoidals, especially because this has a similar effect to the averaging over different object points mentioned earlier. These averages for both $ϵ_{I}^{2}$ and $ϵ_{II}^{2}$ give the simpler (and more representative) results

\bar{ϵ_{I}^{2}} / φ^{2} \approx 4 \sin^{2} (π κ^{2} r_{1}^{2} / 2)

\bar{ϵ_{II}^{2}} / φ^{2} \approx 4 \sin^{2} (π κ^{2} r_{1}^{2} / 2) + 2 \frac{κ^{2} r_{1}^{2}}{C} \cos (π κ^{2} r_{1}^{2} .)

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:

ϵ_{I}^{2} \approx \sum_{m} | a_{m}^{2} | 4 \sin^{2} (\frac{π κ_{m}^{2} r_{1}^{2}}{2}) \approx φ^{2} \frac{Σ_{m} {| a_{m} |}^{2} 4 \sin^{2} (π κ_{m}^{2} r_{1}^{2} / 2)}{Σ_{m} {| a_{m} |}^{2}},

where we used

φ^{2} \approx Σ_{m} | a_{m}^{2} |

. If

κ_{m}^{2} r_{1}^{2} ≪ 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.

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. 262538–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. 1733–55 (1966). [CrossRef]

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

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

12. M. Rimmer, “Analysis of perturbed lens systems,” Appl. Opt. 9533–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]

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

Abstract

1. Introduction

2. Asymptotic propagation estimate based on nominal rays

3. Simple field error estimates in a homogeneous medium

4. Application to imaging systems

5. Concluding remarks

A. MSF-independent rays derivation

B. Invariance of the error estimate in conjugate spaces

C. Perturbation model in angular spectrum domain

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

Funding

Acknowledgments

References

Cited By

Figures (5)

Equations (76)

Optics Express