Effects on the OTF of MSF structures with random variations.

In diamond-machined freeform manufacturing processes, a tool-tip often leaves behind characteristic mid-spatial frequency (MSF) structures on the optical surface. Unwanted movement between the tool-tip and the part results in MSF structures with random variations. Here, we analyze the effects of these MSF structures on the system's optical performance and derive simple analytic estimates for the optical transfer function in terms of the parameters of these structures. These expressions are expected to aid in MSF tolerancing.


Introduction
Precision freeform manufacturing allows the fabrication of surfaces with large numbers of degrees of freedom. However, these processes typically involve the use of subaperture tools, which leave behind surface grooves whose characteristic frequencies are between those of figure and roughness. These errors, known as mid-spatial frequency (MSF) errors, introduce unwanted effects that are difficult to characterize and complicate tolerancing on manufactured parts [1,2].
In prior work, we explored the effects on the optical transfer function (OTF) of wavefront errors, including quadratic aberrations [3] as well as MSF groove structures due to subaperture tools [4]. In the latter case, two types of MSF groove errors common in freeform manufacturing processes were analyzed: parallel straight grooves resulting from diamond milling, and concentric circular (or spiral) grooves resulting from diamond turning. The key theoretical tool proposed there to obtain analytic expressions for the OTF was referred to as the pupil-difference probability density (PDPD), which is defined as the probability density that two random points in the pupil with given vector separation have a given difference in wavefront error. This prior study assumed regular grooves of equal height and width [4]. However, real processes present tool-tip or spindle vibrations, causing the MSF structures to have profiles are not perfectly periodic. Diamond machining usually involves a fixed tool-tip that makes contact with the optical part as the part is moving (rotating or translating). A local coordinate system is used in which the three directions are called the thrust, cutting, and feed directions of the machine, as shown in Fig. 1(a). It is well established that there can be tool-tip vibration in all three of these directions [5][6][7][8][9][10], and so incorporating such effects into the models for the OTF is of central importance. Vibrations in the thrust and cutting directions cause random variations t i in the depth of each groove, following a probability density K T (t i ) with standard deviation µ, as shown in Fig. 1(c). Vibrations in the feed direction, on the other hand, introduce random lateral groove displacements f i following a probability density K F (f i ) with standard deviation σ. Our goal is to provide simple approximate formulas for the OTF of milled and turned surfaces that incorporate these effects.
The outline of this work is as follows: A review of the PDPD is provided in Sec. 2. The analysis of MSF structures including the effects of tool-tip vibrations is given in Sec. 3 for milled surfaces and in Sec. 4 for turned surfaces. For readers more interested in the formulas than the formal derivations, Sec. 5 provides a self-contained summary of the main results. The size of the pupil is 2R, twice the radius. An example of a MSF structure with random variations is also shown (gray). A sample cycle from (b) is shown in (c) with the probability densities K T (red, with uncertainty µ) governing thrust (and cutting) vibrations, and K F (blue, with uncertainty σ) governing feed vibrations.

Review of general PDPD theory
The PDPD, denoted as P(η, ρ), is defined as the probability density that for two random points in the system's pupil whose vector separation is ρ, the difference in aberration function W is equal to some prescribed quantity η. Mathematically, this can be written in the following form: where denotes a definition, δ is the Dirac delta distribution, and O(ρ) is the overlap region of two pupil copies separated by ρ. Note that the area of O(ρ) normalized by its value at ρ = 0 gives precisely the OTF for a perfect unaberrated system. In the case of a system with a circular pupil of radius R (say, in units of NA or spatial frequency), the perfect system's OTF depends only on ρ = |ρ| and is given by [11] OTF perf (ρ) = The usefulness of the PDPD stems from its one-dimensional Fourier relation with the OTF: where the tilde denotes a Fourier transform from error difference η to wavenumber k. As was shown in [3] and [4], working with the PDPD rather than directly with the OTF allows geometrical insights that lead to useful approximations. This is also the case in the treatment that follows.

