GRINCU lens with conicoid iso-indicial surfaces: application for modeling the crystalline lens

Rafael Navarro; Sergio Baquedano; Ana I. Sánchez-Cano

doi:10.1364/OE.435266

1. Introduction

Gradient index (GRIN) lenses are made of inhomogeneous materials, and inside the lens, light propagates through inhomogeneous media [1,2]. Inhomogeneous media are common both in nature and in optical technology. Two paradigmatic examples are air, which always presents some level of inhomogeneity associated with atmospheric turbulence, and the human crystalline lens, which is a classic example of a GRIN lens. They eyes of most animal species contain some type of GRIN optics. One of the first GRIN lens designs was Maxwell’s fisheye dating from 1854 [1]. Advances in theory, design of index distributions, materials, fabrication and applications are important, and GRIN lenses are currently widely used.

Among the wide variety of functions or profiles of refractive index proposed in the literature, we are interested in an onion-type lens. In the continuous case, this lens is made of an infinite number of layers of differential thickness; each layer acts as a thin (infinitesimal) lens. For the differential onion layer to act as a standard optical surface, we want the refractive index n to be constant within that layer. Such layers are iso-indicial surfaces (IISs) with index n. This surface separates the anterior and posterior media with a difference of the refractive index. If the layer, or IIS, is spherical with radius R, then we can compute its optical power as $P = {{d n} / R} = Cd n$, where $C = {1 / R}$ is the curvature and $d n$ is the differential refractive index. This expression from first-order (paraxial) optics suggests that the gradient index (proportional to $d n$) and the curvature C of the IIS are equally important and have a mutual multiplicative effect. The idea of having optical media with two gradients of index and curvature (GRIN/GRCU or GRINCU) was explored in a previous work [3]; the results obtained for the first-order optics approximation confirmed a strong multiplicative effect on the power of the GRINCU lens. Of course, most types of GRIN distributions have an implicit gradient of curvature (the GRCU can be analyzed by computing the curvatures of the IISs), but to the best of our knowledge, this is the first time that the optical effect of the GRCU has been explicitly modeled and studied [3].

In this context, there are two goals of the present work. In the first part, we develop the formulation of the index distribution of a generic GRIN/GRCU lens composed of differential onion layers. To guarantee analytical, treatable expressions, we restrict ourselves to the particular case when the IISs are revolution conicoids, and the gradient of their apical curvature radius is a constant parameter –G [3]. To this end, we generalize and upgrade a formulation developed before to model the crystalline lens [4].

Our second goal is to apply that generic formulation to model the crystalline lens. The model is adaptive [4,5], geometry-invariant [6] or external-surface-following [7] in the sense that the internal distribution automatically adapts to the external geometry.

This type of model is especially useful for modeling the shape of and optical changes in the lens with both age and accommodation. The inhomogeneity of the crystalline lens has long been known [8]; various sources in the extended literature have reported on experimental data [9–13], models [14–16] and finite ray tracing implementations [17–19] of the GRIN human lens.

2. General formulation

Here, we consider a singlet lens with two surfaces, the anterior and posterior surfaces. The signs of the curvature radii of both surfaces are usually used to classify the different types. When both signs are equal (+, +) or (-, -), the lens is a (positive or negative) meniscus, and when the two signs are opposite, the lens is either biconvex (+, -) or biconcave (-, +) lens.

2.1 Meniscus lens

For our formulation, the simplest case is the meniscus. We use cylindrical coordinates $({\omega ,\theta ,z} )$, but in the case of rotational symmetry around the optical axis, the formulation becomes 2-dimensional (2D) $({\omega ,z} )$, where ${\omega ^2} = {x^2} + {y^2}$ is the squared distance to the Z axis and z is the coordinate along the optical axis. In the meniscus, we consider the inner distribution of the refractive index as a function of the normalized distance along the optical axis ${z_0}$ to the anterior surface.

(1)$$n = n({{\raise0.7ex\hbox{${{z_0}}$} \!\mathord{\left/ {\vphantom {{{z_0}} t}} \right.}\!\lower0.7ex\hbox{$t$}}} ).$$

Here, t is the axial thickness of the meniscus lens, and inside the lens, the distance to the anterior surface is $0 \le {z_0} \le t$, so that the function n only depends on a single normalized parameter always ranging from 0 to 1. By definition, IISs are defined by $n = \textrm{constant}$. Then, for each value of ${z_0}$ we can define an IIS with a given value of n. We consider that the first IIS, that is, the IIS for ${z_0} = 0$, coincides with the anterior surface of the lens; the last IIS, the IIS for ${z_0} = t$, coincides with the posterior surface. In a similar way, two IISs separated by a differential distance $d {z_0}$ form a differential meniscus with infinitesimal thickness. We can completely cover the volume of the meniscus lens with an infinite number of infinitesimal menisci, and hence, we can completely define the distribution of the refractive index inside the lens. Below, we consider that both the external and the inner IISs are revolution conicoids (quadratic surfaces) defined by two parameters: their apical curvature radius R and conic constant Q.

When the anterior ${R_a} = R({{z_0} = 0} )$ and posterior ${R_p} = R({{z_0} = t} )$ radii are different, there must be a gradient of curvature radii of the IISs to allow a smooth transition between these two external surfaces. In a previous study [3], a constant gradient -G was introduced so that the curvature radius of the IIS varied in a linear way:

