1. Introduction

Gradient index (GRIN) appears in a wide variety of natural, such as eye optics in most mammals, and artificial optical systems. A paradigmatic example of GRIN optics is the crystalline lens in the human eye. Its inner shell structure, formed by successive layers of tissue, resembles an onion. Each layer can be realized as a thin meniscus lens with a refractive power given by the product of its surface curvature and the refractive index increment relative to the previous layer. In this context, two main types of layers have been proposed in the literature, continuous (differential) [1] and discrete (shell)—although most models consider a continuous distribution refractive index [2–9].

Ray tracing computation in continuous GRIN lens models is typically approached by solving the Euler–Lagrange equation for the ray trajectories [10]. Paraxial ray tracing in particular is essential to determine the first-order properties of an optical system. The cardinal elements (i.e., focii, principal planes, nodal points) and the lens power of a system are the most important parameters. Thus, paraxial ray tracing is the first and most crucial stage in analyzing an optical system, and it has applications in different disciplines ranging from optical design or engineering to clinical practice. Most diagnostic procedures, such as lens prescription, rely on paraxial optics almost exclusively. The ray transfer matrix formalism has the advantage over other methods, such as finite ray-tracing in the paraxial region, to offer an analytical approximation that lowers the computational demand.

Among the different approaches in the literature, Smith and Atchison [11] used a parabolic approximation for the ray trajectory inside the GRIN lens to obtain the equivalent lens power, and Pérez et al. [12] computed a gradient parameter to trace axial and field rays through the lens. More recently, Díaz [13] proposed the powerful and elegant ray transfer (or ABCD) matrix formalism for paraxial ray tracing through different models of the GRIN lens. The two main approaches consider the crystalline lens either a continuous slab [14] or formed by a finite number of shells [15]. Shell-type lenses can be pictured as a usually high number of cemented meniscus lenses, and ray tracing is done via the onion-type layered structure of the lens, not relying on other assumptions such as the parabolic approximation, which is fundamental for the continuous slab method [11,13]. The main drawback of this method is that the computed ABCD matrix is the product of typically hundreds of 2 by 2 matrices.

Our goal is to unify the shell and slab approaches to formulate the ray transfer matrix of onion-type GRIN lenses in the continuous limit– when the layer thickness tends to zero, $\Delta z \rightarrow 0$, and becomes a differential magnitude, $dz$. In these onion-type lenses, either with finite shells or continuous (differential) structure, we can consider each iso-indicial surface in the refractive index distribution as a single refracting optical surface with a known curvature radius and refractive index value. This enables a simple, compact, robust, and general formulation applicable to various GRIN lenses.

The GRIN nature of the crystalline lens has long been known, and even Gullstrand [16] proposed a GRIN lens model a century ago. Experimental studies suggest that the lens cortex is strongly inhomogeneous, with a significant positive gradient towards the center, whereas the nucleus is much more homogeneous [17–20]. Most advanced eye lens models attempt to capture this refractive index distribution. Here we classify models into two groups, slab-type, [2,21] and onion-type [3–9]. The main difference between them is that the onion-type structure is explicit: the analytical expression of the whole family of iso-indicial surfaces, including their curvature radii, is provided. Our method is exceptionally well suited for those types of models, which indeed are among the most recent and advanced. We show below that it is possible to derive an analytical ray transfer (ABCD) matrix for these models, enabling us to perform analytical paraxial ray tracing.

2. Theory

Figure 1 represents a section of the iso-indicial surface distribution of a crystalline lens model [9] as an example of an onion-type biconvex lens. This structure suggests a cemented multiplet formed by a high number of meniscus lenses.

 

Fig. 1. A vertical slice of the iso-indicial surfaces of a 30-year-old GRIN crystalline lens model. The colors indicate refractive index values, G is the curvature gradient parameter, and z0 is the iso-indicial parameter defining each surface.

Download Full Size | PDF

2.1 Ray transfer matrix for onion-type GRIN lenses

Let us formulate the paraxial ray tracing through a GRIN lens structure with a refractive index $n(\omega,\theta, z)$ in cylindrical coordinates. We also apply the standard convention that the light propagates from left to right, that is, towards increasing $z$. We assume that the refractive index depends only on $z$ in the paraxial region, $n\simeq n(z)$. For a layer $k$ we have a iso-indicial surface crossing the axis at a location $z_k$ from the anterior surface. Here the main assumption is that the ray propagation within this layer is described by the corresponding ray transfer matrix that is the product of a translation and a refraction matrix. The translation is equal to the thickness of the layer $\Delta z$, and the refraction occurs at the iso-inidicial surface. Matrices are multiplied in reverse order of the ray path, so as we progress along $z$, each new matrix is added to the left [22]. Thus, for this layer $k$, we have:

