Measurement of gradient index materials by beam deflection, displacement, or mode conversion

Jeremy Teichman

doi:10.1117/1.OE.52.11.112112

30 August 2013 Measurement of gradient index materials by beam deflection, displacement, or mode conversion

Jeremy Teichman

Author Affiliations +

Optical Engineering, Vol. 52, Issue 11, 112112 (August 2013). https://doi.org/10.1117/1.OE.52.11.112112

Abstract

For media with heterogeneous index of refraction, rays follow curved paths. Methods exploiting knowledge of the relationship between refractive index distribution and ray paths to estimate the refractive index distribution based on laser beam deflection, displacement, and mode conversion on passage through the media are described. The work covers axial, radial (cylindrical), and spherical refractive index distributions.

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

Eq. (1)

- \frac{n y^{''}}{1 + y^{' 2}} = \frac{\partial n}{\partial x} y^{'} - \frac{\partial n}{\partial y},

where

n

is the index of refraction, the ray path is defined by

y = y (x)

, and primes indicate differentiation with respect to

x

(e.g., Ref. 1). For axial, radial, and spherical gradients of refractive index, Eq. (1) reduces to an ordinary differential equation. Integration of the resulting equation results in a relationship between the refractive index distribution and the ray paths that can be exploited to estimate the refractive index distribution based on measurements of laser beams passing through the media. The methods presented here do not rely on measurement of the full path of the beams through the media. They use only the beam conditions at the entry and exit points measurable from outside the media. Beam position, direction, and mode shape are considered. If, due to scatter in the media, the entire beam path can be observed and measured, the index all along the beam path relative to that at beam entry can be directly calculated using the local slope of the ray path.

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,² includes examining index reconstruction

• by known ray paths through the full span of a planar or arbitrary medium³^,⁴;
• for planar wave-guides with fixed index faces⁵;
• for thin medium, small displacement, normal beam incidence⁶; and
• for cylindrical fiber preforms.⁷

2. Axial Gradient Across a Slab

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

Eq. (2)

- \frac{y^{''}}{y^{'} (1 + y^{' 2})} = \frac{1}{n} \frac{d n}{d x},

which can be integrated to produce

Eq. (3)

\frac{1 + \frac{1}{y^{' 2}}}{1 + \frac{1}{y_{0}^{' 2}}} = \frac{n^{2} (x)}{n_{0}^{2}},

or using Snell’s law

\frac{1 + \frac{1}{y^{' 2}}}{1 + \frac{1}{{\hat{y}}_{0}^{' 2}}} = \frac{n^{2} (x)}{n_{a}^{2}},

where

n_{0}

is the index of refraction on the entrance face,

n_{a}

is the ambient index, and