Effects of tool-tip vibrations in MSF structures for milled surfaces
The PDPD due to parallel MSF grooves (e.g. resulting from diamond milling) can be approximately expressed in terms of the one-dimensional PDPD of its cross-section [4] by using the scalar pupil displacement variables q n = q ·n and ρ n = ρ ·n, wheren is the unit vector normal to the grooves. The grooves of the nominal MSF structure (without tool-tip vibrations) can be modeled as a piecewise-parabolic function, where h and T are the grooves' nominal peak-to-valley (PV) and width (or period), respectively, as illustrated in Fig. 1(b). It was shown in [4] that such a MSF structure leads to a very simple nominal PDPD (and Fourier transform): where w 0 (h, T, ρ n ) = 8hρ +ρ− , withρ + = (ρ n /T) mod 1,ρ − = (1 − ρ n /T) mod 1 = 1 −ρ + , and rect(u) = 1 for |u| ≤ 1/2 and zero otherwise. The fact that P 1 in Eq. (4) takes the form of a rectangle function can be understood from Fig. 2 where we see that the difference of two displaced versions of W takes on a regular piecewise-linear sawtooth shape (gray lines) constrained in η to a localized region of width w 0 . At this point, it is useful to define the transition function T ρ n as for any f and g. Clearly, T ρ n gives rise to a linear transition between f and g as ρ n passes through the first period. This is useful in Secs. 3 and 4 because the PDPD involves the difference of two copies of W displaced by ρ n . When ρ n < T, this difference includes points from the same period. Although this is unimportant when the MSF structure is periodic (nominal), it is necessary to use T ρ n to account for the case when random variations are introduced to each period independently.

Pistons model for thrust vibrations
Let us start by considering fluctuations in the height of the ridges resulting from variations in thrust (or cutting). We model these fluctuations by letting each parabolic segment be shifted vertically by an amount t i , where each shift is assumed to be independent from all others, following a zero-mean probability density K T (t i ) whose standard deviation µ is assumed to be small compared to h [see Fig. 1(c)]. To begin, we try an overly simplistic model that we call the pistons model, where each complete parabolic segment is rigidly displaced, causing unphysical discontinuities in W as illustrated in Fig. 2(a). The difference of the two oppositely-shifted copies of the error, ∆W = W(q n − ρ n /2) − W(q n + ρ n /2), is then also discontinuous. One can see from Fig. 2 that the rectangular contribution for each straight segment of ∆W is shifted in η by the difference t i − t j . Hence, the mean effect of thrust vibrations is simply given by Equation (6) has the form of a double convolution, so in the Fourier domain it becomes a product: Estimated forms for the probability densities K T are reported in the literature [5,9,10], with models that are based on studies of tool-tip and spindle resonance frequencies. The resulting distribution (and its Fourier transform) governing thrust vibrations is given by where J 0 is the zeroth order Bessel function of the first kind. This leads to the simple expressioñ Note that the difference between the two shifted copies of W in Fig. 2(a) were assumed to involve different parabolic segments, with independent shifts t i and t j . This is a valid assumption as long as ρ n ≥ T. When ρ n < T, the PDPD includes a contribution from the differences of each groove with a displaced version of itself, and hence proportional to P 1 , the nominal PDPD. The weight of this contribution diminishes proportionally to T − ρ n as ρ n increases, and disappears completely once ρ n ≥ T. Hence, we construct the OTF estimate by using the linear transition function in Eq. (5) in the form OTF(k, ρ n ) ≈ OTF perf (ρ)T ρ n (P 1 ,P 1,T ). Figure 3 shows comparisons of this OTF estimate with the average of several numerically computed OTFs of simulated realistic realizations of such surfaces with the corresponding statistics, for several values of the ratio µ/h. In these simulations, each parabolic segment is assigned a random shift t i but the resulting surface is made continuous by extending the lower parabolic segment in each intersection. The plots show that the pistons model predicts well the local maxima of the OTF, but fails even for moderate values of µ/h near the minima (ρ + ≈ 1/2), where the effects of the unphysical discontinuities in W are exacerbated.