We assume that $n'$ is also a function of $z$ in the paraxial region $n'=n'(z)$. The residue $\mathscr {O}(\Delta z)^2$ contains terms of $\Delta z$ to the second power or higher. When the displacement tends to zero, $\Delta z \rightarrow 0$, so does the residue, $\mathscr {O}(\Delta z)^2\rightarrow 0$, so this term can be neglected in a first-order approximation. We can further simplify by rewriting $\Delta n=n'\Delta z$ and approximating the fractions inside the matrix, neglecting again any terms of the order $\mathscr {O}(\Delta z)^2$ that are incorporated in the residue. The resulting matrix is the sum of the identity matrix $\boldsymbol {I}$, a new matrix $\boldsymbol {A}$ that multiplies $\Delta z$, and the residue:

This new matrix $\boldsymbol {A}$ is independent of $\Delta z$. For any $z$:

We approximated this product to a sum by neglecting higher-order powers of $\Delta z$ in all cross terms $\boldsymbol {A}(z_i)\boldsymbol {A}(z_j)(\Delta z)^2$, and substituting $\boldsymbol {IA}(z_k)\Delta z$ with $\boldsymbol {A}(z_k)\Delta z$.

2.1.1 Continuous case

To formulate GRIN lenses in the continuous limit, we replace $\Delta z$ with $dz$ and the sum in Eq. (4) with an integral:

Here, the anterior lens vertex is at $z=0$ and $t$ is the lens thickness. This result is elegant and exciting: the ray transfer matrix $\boldsymbol {M}$ of the whole GRIN distribution can be found by integrating the elements in the differential matrix $\boldsymbol {L}(z_k)$, Eq. (1).

Let us consider the case of a meniscus lens. When the curvature radius is constant, $R(z)=R$, the iso-indicial surfaces are parallel. In this case, the defined integral in element $M_{2,1}$ has a direct solution:

If the refractive index at both surfaces is equal, $n(t)=n(0)$, then $M_{2,1}=0$ and $M_{2,2}=1$. In these conditions, $\boldsymbol {M}$ represents a pure translation, and the GRIN inner structure does not contribute to the lens power.

In the case of a homogeneous biconvex lens, we can picture it as a sort of cemented doublet, having a positive radius $R_a$ for the anterior side of the lens and a negative radius $R_p$ for the posterior. Let us call the thicknesses of the anterior and posterior regions $t_a$ and $t_p$ (with $t=t_a+t_p$). Now we can write the integral in element $M_{21}$ as the sum of the two sections, which results in $M_{21}=\log (n(t_a)/n(0))/R_a + \log (n(t)/n(t_a))/R_p$. With $n(0)=n(t)=n_s$ at the surface, and a refractive index $n_c=n(t_a)$ at the center of the lens, it is straightforward that $M_{21}=\left (1/R_a-1/R_p\right )(\log (n_c) - \log (n_s))$. This expression is equivalent to the power of a thin lens in air where the refractive indices have been replaced with their logarithms. Note that in this case, when the refractive index is the same at the anterior and posterior surface $n(0)= n(t)$, and $M_{22}= 1$.

2.1.2 Surface contribution

Matrix $\boldsymbol {M}$ corresponds to the ray propagation inside the lens where the refractive index $n(z)$ is continuous and has a derivative $n'(z)$. Because of the abrupt change in refractive index, we formulate the refraction at the anterior and posterior surfaces separately. Now the lens matrix can be expressed as the product of three matrices:

Matrices $\boldsymbol {S_a}$ and $\boldsymbol {S_p}$ correspond to pure refraction at the anterior and posterior surfaces (there is no displacement between the surface and the inner structure). Subscripts $a$ and $p$ refer to the radius and refractive index values at the anterior and posterior surfaces ($z=0$ and $z=t$, respectively). The surrounding external media (anterior and posterior sides) have refractive indices $n_{0a}$ and $n_{0p}$.

