## 1.

## Introduction

Recent advances in manufacturing methods for gradient index optical components would benefit from being accompanied by nondestructive metrology techniques that are faster than current precision interferometric methods. Manufacturers require measurement of the spatially varying refractive index profiles in these components to evaluate parts and processes and to locate and orient the index distribution for cutting and forming. This paper describes techniques using geometric optics to reconstruct the refractive index distribution in a sample based on measurement of laser beam displacement, deflection, and, to a lesser degree, mode conversion after passing through samples with axial, radial (cylindrical), or spherical gradients of refractive index. Although one could imagine a number of implementations, the general scheme would likely involve a laser, a series of scattering surfaces such as partially scattering screens or even surface fogging or other temporary surface coatings, and an imaging device. By imposing controlled relative motion of the apparatus and the sample, the laser would pass into the sample at a series of entry points and angles, perhaps measured via images of scattering at the entry surface. Images of scattering of the emergent beam at surfaces on the far side of the sample would allow calculation of the exit location, angle, and possibly the mode shape.

Ray paths confined to a plane, as occur in a plane of refractive index symmetry, obey

## (1)

$$-\frac{n{y}^{\prime \prime}}{1+{y}^{\prime 2}}=\frac{\partial n}{\partial x}{y}^{\prime}-\frac{\partial n}{\partial y},$$Due to similarities in symmetry, radial gradients behave like spherical gradients when viewed normal to the axis of symmetry and like axial gradients when viewed along the axis of symmetry. Therefore, the following three cases suffice to cover the three basic gradient structures: axial viewed normal to the gradient, axial viewed parallel to the gradient, and spherical. The methods are constructed using basic manufacturing shapes for optical blanks (slabs for the axial gradients, rods for the cylindrical gradients, and spheres for the spherical gradients); however, as long as the cut surface shape is known, the methods apply equally to finished components.

Previous work by others, beginning in 1893 with Wiener,^{2} includes examining index reconstruction

## 2.

## Axial Gradient Across a Slab

For this case (see Fig. 1), Eq. (1) reduces to

which can be integrated to produce## (3)

$$\frac{1+\frac{1}{{y}^{\prime 2}}}{1+\frac{1}{{y}_{0}^{\prime 2}}}=\frac{{n}^{2}(x)}{{n}_{0}^{2}},$$## (4)

$$\mathrm{\Delta}y\equiv {y}_{1}-{y}_{0}=\pm {\int}_{0}^{t}\frac{\mathrm{d}x}{\sqrt{\frac{{n}^{2}(x)}{{n}_{0}^{2}}(1+\frac{1}{{y}_{0}^{\prime 2}})-1}}.$$Because the angular deflection and vertical displacement of the beam will not depend on initial vertical position, the mode shape of a beam will be preserved through the slab, and no additional information can be obtained using beam mode shapes. The deflection of the beam ${y}_{1}^{\prime}({y}_{0}^{\prime})$ depends only on the exit index, not on the interior variation. Only the variation of displacement $\mathrm{\Delta}y$ with the entrance slope will provide information about the interior index profile. Without loss of generality, what follows uses only the positive version of Eq. (4) corresponding to ${y}_{0}^{\prime}\ge 0$.

Two basic approaches to solving Eq. (4) could be considered. Particularly in a manufacturing metrology setting where the element has a known target index profile, an iterative approach to solving a discrete version of Eq. (4) would probably be appropriate. In such an approach, the postulated index distribution would be iteratively adjusted, starting with an initial guess given by the target profile to minimize the magnitude of the residual between the measured beam displacement and that predicted by the evaluation of the right-hand side of Eq. (4) using the postulated index distribution.

Alternatively, one could attempt to construct the profile directly. An analytic approach is possible using methods of integral equations under limiting assumptions about the profile. For instance, if the profile is assumed to be monotonic with position [so that $n(x)$ can be inverted to $x(n)$], the following transformation could be applied (without loss of generality assume ${n}_{1}>{n}_{0}$):

## (5)