(2)$$R({{z_0}} )= {R_a} - G{z_0}.$$

Then, it is straightforward to show that $G = {{({{R_a} - {R_p}} )} / t}$ to match the values of both anterior and posterior radii of the meniscus. This means that the GRIN meniscus has two gradients of refractive index GRIN and curvature radius GRCU. Below, we use the acronym GRINCU to refer to this new type of two gradients lens.

2.2 Iso-indicial surfaces

For the particular case of a rotationally-symmetric concentric distribution, we showed in Eq. (3) and Eq. (3) of a previous publication [4] that the conical iso-indicial surfaces (IISs) are given by the expression ${r^2} = {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 f}} \right.}\!\lower0.7ex\hbox{$f$}}\left( {\frac{{{{({z - \Delta } )}^2}}}{{{a^2}}} + \varepsilon \frac{{{\omega^2}}}{{{b^2}}} - O} \right)$ in cylindrical coordinates z, ω$({{\omega^2} = {x^2} + {y^2}} )$; a, b are the conic semiaxes; $\varepsilon ={+} 1$ for ellipses and $\varepsilon ={-} 1$ for hyperbolas; $\Delta = \varepsilon a$ means that the correct (left or right) branch of the conicoid intersects the optical axis at point $({z = {z_0},\omega = 0} )$. The variable r is a normalized (between 0 and 1) non-Euclidean distance to the center; the normalization is done by applying an offset O and a normalization factor f. For each value of r we have a different IIS. This generic expression applies to both anterior (Eq. (3)) and posterior (Eq. (3)) parts of the lens simply by assigning the corresponding values to the parameters. By simple operations, we can separate the standard expression of the conicoid on the left side:

(3)$$\frac{{{{({z - \Delta } )}^2}}}{{{a^2}}} + \varepsilon \frac{{{\omega ^2}}}{{{b^2}}} = f{r^2} + O,$$

On the right side of this equation, $f{r^2} + O$ can take a value between 0 and 1. In the last case, we have the equation of the first (external, anterior) surface; For the internal IISs, the value of $f{r^2} + O$ decreases progressively with the distance to the surface, which means a progressive shift and scaling for the complete family of IISs (see Ref. [4] for details). In general, any IIS intersects the optical axis at a point $({z = {z_0},\omega = 0} )$. In Ref. [4], each IIS, and hence the value of the refractive index, was determined by the variable ${r^2}$: $n = n({{r^2}} )$. However, in the present formulation (Eq. (1)), we changed the parametrization using the normalized variable ${{{z_0}} / t}$ instead of ${r^2}$. Note the difference in nature between z and ${z_0}$: z is the coordinate along the optical axis, whereas ${z_0}$ is the iso-indicial surface parameter. Each IIS intersects the optical axis at point $({z = {z_0},\omega = 0} )$, and thus we can use ${z_0}$ to determine each particular IIS.

To adapt the IIS expression to the present parametrization, we particularize Eq. (3) to the point of intersection with the optical axis $({z = {z_0},\omega = 0} )$, ${\raise0.7ex\hbox{${{{({{z_0} - \Delta } )}^2}}$} \!\mathord{\left/ {\vphantom {{{{({{z_0} - \Delta } )}^2}} {{a^2}}}} \right.}\!\lower0.7ex\hbox{${{a^2}}$}} = f{r^2} + O$, and solve for ${z_0}$:

(4)$${z_0} = \Delta \pm a\sqrt {f{r^2} + O} .$$

Now, we can replace $f{r^2} + O$ with its expression in Eq. (3). In this way, we are replacing ${r^2}$ with the parameter ${z_0}$ in the expression of a generic IIS that intersects the optical axis at $({z = {z_0},\omega = 0} )$:

(5)$$\frac{{{{({z - \Delta } )}^2}}}{{{{({{z_0} - \Delta } )}^2}}} + \varepsilon \frac{{{a^2}{\omega ^2}}}{{{b^2}{{({{z_0} - \Delta } )}^2}}} = 1.$$

Here, it is easy to identify both semiaxes of the IIS as functions of ${z_0}$, where ${a^2}({{z_0}} )= {({{z_0} - \Delta } )^2}$ and ${b^2}({{z_0}} )= {\raise0.7ex\hbox{${{b^2}}$} \!\mathord{\left/ {\vphantom {{{b^2}} {{a^2}}}} \right.}\!\lower0.7ex\hbox{${{a^2}}$}}{({{z_0} - \Delta } )^2}$ with ${a^2}(0 )= {a^2}$ and ${b^2}(0 )= {b^2}$. Note that Eq. (5) represents a complete conicoid, whereas in the meniscus, we only need one half (or a “hemisphere”) of it. For this reason, it is better to use the standard expression of the sagitta z, where the sign of the curvature $C = {1 / R}$ determines the desired (convex or concave) side of the conicoid [20]:

(6)$$z = {z_0} + \frac{{C({{z_0}} ){\omega ^2}}}{{1 + \sqrt {1 - ({1 + Q({{z_0}} )} )C{{({{z_0}} )}^2}{\omega ^2}} }},$$