Realistic model for thrust vibrations
It turns out that an expression can be found based on a more realistic model with continuous grooves, whose form is almost as simple as that of the pistons model. The key is to study the shift s of the intersections of each segment of ∆W. As can be seen from Fig. 4(b), s depends not only on ρ n but also on the shifts t i and t i+1 of two contiguous grooves, as well as on the shift t j of the groove in the displaced replica that overlaps with this junction. The value of s can be found from simple geometric considerations to be Because the slopes of the segments of ∆W do not change with these shifts, no scaling factor of the probability density accompanies them. Since all intersections, whether they correspond to positive or negative ∆W, undergo shifts of this type, the average effect of thrust vibrations is encapsulated by averaging versions of the nominal distributions shifted by s, weighted by the probability densities of the displacements of the grooves: It is convenient to perform a change of variables so that s is one of the variables of integration along with a = t 1 and b = t i+1 . With this, Eq. (11) can be rewritten as Equation (12) can be viewed as a triple convolution, so in the Fourier domain we havẽ Note that Eq. (13) differs from the result for the pistons model, Eq. (7), in that it has three factors of the characteristic functionK T ; furthermore, the arguments of two of these are dependent on ρ n . For K T given by Eq. (8), we have the following simple expression forP 1,T (k, ρ n ): (14) Figure 3 shows comparisons of OTF perf (ρ)T ρ n (P 1 ,P 1,T ) whereP 1,T is given by Eq. (14), with numerically-computed OTFs for several values of the ratio µ/h. Compared with the pistons model's result in Eq. (7), it is evident that Eq. (14) is more accurate, particularly near the local minima of the OTF. Although the form of K T in Eq. (8) is a realistic model for thrust vibrations, it is worth noting that a zero-mean Gaussian form for K T of the same standard deviation, namely predicts very similar results for the OTF. This is because it is the characteristic functions (and not the probability densities) that appear as factors in Eq. (13), and the characteristic functions in Eqs. (8) and (15) largely agree for kµ < 1. Under this assumption,P 1,T can be written as

Feed vibrations model
We now model the effects of feed vibrations by letting each parabolic segment of the groove structure be displaced laterally by an amount f i , as shown in Fig. 5(a). The displacements f i follow a zero-mean probability density K F whose standard deviation σ is assumed to be small compared to the nominal groove width T [see Fig. 1(c)]. These displacements modify the saw-tooth structure of ∆W such that (unlike for thrust vibrations) the slopes of the different line segments in ∆W change, as shown in Fig. 5(b). Therefore, not only do the limits of the contributions to the PDPD change with the fluctuations, but also their relative weights. The lateral displacements of the grooves cause each straight segment of ∆W to provide a different contribution to the PDPD (illustrated by the green rectangle with width w and height α) when compared to the nominal case (gray, with width w 0 and height w −1 0 ). Due to the statistical symmetry of the intersections, we can make the contributions symmetric around the origin in η to facilitate the calculation. The expressions for w and α are then which after averaging lead to the following expression for the PDPD: The Fourier transform of Eq. (19) is given bỹ Estimated forms of K F are similar to those reported for K T [5,9,10]. However, to obtain an analytic result, we assume a zero-mean normal distribution like that in Eq. (15): Appendix A shows how Eq. (20) can then be solved to give the following expression: where erf(u) is the error function and (23) Figure 6 shows the agreement of OTF perf (ρ)T ρ n (P 1 ,P 1,F ) using Eq. (22) with numericallygenerated results for several values of σ/T. The estimate is accurate near the local minima, but there is a noticeable discrepancy near the local maxima (where ρ/T is an integer) due to the appearance of an erroneous cusp, clearly visible in the insets of Figs. 6(a)-6(c). This discrepancy arises from the model's assumption about the relative order of the intersections between parabolic segments in the two replicas of W. This assumption is violated for values of ρ n that fall roughly within σ of each maximum, where the relative order of these intersections can change. [Notice from Fig. 3

