Semi-analytical finite ray-tracing through the quadratic symmetric GRIN lens

Conor Flynn; Alexander V. Goncharov

doi:10.1364/AO.504305

1. INTRODUCTION

There are many uses for GRIN lenses, which is why many of them are of great interest. GRIN lenses are favored in compact imaging systems over conventional curved lenses [1] due to their additional degrees of freedom. It has long been established that GRIN lenses show great potential in the study of the human eye [2–5]. The internal structure of the crystalline lens can be modeled as a GRIN lens and has been the subject of extensive research [6–10]. It has been found that, for some eyes, there is a balancing effect between the corneal surfaces and the crystalline lens equations; in particular, the spherical aberration is positive for the cornea and negative for the crystalline lens [11–15]. The spherical aberration of the lens is seen when light rays entering the lens at different heights are brought to focus at different points on the optical axis. The modeling of the human eye using lenses has been of interest to scientists for more than a century. This was first proposed by Gullstrand’s “No. 1” model in the early 20th century [16]. This model is composed of lenses of different refractive indices with the lens with the highest refractive index being placed at the core and the lenses with the lowest value being placed at the anterior and posterior of the modeled crystalline lens [17,18].

The propagation of light within the GRIN structure is notoriously difficult to model, in view of continuous refraction. We shall consider the paraxial region of the lens, which is the simplest region to analyze. The optical path length (OPL) from point ${z_1}$ to ${z_2}$ can be described by the following function [19]:

(1)$${\rm OPL} = \int_{{z_1}}^{{z_2}} n(z)\;{\rm d}z,$$

where $z$ is the position along the optical axis. When $n$ is a constant or a polynomial of $z$, as above, the integral can be calculated with ease within the system between the limits ${z_1}$ and ${z_2}$. The OPL along the lens axis is an important characteristic of the lens and can be used as a variable constraint when reconstructing the GRIN profile. The distribution of GRIN in lenses found in fish, octopi, and mammals all feature a similar increase in refractive index toward the center of the lens [7,20,21].

In GRIN lenses, the rate of change of refractive index can vary in different ways such as an axial and radial gradients. An axial gradient varies along the principal axis, $n(z)$, or, alternatively, perpendicular to the optical axis, $n(x)$ or $n(y)$, depending on the rotation of the system [22]. This can be represented by the following equation:

n(z) = {n_{00}} + {n_{01}}z + {n_{02}}{z^2} + \cdots ,

where ${n_{\textit{ij}}}$ represents the GRIN coefficients with $i$ the exponent of the $r$ term and $j$ the exponent of the $z$ term, i.e., the term ${n_{01}}$ is the coefficient for ${r^0}{z^1}$. A radial gradient has the maximum index located at the center of the lens, and it decreases radially in all directions toward the surface of the lens [23]. GRIN media of this form are described by

n(r) = {n_{00}} + {n_{10}}r + {n_{20}}{r^2} + \cdots ,

where $r = \sqrt {{x^2} + {y^2}}$. A spherical GRIN lens is the coupling of the axial and radial gradients:

n(z,r) = {n_{00}} + {n_{01}}z + {n_{02}}{z^2} + {n_{10}}r + {n_{20}}{r^2} + \cdots .

GRIN lenses can be composed of two sections, the anterior and the posterior [24,25], or described as one refractive index profile for the entirety of the lens [26]. The GRIN lens maps out a lens with a maximum value for refractive index at the center of the lens, gradually decreasing as it disseminates radially through the system and passes through the iso-indicial contours until reaching the minimum at the surface of the lens [27].

The type of lens proposed in this project is represented by a function of quadratic nature, but with set values for the GRIN lens coefficients. Values were obtained from Table 4 in [17]. The ${n_3}$ and ${n_4}$ values were chosen from the 20U model and simplified; ${n_1}$ was set equal to ${n_4}$ to model a spherical lens. As this method was tested on a quadratic lens, other terms were set to zero. As ${n_2}$ is small, it was assumed zero for this model. The following equation is used to characterize the lens:

(2)$$n(z,h) = {n_0} - {n_1}{h^2} + {n_3}z - {n_4}{z^2},$$

where ${n_i}$, with $i \in \{1,3,4\}$, represents the GRIN lens coefficients, $z$ is the optical axis, and $h$ is the axis perpendicular to this axis, generally represented as the $y$ axis. From the above equation, the refractive index of any part of the lens can be equated. ${n_0}$ is the minimum value of $n$, which is located along the surface of the lens and ${n_{{\max}}}$ is the maximum, which will be the local maximum in the paraxial region, which is the core of the lens. The model used here can be represented as

(3)$$n(z,h) = 1.37 - 0.01{z^2} + 0.04z - 0.01{h^2}.$$

Ray-tracing software may be incorporated into mapping the propagation of light through a medium described by the above equations. The iso-indicial contours of the lens can be mapped out to show shell like structures within the lens, which is represented in Fig. 1.

Fig. 1. General shape of a GRIN lens displaying an illustration of the iso-indicial contours. The GRIN media is given by the equation $n(z,h) = {n_0} - {n_1}{h^2} + {n_3}z - {n_4}{z^2}$, which is the form of the models implemented in this paper.

Download Full Size | PDF

Many characteristics of the lens can be determined from Eq. (2). The radii of the iso-indicial contours can be derived from the above equations and the sag equation below:

(4)$${h^2} = 2R(z)z - (1 + k){z^2},$$

where $k$ is the conic constant and $R(z)$ is the radius of curvature of the lens surface, or iso-indicial contours, which are dependent on the external surface, at point ($z$, $h$). Equation (2) can be written in the above form to compute the necessary parameters. Optimization techniques can be utilized on Eq. (3) to locate the position of ${n_{{\max}}}$ along the optical axis for a ray entering the lens at a known initial height. $\frac{{dn}}{{dz}}$ is found and equated to zero, and solved for $z$, to find the position of and the value of ${n_{{\max}}}$. This is the peak plane [17]. Utilizing Eq. (2), the length (central thickness) of the system in the paraxial region is between ${z_1} = 0$ and ${z_2} = \frac{{{n_3}}}{{{n_4}}}$. The maximum refractive index along the optical axis is located at $z = \frac{{{n_3}}}{{2{n_4}}}$, leading to ${n_{{\max}}} = {n_0} + \frac{{n_3^2}}{{4{n_4}}}$.