where both the curvature $C({{z_0}} )= {1 / {R({{z_0}} )}}$ and the conic constant $Q({{z_0}} )$ are functions of ${z_0}$. On a conicoid surface, the apical curvature radius and the conic constant are $R({{z_0}} )= {{\varepsilon {b^2}({{z_0}} )} / a}({{z_0}} )$ and $Q = {{\varepsilon {b^2}({{z_0}} )} / {{a^2}({{z_0}} )}} - 1$, respectively. Thus,

(7)$$\begin{array}{l} R({{z_0}} )= \varepsilon {\raise0.7ex\hbox{${{b^2}}$} \!\mathord{\left/ {\vphantom {{{b^2}} {{a^2}}}} \right.}\!\lower0.7ex\hbox{${{a^2}}$}}({\Delta - {z_0}} )= R({1 - \varepsilon {\raise0.7ex\hbox{${{z_0}}$} \!\mathord{\left/ {\vphantom {{{z_0}} a}} \right.}\!\lower0.7ex\hbox{$a$}}} )\textrm{ and}\\ Q({{z_0}} )= {\raise0.7ex\hbox{${\varepsilon {b^2}}$} \!\mathord{\left/ {\vphantom {{\varepsilon {b^2}} {{a^2}}}} \right.}\!\lower0.7ex\hbox{${{a^2}}$}} - 1 = Q. \end{array}$$

The sign of $R({{z_0}} )$ is +1 or -1 for the positive and negative sides or hemispheres; its gradient (derivative with respect to ${z_0}$) is $- ({Q + 1} )$ (positive) or $+ ({Q + 1} )$ (negative); $Q({{z_0}} )= Q$ is constant with gradient 0. These gradients correspond to concentric surfaces that only differ in shift and scaling when they maintain the same shape (Q).

Now, let us generalize this formulation to nonconcentric IISs. The simplest way to do this may involve introducing the curvature radius gradient parameter G. Below, we keep $Q({{z_0}} )= Q$ for the sake of simplicity. This means that we restrict ourselves to a meniscus with an external geometry determined by four independent variables (degrees of freedom) t, R_a, R_p and $Q = {Q_a} = {Q_p}$. In the case that ${Q_a} \ne {Q_p}$, one could introduce a constant gradient $- {{({{Q_a} - {Q_p}} )} / t}$ so that $Q({{z_0}} )= {Q_a} - ({{Q_a} - {Q_p}} ){{{z_0}} / t}$. Nevertheless, such generalization to three gradients GRIN/GRCU/GRCC (index/curvature/conic constant) is beyond the scope of the present study and will be the subject of future work.

With the gradient curvature parameter G, the geometry of the internal IIS is determined either by R or Q:

(8a)$$R({{z_0}} )= R({1 - \varepsilon G{\raise0.7ex\hbox{${{z_0}}$} \!\mathord{\left/ {\vphantom {{{z_0}} a}} \right.}\!\lower0.7ex\hbox{$a$}}} )\textrm{ and }Q({{z_0}} )= Q,$$

or by the semiaxes:

(8b)$${a^2}({{z_0}} )= {({\Delta - \varepsilon G{z_0}} )^2}\textrm{ and }{b^2}({{z_0}} )= \frac{{{b^2}}}{{{a^2}}}{({\Delta - \varepsilon G{z_0}} )^2}.$$

Let us recalculate the expression of the conicoid (Eq. (5)) after introducing parameter G to replace ${({{z_0} - \Delta } )^2} = {({\Delta - {z_0}} )^2}$ with ${({\Delta - \varepsilon G{z_0}} )^2}$ using the new expressions of the semiaxes (Eq. (8b)). In addition, we have to recalculate a new $\Delta ^{\prime} = \Delta - ({\varepsilon G - 1} ){z_0}$ to guarantee that the IIS intersects the optical axis at point $({z = {z_0},\omega = 0} )$. Then, we arrive to the equation of the IIS:

(9)$$\frac{{{{({z - \Delta + ({\varepsilon G - 1} ){z_0}} )}^2}}}{{{{({\Delta - \varepsilon G{z_0}} )}^2}}} + \varepsilon \frac{{{a^2}{\omega ^2}}}{{{b^2}{{({\Delta - \varepsilon G{z_0}} )}^2}}} = 1.$$

This is a second-order equation that we can use to solve for ${z_0}$. We can choose between the standard solution of the 2^nd-order equation or its sag expression: it is straightforward to show that ${{{z_0} = \left( { - B + \varepsilon \sqrt {{B^2} - 4AC} } \right)} / {2A}} ={-} {{2C} / {\left( {B + \varepsilon \sqrt {{B^2} - 4AC} } \right)}}$ (where we use $\varepsilon$ instead of ${\pm}$ to guarantee that $0 \le {z_0} \le t$). The sag form is more useful to determine the desired side of conicoid optical surfaces and to avoid indeterminations. After standard operations and simplifications, we arrive at the sag equation:

(10)$${z_0} = \frac{{{\raise0.7ex\hbox{${{z^2}}$} \!\mathord{\left/ {\vphantom {{{z^2}} \Delta }} \right.}\!\lower0.7ex\hbox{$\Delta $}} - 2z + C{\omega ^2}}}{{({1 - \varepsilon G} ){\raise0.7ex\hbox{$z$} \!\mathord{\left/ {\vphantom {z \Delta }} \right.}\!\lower0.7ex\hbox{$\Delta $}} - 1 - \sqrt {{{({{\raise0.7ex\hbox{${\varepsilon Gz}$} \!\mathord{\left/ {\vphantom {{\varepsilon Gz} \Delta }} \right.}\!\lower0.7ex\hbox{$\Delta $}} - 1} )}^2} - ({1 - 2\varepsilon G} )({1 + Q} ){C^2}{\omega ^2}} }}.$$

Thus, ${z_0} = {z_0}({z,\omega } )$ is a function of both coordinates and is defined for any point $({z,\omega } )$ inside the lens. Consequently, $n({{{{z_0}({z,\omega } )} / t}} )$ is a function of both coordinates $({z,\omega } )$, defined inside the meniscus lens. Figure 1 compares two GRINCU menisci with parallel IISs, that is, G = 0 (left), and with a curvature gradient parameter G = 2 (right).

Fig. 1. Example of two GRIN menisci with G = 0 (parallel IIS) and G = 2. The coordinate system is depicted in the first case; the color bar shows the scale of the refractive index values.

Download Full Size | PDF

2.3 General conicoid

Nonrotationally symmetric optical elements such as spherocylindrical lenses are commonly used in visual optics, among other applications. Here, we consider a more general 3-axis conicoid, which is a natural generalization of the rotationally symmetric 2-axis case. The equation for the anterior surface ($f{r^2} + O = 1$) is

(11)$$\frac{{{{({z - \Delta } )}^2}}}{{{a^2}}} + \varepsilon \left( {\frac{{{x^2}}}{{{c^2}}} + \frac{{{y^2}}}{{{d^2}}}} \right) = 1,$$

with semiaxes a, c and d. Compared to Eq. (3), we have simply replaced $\frac{{{\omega ^2}}}{{{b^2}}}$ with $\frac{{{x^2}}}{{{c^2}}} + \frac{{{y^2}}}{{{d^2}}}$. The rotationally symmetric conicoid is determined by two parameters, either $({{a^2},\varepsilon {b^2}} )$ or $({R,Q} )$. Now, we have three semiaxes a, c, d. Let us introduce a parameter D related to the surface astigmatism so that the surface is determined by three alternative parameters $({R,Q,D} )$ where $R = {\raise0.7ex\hbox{${({{R_{\max }} + {R_{\min }}} )}$} \!\mathord{\left/ {\vphantom {{({{R_{\max }} + {R_{\min }}} )} 2}} \right.}\!\lower0.7ex\hbox{$2$}}$ and $D = {\raise0.7ex\hbox{${({{R_{\max }} - {R_{\min }}} )}$} \!\mathord{\left/ {\vphantom {{({{R_{\max }} - {R_{\min }}} )} 2}} \right.}\!\lower0.7ex\hbox{$2$}}$ are the average and the deviation, respectively; equivalently, ${R_{\max }} = R + D$ and ${R_{\min }} = R - D$. For the semiaxes, we obtain ${c^2} = {b^2} + aD$ and ${d^2} = {b^2} - aD$ so that ${b^2} = {\raise0.7ex\hbox{${({{c^2} + {d^2}} )}$} \!\mathord{\left/ {\vphantom {{({{c^2} + {d^2}} )} 2}} \right.}\!\lower0.7ex\hbox{$2$}}$; for the conic constants ${Q_{\max }} = Q + {\raise0.7ex\hbox{$D$} \!\mathord{\left/ {\vphantom {D a}} \right.}\!\lower0.7ex\hbox{$a$}}$ and ${Q_{\min }} = Q - {\raise0.7ex\hbox{$D$} \!\mathord{\left/ {\vphantom {D a}} \right.}\!\lower0.7ex\hbox{$a$}}$, and then $Q = {\raise0.7ex\hbox{${({{Q_{\max }} + {Q_{\min }}} )}$} \!\mathord{\left/ {\vphantom {{({{Q_{\max }} + {Q_{\min }}} )} 2}} \right.}\!\lower0.7ex\hbox{$2$}}$.

Therefore, with these definitions, the parameters of the rotationally symmetric case are equal to the mean of the maximum and minimum values. Finally, the surface can be rotated at any angle $\theta $ around the optical axis. If we apply the corresponding 2D rotation matrix ${\mathbf x^{\prime} = }{{\mathbf U}_{\mathbf \theta }}{\mathbf x}$ to the vector of coordinates ${\mathbf x} = {[{x,y} ]^T}$ we arrive at the expression equivalent to Eq. (3) for the generic sphero-cylindrical IIS.

(12)$$\frac{{{{({Q + 1} )}^2}}}{{{R^2}}}{\left( {z - \frac{R}{{({Q + 1} )}}} \right)^2} + \frac{{({Q + 1} )}}{R}\left( {\frac{{{{({x\cos \theta + y\sin \theta } )}^2}}}{{R + D}} + \frac{{{{({ - x\sin \theta + y\cos \theta } )}^2}}}{{R - D}}} \right) = f{r^2} + O.$$

The geometry of the nonrotationally symmetric IIS is determined by four parameters: R, usually known as the radius of the sphere; D, the radius of the cylinder; $\theta$, the axis (orientation) of the cylinder; and Q, the mean conic constant. To implement this spherocylindrical conicoid geometry, we only have to replace ${\omega ^2}$ with its expression in Eq. (10).

2.4 Biconvex or biconcave lens

The formulation for these lenses is somewhat more complex than that for the meniscus, as the anterior and posterior radii have opposite signs. Not only the external surface but also each ISS has a positive and a negative side or “hemisphere”. Here, it is useful to consider the anterior and posterior sides independently and find the natural central interface of the ISS. This interface is an ellipse or a circumference for the rotationally symmetric case [4]. The axial thickness of the lens is the sum of the anterior and posterior thicknesses $t = {t_a} + {t_p}$, and the anterior and posterior refractive index distributions are functions of the normalized axial distances to their corresponding anterior or posterior external surface:

(13)$$n = \left\{ {\begin{array}{c} {{n_a}({{{1 - {z_{0a}}({z,\omega } )} / {{t_a}}}} )\textrm{for }0 \le {z_0} \le {t_a}}\\ {{n_p}({{{({{z_{0p}}({z,\omega } )- {t_a}} )} / {{t_p}}}} )\textrm{for }{t_a} < {z_0} \le t.} \end{array}} \right.$$

In this formulation, each IIS is composed of two sides or hemispheres with the same biconvex or biconcave pattern as the external surface of the lens. The arguments of ${n_a}({{s_s}} )$ and ${n_p}({{s_p}} )$ are ${{{s_a} = 1 - {z_{0a}}({z,\omega } )} / {{t_a}}}$ and ${{{s_p} = ({{z_{0p}}({z,\omega } )- {t_a}} )} / {{t_p}}}$, respectively. The normalizations and offsets are chosen in such a way that they take values from 1 at the external surface ${s_a}(0 )= {s_p}(t )= 1$ to 0 at the center ${s_a}({{t_a}} )= {s_p}({{t_a}} )= 0$. Note that both ${z_{0a}}$ and ${z_{0p}}$ are given by Eq. (10) but using the specific values for each side; thus, we now use $a_a^2$, $b_a^2$_, ${\varepsilon _a}$ for ${z_{0a}}$ and $a_p^2$, $b_p^2$_, ${\varepsilon _p}$ for ${z_{0p}}$; ${\Delta _a} = {\varepsilon _a}{a_a}$ and ${\Delta _p} = t - {\varepsilon _p}{a_p}$ [4]. Below, we consider the same G for both sides for the sake of simplicity. The complete distribution of the refractive index is obtained by replacing ${z_0}$ with its expression in Eq. (10). However, the curvature gradient has to keep the same sign on both sides, whereas the curvature radius reverses its sign. Therefore, we have to reverse the expression of the radius in Eq. (8a) by replacing ${z_0}$ with $({t - {z_0}} )$ so that ${R_a}({{z_0}} )= {R_a}({1 - {\raise0.7ex\hbox{${{\varepsilon_a}G{z_0}}$} \!\mathord{\left/ {\vphantom {{{\varepsilon_a}G{z_0}} {{a_a}}}} \right.}\!\lower0.7ex\hbox{${{a_a}}$}}} )$ and ${R_p}({{z_0}} )= {R_p}({1 - {\raise0.7ex\hbox{${{\varepsilon_p}G({t - {z_0}} )}$} \!\mathord{\left/ {\vphantom {{{\varepsilon_p}G({t - {z_0}} )} {{a_p}}}} \right.}\!\lower0.7ex\hbox{${{a_p}}$}}} )$. Then, we can directly use Eq. (10) for ${z_{0a}}$ and compute the corresponding expression for ${z_{0p}}$ so that

(14a)$${z_{0a}} = \frac{{({1 - {\varepsilon_a}G} )z - {\varepsilon _a}{a_a} + {\varepsilon _a}\sqrt {{{({Gz - {a_a}} )}^2} + ({2{\varepsilon_a}G - 1} ){\varepsilon _a}\frac{{{a_a}^2}}{{{b_a}^2}}{\omega ^2}} }}{{({1 - 2{\varepsilon_a}G} )}}$$

and

(14b)$${z_{0p}} = \frac{{({1 - {\varepsilon_p}G} )z - {\Delta _p} - {\varepsilon _p}Gt - {\varepsilon _p}\sqrt {{{({{\Delta_p} - {\varepsilon_p}G({z - t} )} )}^2} - ({1 - 2{\varepsilon_p}G} ){\varepsilon _p}\frac{{{a_p}^2{\omega ^2}}}{{{b_p}^2}}} }}{{({1 - 2{\varepsilon_p}G} )}}.$$

The anterior and posterior sides intersect the optical axis at $({{z_0}_a,\omega } )$ and $({{z_0}_p,\omega } )$, respectively, and they join each other at a given point of intersection $({{z_i},{\omega_i}} )$. In the rotationally symmetric case, the intersection of the two halves of the IIS is a circumference of radius ${\omega _i}$. By definition, both sides of the IIS share the same refractive index value. Therefore, the intersection occurs when the arguments are equal ${s_a} = {s_p}$:

(15)$$\frac{{{z_0}_p({{z_i},{\omega_i}} )- {t_a}}}{{{t_p}}} = \frac{{{t_a} - {z_0}_a({{z_i},{\omega_i}} )}}{{{t_a}}}.$$

The set of points $({{z_i},{\omega_i}} )$, i.e., the locus formed by the intersections $({{z_i},{\omega_i}} )$ (or crossings) between the two hemispheres of every IIS forms a central interface surface. The intersection of this surface with the optical axis occurs at ${z_0}_a = {z_0}_p = {t_a}$. In this sense, as suggested by Eq. (11), the GRIN lens is a sort of cemented doublet with continuity (no change) of refractive index at the interface. Note that in the general nonrotationally symmetric case, ${\omega _i} = {\omega _i}({{x_i},{y_i}} )$. Then, we can impose ${\omega _{ia}} = {\omega _{ip}}$ as a sufficient condition to keep Eq. (15) in its 2D form. It is straightforward to show that this is achieved by imposing ${\theta _a} = {\theta _p}$ and ${D_p} = {\raise0.7ex\hbox{${{R_p}}$} \!\mathord{\left/ {\vphantom {{{R_p}} {{R_a}}}} \right.}\!\lower0.7ex\hbox{${{R_a}}$}}{D_a}$.

Now, we can obtain the equation for the central interface surface by replacing ${z_0}_a$ and ${z_0}_p$ with their expressions (Eq. (10)). However, it is straightforward to show that it is a fourth-order equation instead of a second-order conicoid as in the former case of concentric IISs [4].

We implemented two alternative analytical [21] and numerical solutions for the quartic equation. The detailed formulation of the analytical solution can be found in [22], and we include our custom implementation as we show in Code 1 (Ref. [23]). Nevertheless, we implemented a more efficient solution for practical numerical implementation based on transforming Eq. (15) in the associated inequalities for the arguments ${s_a} = \frac{{{t_a} - {z_0}_a({{z_i},{\omega_i}} )}}{{{t_a}}}$ and ${s_p} = \frac{{{z_0}_p({{z_i},{\omega_i}} )- {t_a}}}{{{t_p}}}$. For any point inside the lens, we have three possibilities: (1) ${s_a} = {s_p}$ means that we are at the interface; (2) ${s_a} > {s_p}$ means that the point belongs to the anterior part; and (3) ${s_a} < {s_p}$ means that the point belongs to the posterior part. In addition, a point is inside the lens if both arguments are less than or equal to 1: ${s_a} \le 1$ and ${s_p} \le 1$.

3. GRINCU model of the crystalline lens

In this section, we propose a new general (adaptive [4] or geometrically invariant [6] or external surface-following [7]) model of the crystalline lens based on the above formulation. Basically, we reformulate a previous model [4,24], particularizing the general expression in Eq. (11) for that particular case. That index distribution model was based on a power law, with a noninteger degree (exponent) 2p; the index distribution for the anterior side ${n_a}$ was monotonically increasing between a minimum value at the anterior surface, and conversely, ${n_p}$ was monotonically decreasing between a minimum value at the external surfaces ${n_s} = {n_a}(0 )= {n_p}(0 )$ and a maximum value at the central interface surface ${n_c} = {n_a}(1 )= {n_p}(1 )$.

(16)$$\begin{array}{l} {n_a} = {n_c} - ({{n_c} - {n_s}} ){({{{1 - {z_0}({z,\omega } )} / {{t_a}}}} )^{2p}}\textrm{ for }0 \le {z_0} \le {t_a}\textrm{ and}\\ {n_p} = {n_c} - ({{n_c} - {n_s}} ){({{{({{z_0}({z,\omega } )- {t_a}} )} / {{t_p}}}} )^{2p}}\textrm{ for }{t_a} < {z_0} \le t. \end{array}$$

The generic expression of the IISs ${z_0}({z,\omega } )$ is given in Eq. (10), and the central interface can be obtained either in an analytical or numerical way, as explained above.

3.1 Numerical implementation

We implemented two versions of the model in MATLAB (The MathWorks, Inc.) scripts, which permitted us to validate the formulation and was a necessary step toward future implementation of finite ray tracing through the lens model. For the general 3-dimensional (3D) model, we computed both $n({x,y,z} )$ and its gradient $\nabla n({x,y,z} )= ({{{\partial n} / {\partial x,}}{{\partial n} / {\partial y,}}{{\partial n} / {\partial z}}} )$, that is, the three partial derivatives in a 3D sampling grid. In the present formulation, we mainly used cylindrical coordinates so that we applied the chain rule to compute the partial derivatives with respect to the Cartesian coordinates from the partial derivatives with respect to the cylindrical coordinates. Since the calculations of these derivatives are standard and the final expressions are long, they are not included here in the text but are included in the code (in MATLAB notation).

We used a sampling resolution of 0.05 mm. This means that for typical lens dimensions (9 mm diameter and 4 mm thickness), we have 2,592,000 sampling points (with finer resolutions, the computation slows down dramatically in normal PC computers with the risk of memory overflow). In this 3D implementation, we avoided solving the quartic equation in 3D because it strongly increases the computing time. Instead, we included the 2D (rotationally symmetric case) analytical solution in the supplementary script because it is not difficult to generalize it to 3D and then obtain the complete central interface surface. Nevertheless, we applied the numerical solution (much faster) based on the inequalities to the arguments (as explained above) to classify every point among four different options:

- if either ${s_a} > 1$ or ${s_p} > 1$, then the point is outside the lens; otherwise,
- if ${s_a} = {s_p}$ the point belongs to the central interface; otherwise,
- if ${s_a} > {s_p}$ the point belongs to the anterior part; otherwise
- if ${s_a} < {s_p}$ the point belongs to the posterior part.

The script for the computation of n and $\nabla n$ in 3D is provided as a supplementary file. It also includes the approximated computation of the paraxial power [3]. Additionally, we wrote a 2D version for the rotationally symmetric case, which enables us to work with much finer resolution (not included in the supplementary material).

3.2 Results

The resulting models for a crystalline lens with different curvature gradients G, ranging from 0 to 3, are compared in Fig. 2 for the rotationally symmetric case. The upper row shows color maps of the refractive index distributions; the black line represents the central interface computed by the analytical method. The lower row shows IISs. The parameters used to compute these models are adapted from Ref. 19 for the unaccommodated lens of a 30-year-old: R_a = 10.96 mm, Q_a = -3; R_p = -5.51 mm; t = 3.638 mm, t_a = 0.6t; n_s = 1.3726, n_c = 1.4301, p = 3.741.

Fig. 2. Refractive index distribution (upper row) of lens models with different curvature gradients, including the central interface. The contours in the bottom row represent IISs.

Download Full Size | PDF

Figure 3 shows a transverse slice for the plane z = t_a for a nonrotationally asymmetric case (with G = 2) with D_a = 1.096 and θ = 30°. Note the elliptical shape with the meridian of maximum radius oriented along the 30° meridian.

Fig. 3. Transversal slice corresponding to the plane z = t_a of the refractive index distribution for the case G = 2.

Download Full Size | PDF

The paraxial power of the lens models is represented as a function of the curvature gradient parameter G in Fig. 4. This approximated power was computed by the method described in Ref. 3 (it is included in the last part of the supplementary MATLAB script) using n_a = 1.3374 for the surrounding aqueous and vitreous humors.

Fig. 4. Refractive power of the lens models as a function of the curvature gradient G.

Download Full Size | PDF

Finally, the three partial derivatives of the gradient $\nabla n = ({{{\partial n} / {\partial x,}}{{\partial n} / {\partial y,}}{{\partial n} / {\partial z}}} )$ are shown in Fig. 5 for G = 2. The computation of the gradient is included here, as it is essential in ray tracing computations. Note the discontinuity in the three partial derivatives at the central interface, where the gradient reverses it sign from positive to negative values. This discontinuity could be avoided in different ways, as for instance using higher order models for the ISSs [6]. However experimental evidence suggests that such discontinuity could be present in real lenses [10].

Fig. 5. Three partial derivatives of the refractive index n (for G = 2).

Download Full Size | PDF

4. Discussion and conclusions

A generalization of a previous adaptive or external-surface-following type of GRIN lens [4] has been proposed. The first and most important new feature was to include a curvature radius gradient parameter –G of the IIS, which allows us direct and explicit control in the formulation. In this way, we formulated a new type of lens (GRINCU) characterized by two gradients of refractive index (GRIN) and of the curvature of the IIS (GRCU). We believe that the multiplicative effect of the two gradients has a high impact on the optical performance of this type of lens. In fact, we found a mutual multiplicative effect of both gradients on the lens power [3]. As a result, in the case of biconvex (or biconcave) lenses, the second new feature is that the central interface surface between the two anterior and posterior sides of the lens is no longer a plane or a conicoid but a fourth-order surface to guarantee the continuity of the refractive index distribution at the interface. We implemented two alternative analytical and numerical solutions to solve the fourth-order equation to find the surface. As in the former model [4], this does not guarantee continuity in the gradient $\nabla n$. In fact, we can clearly see these discontinuities (sign reversals) in the partial derivatives toward the peripheral parts of the interface (Fig. 5). Therefore, one must be cautious when implementing ray tracing through these GRINCU lenses. One of the simplest solutions is to consider the anterior and posterior parts separately (as a sort of cemented doublet) for ray tracing [4]. Bahrami and Goncharov [6] added a third-order term to the second-order conicoid IIS, which was adjusted to avoid discontinuity in the derivatives.

The third relevant new feature is formulation and implementation of a general nonrotationally symmetric case to include the case of toric-like (spherocylindrical) lenses with two different curvatures along the principal meridians. This is an important feature when modeling biological lenses, such as the crystalline lens, and potentially important as well, for example, when designing intraocular lenses.

The last part of this study consisted of applying this formulation to implement a GRINCU nonrotationally symmetric model of the crystalline lens. We confirmed that the refractive power (diopters, D) increased as we increased the curvature radius parameter G. The power increased by approximately 13.8% when changing from G = 0 (24.25 D) to G = 3 (27.6 D).

To the best of our knowledge, this is the first case of GRINCU lens formulation with an explicit control parameter of the inner curvature gradient. This is the main novelty of the present model. Sheil and Goncharov [6] introduced a new parameter m in their AVOCADO lens model. This parameter was used to control a nonlinear scaling of both curvature radii R and conic constant (shape factor) Q of the IIS and hence the index gradient along the radial variable ω. The main difference is that the gradient $\nabla R ={-} G$ is constant in our model and Q = constant (versus two different nonlinear gradients for both R and Q in Ref. 6). Our goal was to obtain the simplest mathematical formulation, that is, second-order conicoid surfaces: linear scaling (constant gradient) of R and constant Q. There are two main reasons for choosing the simplest approach compatible with experimental data and evidence.

On the one hand, simple (linear or second-order) models strongly facilitate both geometrical and physical (optical) interpretation. On the other hand, the experimental measurement of the refractive index distribution of the lens is extremely difficult. A variety of non-destructive methods have been developed to measure the distribution of refractive index within the crystalline lens. They are based on optical fibers [9], magnetic resonance imaging (MRI) [10], optical coherence tomography (OCT) [12], X-ray interferometry [13] or laser ray tracing [25]. In most cases numerical models are needed for understanding the raw data obtained with the different techniques. For this reason we need to use the models with the least sophistication (maximum simplicity) to analyze those data. The same principle applies to the design of GRIN lenses. In this sense, we believe that further sophistication should be added only when required. For example, experimental evidence suggests that the human lens may exhibit toricity to different extents on its surfaces [12], which motivated us to generalize the results to nonrotationally symmetric conicoids (Subsection 2.3). Further work will include adding a third gradient of conic constant Q, implementing finite ray tracing, and applying the GRINCU lens to develop accommodating age-dependent optical models of the human eye.

Funding

Ministerio de Ciencia, Innovación y Universidades (PID2019-107058RB-I00).

Disclosures

The authors declare no conflicts of interest.

Data availability

No data were generated or analyzed in the presented research.

References

1. E. W. Marchand, Gradient Index Optics (Academic, 1978).

2. C. Gomez-Reino, M.V. Perez, and C. Bao, Gradient-Index Optics: Fundamentals and Applications (Springer-Verlag, 2002). [CrossRef]

3. R. Navarro and N. López-Gil, “Impact of internal curvature gradient on the power and accommodation of the crystalline lens,” Optica 4(3), 334–340 (2017). [CrossRef]

4. 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(8), 2175–2185 (2007). [CrossRef]

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

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

7. Q. Li and F. Fang, “Physiology-like crystalline lens modelling for children,” Opt. Express 28(18), 27155–27180 (2020). [CrossRef]

8. A. Gullstrand, “Appendix II,” in Handbuch der Physiologischen Optik, Helmholtz, 3rd ed. English translation edited by J. P. Southall (Optical Society of America, 1962) Vol. 1, 351–352.

9. B. K. Pierscionek, “Refractive index contours in the human lens,” Exp. Eye Res. 64(6), 887–893 (1997). [CrossRef]

10. C. Jones, D. Atchison, R. Meder, and J. Pope, “Refractive index distribution and optical properties of the isolated human lens measured using magnetic resonance imaging (MRI),” Vision Res. 45(18), 2352–2366 (2005). [CrossRef]

11. S. Kasthurirangan, E. L. Markwell, D. A. Atchison, and J. M. Pope, “In vivo study of changes in refractive index distribution in the human crystalline lens with age and accommodation,” Invest. Ophthalmol. Vis. Sci. 49(6), 2531–2540 (2008). [CrossRef]

12. A. de Castro, S. Ortiz, E. Gambra, D. Siedlecki, and S. Marcos, “Three-dimensional reconstruction of the crystalline lens gradient index distribution from OCT imaging,” Opt. Express 18(21), 21905–21917 (2010). [CrossRef]

13. B. Pierscionek, M. Bahrami, M. Hoshino, K. Uesugi, J. Regini, and N. Yagi, “The eye lens: age-related trends and individual variations in refractive index and shape parameters,” Oncotarget 6(31), 30532–30544 (2015). [CrossRef]

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

15. J. A. Díaz, C. Pizarro, and J. Arasa, “Single dispersive gradient-index profile for the aging human lens,” J. Opt. Soc. Am. A 25(1), 250–261 (2008). [CrossRef]

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

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

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

19. R. Navarro, “Adaptive model of the aging emmetropic eye and its changes with accommodation,” J. Vision 14(13), 21 (2014). [CrossRef]

20. W. Smith, Modern Optical Engineering (McGraw Hill, 2007).

21. M. Abramowitz and I. A. Stegun, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, (Dover, 1972), chap. 3.

22. E. W. Weisstein, “Quartic Equation,” in MathWorld–A Wolfram Web Resource. https://mathworld.wolfram.com/QuarticEquation.html.

23. R. Navarro and S. Baquedano, “Numerical implementation of GRINCU model of the crystalline lens,” figshare (2021), https://doi.org/10.6084/m9.figshare.14865105

24. R. Navarro, F. Palos, and L. González, “Adaptive model of the gradient index of the human lens, II: optics of the accommodating aging lens,” J. Opt. Soc. Am. A 24(9), 2911 (2007). [CrossRef]

25. C. Qiu, B. Maceo Heilman, J. Kaipio, P. Donaldson, and E Vaghefi, “Fully automated laser ray tracing system to measure changes in the crystalline lens GRIN profile,” Biomed. Opt. Express 8(11), 4947–4964 (2017). [CrossRef]

GRINCU lens with conicoid iso-indicial surfaces: application for modeling the crystalline lens

Abstract

1. Introduction

2. General formulation

2.1 Meniscus lens

2.2 Iso-indicial surfaces

2.3 General conicoid

2.4 Biconvex or biconcave lens

3. GRINCU model of the crystalline lens

3.1 Numerical implementation

3.2 Results

4. Discussion and conclusions

Funding

Disclosures

Data availability

References

Supplementary Material (1)

Data availability

Cited By

Figures (5)

Equations (18)

Optics Express