(d) that a similar behavior, albeit less important, was observed for thrust vibration.]
Although a more sophisticated model could be formulated to account for these effects, it would be considerably more complicated, so we propose instead a simple sinusoidal model, based on the observed shape of the numerically-generated OTFs shown in Figs. 6(a)-6(c), namelỹ where C 0 (σ) and C 1 (σ) are chosen so that this sinusoidal coincides with Eq. (22) in regions away from the local maxima. We choose as fitting points the local minima (ρ + =ρ − = 1/2), where Eq. (22) is accurate, and points slightly away from the local maxima (ρ ± = 1/8). The full expressions are given in Appendix A. It turns out, however, that for kh < 3/2 the following simple approximations to those expressions can be used where a 0 (kh)/2 and a 1 (kh), derived in [4], are the coefficients for the Fourier form ofP 1 ; note that C 0 (0) = a 0 /2 and C 1 (0) = a 1 . The full expressions of a 0 and a 1 are given in Appendix B. However, for kh < 3/2, it is sufficient to use their following truncated Taylor series: (26) Figure 6 shows comparisons of the resulting estimate OTF ≈ OTF perf (ρ)T ρ n (P 1 ,P 1,F ) withP 1,F given by Eq. (24), with numerical simulations. It is evident that this simple model works well even at the local maxima, where Eq. (22) fails.

Combination of thrust and feed vibrations
The combination of the two effects described so far, without the sinusoidal fitting technique discussed in Sec. 3.3, is given by the product of Eqs. (16) and (22): The effects of thrust vibrations encapsulated in the exponential in Eq. (27) can also be incorporated in the sinusoidal model described after Eq. (23) in order to arrive at a new simple sinusoidal model for the combined effects of thrust and feed vibrations: with G 0 (µ, σ) and G 1 (µ, σ) given explicitly in Appendix A. For kh < 3/2, the expressions for G 0 and G 1 are approximately given by with C 0 (σ) and C 1 (σ) given by their approximations in Eq. (25). The resulting estimate for the OTF, given by OTF perf (ρ)T ρ n (P 1 ,P 1,C ), is compared to numerical simulations in Fig. 7.

Standard deviation estimate
Tolerancing requires not only an estimate for the mean of the OTF but also for its spread. Based on numerical simulations we propose the following expression for the standard deviation ofP 1,C in terms of the dimensionless combinations kh (nominal PV in waves), 2R/T (number of cycles), σ/T (relative feed vibration strength), and µ/h (relative thrust vibration strength): (30) Figure 8 shows the resulting estimate for the standard deviation for the OTF, namely ±∆OTF(k, ρ n ) ≈ OTF perf (ρ)[T ρ n (P 1 ,P 1,C ± ∆P 1,C ) − T ρ n (P 1 ,P 1,C )] = ±OTF perf (ρ)T ρ n (0, ∆P 1,C ), compared with numerically generated OTF standard deviations assuming normal distributions for K T and K F . Of course, OTF perf is an upperbound to OTF + ∆OTF, regardless of Eq. (30).

Effects of tool-tip vibrations in MSF structures for turned surfaces
We now extend the model to surfaces with concentric grooves such as those approximating the MSF geometry on diamond-turned surfaces, where the description of the groove shapes given earlier applies to the radial cross section, so that the relevant one-dimensional parameter is ρ = |ρ|. As shown in [4], the PDPD for such groove structures can be decomposed in two contributions: one due to the parts of the overlap region where the grooves from two shifted copies of the error align roughly, and one where they do not (as shown in Appendix C.). The latter usually dominates; this gives rise to a strong attenuation of the oscillations in the OTF and makes the effect of variations along the grooves (which break the rotational symmetry) negligible.

Effects of thrust and feed vibrations
As shown in [4], the (Fourier-transformed) PDPD of a nominal groove structure with turned geometry,P 2 , can be calculated from the corresponding quantity for the radial cross-section,P 1 , by expressing the latter as a Fourier series in ρ and multiplying each term by an appropriate factor. SinceP 1,C in Eq. (28) is already a truncated Fourier series, this transformation is straightforward, and (as shown in Appendix B. with a 0 → 2G 0 and a 1 → G 1 ) it leads tõ which combines the effects of thrust and feed vibrations in the same manner that Eq. (27) does.
Here κ ≈ 0.597 and φ 0 ≈ π/5 and C is the Fresnel cosine integral. The resulting OTF estimate is given by OTF perf (ρ)T ρ (P 2 ,P 2,C ); this quantity, withK T given by Eq. (15), is compared with numerical simulations in Fig. 9. Note that the parts proportional to G 0 and G 1 in Eq. (31) correspond to the average value (baseline) and the oscillations of the OTF, respectively. Since the OTF oscillations for a turned surface are small for low values of kh, it is possible to further approximate G 1 (kh) ≈ 0 in Eq. (31) for an even simpler expression forP 2,C .

Summary of results
The OTF estimates derived in this work can be summarized by the following simple expression: where x = ρ ·n for a milled surface with grooves aligned perpendicularly to the unit vectorn and x = ρ = |ρ| for a turned surface, OTF perf (ρ) is given in Eq. (2), G 0 and G 1 are given by B(ρ) 1 and B(ρ) 0.07k 2 h 2 exp(−16ρ 2 /T 2 ) + 1 for milled and turned surfaces, respectively, and the oscillatory function C is defined as where A, B and F are given in Appendix C. Recall that k is the wavenumber, h is the nominal PV of the wavefront error, T is the nominal groove width, and µ and σ are the standard deviations of the the displacements caused by thrust/cutting and feed vibrations, respectively (see Fig. 1). The corresponding estimates for the standard deviations of the OTFs are given by It is important to remember that OTF perf is a strict upperbound to OTF + ∆OTF.

Concluding remarks
Diamond-machined freeform optical surfaces are inevitably affected by the presence of MSF structures, which cause degradation in optical quality and performance. In this work, we extend the analysis of these errors initiated in [4] by incorporating the effects of vibrations in the thrust, cutting and feed directions, resulting from unwanted relative movement between the diamond tool-tip and the optical element under fabrication. The study was based on a mathematical representation referred to as the PDPD, which allows a geometric/visual interpretation of the effects of these errors. By working in the PDPD domain, we show how thrust and feed vibrations modify the nominal PDPD through several probability integrals. Although most of the examples presented assumed that the vibrations follow a normal distribution, the theory is general for any probability distribution. Parallel and concentric MSF groove patterns were considered, and expressions for the corresponding mean OTF and spread were given. It should be noted, though, that several approximations were assumed when modeling MSF structures with random variations in both milled and turned geometries. For example, in the consideration of milled surfaces, Eq. (28) does not account for variations along the grooves' direction; these effects can be incorporated by multiplying the result by the Fourier transform of a PDPD describing those variations. Furthermore, for turned surfaces, Eq. (31) does not account for variations that break the rotational symmetry of turned MSF structures; the effect of these variations are negligible.
The results serve to further complement previous tolerancing tools regarding MSF structures found on freeform optical surfaces. Since all realistic diamond-machining processes involve tool-tip vibrations MSF structures, the expressions presented here provide a means to further predict the degradation of the OTF beyond the nominally periodic MSF structures.

Appendix A: Milled feed vibration derivation
We begin by defining the characteristic functionsK F as With this substitution, Eq. (20) becomes a sextuple integral: After making a change of variables (with unit Jacobian) from f i , f i+1 to q, r according to , the integrals in f j and p j can be solved analytically to yield T 2 qr dqdr dp i dp i+1 . (40) The double integral in the square brackets can also be solved analytically, so that Eq. (40) becomesP By assuming the Gaussian form for K F in Eq. (15), Eq. (41) can be solved to give Eq. (22). The exact forms of the coefficients C 0 and C 1 in Eq. (24) are given by Re exp − T 2 256σ 2 erf γ + for γ 0 γ + (T/2) = γ − (T/2). The corresponding forms for G 0 and G 1 in Eq. (28) are Re exp − T 2 256σ 2 erf γ + (45) Note that, for µ = 0 and σ = 0, a m were found in [4] where erfi(u) is the imaginary error function. Equations (26) give the truncated (to second order) Taylor series of a 0 (kh) and a 1 (kh), which are a valid approximation for kh 3/2.