The aim of this paper is to derive a differential equation that can describe the light path in a medium with varying refractive index. Using the paraxial ray-tracing equations, a differential equation can be formed to describe the light trajectory close to the $z$ axis, the paraxial region. The aim is to find a method to extend the validity of this equation for light entering the lens at finite heights, outside the paraxial region.

This can then be compared with an optical design software, such as Zemax OpticStudio, that can be used to model the lens. GRIN5 is a surface type available here that can be utilized for ray-tracing through a medium described by Eq. (2) and can compute the effective focal length (EFL) of the designed lens as a reference.

2. PARAXIAL DIFFERENTIAL EQUATION

A. Deriving an Equation

Ray-tracing permits the study of an optical system by calculating different properties of a lens, such as the EFL or back focal distance (BFD), using geometrical optics [28]. Paraxial ray-tracing involves rays near the optical axis. The paraxial ray-tracing equations can be stated as follows:

(5)$${n_{j + 1}}{u_{j + 1}} = {n_j}{u_j} - {h_j}{F_j},$$

(6)$${F_j} = \frac{{{n_{j + 1}} - {n_j}}}{{{R_j}}},$$

(7)$${h_{j + 1}} = {h_j} + {d_{j + 1}}{u_{j + 1}},$$

where ${n_j}$ is the refractive index of the ${j}$th surface, ${u_j}$ is the convergence angle after the ${j}$th surface, ${h_j}$ is the height that the ray enters the ${(j + 1)}$th surface of the lens, ${F_j}$ is the power of the ${j}$th lens surface, ${R_j}$ is the radius of curvature of the ${j}$th lens surface, or iso-indicial contour, and ${d_j}$ is the distance between the ${j}$th and ${(j + 1)}$th surfaces, and it can be written as $\Delta z$ to show the change along the optical axis.

Due to the continuous nature of the quadratic GRIN profile, a differential equation is more suitable to ray-trace light as it propagates through the GRIN media. Iterative methods are better suited for a large number of iso-indicial contours, such as shell models [29]. They may be a better representation of real GRIN profiles, but a large number of parameters need to be tracked at each shell, such as $n,u,h,F$, and $d$. A continuous GRIN profile is easier to analyze with a differential ray equation. The paraxial ray-tracing equations, Eqs. (5)–(7), are used for such a lens. This set of equations will be utilized to form a differential equation that will track the change of height within a medium with a varying refractive index. Numerical methods are still used for GRIN lenses as there is no general paraxial ray-tracing equation with analytical solutions available [30]. This makes it sufficiently more difficult to express the EFL and BFD and thus, in turn, the power of the lens. One parameter of the lens that effects the refractive power, and in turn the EFL, is the radii of curvature of the iso-indicial contours of the lens [15].

Combining Eqs. (5) and (6), a first-order differential equation can be formed with $n$ and $u$. Using Eq. (7), this can then be changed from a first-order to a second-order differential with respect to $z$,

(8)$$\begin{split}{n_{j + 1}}{u_{j + 1}} & = {n_j}{u_j} - {h_j}{F_j} \\[-3pt] & = {n_j}{u_j} - {h_j}\left(\frac{{{n_{j + 1}} - {n_j}}}{{{R_j}}}\right) \\[-3pt] & = {n_j}{u_j} - {h_j}\frac{{\Delta n}}{{{R_j}}}, \\[-3pt] {n_{j + 1}}{u_{j + 1}} - {n_j}{u_j} & = - {h_j}\frac{{\Delta n}}{{{R_j}}}, \\[-3pt] \Delta (nu) & = - {h_j}\frac{{\Delta n}}{{{R_j}}}, \\[-3pt] \frac{{\Delta (nu)}}{{\Delta z}} & = - \frac{{{h_j}}}{{{R_j}}}\frac{{\Delta n}}{{\Delta z}}, \\[-3pt] \mathop {\lim}\limits_{\Delta z \to 0} \left(\frac{{\Delta (nu)}}{{\Delta z}}\right.& = - \left.\frac{{{h_j}}}{{{R_j}}}\frac{{\Delta n}}{{\Delta z}}\right), \\[-3pt] \frac{d}{{dz}}(nu) & = - \frac{h}{R}\frac{{dn}}{{dz}}, \\[-3pt] n\frac{{du}}{{dz}} + u\frac{{dn}}{{dz}} & = - \frac{h}{R}\frac{{dn}}{{dz}}.\end{split}$$

From Eq. (6), it can also be noted that

\begin{split}F &= \frac{{\Delta n}}{R}, \\ \frac{F}{{\Delta z}} &= \frac{1}{R}\frac{{\Delta n}}{{\Delta z}}, \\ \mathop {\lim}\limits_{\Delta z \to 0} \left(\frac{F}{{\Delta z}}\right. &= \left.\frac{1}{R}\frac{{\Delta n}}{{\Delta z}}\right), \\ \frac{F}{{dz}}& = \frac{1}{R}\frac{{dn}}{{dz}}.\end{split}

From Eq. (7),

(9)$$\begin{split}\frac{{dh}}{{dz}} &= u, \\ \frac{{{d^2}h}}{{d{z^2}}} &= \frac{{du}}{{dz}}.\end{split}$$

Using results from Eqs. (8) and (9), the following equation can be formed:

(10)$$n\frac{{{d^2}h}}{{d{z^2}}} + \frac{{dn}}{{dz}}\frac{{dh}}{{dz}} + \frac{h}{R}\frac{{dn}}{{dz}} = 0.$$

This will be labeled as the “general paraxial equation” and is used to describe the way in which light rays pass through the GRIN lens.

From the above expressions, using Fermat’s principle on a vertical slab with the refractive index varying with height, $n = n(h)$, only, the relationship below can be extrapolated:

(11)$$\frac{{n^\prime}}{R}h = 2{n_1}h,$$

where $n^\prime $ is equal to the first derivative of $n(z)$ with respect to $z$ and ${n_1}$ is the GRIN coefficient used previously, which is the term used for the $h$ squared expression in Eq. (2). Additional details for this can be found in Appendix A.

Using Eqs. (10) and (11), the final form of the paraxial differential equation, the “reduced paraxial equation,” can be found and expressed as

(12)$$n\frac{{{d^2}h}}{{d{z^2}}} + \frac{{dn}}{{dz}}\frac{{dh}}{{dz}} + 2{n_1}h = 0.$$

This equation can be expressed as a form of the ray equation [31] within the paraxial region, which takes the following form:

\frac{d}{{ds}}\left(n(z,h)\frac{{dh}}{{ds}}\right) = \frac{{dn}}{{dh}}.

Equation (12) can also be compared to a rearranged equation proposed by Beliakov and Chan [32] and Liu [1],

\frac{{nh^{\prime \prime}}}{{1 + {{(h^\prime)}^2}}} + \frac{{\partial n}}{{\partial z}}h^\prime - \frac{{\partial n}}{{\partial h}} = 0.

Implementing this equation in the paraxial region will set $h^\prime (z)$ to be small; thus, ${(h^\prime (z))^2} \ll 1$ and, hence, can be neglected.

B. Solving the Equation

Despite $h(z)$ having an analytic solution, it would prove difficult to solve without using numerical techniques and computer software. Even with the complicated analytical form found, it is suitable for ray-tracing software such as Zemax. GRIN5 is available in Zemax OpticStudio optical software, which can be used to test any analytical results. This work will help toward understanding the focusing effect of the crystalline lens of the eye [33,34].

It can be assumed that the value of $h(z)$ is small in order to take the incident rays entering in the paraxial region into account. Wolfram Mathematica 9.0 is utilized to solve this differential equation (no additional packages or libraries are required for the methods outlined here). It can be observed that the power law equation, Eq. (2), does not provide a solution due to the $h{(z)^2}$ term. To counter this, a new expression is proposed. By setting $h(z) = {h_0} - {h^*}(z)$, where ${h_0}$ is the initial height that the ray enters the lens and ${h^*}(z)$ is the height that the ray decreases as it propagates through the system, the solution will now be affected by the initial conditions. The new power law is now represented by

n(z,h(z)) = {n_0} - {n_1}{({h_0} - {h^*}(z))^2} + {n_3}z - {n_4}{z^2}.

For the paraxial region, ${h^*}(z) \ll {h_0}$ so ${h^*}(z)$ can be neglected, and $n$ can be expressed as

(13)$$n(z,h(z)) = {n_0} - {n_1}h_0^2 + {n_3}z - {n_4}{z^2}.$$

Solving Eq. (13) with the above refractive index profile, the solution is

(14)$$\begin{split}{h(z) =}&{{C_1}{P_{\frac{{\sqrt {8{n_1} + {n_4}} - \sqrt {{n_4}}}}{{2\sqrt {{n_4}}}}}}\left({\frac{{2{n_4}z - {n_3}}}{{\sqrt {- 4{n_1}{n_4}h_0^2 + n_3^2 + 4{n_0}{n_4}}}}} \right) }\\[-3pt] &\quad+{ {C_2}{Q_{\frac{{\sqrt {8{n_1} + {n_4}} - \sqrt {{n_4}}}}{{2\sqrt {{n_4}}}}}}\!\left({\frac{{2{n_4}z - {n_3}}}{{\sqrt {- 4{n_1}{n_4}h_0^2 + n_3^2 + 4{n_0}{n_4}}}}} \right)\!,}\end{split}$$

where ${P_n}(x)$ and ${Q_n}(x)$ are Legendre polynomials. This expression will be labeled as the “Legendre paraxial solution.” The constant terms ${C_1}$ and ${C_2}$ can be found through examination of the initial conditions. Setting the initial height to the coordinates at which the rays enter the lens, $({z_0},h({z_0}))$, to a known value, say ${h_0}$, and the derivative of the function of height, $h^\prime (z)$, which is the expression for the angle the trajectory of light is traveling with respect to the optical axis, gives the initial conditions,

(15)$$\begin{split}h({z_0}) &= {h_0},\\[-3pt]h^\prime ({z_0}) &= {u_0},\end{split}$$

where ${u_0}$ is the angle of incidence at ${z_0}$. A plot of the light path through the paraxial region of a GRIN lens is shown below in Fig. 2. A non-zero initial height is required to calculate EFL by the following equation:

{\rm EFL} = \frac{{- {h_0}}}{{{n_0}h^\prime ({z_2})}},

where ${z_2}$ is the point along the optical axis where the light departs from the GRIN system.

Fig. 2. Plot showing the change of height of the light as it propagates through the lens modeled by Eq. (3) with an initial height, ${h_0} = 0.0001$, within the paraxial region and moving parallel to the optical axis, ${u_0} = 0$, prior to entering the GRIN lens. This plot is given by Eq. (14) after solving the initial conditions in Eq. (15).

Download Full Size | PDF

3. FINITE DIFFERENTIAL EQUATION

A. Deriving a Non-Paraxial Equation

The differential equation will be derived from the ray equation as opposed to the paraxial ray-tracing equations. The equation that shows the path light travels can be used to derive the desired equation as follows:

{(ds)^2} = {(dx)^2} + {(dy)^2} + {(dz)^2},

but using only two of the three planes by setting $dx = 0$ and renaming the vertical axis as $y = h(z)$, it can be stated that

ds = \sqrt {1 + {{(h^\prime (z))}^2}} dz,