For simplicity, we assume that the anterior and posterior surrounding media have equal refractive indices, $n_{0a}=n_{0p}=n_0$, and so do the anterior and posterior lens surfaces, $n_a=n_p=n_s$. The matrix product yields:

The main theoretical result of this work is the ray transfer matrix of a biconvex GRIN lens comprised of infinitesimal shells or layers –each characterized by its refractive index and curvature radius–. The most relevant element in this matrix is the focal length’s inverse, $L_{21}=1/f'$, which determines the lens power. The first term of $L_{21}$ mirrors the power of a standard thin lens with a refractive index equal to the surface index value $n_s$. The second term is associated with the power of a thick lens. And the third term containing the integral is the contribution of the GRIN inner structure to the lens power. It is a simple additive law. The remaining elements of $\boldsymbol {L}$ parallel those of a homogeneous lens, $\boldsymbol {H}=\boldsymbol {S_p T S_a}$, where the GRIN contribution is replaced with a translation matrix $\boldsymbol {T}$. To convert the ray transfer matrix of a homogeneous lens $\boldsymbol {H}$ into a GRIN lens matrix $\boldsymbol {L}$, we only need to add the contribution of the GRIN structure (integral term) to $H_{21}$.

This theoretical result is significant as it provides a straightforward and intuitive formulation. When $n(z)$ is continuous and differentiable and $n'(z)$ exists in the interval, the integral can be solved numerically, for instance, by applying the trapezoidal rule. In some cases (see below), it can be solved either analytically or with a close analytical approximation.

Once we have the characteristic $ABCD$ matrix, we can compute physical magnitudes such as the power, $\phi _L$, the distance from the anterior lens vertex to the principal object plane, $d_H$, and the distance from the posterior vertex to the principal image plane, $d'_H$ [22].

3. Application to the crystalline lens

The human lens is a paradigmatic example of this onion-type GRIN distribution. Most experimental measurements [17–20] and models [4–6,8,21] suggest that the axial distribution of the anterior and posterior regions follows a power law. One of the most recent models is the GRINCU lens [9]. In this model, the anterior and posterior refractive index distributions are given by Eq. (13) in [9], expressed in terms of the parameter $z_0=z_0(\omega, \theta, z)$ that describes the iso-indicial surfaces for different values of $n$. In the paraxial region, we can assume that the index distribution depends only on $z$, so $z_0=z_0(0,0,z)$:

The GRINCU lens is characterized by having a curvature radius gradient parameter $G$ in addition to the refractive index gradients. This model assumes the simplest mathematical formulation, which is to use a constant G value. The apical curvature radii of the iso-indicial surfaces become:

3.1 Numerical Implementation

Equations (13) and (12) are explicit expressions of $R(z)$ and $n(z)$ in terms of $z$, and we can easily compute $n'(z)=dn/dz$ from Eq. (12). These three functions make up the integral in Eqs. (9) and (10), which express the GRIN contribution to the lens power as an explicit function of $z$. To compute the integral, we split it into the anterior and posterior contributions and multiply by the refractive index of the surrounding medium $n_0$.

Table 1. Fixed variable values used for the integration in Fig. (2)

View Table

Figure 2 shows the refractive index axial distribution $n(z)$ and the contribution of each layer to the total lens power, as an approximate differential $dP(z)/dz$, for different values of the gradient parameter $G$. In this example (a 30 years old GRINCU lens model), we used a finite value of $dz=36\mu m$. The integration of $dP(z)$ has no analytical solution in general, so a numerical integration was implemented based on the trapezoidal rule. The resulting total lens powers $P$ for this example are included in the legend. The convergence obtained with the trapezoidal rule is shown in Fig. 3 and discussed below.

 

Fig. 2. The convex-up solid line shows the variation of $n$ along $z$ in the paraxial region. The concave-up lines represent $dP(z)$ for different G values. See Table (1)

Download Full Size | PDF

 

Fig. 3. Power of the lens model against the number of steps for two methods: matrix product and numerical integration. Percent difference with the last point ($N = 33,000$) is shown for three values of $N$.

Download Full Size | PDF

3.2 Exact integration

A particular case of interest is that of $G=0$, a homogeneous biconvex lens. In the absence of a curvature gradient, the iso-indicial surfaces are parallel in the anterior and posterior regions and Eq. (14) has an exact analytical solution: $R_a$ and $R_p$ can be factored out of the integrals, which are resolved into logarithms:

This expression resembles that of a thin biconvex lens in which the difference of refractive indexes is replaced by the difference of their respective logarithms (multiplied by the surface index).

3.3 Analytical ray transfer matrix

The integrals in Eq. (14) cannot be solved analytically; however, we worked out a highly accurate second-order approximation in terms of hypergeometric functions. We discuss the validity and accuracy of this approximation below, and the step-by-step mathematical derivation is provided in the Supplemental Material. The approximated expressions obtained for the integrals are:

4. Results and discussion

The mathematical ray transfer matrix formulation derived here is applied to analyze the effect of the curvature radius gradient parameters $G$ and $Q$ on the lens refractive power. The results are plotted in Figs. 4 and 5, respectively. We observe similar non-linear dependencies of the power on both parameters but with opposite trends: Increasing $G$ decreases power, whereas increasing $Q$ increases power, which means that ellipsoid surfaces ($Q > 0$) deliver more power than hyperbolic surfaces ($Q < 0$). The reduction of GRIN lens power with increasingly negative $Q$ values has been previously reported in the literature [5], but because the paraxial power is unaffected by conic constants in homogeneous lenses, the authors found the result non-intuitive. This new formulation provides a relatively simple and direct explanation to this effect, as the curvature radius gradient is proportional to $Q+1$, as described in Eq. (13). We believe further study is needed to fully understand the relationship between the conic constant and the paraxial lens power.

 

Fig. 4. The lines represent $P(G)$ calculated with different methods, the anterior and posterior conic constants are $Q_a=-4, Q_p=-3$.

Download Full Size | PDF

4.1 Methods assessment

We have compared the numerical implementation of the present method with the original direct matrix product in terms of convergence. In our implementation we used the matrices $\boldsymbol {L_k}(z_k)$ of Eq. (2) for our matrix product, with $k = 1, 2,{\ldots }N$ (unlike Giovanzana et al., [15] who used a normalized version of these matrices). The results of both methods are compared in Fig.3. The horizontal axis represents $N$, the number of integration steps or matrices (shells) used to compute power. The vertical axis is the total lens power versus the number of steps used. The dashed blue line represents the numerical integration method, and the continuous red line corresponds to the matrix product. The two methods provide practically the same result, with a difference of 0.06 D (0.2 %).

However, convergence is considerably faster with the integration method, which is strongly enhanced by the trapezoidal rule [24]. With only ten steps, the error relative to the convergence value is around $+ 0.5 D$. In contrast, the matrix product provides a much higher error, nearly $-8 D$ (consistent with the trend reported in Fig. 1 in Ref. [15]). With the numerical integration method, we reach a reasonable convergence with less than 50 steps; with 100, the error is negligible, and the result is stable. Conversely, the matrix product requires a much higher number of steps to converge and is only stable after around 10,000 steps. One possible explanation is that there is no equivalent to the trapezoidal rule for matrix multiplication.

Figures 4 and 5 summarize the results of different refractive power matrix computation methods compared with the thin lens approximation (dotted red line) used by several authors before [8,9,25]. The thin lens approximation involves a more straightforward formulation and does not require matrix formalism. The methods proposed in this work are labeled as analytical (continuous black lines), numerical (dashed green lines), and exact (red asterisk) for G = 0. The agreement between them is almost perfect. The matrix product (dashed blue line) shows a nearly constant very small underestimation of 0.06 D relative to the numerical integration. Our analytical and numerical methods give powers of about 0.6 D less than a thin lens.

 

Fig. 5. The lines represent $P(Q_a)$ calculated with different methods. The posterior conic constant is $Q_p=Q_a+1$

Download Full Size | PDF

For comparison purposes, we also implemented the ABCD matrix slab method (point dashed purple lines) based on the parabolic approximation [13]. This method slightly underestimates power compared with the matrix product and the differential shell method proposed here. In contrast with the other methods, the line is not generally parallel (see Figs. 4 and 5). There is a closer match for more negative $Q$ values (minimum $0.67\%$ difference) and higher $G$ values (minimum $0.40\%$ difference). On the contrary, the more significant difference of $1.4\%$ corresponds to non-negative $Q$ values (see Fig. 5). Nevertheless, the differences are consistently below $0.4 D$. The agreement between these two methods is reasonable, especially for the typical negative conic constant values used in lens models.