$$f(\zeta )\equiv \frac{\mathrm{\Delta}y[{y}_{0}^{\prime}(\zeta )]}{\sqrt{\zeta +{n}_{0}^{2}}}={\int}_{0}^{t}\frac{\mathrm{d}x}{\sqrt{{n}^{2}(x)-\zeta -{n}_{0}^{2}}}\phantom{\rule{0ex}{0ex}}={\int}_{0}^{{n}_{1}^{2}-{n}_{0}^{2}}\frac{g(\eta )\mathrm{d}\eta}{\sqrt{\eta -\zeta}},$$Equation (5) theoretically has an analytic solution,^{8}

## (6)

$$g(\eta )=\frac{-1}{\sqrt{8\pi}{\mathrm{\Gamma}}^{2}\left(\frac{3}{4}\right){\eta}^{\frac{1}{4}}}\frac{d}{d\eta}\left[{\int}_{\eta}^{{n}_{1}^{2}-{n}_{0}^{2}}\frac{\mathrm{d}v}{\sqrt{v-\eta}}{\int}_{0}^{v}\frac{f(\zeta )\mathrm{d}\zeta}{{\zeta}^{\frac{1}{4}}{(v-\zeta )}^{\frac{1}{4}}}\right].$$Unfortunately, note that $\zeta <0$ and the innermost integral must take positive values from 0 to ${n}_{1}^{2}-{n}_{0}^{2}$. Alternatively, for inverting a finite set of measurements ($m=1\dots M$), a piecewise construction of the index using any explicitly integrable functional form for the segments would yield the set of $M$ nonlinear equations

## (7)

$$f({\zeta}_{m})=\sum _{k=1}^{K}{\int}_{\frac{(k-1)}{K}t}^{\frac{k}{K}t}\frac{\mathrm{d}x}{\sqrt{{n}_{k}^{2}(x)-{n}_{0}^{2}-{\zeta}_{m}}}.$$By choosing the number of segments such that the total number of unknowns is equal to or less than the number of measurements, the system could be solved using standard methods. This second method closely resembles the first, although in the first method, numerical integration and any reasonably parameterized form of index profile could be used. With the first method, care must be taken to ensure a well-posed problem because an index profile form with more parameters than the number of measurements would be under-constrained and therefore not have a unique solution.

## 3.

## Axial Gradient Along a Slab or Cross-Section of a Radial Gradient

For this case (see Fig. 2), Eq. (1) reduces to

Examination of Eq. (8) reveals that the ray path is dependent only on the local logarithmic derivative of the index of refraction. A boundary condition on the index must be known to solve for the absolute index of refraction. Integration of Eq. (8) yields

Therefore, the measurement of the entrance and exit slopes of the beam directly gives the ratio of indices at the beam entry and exit points. The measurement of the slope could, e.g., be done with two partly scattering screens at different known distances from the beam exit. Equation (9) expresses the same invariant as Eq. (3) in a rotated reference frame, where the iso-indicial surfaces are now normal to the $y$-axis. In principle, Eq. (9) with data from a sweep over the full height of the slab would enable what amounts to numerical integration of the index (assuming as a boundary condition a region of known index).

Rearranging Eq. (9) and integrating yields

## (10)

$$t={\int}_{{y}_{0}}^{{y}_{1}}\frac{\mathrm{d}y}{\sqrt{(1+{y}_{0}^{\prime 2})\frac{{n}^{2}(y)}{{n}_{0}^{2}}-1}},$$In the ray-optic approximation, mode conversion is like a distributed or convolved version of beam displacement. Instead of a ray, we now have an intensity distribution at each face of the slab, ${I}_{0}(y)$, ${I}_{1}(y)$. For the axial case (the cylindrical case is more complicated due to off-axis rays from finite beam width and is beyond the scope of the current work), assuming the index of refraction and intensity distribution are continuous and smooth and that the incident beam has parallel rays, we have

where $d({y}_{0})\equiv {y}_{1}\u2013{y}_{0}$ is the displacement of a ray and depends on the index distribution and the incidence angle. Equation (11) becomes singular if there is a focal point at ${y}_{1}$. In this manner, intensity measurement is a proxy for displacement measurement, only making displacement measurement more complicated. However, it may be that either the beams used are inherently wide compared with the measurement precision, meaning one could not escape making mode-shape measurements, or that due to precision differences in displacement measurement and intensity measurement, there would still be an advantage to the latter. In any case, Eq. (11) could be useful to predict mode conversion effects based on displacement. Making multiple measurements per beam position (mode characteristics, displacement, and angular deflection) could allow some reduction in measurement-driven errors. For normally incident Gaussian beams with small displacements (single parameter fit to index distribution over any given beam path), published calculations show detailed wave-optic predictions of mode conversion.^{6}

## 3.1.

### Index Profile Computation

## 3.1.1.

#### Deflection-only method

It may be useful to reframe Eq. (9) in its analytically distilled form to understand how it might lead to computation of index profile. Equation (9) evaluated at the exit side of the slab

has measured quantities on the left-hand side and values of the unknown index profile on the right-hand side. Taking the natural logarithm of both sides and calling and yieldsPresuming that $(d{y}_{1}/d{y}_{0})>0$ (which is satisfied if there is no ray crossing the medium) and the index is known either along an interval [${y}_{0}$, ${y}_{0}+d({y}_{0})$] or at a point where $(dn/dy)=0$, the index could be calculated everywhere [calculating each successive value of $\zeta $ from previously known values, where the value in Eq. (13) that would be previously known depends on the sign of $d$]. Equation (13) could be approached analytically (see Appendix), but, in practice, the progressive sweep described here or any finite element method suited to linear relationships like Eq. (13) is more likely to be useful.

## 3.1.2.

#### Displacement only and mixed methods

Equations (10) and (12), together with measurements of $d$ and ${y}_{1}^{\prime}$, form a basis for index determination. Because Eq. (10) is not analytically integrable for all $n(y)$, index profile estimation requires some numerical approximation. The approaches to approximating solutions to Eq. (10) can take the form of assuming properties, such as piecewise linearity, of the index profile or properties of the ray paths. These approaches are joined by assumptions about small displacements relative to the index gradient or thin slabs.

**Approximation A:** For the special case where the incident ray is almost normal to the slab surface, it is useful to rewrite Eq. (10) as

## (14)

$$t={\int}_{0}^{d({y}_{0})}\frac{\mathrm{d}\mathrm{\Delta}y}{\sqrt{[\sqrt{1+{y}_{0}^{\prime 2}}\frac{n(y)}{{n}_{0}}-1][\sqrt{1+{y}_{0}^{\prime 2}}\frac{n(y)}{{n}_{0}}+1]}}.$$For $[n(y)-{n}_{0}]/{n}_{0}\equiv \tilde{n}(y)\ll 1$ and ${y}_{0}^{\prime}\ll 1$, Eq. (10) can be approximated by

## (15)

$$t\approx {\int}_{0}^{d({y}_{0})}\frac{\mathrm{d}\mathrm{\Delta}y}{\sqrt{{y}_{0}^{\prime 2}+2\tilde{n}(y)}}.$$### Approach 1: Linear index over the displacement of any given ray

If fractional index changes are small over length scales of the order of the slab thickness, rays deflect only a small amount and therefore traverse indices of the order of the index at the entry point. Assume the index profile encountered by any given ray can be well approximated by a linear function:

Plugging this form into Eq. (10) yields

## (17)

$$t={\int}_{0}^{d({y}_{0})}\frac{\mathrm{d}\mathrm{\Delta}y}{\sqrt{(1+{y}_{0}^{\prime 2})({a}^{2}\mathrm{\Delta}{y}^{2}+2a\mathrm{\Delta}y+1)-1}}.$$Integration yields [Wolfram’s online Integrator Engine (http://integrals.wolfram.com/index.jsp)]

## (18)

$$(1+\frac{{y}_{0}^{\prime}}{\sqrt{1+{y}_{0}^{\prime 2}}}){e}^{at\sqrt{1+{y}_{0}^{\prime 2}}}=(1+ad)+\sqrt{{(1+ad)}^{2}-\frac{1}{1+{y}_{0}^{\prime 2}}},$$Equation (17) can be reduced to Eq. (15) under the assumptions of Approximation A, noting that $\tilde{n}=a\mathrm{\Delta}y$ for this case, and then integrated and solved for the local index gradient

As with Eq. (18), this overconstrains the independently measured ray exit slope (unless the true index profile is linear). Under these approximations of a linear index profile with small angle of incidence and small index change over the ray path, the ray path internal to the slab will be given by

so these approximations are equivalent to assuming parabolic ray paths as proposed by Barnard and Ahlborn.^{9}

### Approach 2: Quadratic index over the displacement of any given ray

The quadratic approximation has the advantage that the number of unknown parameters matches the number of measurements:

I am unaware of an explicit solution of Eq. (10) emerging from this form; however, under Approximation A [Eq. (15)] with

## (21)

$$t={\int}_{0}^{d({y}_{0})}\frac{\mathrm{d}\mathrm{\Delta}y}{\sqrt{{y}_{0}^{\prime 2}+2(a\mathrm{\Delta}y+b\mathrm{\Delta}{y}^{2})}},$$## (22)

$$t\sqrt{2b}=\mathrm{ln}\text{\hspace{0.17em}}\frac{a+2bd+\sqrt{2b}\sqrt{{y}_{0}^{\prime 2}+2(ad+b{d}^{2})}}{a+{y}_{0}^{\prime}\sqrt{2b}}.$$Equation (22) could be solved numerically, together with the Approximation A equivalent of Eq. (12),

### Other approximations

While this type of approach could theoretically be extended to higher-order polynomial approximations to the local index profile (representing Taylor series expansions), there would not be enough measurements to evaluate the parameters of such a framework. An alternative would be to look at different expansions of the local index profile, particularly those admitting explicit integration of Eq. (10). A one-parameter example of such would be an exponential profile; however, a two parameter alternative would allow the exit slope and displacement measurements to be both used without over-constraining the problem. Published work includes a list of explicitly integrable axial index distributions (e.g., Ref. 5).

## 3.2.

### Refraction at Interfaces

This analysis has so far given results in terms of the ray slopes just inside the surface. Discrete refraction at the surface in accordance with Snell’s law,

## (24)

$${y}_{i}^{\prime}=\frac{\frac{{n}_{a}}{{n}_{i}}{\hat{y}}_{1}^{\prime}}{\sqrt{1+(1-\frac{{n}_{a}^{2}}{{n}_{i}^{2}}){\hat{y}}_{i}^{\prime 2}}},$$## (25)

$${n}_{1}^{2}={n}_{0}^{2}+{n}_{a}^{2}(\frac{{\hat{y}}_{1}^{\prime 2}}{1+{\hat{y}}_{1}^{\prime 2}}-\frac{{\hat{y}}_{0}^{\prime 2}}{1+{\hat{y}}_{0}^{\prime 2}}),$$## (26)

$$t={\int}_{{y}_{0}}^{{y}_{1}}\frac{\mathrm{d}y}{\sqrt{\left(\frac{1+{\hat{y}}_{0}^{\prime 2}}{{n}_{0}^{2}+({n}_{0}^{2}-{n}_{a}^{2}){\hat{y}}_{0}^{\prime 2}}\right){n}^{2}-1}}.$$Note that for normal incidence, there is no refraction at the entrance surface.

## 4.

## Radial or Spherical Gradients

Due to symmetry, the rays in a planar circular slice through a radially symmetric cylinder will remain in the plane of the slice. In a concentric spherical medium (see Fig. 3), ray paths will all remain in the plane they share with the center of the sphere and they will be identical to those in a cylindrical medium with the same radial index dependence.

For this case, Eq. (1) reduces to

## (27)

$$\frac{d\text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}n}{dr}=\frac{r{r}^{\prime \prime}-2{r}^{\prime 2}-{r}^{2}}{r({r}^{2}+{r}^{\prime 2})},$$## (28)

$$\frac{dr}{d\theta}=\pm \sqrt{\frac{{n}^{2}(r){r}^{4}}{{n}_{0}^{2}{R}^{4}}({R}^{2}+{r}_{0}^{\prime 2})-{r}^{2}},$$## (29)

$$\frac{d\widehat{r}}{d\theta}=\pm \sqrt{\frac{{\widehat{n}}^{2}(\widehat{r}){\widehat{r}}^{4}}{{\widehat{n}}_{0}^{2}}(1+{\hat{r}}_{0}^{\prime 2})-{\widehat{r}}^{2}}.$$Hence, for the radially symmetric medium as with the planar media, the local ray slope is dependent only on the boundary condition and the local index of refraction, not the intermediate path. The sign ambiguity emerges from the symmetry of the inbound and outbound paths in the medium.

Previous work on radially symmetric media limited investigation to azimuthally complete media, such as fiber-optic preforms. For manufacture of spherical shell sections, such as hemispherical lens components, exposed surfaces provide an opportunity to use the deflection technique embodied in Eq. (29) to directly measure the refractive index profiles.

For azimuthally complete media (full spheres or cylinders), the only accessible surface has a uniform index of refraction, and therefore nothing can be learned from the angle of the emerging beam relative to the surface (deflection). Displacement in the radial realm measures the difference in azimuth between the entry and exit points of the beam. To integrate Eq. (29) into the displacement relation, one must first identify the innermost transit point of the beam, ${\widehat{r}}^{*}$, where symmetry dictates $d\widehat{r}/d\theta =0$. This satisfies the relation

## (30)

$${\widehat{r}}^{*}=\frac{{\widehat{n}}_{0}}{{\widehat{n}}^{*}}\frac{1}{\sqrt{1+{\widehat{r}}_{0}^{\prime 2}}},$$In these terms, Eqs. (29) and (30) can be written as

## (31)

$$\frac{1}{\widehat{r}}\frac{d\widehat{r}}{d\theta}=\pm \sqrt{{\widehat{n}}^{2}{\widehat{r}}^{2}\text{\hspace{0.17em}}{\mathrm{csc}}^{2}\text{\hspace{0.17em}}\widehat{\varphi}-1}$$## (32)

$${\widehat{n}}^{*}{\widehat{r}}^{*}={\widehat{n}}_{0}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\varphi =\mathrm{sin}\text{\hspace{0.17em}}\widehat{\varphi},$$In terms of this angle, the invariant of Eq. (31) can be written as

## (33)

$${\mathrm{sin}}^{2}\text{\hspace{0.17em}}\widehat{\varphi}={\widehat{n}}^{2}{\widehat{r}}^{2}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\psi ,$$Integration of Eq. (31) from the outer boundary to ${\widehat{r}}^{*}$ and back yields

## (34)

$$\mathrm{\Delta}\theta =\pm 2{\int}_{{\widehat{r}}^{*}}^{1}\frac{\mathrm{d}\widehat{r}}{\widehat{r}\sqrt{{\widehat{n}}^{2}{\widehat{r}}^{2}\text{\hspace{0.17em}}{\mathrm{csc}}^{2}\text{\hspace{0.17em}}\widehat{\varphi}-1}}.$$Index reconstruction proceeds as follows. Let

Suppose $\zeta (\widehat{r})$ is a one-to-one mapping. Substitution of Eq. (35) into Eq. (34) yields

## (36)

$$\mathrm{\Delta}\theta (\widehat{\phi})=\pm 2\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\widehat{\varphi}{\int}_{\mathrm{sin}\text{\hspace{0.17em}}\widehat{\phi}}^{{\widehat{n}}_{0}}\frac{{f}^{\prime}(\zeta )\mathrm{d}\zeta}{\sqrt{{\zeta}^{2}-{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\widehat{\varphi}}},$$In a practical scanning system, a beam parallel to the $x$-axis will be translated along the $y$-axis. Letting the position of the beam be

and considering the negative instance of Eq. (36) (e.g., rays above the $x$-axis traveling in the positive $x$-direction), multiplication and integration in a manner similar to that of an Abel integral transform yields## (39)

$${\int}_{\zeta}^{1}\frac{\mathrm{\Delta}\theta (\widehat{y})\mathrm{d}\widehat{y}}{\sqrt{{\widehat{y}}^{2}-{\zeta}^{2}}}=-{\int}_{\zeta}^{1}{\int}_{\widehat{y}}^{{\widehat{n}}_{0}}\frac{2\widehat{y}{f}^{\prime}(\tau )}{\sqrt{{\widehat{y}}^{2}-{\zeta}^{2}}\sqrt{{\tau}^{2}-{\widehat{y}}^{2}}}\mathrm{d}\tau \mathrm{d}\widehat{y}.$$Switching the order of integration and presuming ${\widehat{n}}_{0}>1$,

Applying the boundary condition $f({\widehat{n}}_{0})=0$ and rearranging yields

## (40)

$$f(\zeta )=\frac{1}{\pi}{\int}_{\zeta}^{1}\frac{\mathrm{\Delta}\theta (\widehat{y})\mathrm{d}\widehat{y}}{\sqrt{{\widehat{y}}^{2}-{\zeta}^{2}}}+\frac{1}{2}f(1)\phantom{\rule{0ex}{0ex}}-\frac{1}{\pi}{\int}_{1}^{{\widehat{n}}_{0}}{f}^{\prime}(\tau ){\mathrm{tan}}^{-1}\left(\frac{{\tau}^{2}+{\zeta}^{2}-2}{2\sqrt{1-{\zeta}^{2}}\sqrt{{\tau}^{2}-1}}\right)\mathrm{d}\tau .$$Equation (41) indicates that from the tangential ray at $\widehat{y}=1$ to the normal ray at $\widehat{y}=0$, each ray penetrates progressively deeper into the medium so that reconstruction of the refractive index profile at position $\zeta =\widehat{y}$ in Eq. (40) depends only on the measured displacement of rays corresponding to greater values of $\widehat{y}$. However, finite surface refraction when the medium has an index of refraction greater than ambient prevents progressive interrogation up to the surface. From Eq. (32), even a ray tangent to the surface ($\mathrm{sin}\text{\hspace{0.17em}}\widehat{\varphi}=1$) penetrates to $\zeta =1$. Thus, for every incident ray, the effect of the medium from the boundary to $\zeta =1$ manifests only in an aggregated fashion. As a result, unless the ambient index of refraction matches or exceeds the surface index of refraction, displacement results will not allow reconstruction of any part of the index profile. Each of the innumerable outer profiles (surface to $\zeta =1$) satisfying the observed aggregate behavior would lead to a different reconstructed inner profile [For the lens design problem, which is the inverse of the metrology problem described herein (given a set of desired ray displacements, what index would be needed), Eq. (40) implies that the refractive index profile from the boundary to $\zeta =1$ could be chosen arbitrarily and the profile from $\zeta =1$ to the center would follow]. For a medium immersed in index matched fluid, the last two terms in Eq. (40) vanish, and the full index profile can be reconstructed.

Returning to azimuthally incomplete media, which would allow beam exit from boundaries of varying refractive index (e.g., a lens cut from a spherically symmetric distribution or a hemisphere or hemispherical shell), Eq. (33) allows direct index reconstruction. If the boundary of the medium is along a radial line (normal to the iso-indicial surfaces), the external exit angle relative to the surface normal $\widehat{\alpha}$ gives the following relationship for deriving the index at the exit point:

## (42)

$$\widehat{n}=\sqrt{{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\widehat{\alpha}+\frac{{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\widehat{\varphi}}{{\widehat{r}}^{2}}}.$$This would allow the nondestructive evaluation by the deflection method of hemispherical lens preforms or the like.

Mode conversion for the spherical case is more complicated because finite beam width leads to rays in different planes of refractive index symmetry. Treatment of such is beyond the scope of this work.

## 5.

## Accuracy

Any implementation of the methods described here will have two sources of error: projection error and measurement error. Relying on a finite number of discrete measurements, the methods will produce reconstructed index distributions characterized by a finite number of parameters. Calculation of the index distribution from the measurements amounts to a projection of the true distribution onto the more limited function space of the finite parameter representation. The first source of error is the difference between the true distribution and its projection, and, for each method described here, the error will depend on the set of independent measurements taken, the form of the true distribution, and the finite parameter representation. For instance, for a continuous distribution well sampled by the measurements and a piecewise-linear representation, the projection error should be

where $\delta $ is the size of the linear segments and ${n}^{\prime \prime}$ is the second derivative of the index profile.Measurement error comprises random measurement error reducible through redundant measurement as well as systematic measurement error. Redundant measurement could take the form of displacement and deflection techniques used together. The magnitude of the measurement errors will depend on the particular measurement technique and setup, but the sensitivity of the reconstructed profile to these errors can at least be estimated from the equations used for reconstruction. However, the integral equations do not readily lend themselves to analytic evaluation of the sensitivity derivatives, so the sensitivity would, in general, need to be evaluated numerically. Sensitivity will also depend on the index profile. These complexities make it difficult to generically characterize the sensitivity of these methods. For a specific example, Lin et al. examine sensitivity of an implementation of the deflection method for axial gradients along a slab.^{11}

## 6.

## Conclusion

For media with axial, cylindrical, or spherical distributions of refractive index, a ray slope invariant [Eqs. (3), (9), and (33)] allows the use of ray deflection as a means of measuring the index of refraction at the exit point as a function of the index of refraction at the entrance. This means of measurement is only useful when surfaces of varying index are exposed, in which case approximate reconstructions of the index profile are possible with a finite number of measurements. Ray displacement and beam mode conversion also capture information about the index distribution and can be used to reconstruct or help reconstruct the index distribution. For azimuthally complete spherical gradients or radial gradients viewed normal to the axis of symmetry, index reconstruction is only possible in an index-matching fluid. For the other cases considered, numerical approaches are required, but index reconstruction is possible. Table 1 summarizes these results. The analysis here provides a survey of the derivations and proposed solution methods for these approaches with an eye toward metrology applications for gradient index lens manufacture.

## Table 1

Summary of relations in uniform terms with ψ indicating the angle between the ray and the index gradient.

Displacement | Deflection | Mode shape | |
---|---|---|---|

Axial across slab | Δy=±∫0tdx/n^2(x)csc2 ψ^0−1 | Invariant | Invariant |

Axial along slab | t=∫y0y1dy/n^2(y)/[n^02−cos2 ψ^0]−1 | n^12=n^02+cos2 ψ^1−cos2 ψ^0 | I0(y0)/I1(y1)=∂y1/∂y0 |

Spherical | Δθ=±2∫r^*1dr^/r^n^2r^2 csc2 ψ^−1, sin2 ψ^=n^(r^*)2r^*2 | sin2 ψ^=n^2r^2 sin2 ψ | Out-of-scope |

## Appendices

## Appendix:

### Solution to Beam Deflection Angle Method

Beginning with Eq. (13), the general problem consists of finding $\zeta (y)$ satisfying

given $g(y)$ and $d(y)$. Equation (44) is like a continuous version of a difference equation. In that vein, the solution can be separated into a homogeneous and a particular solution, where the homogeneous solution satisfiesSuppose there are two solutions to Eq. (44), ${\zeta}_{1}(y)$ and ${\zeta}_{2}(y)$. Taking the difference between the two versions of Eq. (44), one with each solution plugged in,

where $\mathrm{\Delta}\zeta (y)\equiv {\zeta}_{1}(y)-{\zeta}_{2}(y)$. Thus, any solutions to Eq. (44) differ by a homogeneous solution. Therefore, finding any particular solution and adding it to the general homogeneous solution will provide a complete solution to Eq. (44). By inspection, the general solution to Eq. (45) (guaranteeing only continuity) is any periodic function transformed over successive periods by the known function $d(y)$ as long as $(d/dy)[y+d(y)]>0$. Let $h(y)\equiv y+d(y)$. Taking an arbitrary starting point ${y}_{a}$ and defining an arbitrary periodic function $\varphi (y)$ over [${y}_{a}$, ${y}_{a}+d({y}_{a})$],## (46)

$${\zeta}_{h}=\{\begin{array}{cc}\begin{array}{l}\varphi (y)\\ \varphi [{h}^{-1}(y)]\\ \varphi \{{h}^{-1}[{h}^{-1}(y)]\}\\ \vdots \end{array}& \begin{array}{l}{y}_{a}\le y<h({y}_{a})\\ h({y}_{a})\le y<h[h({y}_{a})]\\ \vdots \\ \vdots \end{array}\end{array}.$$Because of the degrees of freedom inherent in the arbitrary periodic function, the boundary condition necessary to determine a solution uniquely must define the value of $\zeta (y)$over a full interval [${y}_{a}$, ${y}_{a}+d({y}_{a})$]. Such a boundary condition would reduce the solution of Eq. (44) immediately to a recursive form like Eq. (46).

A particular solution can be identified using a Green’s function approach. Let

where $\delta (y)$ is the Dirac delta function. LetTo construct ${\zeta}_{p}(y)$, integrate Eq. (47) times the right-hand side of Eq. (44),

## (49)

$${\int}_{-\infty}^{\infty}g(\tilde{y})(u(h(y),\tilde{y})-u(y,\tilde{y}))\mathrm{d}\tilde{y}={\int}_{-\infty}^{\infty}g(\tilde{y})\delta (y-\tilde{y})\mathrm{d}\tilde{y}\phantom{\rule{0ex}{0ex}}=g(y).$$Comparing the terms of Eqs. (49) and (44) reveals that

Taking the Green’s function from Eq. (48),

## (51)

$${\zeta}_{p}(y)={\int}_{-\infty}^{\infty}g(\tilde{y})\sum _{n=1}^{\infty}\pm \delta [{h}^{\mp n}(y)-\tilde{y}]\mathrm{d}\tilde{y}=\pm \sum _{n=1}^{\infty}g[{h}^{\mp n}(y)].$$If either $g[{h}^{n}(y)]$ or $g[{h}^{-n}(y)]$ is appropriately convergent, then a solution has been identified. If, as may often be the case, on at least one end of the slab, the index becomes homogeneous, then

and the convergence of Eq. (51) is possible. The index profile follows fromThe solution method above guarantees that, for smooth $\varphi (y)$, $\zeta (y)$ will be at least as smooth as $g(y)$ within each interval [${h}^{n}({y}_{a})$, ${h}^{n+1}({y}_{a})$]. Guaranteeing smoothness across the interval boundaries would further constrain the solution.

## Acknowledgments

This work was conducted by the Institute for Defense Analyses under contract W91WAW-11-C-0003 for the Defense Advanced Research Projects Agency Manufacturable Gradient Index Optics Program. The views, opinions, and findings should not be construed as representing the official position of either the Department of Defense or the sponsoring organization. Approved for public release, distribution unlimited.

## References

E. W. Marchand, Gradient Index Optics, Academic Press, New York (1978).Google Scholar

O. Wiener, “Darstellung gekrümmter Lichtstrahlen und Verwerthung derselben zur Untersuchung von Diffusion und Wärmeleitung,” Ann. Phys. 285(5), 105–149 (1893).ANPYA20003-3804http://dx.doi.org/10.1002/(ISSN)1521-3889Google Scholar

G. BeliakovD. Y. C. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. I. Planar systems,” Appl. Opt. 36(22), 5303–5309 (1997).APOPAI0003-6935http://dx.doi.org/10.1364/AO.36.005303Google Scholar

G. BeliakovD. Y. C. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. II. Three-dimensional systems with symmetry,” Appl. Opt. 37(22), 5106–5111 (1998).APOPAI0003-6935http://dx.doi.org/10.1364/AO.37.005106Google Scholar

G. BeliakovS. Buckley, “Reconstruction of the refractive index profile of planar waveguides using ray tracing analysis,” in Proc. Optical Networking: Technologies, Traffic Engineering, and Management, Melbourne, Australia (2003).Google Scholar

A. MandelisB. S. H. Royce, “Fundamental-mode laser-beam propagation in optically inhomogeneous electrochemical media with chemical species concentration gradients,” Appl. Opt. 23(17), 2892–2901 (1984).APOPAI0003-6935http://dx.doi.org/10.1364/AO.23.002892Google Scholar

P. L. Chu, “Nondestructive measurement of index profile of an optical-fibre preform,” Electron. Lett. 13(24), 736–738 (1977).ELLEAK0013-5194http://dx.doi.org/10.1049/el:19770520Google Scholar

A. D. Polyanin, Handbook of Integral Equations, CRC Press, Washington, DC (1998).Google Scholar

A. J. BarnardB. Ahlborn, “Measurement of refractive index gradients by deflection of a laser beam,” Am. J. Phys. 43(7), 573–574 (1975).AJPIAS0002-9505http://dx.doi.org/10.1119/1.9769Google Scholar

R. K. Luneburg, Mathematical Theory of Optics, University of California Press, Berkeley (1964).Google Scholar

D. Linet al., “One-dimensional gradient-index metrology based on ray slope measurements using a bootstrap algorithm,” Opt. Eng. 52(11), 112108 (2013).OPEGAR0091-3286http://dx.doi.org/10.1117/1.OE.52.11.112108Google Scholar

## Biography

**Jeremy A. Teichman** received a BS degree in mechanical engineering and in computer science from Yale University in 1996, and SM and PhD degrees in mechanical engineering (fluid and thermal sciences) from the Massachusetts Institute of Technology in 1998 and 2002, respectively. Since 2002, he has been a research staff member at the Institute for Defense Analyses in Alexandria, Virginia. His research interests include statistical and mathematical modeling of physical systems, treatment of uncertainty in measurement systems, fluid dynamics, and heat transfer.