where $h^\prime (z) = \frac{{dh(z)}}{{dz}}$. Multiplying both sides by $n$ and taking the integral of each side, Fermat’s principle can be seen on the left hand side of

\int_P^Q n{\rm d}s = \int_{{z_P}}^{{z_Q}} n(z,h(z))\sqrt {1 + {{(h^\prime (z))}^2}} {\rm d}z,

where ${z_P}$ and ${z_Q}$ are the $z$ COORDINATES at points $P$ and $Q$, respectively.

Using Fermat’s Principle, the integrand of the right hand side [35,36],

f(z,h(z),h^\prime (z)) = n(z,h(z))\sqrt {1 + {{(h^\prime (z))}^2}} ,

can be used with the Euler equation [35,36],

\frac{d}{{dz}}({f_{h^\prime (z)}}) = {f_{h(z)}},

where ${f_x}$ represents the derivative $\frac{{df}}{{dx}}$, such that the above identity can be written as

\frac{d}{{dz}}\left(\frac{{df}}{{dh^\prime (z)}}\right) = \frac{{df}}{{dh(z)}}.

Calculating ${f_{h^\prime (z)}}$ and ${f_{h(z)}}$, the following expressions are obtained:

\begin{split}{f_{h^\prime (z)}} &= \frac{{n(z,h(z))h^\prime (z)}}{{\sqrt {1 + {{(h^\prime (z))}^2}}}}, \\ {f_{h(z)}} &= \frac{{dn}}{{dh(z)}}\frac{{ds}}{{dz}}.\end{split}

This can then be shown to derive the following equation:

\frac{{n\frac{{{d^2}h(z)}}{{d{z^2}}}}}{{\sqrt {1 + {{(h^\prime (z))}^2}}}} + \frac{{dn}}{{dz}}\frac{{dh(z)}}{{dz}} - \frac{{dn}}{{dh(z)}} = 0,

which can be reduced, using a Taylor series, to

(16)$$\frac{{n\frac{{{d^2}h(z)}}{{d{z^2}}}}}{{1 + \frac{1}{2}{{(h^\prime (z))}^2}}} + \frac{{dn}}{{dz}}\frac{{dh(z)}}{{dz}} + 2{n_1}h(z) = 0.$$

This can be referred to as the “non-paraxial equation.”

B. Frobenius Method

A power series can be used to derive a recurrence relation in order to solve the differential equation. This will allow one to use a power series to map the propagation of light through GRIN media without the need of Legendre polynomials, as in Eq. (14). The Frobenius method is used to solve the non-paraxial equation, Eq. (16). Expanding the function $h(z)$ as a summation,

h(z) = \sum\limits_{m = 0}^\infty {a_m}{z^{m + r}},

with the polynomial coefficients represented by ${a_m}$, and applying this to the first and second derivatives,

\begin{split}h^\prime (z) &= \sum\limits_{m = 0}^\infty (m + r){a_m}{z^{m + r - 1}}, \\ h^{\prime \prime} (z) &= \sum\limits_{m = 0}^\infty (m + r - 1)(m + r){a_m}{z^{m + r - 2}}.\end{split}

In order to reduce this equation further, the expression $h^\prime (z)$ in the denominator of the first term can be represented by $h^\prime (z) = \beta + \alpha z$, where $\beta$ is the intercept on the $h$ axis and $\alpha$ is the slope of the line $h^\prime (z)$. This will allow one to map the rate of change of height of the ray as a linear line. Higher order terms can be included if required for lenses with higher power. These terms are calculated by using Eq. (14) at finite heights to obtain the expression of a linear line to represent the slope of height with respect to the optical axis.

Solving the recursion relation for $r$ leads to two values, $r = 0$ and $r = 1$, leading to the expression that will be labeled the “Frobenius solution,”

(17)$$h(z) = \sum\limits_{m = 0}^\infty {a_m}{z^m} + \sum\limits_{m = 0}^\infty {b_m}{z^{m + 1}}.$$

Expanding ${a_m}$ and ${b_m}$ in terms of ${a_0}$ and ${b_0}$, respectively, which are analogs to the constants ${C_1}$ and ${C_2}$, was utilized in the paraxial approach, and they are derived from the initial conditions. The recurrence relations for ${a_m}$ and ${b_m}$ can be described as follows:

\begin{split}{a_0} &= {a_0},\quad {a_1} = 0,\quad {a_2} = - \frac{{{n_1}\left({{\beta ^2} + 1} \right)}}{{{n_0} - h_0^2{n_1}}}{a_0}, \\ {a_3}& = \frac{{- 1}}{{3\left({{n_0} - h_0^2{n_1}} \right)}}({n_3}({\beta ^2} + 2){a_2} + 2{n_1}\alpha \beta {a_0}), \\ {a_m} &= \frac{{- 1}}{{(m - 1)m\left({{n_0} - h_0^2{n_1}} \right)}} \\ &\times({a_{m - 1}}((m - 1){n_3}((m - 2) + ({\beta ^2} + 1))) \\ &\quad+ {a_{m - 2}}((m - 2)\left({2{n_3}\alpha \beta - 2{n_4}({\beta ^2} + 1)} \right) \\ &\quad- (m - 2)(m - 3){n_4} + 2{n_1}\left({{\beta ^2} + 1} \right)\!\big) \\ &\quad+ {a_{m - 3}}((m - 3)\left({{n_3}{\alpha ^2} - 4{n_4}\alpha \beta} \right) + 4{n_1}\alpha \beta) \\&\quad + 2{a_{m - 4}}({n_1}{\alpha ^2} - (m - 4){n_4}{\alpha ^2})),\end{split}

\begin{split}{b_0} &= {b_0},\quad {b_1} = - \frac{{{n_3}\left({{\beta ^2} + 1} \right)}}{{2\left({{n_0} - h_0^2{n_1}} \right)}}{b_0}, \\ {b_2} &= - \frac{{{n_3}\left({{\beta ^2} + 2} \right){b_1} + (({\beta ^2} + 1)({n_1} - {n_4}) + {n_3}\alpha \beta){b_0}}}{{3\left({{n_0} - h_0^2{n_1}} \right)}}, \\ {b_3}& = \frac{{- 1}}{{12\left({{n_0} - h_0^2{n_1}} \right)}}(3{n_3}({\beta ^2} + 3){b_2} + \left({4\alpha \beta ({n_1} - {n_4}) + {n_3}{\alpha ^2}} \right){b_0} \\ &\quad+2\left({- {n_4} + 2\alpha \beta {n_3} + ({\beta ^2} + 1)(n1 - 2{n_4})} \right){b_1}), \\ {b_m} &= \frac{{- 1}}{{(m + 1)m\left({{n_0} - h_0^2{n_1}} \right)}} ({b_{m - 1}}(m{n_3}(m + {\beta ^2})) \\&\quad + {b_{m - 2}}((m - 1)\left({2{n_3}\alpha \beta - 2{n_4}({\beta ^2} + 1)} \right) - (m - 2)(m - 1){n_4} \\ &\quad+2{n_1}({\beta ^2} + 1)) + {b_{m - 3}}((m - 2)({n_3}{\beta ^2} - 4{n_4}\alpha \beta) + 4{n_1}\alpha \beta) \\ &\quad+ {b_{m - 4}}(2{\beta ^2}({{n_1} - (m - 3){n_4}} ))).\end{split}

Making use of the same initial conditions as stated in Eq. (15), Eq. (17) can be expressed as

(18)$$h(z) = \sum\limits_{m = 0}^\infty {a_m}{z^m},$$

for rays entering the lens parallel to the optical axis within the paraxial region, otherwise written as ${u_0} = 0$, namely $h^\prime (z) = 0$, leading to ${b_0} = 0$ and thus ${b_m} = 0$ for all values of $m$.

Moore made use of a single recursion relationship for the paraxial region [37], and as seen in Eq. (18), this is sufficient for the paraxial region. By introducing a second summation, this allows ray-tracing for finite-heights. Terms ${a_m}$ and ${b_m}$ are also affected by the initial height of the rays to allow for study of non-paraxial light propagation. A total of 61 terms are used for each of the summations stated above for Eq. (17) to test these methods.

An approximate expression for the EFL can be found in Appendix B.

The exit coordinate ${z_2}$ is required for calculating the BFDs of the optical system at varying initial heights. The point at which the light ray emerges from the lens is computed from

h(z) = h(d) + \frac{{\tan (u(z))d}}{2} - \frac{{\tan (u(z)){z^2}}}{{2d}},

where $d$ is the thickness of the lens and $\tan (u(z))$ can be equated to the negative of the rate of change of height with respect to $z$, with $u(z)$ representing the angle of incidence at point $z$. The above expression can be equated to

h(z) = \sqrt {\frac{{- {n_4}{z^2} + {n_3}z}}{{{n_1}}}} ,

in order to locate the exit position of the light ray.

Fig. 3. EFL (blue) and BFD (red) calculated using Eq. (17) up to ${a_m}$ for the refractive index profile $n(z,h) = 1.37 - 0.01{h^2} + 0.04z - 0.01{z^2}$, where $m$ is the number of terms used to generate the power series.

Download Full Size | PDF

Fig. 4. Shape of a quadratic elliptical GRIN lens displaying an illustration of the iso-indicial contours.

Download Full Size | PDF

4. ELLIPTICAL LENS

The methods above are tested on a quadratic lens with the identical values for ${n_0}$ and ${n_{{\max}}}$, while being represented by a lens with the same semi-diameter, but half the thickness. By making this change, it is ensured that the lens is no longer spherical and thus these results can be evaluated on an aspheric GRIN lens, testing the lens for a different conic constant. These lens characteristics can be altered by manipulating the values of ${n_1}$, ${n_3}$, and ${n_4}$ with the use of the following geometry-based formulas to achieve the desired lens:

\begin{split}d &= \frac{{{n_3}}}{{{n_4}}}, \\ D &= \sqrt {\frac{{n_3^2}}{{{n_1}{n_4}}}} , \\ {n_{{\max}}} & = {n_0} + \frac{{n_3^2}}{{4{n_4}}},\end{split}

where $d$ is the thickness of the lens along the optical axis and $D$ is the diameter of the lens. Solving the above parameters when $d = 2$, $D = 4$, ${n_{{\max}}} = 1.41$, and ${n_0} = 1.37$, the desired lens can be represented by the following refractive profile:

n(z,h) = 1.37 - 0.01{h^2} + 0.08z - 0.04{z^2}.

, (The plot for this is shown in Fig. 4.)

The radius and the conics of this lens are given by

R = \frac{{{n_3}}}{{2{n_1}}}, \\ k = \frac{{{n_4}}}{{{n_1}}} - 1.

The above relationships are obtained for the normalized power law for the external surface, ${h^2} = \frac{{n3}}{{n1}}z - \frac{{n4}}{{n1}}{z^2}$, and compared to Eq. (4). This quadratic lens can now be generated on Zemax OpticStudio using GRIN5 to test the methods on an elliptical lens.

As this lens has the same ${n_0}$ and ${n_{{\max}}}$, but with only half the thickness, higher order terms are needed to map the trajectory for $h^\prime (z)$. It is proposed that the polynomial $h^\prime (z) = \beta + \alpha z + \gamma {z^3}$ can be used to obtain the power series for the elliptical lens. The Frobenius method is carried out again with a distinct set ${a_m}$ and ${b_m}$ terms required to deal with the now cubic nature of $h^\prime (z)$.

5. RESULTS

The EFL is calculated using the Legendre and Frobenius solutions outlined above. The quadratic GRIN model, Eq. (3), is simulated using GRIN5 in Zemax OpticStudio and is used as a comparison, where the step size $\Delta t$ is set as $1.0 \times {10^{- 5}}$. These results can be seen in Table 1. These values are calculated within the paraxial region of the lens. An initial height of 0.0001 mm was proposed for this and from tests carried out as it is an arbitrarily small, finite number. Results are accurate up to eight significant figures, in order to show the difference in nm.

Table 1. EFL Calculated for Each Method for the Lens Represented by $n(z,h) = 1.37 - 0.01{h^2} \;+\; 0.04z - 0.01{z^2}$ for the Spherical and $n(z,h) \;=\; 1.37 \,-\, 0.04{h^2} \,+\, 0.08z \,-\def\LDeqbreak{} 0.01{z^2}$ for the Elliptical Lens

View Table

The calculated EFL for the spherical lens is less than $\frac{\lambda}{{100}}$ from the accepted value of Zemax, while the elliptical lens has an answer within $\frac{\lambda}{{25}}$, where $\lambda = 550\;{\rm nm} $.

The f-number, $f/\#$, of the spherical lens can now be calculated using the following formula:

f/\# = \frac{{\rm EFL}}{D},

where $D$ is the diameter of the lens. Here it can be seen that the $f/\#$ is equal to approximately 3.247. The depth of focus (DOF) formula is

\Delta z = \pm 2\lambda {(f/\#)^2},

where $\Delta z$ is the change in DOF and $\lambda$ is the wavelength 550 nm. The result of this leads to a $\Delta z$ value of 0.01597 mm. Looking at the LSA plot, a distance of $\Delta z$ corresponds to an initial height of 0.12 mm, which represents the paraxial region’s maximum height.

The BFDs of the lens for light rays entering the lens at ${z_1}$, and emerging at ${z_2}$, from finite heights are calculated to study the longitudinal spherical aberration (LSA) of the lens. The BFD is calculated using the following expression:

{\rm BFD} = - \frac{{h({z_2})}}{{h^\prime ({z_2})}} - \left({\frac{{{n_3}}}{{{n_4}}} - {z_2}} \right).

Figures 5 and 6 show the trajectory of light rays of varying heights passing through a spherical and elliptical GRIN optical system, respectively. From this, the principal surfaces of the lens can be located. The LSA of both proposed lenses can be seen in Fig. 7. Figure 8 shows the difference between the BFD calculated using the Frobenius and those from ray-tracing software. These methods provide an accurate, better than $\frac{\lambda}{{14}}$, estimate of the ray trajectory within the GRIN, allowing one to predict the ray exit height and exit angle with respect to the optical axis and thus giving the information about lens aberrations, such as spherical aberration and coma.

Fig. 5. Plot showing the trajectory of light propagating through the spherical GRIN lens, with an example of iso-indicial contours present. The primary and secondary principal surfaces are shown in blue and red, respectively. It can be seen that the principal points coincide. A density plot bar is located on the right to show the refractive index at each point of the lens.

Download Full Size | PDF

Fig. 6. Plot showing the trajectory of light propagating through the elliptical GRIN lens, with an example of iso-indicial contours present. The primary and secondary principal surfaces are shown in blue and red, respectively. It can be seen that the principal points do not coincide. A density plot bar is located on the right to show the refractive index at each point of the lens.

Download Full Size | PDF

Fig. 7. LSA plots for the two optical systems outlined here. For the spherical lens, Zemax OpticStudio is shown in black and the Frobenius method in red, and for the elliptical lens, Zemax OpticStudio data are in gray and the Frobenius method in cyan.

Download Full Size | PDF

Fig. 8. Plot of the difference between Zemax OpticStudio and the Frobenius method for the spherical GRIN lens in blue and the elliptical lens in red.

Download Full Size | PDF

The limitations of this model are examined. As the radius of curvature of the external shape, and thus the iso-indicial contours, the conic and the maximum refractive index are dependent on the GRIN coefficients ${n_1}$, ${n_3}$, and ${n_4}$ for a quadratic lens, and there will be no restrictions here. The methods developed here are shown using the lens immersed in a medium that is normalized to the external surface of the lens, but they can be used for a lens in air, or any other medium. The lens is generally studied within a normalized medium, but this can restrict the method [38]. Snell’s law will need to be used to calculate the new angle of incidence, ${u_0} = - {\rm arctan} ({\rm arcsin} (\frac{{{n_m}}}{{{n_0}}}\sin {\theta _m}))$, where ${n_m}$ is the refractive index of the medium the lens is immersed in and ${\theta _m}$ is the angle of convergence within this medium. Once this is established, then this method holds true for a symmetric GRIN lens in any medium, not just a normalized medium.

6. CONCLUSION

The proposed method for semi-analytical ray-tracing has been presented and used in two examples of symmetric GRIN lenses with different conics (spherical and elliptical). The light rays propagating through a GRIN system following the power law of Eq. (2) can be modeled using the two derived expressions in this paper. Equation (14), the Legendre solution to Eq. (12), will provide information for the paraxial region, but it does not meet the required accuracy outside of this. Equation (17), the Frobenius solution to Eq. (16), allows one to ray-trace through the GRIN lens at finite initial heights. These equations will allow one to accurately measure the EFL and BFDs of the lens for rays at finite heights in absence of ray-tracing software, while also accurately estimating the power of the GRIN lens at any annular zone and the LSA.

Due to the nature of the differential equation, and the proposed solution in this research, it would be possible to reconstruct the refractive index profile of a GRIN lens. The inverse problem for GRIN lenses would be to observe a predetermined path where it is known where light enters and emerges from the lens and, from there, to calculate the GRIN power law of the lens, such as Eq. (2). This is an area of interest currently [39–44]. The techniques proposed here would allow for this as it has the potential to calculate the refractive index along the ray path. In principal, this will allow a user to reconstruct a quadratic GRIN lens of a given geometry and optical power.

The semi-analytical nature of our method is due the inherent residual error as a result of the number of finite terms used in the Frobenius polynomial. However, at every step in the method, no numerical iterations were used to obtain the required parameters. Therefore, the proposed method is superior to any numerical methods when solving inverse problems.

APPENDIX A

Using Fermat’s principle on a vertical slab of refractive index $n(h) = {n_0} - {n_1}{h^2}$, surrounded by air, we see the following relationship:

n(h)dz + \sqrt {{h^2} + {f^2}} = {n_0}dz + f,

where ${n_1}$ is the GRIN coefficient from Eq. (2) that is responsible for the rate of change of refractive index perpendicular to the optical axis. This leads to

\begin{split}{h^2} &= 2f({n_0} - n(h))dz, \\[-3pt] \frac{F}{{dz}} &= \frac{{2({n_0} - n(h))}}{{{h^2}}},\end{split}

where

\begin{split}\frac{F}{{dz}} &= \frac{{n^\prime}}{R}, \\[-3pt] \frac{{n^\prime}}{R} &= \frac{{2({n_0} - n(h))}}{{{h^2}}},\end{split}

such that

\frac{{n^\prime}}{R}h = 2{n_1}h.

APPENDIX B

An approximate equation was developed for the effective focal length from the Frobenius method when used for the paraxial region. The equation for the EFL is ${-}\frac{{{h_0}}}{{{n_0}h^\prime ({z_2})}}$. For the paraxial region using the recurrence relations, the values of the constants are as follows:

\begin{split}{a_0} &= {h_0}, \\ {b_0}&= 0.\end{split}

Using the Frobenius method, a polynomial for $h(z)$ is found using 19 terms as this is sufficient as per Fig. 3, and then this is differentiated with respect to $z$. For the paraxial zone, the exit point can be approximated by the thickness of the lens. This means that ${z_2} = \frac{{{n_3}}}{{{n_4}}}$. The polynomial is reduced by removing negligible terms that do not have an impact on the EFL, and the approximate solution is presented here. This can be used for a spherical or elliptical lens,

\begin{split}{\rm EFL} &= (315n_0^3n_4^7)/({n_1}{n_3}(630n_0^3n_4^6 \\&\quad-210n_0^2{n_1}n_3^2n_4^4 + 21{n_0}{n_1}n_3^4n_4^2({n_1} + 2{n_4}) \\ &\quad-{n_1}n_3^6(n_1^2 + 8{n_1}{n_4} + 9n_4^2))).\end{split}

The above formula gives a difference of 18 nm and 3 nm for the spherical and elliptical lens, respectively.

Funding

Hardiman Research Scholarship at the University of Galway, Ireland; Irish Research eLibrary.

Acknowledgment

Open access funding provided by Irish Research eLibrary. This research was also supported by the Hardiman Research Scholarship at the University of Galway, Ireland.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. W. Liu, H. Hu, F. Liu, and H. Zhao, “Manipulating light trace in a gradient-refractive-index media: a Lagrangian optics method,” Opt. Express 27, 4714–4726 (2019). [CrossRef]

2. R. Navarro, “The optical design of the human eye: a critical review,”J. Optometry 2, 3–18 (2009). [CrossRef]

3. R. Navarro, F. Palos, and L. González, “Adaptive model of the gradient index of the human lens. I. Formulation and model of aging ex vivo lenses,” J. Opt. Soc. Am. A 24, 2175–2185 (2007). [CrossRef]

4. A. Jaimes-Nájera, J. E. Gómez-Correa, V. Coello, B. K. Pierscionek, and S. Chávez-Cerda, “Single function crystalline lens capable of mimicking ciliary body accommodation,” Biomed. Opt. Express 11, 3699–3716 (2020). [CrossRef]

5. D. A. Atchison and L. N. Thibos, “Optical models of the human eye,” Clin. Exp. Optometry 99, 99–106 (2016). [CrossRef]

6. R. G. Zainullin, A. B. Kravtsov, and E. P. Shaitor, “The crystalline lens as a Luneburg lens,” Biofizica 19, 913–915 (1974).

7. A. V. Goncharov, “Vision in nature through GRIN media: smart optical design,” Proc. SPIE 11481, 1148103 (2020). [CrossRef]

8. J. W. Blaker, “Toward an adaptive model of the human eye,” J. Opt. Soc. Am. 70, 220–223 (1980). [CrossRef]

9. H. Liou and N. Brennan, “Anatomically accurate, finite model eye for optical modeling,” J. Opt. Soc. Am. A 14, 1684–1695 (1997). [CrossRef]

10. G. Smith, D. A. Atchison, and B. K. Pierscionek, “Modeling the power of the aging human eye,” J. Opt. Soc. Am. A 9, 2111–2117 (1992). [CrossRef]

11. J. He, J. Gwiazda, F. Thorn, and R. Held, “Wave-front aberrations in the anterior corneal surface and the whole eye,” J. Opt. Soc. Am. A 20, 1155–1163 (2003). [CrossRef]

12. A. Popiolek-Masajada and H. T. Kasprzak, “A new schematic eye model incorporating accommodation,” Optom. Vis. Sci. 76, 720–727 (1999). [CrossRef]

13. G. Smith, P. Bedggood, R. Ashman, M. Daaboul, and A. Metha, “Exploring ocular aberrations with a schematic human eye model,” Optom. Vis. Sci. 85, 330–340 (2008). [CrossRef]

14. G. Smith, “The optical properties of the crystalline lens and their significance,” Clin. Exp. Optom. 86, 3–18 (2003). [CrossRef]

15. C. J. Sheil and A. V. Goncharov, “Accommodating volume-constant age-dependent optical (AVOCADO) model of the crystalline GRIN lens,” Biomed. Opt. Express 7, 1985–1999 (2016). [CrossRef]

16. A. Gullstrand, “Appendix II,” in Handbuch der Physiologischen Optik, J. P. Southhall, ed., 3rd ed. (Optical Society of America, 1924), Vol. 1, pp. 351–352.

17. A. V. Goncharov and C. Dainty, “Wide-field schematic eye models with gradient-index lens,” J. Opt. Soc. Am. A 24, 2157–2174 (2007). [CrossRef]

18. H. T. Kasprzak, “New approximation for the whole profile of the human crystalline lens,” Ophthalmic Physiolog. Opt. 20, 31–43 (2000). [CrossRef]

19. T. Sakamto, “Ray trace algorithms for GRIN media,” Appl. Opt. 26, 2943–2946 (1987). [CrossRef]

20. W. S. Jagger and P. J. Sands, “A wide-angle gradient index optical model of the crystalline lens and eye of the octopus,” Vis. Res. 39, 2841–2852 (1999). [CrossRef]

21. W. S. Jagger and P. J. Sands, “A wide-angle gradient index optical model of the crystalline lens and eye of the rainbow trout,” Vis. Res. 39, 2841–2852 (1996).

22. B. K. Pierscionek and J. W. Regini, “The gradient index lens of the eye: an optobiological synchrony,” Prog. Retinal Eye Res. 31, 332–349 (2012). [CrossRef]

23. D. T. Moore, “Gradient-index optics: a review,” Appl. Opt. 19, 1035–1038 (1980). [CrossRef]

24. J. J. Esteve-Taboada, R. Montés-Micó, and T. Ferrer-Blasco, “Schematic eye models to mimic the behavior of the accommodating human eye,” J. Cataract Refractive Surg. 44, 627–641 (2018). [CrossRef]

25. J. E. Gómez-Correa, S. E. Balderas-Mata, B. K. Pierscionek, and S. Chávez-Cerda, “Composite modified Luneburg model of human eye lens,” Opt. Lett. 40, 3990–3993 (2015). [CrossRef]

26. Y. Huang and D. T. Moore, “Human eye modeling using a single equation of gradient index crystalline lens for relaxed and accommodated states,” in International Optical Design, Technical digest (CD) (Optica Publishing Group, 2006), paper MD1.

27. A. Popiolek-Masajada and H. Kasprzak, “Model of the optical system of the human eye during accommodation,” Ophthalmic Physiolog. Opt. 22, 201–208 (2002). [CrossRef]

28. A. Walther, The Ray and Wave Theory of Lenses (Cambridge University, 1995).

29. M. Bahrami and A. V. Goncharov, “Geometry-invariant GRIN lens: finite ray tracing,” Opt. Express 22, 27797–27810 (2014). [CrossRef]

30. M. Bahrami and A. V. Goncharov, “Geometry-invariant gradient refractive index lens: analytical ray tracing,” J. Biomed. Opt. 17, 055001 (2012). [CrossRef]

31. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 1999), pp. 129–131.

32. G. Beliakov and D. Y. C. Chan, “Analysis of inhomogeneous optical systems by the use of ray tracing. I. Planar systems,” Appl. Opt. 36, 5303–5309 (1997). [CrossRef]

33. F. Manns and A. Ho, “Paraxial equivalent of the gradient-index lens of the human eye,” Biomed. Opt. Express 13, 5131–5150 (2022). [CrossRef]

34. R. Navarro, V. Lockett-Ruiz, and J. L. López, “Analytical ray transfer matrix for the crystalline lens,” Biomed. Opt. Express 13, 5836–5848 (2022). [CrossRef]

35. V. Lakshimarayanan, A. K. Ghatsk, and K. Thyagarajan, Lagrangian Optics (Kluwer Academic, 2002).

36. E. Miersemann, Calculus of Variations Lecture Notes (Leipzig University, 2012).

37. D. T. Moore, “Design of singlets with continuously varying indices of refraction,” J. Opt. Soc. Am. 61, 886–894 (1971). [CrossRef]

38. J. E. Gómez-Correa, “Geometrical-light-propagation in non-normalized symmetric gradient-index media,” Opt. Express 30, 33896–33910 (2022). [CrossRef]

39. S. Nemoto and T. Makimoto, “Refractive-index distribution for a prescribed ray path,” J. Opt. Soc. Am. 69, 450–454 (1979). [CrossRef]

40. A. V. Goncharov, M. Nowakowski, M. T. Sheehan, and C. Dainty, “Reconstruction of the optical system of the human eye with reverse ray-tracing,” Opt. Express 16, 1692–1703 (2008). [CrossRef]

41. H. Ohno and T. Usui, “Neural network gradient-index mapping,” OSA Contin. 4, 2543–2551 (2021). [CrossRef]

42. W. Liu, “Ray tracing in concentric gradient-index media: optical Binet equation,” J. Opt. Soc. Am. A 39, 1025–1033 (2022). [CrossRef]

43. H. Ohno and T. Usui, “Design method of gradient-index optics in homogeneous medium using neural network gradient-index mapping,” Opt. Rev. 30, 355–360 (2023). [CrossRef]

44. J. E. Gómez-Correa, A. L. Padilla-Ortiz, J. P. Trevino, A. Jaimes-Nájera, N. Lozano-Crisóstomo, A. Cornejo-Rodriguez, and S. Chávez-Cerda, “Symmetric gradient-index media reconstruction,” Opt. Express 31, 29196–29212 (2023). [CrossRef]

	EFL (mm)
	Sphere	Elliptical
Zemax	12.988489	25.239333
Legendre	12.988492	25.239355
Frobenius	12.988492	25.239355

Semi-analytical finite ray-tracing through the quadratic symmetric GRIN lens

Abstract

1. INTRODUCTION

2. PARAXIAL DIFFERENTIAL EQUATION

A. Deriving an Equation

B. Solving the Equation

3. FINITE DIFFERENTIAL EQUATION

A. Deriving a Non-Paraxial Equation

B. Frobenius Method

4. ELLIPTICAL LENS

5. RESULTS

6. CONCLUSION

APPENDIX A

APPENDIX B

Funding

Acknowledgment

Disclosures

Data availability

REFERENCES

Data availability

Cited By

Figures (8)

Tables (1)

Equations (52)

Applied Optics