The GRIN matrices that we obtained for the matrix product $\boldsymbol {M_p}$, our numerical integration $\boldsymbol {M_n}$, and the slab $\boldsymbol {M_s}$ methods are:

We see a close agreement in element $M_{21}$, which is proportional to the lens power, but slight differences in the other three elements. In our approximation, the diagonal elements are precisely 1, and the element $M_{12} = -t$ equals the lens thickness. In the other two matrices, these elements differ from ours. The most significant deviation seems to correspond to the "equivalent" thickness, element $M_{12}$, which is nearly $2.6\%$ thinner in the matrix product. We conclude that, excluding the thin lens approximation (dotted red lines in Figs. 4 and 5), all the ray transfer matrices compared here provide a reasonable agreement. The ABCD matrices are nominally unitary and thus their determinants should be equal to one. This is true for $M_p$ and $M_s$, but due to the GRIN contribution to the homogeneous lens described by the result in Eq. (10), the determinant of $M_n$ is slightly greater than one.

4.2 Application to other GRIN models

As mentioned above, the theory and numerical implementation of this matrix formulation can be general not only for the human lens but for other GRIN lenses, provided that:

To demonstrate the generality of our method, we apply it to the AVOCADO crystalline lens model [8], one of the most advanced in the literature. Another important reason we chose the AVOCADO model is that the authors provide explicit results of lens power computation and compare the ray tracing with the thin lens formula in Fig. 3 of Ref. [8]. We have reproduced their results employing our method, as shown in Fig. 6. Comparing our plots and values with theirs, we see that they are almost identical; our model replicated both the thin lens and ray tracing lens powers with high accuracy.

 

Fig. 6. Lens power vs. $m$ for different values of $p$. The solid lines are calculated with the thin lens approximation of the integral, and the dotted lines correspond to the numerical calculation.

Download Full Size | PDF

5. Summary and conclusions

The main result of this work is theoretical: the ray transfer matrix of the GRIN media is similar to that of a single layer, where matrix elements $A, B, C$, and $D$ are replaced with their corresponding integrals. This substantial simplification enables an affordable and straightforward computation of power and other first-order properties of the optical system.

Furthermore, the result of Eq. (10) involves an even more crucial simplification. We can rewrite the ray transfer matrix of a GRIN lens as

Equation (19) shows that the lens transfer matrix is the sum of the ray transfer matrix of a standard homogeneous lens and a single-element matrix $\boldsymbol {G}$ containing the contribution of the GRIN inner structure. All elements of $\boldsymbol {G}$ are equal to zero except $G_{21}$. This single-element matrix is even simpler than a pure refraction matrix. Therefore the initially complex product of a high number of matrices becomes a much simpler computation with straightforward additive rules and integrals, and the triple matrix product in Eq. (7) becomes the sum of matrices in Eq. (18). The apparent complexity of the GRIN paraxial ray tracing has motivated many authors to use the idea of an equivalent refractive index [26–28], a theoretical equivalent refractive index $n_{e}$ so that the corresponding homogeneous lens matrix is equal to that of the lens: $\boldsymbol {L}=\boldsymbol {H'}$. But after Eq. (18), it is obvious that $\boldsymbol {H'}=\boldsymbol {H} + \boldsymbol {G}$. This simple additive rule for the GRIN contribution renders the use of an equivalent refractive index unnecessary.

Another exciting and general result is that when the curvature radius $R$ is constant, the GRIN does not contribute to the refractive power of the lens. This only happens in meniscus lenses with equal curvature radii. In biconvex lenses, however, the radii of the anterior and posterior parts have opposite signs, so there is a GRIN contribution to the refractive power—that of a conventional thin lens with equivalent refractive power but replacing the refractive index difference with the difference of their logarithms (see Eq. (15)).

We found that the curvature radius gradient, $G$, and the conic constant, $Q$, strongly impact the lens power in crystalline lens models. Negative values of $Q$ (hyperbolas) give less power, and as $Q$ becomes more positive, power increases. We have the opposite effect regarding the gradient parameter $G$ (Fig. 4): as $G$ increases, the lens power decreases.

Finally, we believe that the theory and methods developed here provide a powerful, entirely general, and practical (robust and easy to use) tool for paraxial ray tracing, applicable to both GRIN lens models and the analysis of experimental data. In particular, the analytical approximation is especially well suited to assess and model age-related or accommodative changes in the crystalline lens, which will be the subject of future work.