y_{0}^{'}

and

{\hat{y}}_{0}^{'}

are the internal and external entrance slopes, respectively. [This paper adopts the notation that external slopes (related by Snell’s law to internal slopes) are represented by a “^”.] Equation (3) indicates that the local slope of a ray path with given entrance conditions will depend only on the local index of refraction and the index and slope at the entry face. Equation (3) is equivalent to Snell’s law invariant along any ray,

n^{2} \sin^{2} θ = constant

, where

θ

is the angle between the ray and the

x

-axis and where the surface normal in Snell’s law has been replaced by the refractive index gradient, which is normal to the iso-indicial surfaces. Measurement of the exit slope would allow direct calculation of the index on the exit face. Further integration of Eq. (3) across the thickness of the slab yields

Eq. (4)

Δ y \equiv y_{1} - y_{0} = \pm \int_{0}^{t} \frac{d x}{\sqrt{\frac{n^{2} (x)}{n_{0}^{2}} (1 + \frac{1}{y_{0}^{' 2}}) - 1}} .

Fig. 1

Horizontal axial gradient schematic.

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}^{'} (y_{0}^{'})$ depends only on the exit index, not on the interior variation. Only the variation of displacement $Δ 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}^{'} \geq 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}$ ):

ζ \equiv \frac{- n_{0}^{2}}{1 + y_{0}^{' 2}}

η \equiv n^{2} (x) - n_{0}^{2}

Eq. (5)

f (ζ) \equiv \frac{Δ y [y_{0}^{'} (ζ)]}{\sqrt{ζ + n_{0}^{2}}} = \int_{0}^{t} \frac{d x}{\sqrt{n^{2} (x) - ζ - n_{0}^{2}}} = \int_{0}^{n_{1}^{2} - n_{0}^{2}} \frac{g (η) d η}{\sqrt{η - ζ}},

where

g (η) \equiv \frac{d x}{d η} .

n (x)

can be reconstructed if

g (η)

is known by inverting

x (n) = \int_{0}^{n} 2 n g (n^{2} - n_{0}^{2}) d n .

Equation (5) theoretically has an analytic solution,⁸

Eq. (6)

g (η) = \frac{- 1}{\sqrt{8 π} Γ^{2} (\frac{3}{4}) η^{\frac{1}{4}}} \frac{d}{d η} [\int_{η}^{n_{1}^{2} - n_{0}^{2}} \frac{d v}{\sqrt{v - η}} \int_{0}^{v} \frac{f (ζ) d ζ}{ζ^{\frac{1}{4}} {(v - ζ)}^{\frac{1}{4}}}] .

Unfortunately, note that $ζ < 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

Eq. (7)

f (ζ_{m}) = \sum_{k = 1}^{K} \int_{\frac{(k - 1)}{K} t}^{\frac{k}{K} t} \frac{d x}{\sqrt{n_{k}^{2} (x) - n_{0}^{2} - ζ_{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

Eq. (8)

n y^{''} = \frac{d n}{d y} (1 + y^{' 2}) .

Fig. 2

Vertical axial gradient schematic.

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

Eq. (9)

\frac{1 + y^{' 2}}{1 + y_{0}^{' 2}} = \frac{n^{2} (y)}{n_{0}^{2}} .

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

Eq. (10)

t = \int_{y_{0}}^{y_{1}} \frac{d y}{\sqrt{(1 + y_{0}^{' 2}) \frac{n^{2} (y)}{n_{0}^{2}} - 1}},

presuming the rays are monotonic in

y

and, without loss of generality, that

y_{1} > y_{0}

. Measurement of the displacement (

y_{0}

and

y_{1}

) gives a second relationship with information about an integrated function of index along the ray path. Measurement of exit angle and displacement—two scalar quantities—does not provide enough information to determine the function

n (y)

along the ray path between

y_{0}

and

y_{1}

but could be used to evaluate two parameters of a local approximation to an index profile.

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

Eq. (11)

I_{1} [y_{0} + d (y_{0})] = \frac{I_{0} (y_{0})}{1 + d^{'} (y_{0})},

where

d (y_{0}) \equiv y_{1} - 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.⁶

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

Eq. (12)

\frac{1 + y_{1}^{' 2}}{1 + y_{0}^{' 2}} = \frac{n_{1}^{2}}{n_{0}^{2}}

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

g (y_{0}) \equiv \ln (\frac{1 + y_{1}^{' 2}}{1 + y_{0}^{' 2}})

and

ζ (y) \equiv 2 \ln n (y)

yields

Eq. (13)

g (y_{0}) = ζ [y_{0} + d (y_{0})] - ζ (y_{0}) .

Presuming 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 $(d n / d y) = 0$ , the index could be calculated everywhere [calculating each successive value of $ζ$ 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}^{'}$ , 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

Eq. (14)

t = \int_{0}^{d (y_{0})} \frac{d Δ y}{\sqrt{[\sqrt{1 + y_{0}^{' 2}} \frac{n (y)}{n_{0}} - 1] [\sqrt{1 + y_{0}^{' 2}} \frac{n (y)}{n_{0}} + 1]}} .

For $[n (y) - n_{0}] / n_{0} \equiv \tilde{n} (y) ≪ 1$ and $y_{0}^{'} ≪ 1$ , Eq. (10) can be approximated by

Eq. (15)

t \approx \int_{0}^{d (y_{0})} \frac{d Δ y}{\sqrt{y_{0}^{' 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:

Eq. (16)

n (y) = a n_{0} (y - y_{0}) + n_{0} .

Plugging this form into Eq. (10) yields

Eq. (17)

t = \int_{0}^{d (y_{0})} \frac{d Δ y}{\sqrt{(1 + y_{0}^{' 2}) (a^{2} Δ y^{2} + 2 a Δ y + 1) - 1}} .

Integration yields [Wolfram’s online Integrator Engine ( http://integrals.wolfram.com/index.jsp)]

Eq. (18)

(1 + \frac{y_{0}^{'}}{\sqrt{1 + y_{0}^{' 2}}}) e^{a t \sqrt{1 + y_{0}^{' 2}}} = (1 + a d) + \sqrt{{(1 + a d)}^{2} - \frac{1}{1 + y_{0}^{' 2}}},

which cannot be solved explicitly but can readily be solved numerically for

a

. This solution would yield a piecewise linear solution to the index, which would in general not be precisely compliant with Eq. (12).

Equation (17) can be reduced to Eq. (15) under the assumptions of Approximation A, noting that $\tilde{n} = a Δ y$ for this case, and then integrated and solved for the local index gradient

Eq. (19)

a = \frac{2}{t} (\frac{d}{t} - y_{0}^{'}) .

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

y = \frac{1}{2} a x^{2} + y_{0}^{'} x + y_{0},

so these approximations are equivalent to assuming parabolic ray paths as proposed by Barnard and Ahlborn.⁹

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:

Eq. (20)

n (y) = n_{0} a (y - y_{0}) + n_{0} b {(y - y_{0})}^{2} + n_{0} .

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

\tilde{n} = a Δ y + b Δ y^{2},

Eq. (21)

t = \int_{0}^{d (y_{0})} \frac{d Δ y}{\sqrt{y_{0}^{' 2} + 2 (a Δ y + b Δ y^{2})}},

which yields (Wolfram Integrator)

Eq. (22)

t \sqrt{2 b} = \ln \frac{a + 2 b d + \sqrt{2 b} \sqrt{y_{0}^{' 2} + 2 (a d + b d^{2})}}{a + y_{0}^{'} \sqrt{2 b}} .

Equation (22) could be solved numerically, together with the Approximation A equivalent of Eq. (12),

Eq. (23)

y_{1}^{' 2} = y_{0}^{' 2} + 2 (a d + b d^{2}) .

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,

n_{1} \sin θ_{1} = n_{2} \sin θ_{2}

relates the internal entrance and exit slopes,

y_{0}^{'}

and

y_{1}^{'}

to the measurable slopes just outside the slab surface,

{\hat{y}}_{0}^{'}

and

{\hat{y}}_{1}^{'}

. Rewriting Snell’s law in terms of slopes leads to

Eq. (24)

y_{i}^{'} = \frac{\frac{n_{a}}{n_{i}} {\hat{y}}_{1}^{'}}{\sqrt{1 + (1 - \frac{n_{a}^{2}}{n_{i}^{2}}) {\hat{y}}_{i}^{' 2}}},

where

n_{a}

represents the ambient index of refraction and

i

is either 0 or 1. Combining Eqs. (24) and (12) leads to

Eq. (25)

n_{1}^{2} = n_{0}^{2} + n_{a}^{2} (\frac{{\hat{y}}_{1}^{' 2}}{1 + {\hat{y}}_{1}^{' 2}} - \frac{{\hat{y}}_{0}^{' 2}}{1 + {\hat{y}}_{0}^{' 2}}),

which is equivalent to the deflection relation, Eq. (12), directly in terms of measurable slopes. The similar equivalent of Eq. (10), the displacement relation, is

Eq. (26)

t = \int_{y_{0}}^{y_{1}} \frac{d y}{\sqrt{(\frac{1 + {\hat{y}}_{0}^{' 2}}{n_{0}^{2} + (n_{0}^{2} - n_{a}^{2}) {\hat{y}}_{0}^{' 2}}) 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.

Fig. 3

Spherical gradient schematic.

For this case, Eq. (1) reduces to

Eq. (27)

\frac{d \ln n}{d r} = \frac{r r^{''} - 2 r^{' 2} - r^{2}}{r (r^{2} + r^{' 2})},

where primes indicate differentiation with respect to

θ

along a ray. Integration of Eq. (27) yields

Eq. (28)

\frac{d r}{d θ} = \pm \sqrt{\frac{n^{2} (r) r^{4}}{n_{0}^{2} R^{4}} (R^{2} + r_{0}^{' 2}) - r^{2}},

where

R

is the outer radius of the element,

n_{0}

is the index on the outer boundary, and

r_{0}^{'}

is

{d r / d θ |}_{r = R}

. In dimensionless terms, letting

\hat{n} \equiv \frac{n}{n_{a}}, \hat{r} \equiv \frac{r}{R},

where

n_{a}

is the ambient index of refraction,

Eq. (29)

\frac{d \hat{r}}{d θ} = \pm \sqrt{\frac{{\hat{n}}^{2} (\hat{r}) {\hat{r}}^{4}}{{\hat{n}}_{0}^{2}} (1 + {\hat{r}}_{0}^{' 2}) - {\hat{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, ${\hat{r}}^{*}$ , where symmetry dictates $d \hat{r} / d θ = 0$ . This satisfies the relation

Eq. (30)

{\hat{r}}^{*} = \frac{{\hat{n}}_{0}}{{\hat{n}}^{*}} \frac{1}{\sqrt{1 + {\hat{r}}_{0}^{' 2}}},

where

{\hat{n}}^{*}

is the relative index at

{\hat{r}}^{*}

. In terms of angle of incidence

\hat{ϕ}

, where the initial internal angle

ϕ

after refraction is given by Snell’s law,

{\hat{n}}_{0} \sin ϕ = \sin \hat{ϕ}, \cot ϕ = {\frac{d \hat{r}}{d θ} |}_{\hat{r} = 1} .

In these terms, Eqs. (29) and (30) can be written as

Eq. (31)

\frac{1}{\hat{r}} \frac{d \hat{r}}{d θ} = \pm \sqrt{{\hat{n}}^{2} {\hat{r}}^{2} \csc^{2} \hat{ϕ} - 1}

and

Eq. (32)

{\hat{n}}^{*} {\hat{r}}^{*} = {\hat{n}}_{0} \sin ϕ = \sin \hat{ϕ},

indicating that the behavior is scale invariant. The angle of the ray with the local gradient,

ψ

, is given by

\cot ψ = \frac{1}{\hat{r}} \frac{d \hat{r}}{d θ} .

In terms of this angle, the invariant of Eq. (31) can be written as

Eq. (33)

\sin^{2} \hat{ϕ} = {\hat{n}}^{2} {\hat{r}}^{2} \sin^{2} ψ,

equivalent to that in Ref. 10 or 1.

Integration of Eq. (31) from the outer boundary to ${\hat{r}}^{*}$ and back yields

Eq. (34)

Δ θ = \pm 2 \int_{{\hat{r}}^{*}}^{1} \frac{d \hat{r}}{\hat{r} \sqrt{{\hat{n}}^{2} {\hat{r}}^{2} \csc^{2} \hat{ϕ} - 1}} .

Index reconstruction proceeds as follows. Let

Eq. (35)

ζ (\hat{r}) \equiv \hat{n} (\hat{r}) \hat{r} .

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

Eq. (36)

Δ θ (\hat{φ}) = \pm 2 \sin \hat{ϕ} \int_{\sin \hat{φ}}^{{\hat{n}}_{0}} \frac{f^{'} (ζ) d ζ}{\sqrt{ζ^{2} - \sin^{2} \hat{ϕ}}},

where

Eq. (37)

f (ζ) \equiv \ln \hat{r} (ζ) .

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

Eq. (38)

\hat{y} \equiv \frac{y}{R} = \sin \hat{ϕ}

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

Eq. (39)

\int_{ζ}^{1} \frac{Δ θ (\hat{y}) d \hat{y}}{\sqrt{{\hat{y}}^{2} - ζ^{2}}} = - \int_{ζ}^{1} \int_{\hat{y}}^{{\hat{n}}_{0}} \frac{2 \hat{y} f^{'} (τ)}{\sqrt{{\hat{y}}^{2} - ζ^{2}} \sqrt{τ^{2} - {\hat{y}}^{2}}} d τ d \hat{y} .

Switching the order of integration and presuming ${\hat{n}}_{0} > 1$ ,

- \int_{ζ}^{1} \frac{Δ θ (\hat{y}) d \hat{y}}{\sqrt{{\hat{y}}^{2} - ζ^{2}}} = \int_{ζ}^{1} \int_{ζ}^{τ} \frac{2 \hat{y} f^{'} (τ)}{\sqrt{{\hat{y}}^{2} - ζ^{2}} \sqrt{τ^{2} - {\hat{y}}^{2}}} d \hat{y} d τ + \int_{1}^{{\hat{n}}_{0}} \int_{ζ}^{1} \frac{2 \hat{y} f^{'} (τ)}{\sqrt{{\hat{y}}^{2} - ζ^{2}} \sqrt{τ^{2} - {\hat{y}}^{2}}} d \hat{y} d τ = \int_{ζ}^{1} f^{'} (τ) {[- \tan^{- 1} (\frac{τ^{2} + ζ^{2} - 2 {\hat{y}}^{2}}{2 \sqrt{{\hat{y}}^{2} - ζ^{2}} \sqrt{τ^{2} - {\hat{y}}^{2}}})]}_{\hat{y} = ζ}^{τ} d τ + \int_{1}^{{\hat{n}}_{0}} f^{'} (τ) {[- \tan^{- 1} (\frac{τ^{2} + ζ^{2} - 2 {\hat{y}}^{2}}{2 \sqrt{{\hat{y}}^{2} - ζ^{2}} \sqrt{τ^{2} - {\hat{y}}^{2}}})]}_{\hat{y} = ζ}^{1} d τ .

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

Eq. (40)

f (ζ) = \frac{1}{π} \int_{ζ}^{1} \frac{Δ θ (\hat{y}) d \hat{y}}{\sqrt{{\hat{y}}^{2} - ζ^{2}}} + \frac{1}{2} f (1) - \frac{1}{π} \int_{1}^{{\hat{n}}_{0}} f^{'} (τ) \tan^{- 1} (\frac{τ^{2} + ζ^{2} - 2}{2 \sqrt{1 - ζ^{2}} \sqrt{τ^{2} - 1}}) d τ .

From Eqs. (32) and (38),

Eq. (41)

ζ^{*} \equiv {\hat{n}}^{*} {\hat{r}}^{*} = \hat{y} .

Equation (41) indicates that from the tangential ray at $\hat{y} = 1$ to the normal ray at $\hat{y} = 0$ , each ray penetrates progressively deeper into the medium so that reconstruction of the refractive index profile at position $ζ = \hat{y}$ in Eq. (40) depends only on the measured displacement of rays corresponding to greater values of $\hat{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 ( $\sin \hat{ϕ} = 1$ ) penetrates to $ζ = 1$ . Thus, for every incident ray, the effect of the medium from the boundary to $ζ = 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 $ζ = 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 $ζ = 1$ could be chosen arbitrarily and the profile from $ζ = 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 $\hat{α}$ gives the following relationship for deriving the index at the exit point:

Eq. (42)

\hat{n} = \sqrt{\sin^{2} \hat{α} + \frac{\sin^{2} \hat{ϕ}}{{\hat{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

Eq. (43)

ε_{projection} \sim δ^{2} n^{''},

where

δ

is the size of the linear segments and

n^{''}

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

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 = \pm \int_{0}^{t} d x / \sqrt{{\hat{n}}^{2} (x) \csc^{2} {\hat{ψ}}_{0} - 1}$	Invariant	Invariant
Axial along slab	$t = \int_{y_{0}}^{y_{1}} d y / \sqrt{{\hat{n}}^{2} (y) / [{\hat{n}}_{0}^{2} - \cos^{2} {\hat{ψ}}_{0}] - 1}$	${\hat{n}}_{1}^{2} = {\hat{n}}_{0}^{2} + \cos^{2} {\hat{ψ}}_{1} - \cos^{2} {\hat{ψ}}_{0}$	$I_{0} (y_{0}) / I_{1} (y_{1}) = \partial y_{1} / \partial y_{0}$
Spherical	$Δ θ = \pm 2 \int_{{\hat{r}}^{}}^{1} d \hat{r} / \hat{r} \sqrt{{\hat{n}}^{2} {\hat{r}}^{2} \csc^{2} \hat{ψ} - 1}$ , $\sin^{2} \hat{ψ} = \hat{n} {({\hat{r}}^{})}^{2} {\hat{r}}^{* 2}$	$\sin^{2} \hat{ψ} = {\hat{n}}^{2} {\hat{r}}^{2} \sin^{2} ψ$	Out-of-scope

Appendices

Appendix:

Solution to Beam Deflection Angle Method

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

Eq. (44)

ζ [y + d (y)] - ζ (y) = g (y),

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,

ζ (y) = ζ_{h} (y) + ζ_{p} (y),

where the homogeneous solution satisfies

Eq. (45)

ζ_{h} [y + d (y)] - ζ_{h} (y) = 0 .

Suppose there are two solutions to Eq. (44), $ζ_{1} (y)$ and $ζ_{2} (y)$ . Taking the difference between the two versions of Eq. (44), one with each solution plugged in,

Δ ζ [y + d (y)] - Δ ζ (y) = 0,

where

Δ ζ (y) \equiv ζ_{1} (y) - ζ_{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 / d y) [y + d (y)] > 0

. Let

h (y) \equiv y + d (y)

. Taking an arbitrary starting point

y_{a}

and defining an arbitrary periodic function

ϕ (y)

over [

y_{a}

,

y_{a} + d (y_{a})

],

Eq. (46)

ζ_{h} = {\begin{matrix} \begin{array}{l} ϕ (y) \\ ϕ [h^{- 1} (y)] \\ ϕ {h^{- 1} [h^{- 1} (y)]} \\ ⋮ \end{array} & \begin{array}{l} y_{a} \leq y < h (y_{a}) \\ h (y_{a}) \leq y < h [h (y_{a})] \\ ⋮ \\ ⋮ \end{array} \end{matrix} .

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 $ζ (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

Eq. (47)

u [h (y), \tilde{y}] - u (y, \tilde{y}) = δ (y - \tilde{y}),

where

δ (y)

is the Dirac delta function. Let

h^{n} (y) \equiv \binom{\underset{︸}{h (h {h [\dots h (y) \dots]})}}{n nested compositions}

and likewise for

h^{- n} (y)

as

n

nested compositions of

h^{- 1} (y)

. Two particular solutions to Eq. (47), differing, as required, by a function of the form of Eq. (46), are

Eq. (48)

u (y, \tilde{y}) = \sum_{n = 1}^{\infty} \pm δ [h^{\mp n} (y) - \tilde{y}],

both of which are forms of Dirac combs (picket fence) on a half space.

To construct $ζ_{p} (y)$ , integrate Eq. (47) times the right-hand side of Eq. (44),

Eq. (49)

\int_{- \infty}^{\infty} g (\tilde{y}) (u (h (y), \tilde{y}) - u (y, \tilde{y})) d \tilde{y} = \int_{- \infty}^{\infty} g (\tilde{y}) δ (y - \tilde{y}) d \tilde{y} = g (y) .

Comparing the terms of Eqs. (49) and (44) reveals that

Eq. (50)

ζ_{p} (y) = \int_{- \infty}^{\infty} g (\tilde{y}) u (y, \tilde{y}) d \tilde{y} .

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

Eq. (51)

ζ_{p} (y) = \int_{- \infty}^{\infty} g (\tilde{y}) \sum_{n = 1}^{\infty} \pm δ [h^{\mp n} (y) - \tilde{y}] 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

Eq. (52)

\lim_{n \to \infty} g [h^{\pm n} (y)] = 0,

and the convergence of Eq. (51) is possible. The index profile follows from

Eq. (53)

n (y) = e^{ζ (y) / 2} .

The solution method above guarantees that, for smooth $ϕ (y)$ , $ζ (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

1.

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

2.

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

3.

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). http://dx.doi.org/10.1364/AO.36.005303 APOPAI 0003-6935 Google Scholar

4.

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). http://dx.doi.org/10.1364/AO.37.005106 APOPAI 0003-6935 Google Scholar

5.

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, (2003). Google Scholar

6.

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). http://dx.doi.org/10.1364/AO.23.002892 APOPAI 0003-6935 Google Scholar

7.

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

8.

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

9.

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

10.

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

11.

D. Linet al., “One-dimensional gradient-index metrology based on ray slope measurements using a bootstrap algorithm,” Opt. Eng., 52 (11), 112108 (2013). http://dx.doi.org/10.1117/1.OE.52.11.112108 OPEGAR 0091-3286 Google 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.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Jeremy Teichman "Measurement of gradient index materials by beam deflection, displacement, or mode conversion," Optical Engineering 52(11), 112112 (30 August 2013). https://doi.org/10.1117/1.OE.52.11.112112

Published: 30 August 2013

Access the abstract

JOURNAL ARTICLE
8 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY