## Abstract

Ray tracing in gradient-index (GRIN) media has been traditionally performed either by using the analytical or numerical solutions to the Eikonal equation or by creating a layered medium where Snell’s law is calculated in each layer. In this paper, an exact general method to perform ray tracing in GRIN media is presented based on the invariants of the system as stated by Fermat’s principle when the media presents symmetries. Its advantage, compared with other methods reported in the literature, relies on its easy implementation. Besides the GRIN distribution and the initial conditions of the incident ray, once the invariants of the system are stated the resulting math is simple to solve and interpret. To benchmark the algorithm, ray tracing in typical cases of GRIN media is calculated, finding minimal discrepancies between the analytical solutions and our simulations. The used media are axial refractive index and parabolic index fiber and lenses with spherical gradient-index symmetry, such as: Luneburg’s, Gutman’s, generalized Maxwell’s Fish-eye, Eaton’s, and concentrator lenses. Our method can be further applied to distributions with symmetries associated with other common curvilinear orthogonal coordinate systems, in particular to those associated to the separability of the Helmholtz equation that would allow us to investigate wave optics in these GRIN media with the associated geometries.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Inhomogeneous materials possessing a gradient in their refractive index are known as gradient-index (GRIN) media [1,2]. GRIN media are commonly investigated considering certain geometries of its isoindicial surfaces (surfaces of constant refractive index). For instance, investigations have been done mainly in inhomogenoeous media whose gradient can be axial, cylindrical and spherical. In the first case, the isoindicial surfaces are planes, for the second case they are concentric cylinders about the optical axis, whereas for the third kind the isoindicial surfaces are concentric spheres.

The propagation of light through a medium in geometrical optics is ruled by the ray equation [3]. For a GRIN medium, rays travel through non-linear paths, while for the simplest case, i.e., homogeneous media, rays describe straight lines [3]. It is possible to find closed analytical solutions to the ray equation in some particular cases [4–10], but for the general case it is required to use numerical methods [11–13]. These kind of solutions require discretizing the domain and then, to implement an iterative calculation of the ray paths. Different integration techniques might involve numerical methods like finite element or quadratures, where segments of rays are calculated between interfaces [14]. Some other numerical methods consider a discretization (grid) of the continuous medium and the path is then constructed by calculating Snell’s law at each interface, or point of the grid [15]. A similar approach is used to compute rays under the paraxial approximation, dividing the domain into thin vertical slices, where the refractive index of each slice is replaced by a quadratic medium, in which the ray is traced by means of an ABCD matrix [16]. In 1988, an interesting ray-tracing method for acoustic surfaces with smooth large-scale surface irregularities was proposed by Krylov [17]. The irregularities were in the form of axisymmetrical recesses, or protuberances, and cylinders of variable curvature. This method uses three equations to calculate the propagation of the rays: two differential equations and a generalized Snell’s law as well for inhomogeneous media with spherical symmetry [18,19]. However, in some cases, depending on the GRIN distribution, obtaining such a solution results rather complicated and it is often necessary to make several approximations to find a solution for the ray trajectory [20].

Most remarkably, it should be noted that the generalized Snell’s law for inhomogeneous media with spherical symmetry involves a conserved quantity, which comes from Fermat’s variational principle. Conserved quantities are useful in solving several problems, since their complexity is usually reduced. Now, it can be useful to recall that a symmetry present in the Lagrangian of a system has an associated conservation law, as Noether’s theorem states [21]. The Hamiltonian and Lagrangian formalisms have been developed and implemented in different versions along the years being the first approaches applied to electron beams [3], in more recent investigations the analogy has been developed [22–25], and ray tracing algorithms have been implemented [26–28].

In this paper, we show that the conserved quantities induced by the symmetries of the system, henceforth Fermat’s ray invariants, allow to establish a simple exact numerical method for ray tracing. Due to the elegance of the principle on which it is based, the method can be easily implemented. A remarkable advantage of this method is that it is not necessary to find an explicit solution to the ray equation, only the GRIN distribution, and the corresponding invariants, are required. This represents an advantage over the aforementioned methods which work only for specific cases. Another important advantage over most of the non-exact numerical methods is that the discretization is not tied to the symmetry of the refractive index. To show its power, the algorithm is implemented in different important cases with rectangular, cylindrical and spherical symmetry. The cases included are axial refractive index lenses, parabolic index fiber and Luneburg’s, Gutman’s, Generalized Maxwell’s Fish-eye, Eaton’s, and Concentrator lenses. In the examples provided, the ray paths are not calculated by integration, but by calculating the direction of the rays, at every step of the path, as dictated by Fermat’s ray invariants. In every case, the invariant is calculated at the beginning of the path through the medium, and then used to compute the rest of the path. This method is simple, elegant and does not require any approximation, except for the computational step. Even though analytical solutions can be obtained for GRIN distributions with spherical or cylindrical symmetry [29], in general it is not an easy task.

## 2. Fermat’s ray invariants

Variational Calculus provides a very nice and elegant description of fundamental laws in Physics. When studying a physical system it involves finding the extremum of an integral of a function that depends on one or several independent variables, dependent variables and their derivatives. This function in Physics is commonly referred to as Lagrangian. The requirement of finding the extremum of the integral yields to the Euler-Lagrange equations of the system [30]. For instance, in Classical Mechanics such quantity is the action, while in the calculation of geodesics it is the distance between two given points on a given surface, whereas in Geometrical Optics it is the optical path length of a given light ray path in a medium.

A physical system can be described by a respective Lagrangian that carries all the dynamical information of the system. For the simplest case of Classical Mechanics, in which there is only one single independent variable, the time, the Lagrangian is written in terms of the kinetic and potential energies of the particular system and the extremum is determinend by the Hamilton’s Principle of least action. For geodesics it is written in terms of the element of length of the corresponding space and the extremum determines the shortest path on a given 3D surface. And for the optical case the Lagrangian is written in terms of both the refractive index and the element of length in a given coordinate system. The final form of the Lagrangian will be determined by the 3D coordinate system used.

Assuming that the refractive index of the medium is inhomogeneous, its value depending on the position, its gradient will be non-zero. This kind of medium is referred to as a medium with gradient refractive index (GRIN). In this case, the extremum will be given by the path that makes the ray of light traveling in the minimum time. Fermat’s principle for a path $\mathbf {r}(s)$ taken by a ray of light passing from A to B in three-dimensional space can be expressed as

where $C$ represents the curve joining the points A and B and $\textrm{d}s$ is the element of length along the curve. Considering generalized coordinates associated to an arbitrary curvilinear orthogonal coordinate system, $\mathbf {s}=\mathbf {s}(q_{1}, q_{2}, q_{3})$, and using the property of invariance under reparameterisation, $(q_{1}(\tau ), q_{2}(\tau ), q_{3}(\tau ))$, Eq. (1) can be rewritten as [31]The solution of this problem must satisfy the Euler-Lagrange equations, which in this case are given by

From the Euler-Lagrange Eq. (4), it is possible to obtain two important invariants. If the Lagrangian (3) does not depend on a given generalized coordinate $q_i$, Eq. (4) yields

where $K_1$ is a constant related to a conserved quantity that we will call from now on a Fermat’s ray invariant.There is another invariant obtained by taking the total derivative of the Lagrangian $\mathcal {L}$ with respect to $\tau$ and substituting in the corresponding Euler-Lagrange equation, namely

When the Lagrangian does not depend on the independent variable $\tau$ we immediately get

Equations (5) and (7) are the invariants that we will use to study ray propagation in GRIN media.

In the following sections, we will apply the above invariants to axial, radial cylindrical, and spherical GRIN geometries. We will start considering the spherically symmetric case in great detail for a variety of GRIN media.

## 3. Physical ray tracing in spherical lenses

The spherical symmetric case is the most general and allows for a clear explicit definition of all the parameters involved in our description. So we will start our discussion with this geometry. Let us consider a spherical GRIN lens of radius $R$, with a GRIN described by a function $n\left (r\right )$, where $r=\sqrt {x^2+y^2+z^2}$ and $x$, $y$, and $z$ are Cartesian coordinates. This function represents a GRIN distribution where its maximum value is located at the center of the lens, and decreases from it to its minimum, on the lens surface.

For media possessing radial symmetry, the Optical Lagrangian in spherical coordinates $\left (r,\theta ,\phi \right )$ obtained from Eq. (3). It can be shown that the rays are confined to a plane [18]; by choosing the plane to be the equatorial plane $x-y$, for which $\theta =\pi /2$, the Lagrangian reduces to

where $\phi$ is the dependent variable, $r$ the independent one and $\phi _{r}=\textrm {d}\phi /\textrm {d}r$. Observe that the Lagrangian does not dependend on $\phi$, then by Eq. (5) we get where $K$ is a constant. Substituting the Lagrangian and using the geometry shown in Fig. 1, after some algebra we find that where $\varphi$ is the angle between the position vector $\textbf {r}$ of the ray path, and its unit tangent vector $\textbf {s}$, see Fig. 1.Equation (10) is known as the generalized Snell’s law for inhomogeneous media with spherical symmetry that we identify as a Fermat’s ray invariant, this relation is satisfied by each ray inside the lens. In other words, each ray possesses a definite value of $K$ determined by the entrance point and its value does not change throughout its propagation [29]. Also, this equation implies that the ray paths in a medium defined by $n=n(r)$ are planar curves contained in a plane that passes through the origin [3,18].

Without loss of generality, consider the following boundary conditions, the maximum radius of the lens is $R=1$, and $n\left (R\right )=1$, which are analogous to consider a normalized coordinate $r$ and to consider the value of the GRIN at the boundary of the lens to be normalized as well with respect to the surrounding refractive index in which the GRIN is embedded [34]. If we consider these boundary conditions, it is easy to calculate the variation range of $K$, due to the fact that the maximum point of incidence can be obtained when $r=1$, $n=1$ and $\varphi =\frac {\pi }{2}$. At this point, from Eq. (10), we have that $K=1$ and the minimum point is when we are on the axis of propagation, where $\varphi =0$ and in consequence $K=0$, then the range of $K$ is

#### 3.1 Algorithm

The values of $K$, for each ray, can be calculated on the lens surface ($R=1$) and the parameters necessary to calculate them are shown in Fig. 1. From this Figure, it is possible to see that the incident ray is originated at the point $P_{0}=\left (x_{0},y_{0}\right )$, while $P_{i}=\left (x_{i},y_{i}\right )$ indicates the point where the incident ray enters at the spherical lens, whereas the point $P_{c}=\left (x_{c},y_{c}\right )$ is located at the center of the lens and at the same time is the origin of the coordinate system, where $x_{c}=0$, and $y_{c}=0$. Using these three points, the value of $K$ for the most general case can be calculated as

The importance of calculating the value of $K$ is that if it is known at the moment the ray enters the lens, it is possible to know the path of the entire ray inside the lens. This is because the value of $K$ remains constant over the entire path of the ray. For this reason, the principal objective of the algorithm is to preserve the value of $K$ along the ray trajectory. To achieve this objective, it is necessary to calculate the value of $\varphi$ at each instant of the ray propagation. From Eq. (10) and Eq. (12), the value of $\varphi$ along the ray trajectory is given by

Consider a ray coming from the point $P_0$ that impinges at the boundary of the lens ($r=1$), at $P_{i}=\left (x_{i},y_{i}\right )$ with an angle $\alpha _0$, as shown in Fig. 2. From Eq. (13), the value of the angle $\varphi$ at $P_i$ is given by $\varphi \left (r=1\right )=\sin ^{-1}{(K)}$.

Now let us trace, from the point $P_i$, a straight line of length $\Delta \mathbf {s}$ with angle $\gamma =\pi -\mu _0-\varphi$ formed between $\mathbf {s}$ and the horizontal, (see Fig. 1 and 2). Then, the ray after this step will reach the point $P_{i+\Delta }$ given by

Notice that depending on the GRIN of the lens (distribution of refractive index in the shells), the rays can propagate around the origin, traveling backwards or spiraling towards the center of the lens as will be shown below in the examples. This is an advantage over other ray tracing algorithms whose steps are such that the rays usually only traveling in the forward direction and are incapable of getting the rays traveling backwards.

We have just elaborated the algorithm of our work; once we know the position of the point $P_{i+\Delta }$, by means of Eq. (13) we are able to calculate the angle of the new step and trace another straight line in the direction given by the corresponding angle $\gamma$.

The stages of the algorithm are presented in Fig. 3. We call this method Physical Ray Tracing (PRT), since it is based on exploiting Fermat’s ray invariants.

#### 3.2 Numerical implementation for different spherical GRIN lenses

To test the algorithm presented in Fig. 3, a ray tracing through of a set of GRIN lenses with spherical symmetry has been performed. The GRIN distribution of each of these lenses is shown in Table 1. To clearly observe the variation of the refractive index on the propagation axis, the set of lenses have been divided into three subsets: The first subset consists of the Gutman lens [35], the second subset contains both the Eaton [36] and Concentrator lens [37], and the third subset consists of the Generalized Maxwell’s Fish-Eye lens [38,39]. We emphasize that the Luneburg lens [29] is included in the first subset, because the Luneburg lens is a special case of the Gutman lens for $f=1$. Besides, the Maxwell’s Fish-Eye lens [40] is included in the third subset, because it is a special case of the Generalized Maxwell’s Fish-Eye lens when $m=1$.

As discussed above, in spherical GRIN lenses the rays are confined to a plane, thus we will use the the $x-y$ plane. The $x-$axis will be the propagation axis and the $y-$axis the vertical axis. It is worthy to note that, once the ray tracing is performed in that plane, due to the spherical symmetry of the GRIN, any rotation with respect to any axis passing through the origin of the lens (origin of the coordinate system), will yield an actual ray path, which permits to generate a three-dimensional ray tracing [41].

In the first and second subsets mentioned above, the value of $K$ is given by [42]:

since the incident rays come from infinity, that is, $P_{0}=\left (\infty ,0\right )$, and $P_{i}=\left (x_{i},y_{i}\right )$.Let us now focus on the first subset. From Fig. 4(a), it is possible to observe that rays coming from infinity into a Gutman lens with $f=0.75$ are focused inside the lens, precisely at $x=0.75$. On the other hand, $f=1$ represents a special case: the Gutman lens becomes a Luneburg lens; the incident rays are focused on the lens surface, as shown in Fig. 4(b). It is remarkable to note that Gutman’s solution for $f>1$ has spherical aberration [43,44], i.e. the lens is not stigmatic, as shown in Fig. 4(c), where $f=1.25$.

On the other hand, for the second subset, the Eaton lens has the remarkable ability to bend the ray path to such a degree that it turns around [45] as it is shown in Fig. 4(d). Besides, the Concentrator lens, which also belongs to the second subset, is such that all the light rays are “absorbed” towards its center [37], as shown in Fig. 4(e). Notice that, as mentioned above, our algorithm is able to trace rays that make loops.

For the third subset, the value of $K$ can be written as

since the incident ray parameters are given by $P_{0}=\left (x_{0},y_{0}\right )$, and $P_{i}=\left (R,0\right )$, i.e., all the incident rays enter the lens at the point $P_{i}=\left (R=-1,0\right )$, see Fig. 1.From Table 1 it is possible to observe that the Generalized Maxwell’s Fish-Eye lens expression depends on the parameter $m$. If $m=1$, the traditional Maxwell’s Fish-Eye lens equation is recovered (see Table 1). The ray tracing in this lens is shown in Fig. 4(f), in which it can be noticed that the incident rays, which come from the surface of the lens, focus at the opposite point on the surface, as the theory predicts [40]. A special case is given for $m <1$, in which two foci are generated in the plane of the rays but it is in fact a ring-focus on the surface of the spherical lens. As $m$ decreases, the foci separate, as it is shown in Figs. 4(g), (h) and (i).

#### 3.3 Benchmarking the PRT method

Now, to verify the validity and precision of our physical ray tracing algorithm, we follow the common procedure of applying it to a problem for which there exists a known analytical solution and comparing the numerical method. We use the well studied problem of the Luneburg lens and compare it with the analytical solution given by J. A. Grzesik [46],

From Fig. 5, it is possible to observe that our results are in good agreement with the analytical solution. As expected, the numerical error is step dependent. However, it can also be observed that in the worst case ($\Delta \mathbf {s}=0.1$), the error is of $10^{-1}$ order. Figure 5(b) shows how as the step size decreases, the relative error also decreases with a comparable order of magnitude. Such degree of agreement is not surprising, since our method is a seminumerical method, whose accuracy depends on the value of $\Delta \mathbf {s}$.

After comparing several features of the RK4 method with the PRT we obtained the following results:

- 1. The RK4 method is restricted to the computational grid defined to apply the method, whereas the PRT requires only the function $n(r)$.
- 2. The PRT is perfectly capable of matching the RK4 precision albeit with higher computation times.
- 3. The computation time was measured in Matlab and the running times were always of the same order for each step size up to $10^{-4}$. For smaller steps, the RK4 method had a lower computation time than that of the PRT.
- 4. The accuracy of both methods is comparable. However, the RK4 method is more accurate up to a certain limit. The PRT method can achieve the same accuracy at the expense of more computational time. The accuracy of both methods was computed and shown in Fig. 5(c).
- 5. The PRT method shows stable results across simulations and computational steps. Although the RK4 method is stable along the ray path, it presents a divergent behavior beyond the point where the ray exits the GRIN medium. The PRT does not exhibit this problem.
- 6. The PRT method continues working outside the GRIN medium while the RK4 method has no validity there. For the RK4 method, the calculations brake down outside the medium, right after the interphase. The PRT does not have this drawback because the equation for the invariant $K$ remains valid outside the GRIN medium and the program yields straight rays outside the medium as expected, even if $n>1$ outside the GRIN.
- 7. For a step of $10^{-4}$ the PRT method performs with a $10^{-2}$ error, which is sufficient accuracy for experimental purposes.
- 8. The PRT has the remarkable quality of functioning with discontinuous media (with step variation) which is impossible for a typical numerical method, in a single stroke, i.e., the program would have to be run between interphases.

The RK4 method was selected among others for its power and universality. Different step sizes yield a higher accuracy for the RK4 method and comparable computation time, however, the PRT has a main advantage over standard numerical methods since the computations are intimately related to the physics, namely the Fermat’s and Noether’s theorems which are the subject of this work. The evidence shown is important to assert that the PRT method represents an advantage over other common numerical methods which solve the ray equation.

## 4. Physical ray tracing in axial refractive index lenses

Now let us consider axial GRIN media, whose isoindicial surfaces are planes that are parallel to the optical axis. In this case, the GRIN is represented by a function $n=n\left (y\right )$. The Optical Lagrangian for this axial medium is given by

where $y$ is the dependent variable, $x$ the independent variable and $y_{x}=\textrm {d}y/\textrm {d}x$. Noticing that the Lagrangian does not depend explicitly on $x$, the corresponding invariant is determined using Eq. (7), namely After substitution of the Lagrangian we getOn the other hand, from Fig. 6 it is possible to deduce that

Combining the last two equations yields

Equation (22) is known as the generalized Snell’s law for an inhomogeneous medium described by $n\left (y\right )$, where $K_L$ is the corresponding Fermat’s ray invariant.

The treatment of this case is analogous to that of Section 3.; for each ray, $K_{L}$ will be calculated on the surface of the medium, while $\varphi _{L}$ will be obtained through Eq. (22). To do so, let us consider that the ray starts from the point $P_ {0}=\left (x_ {0}, y_{0}\right )$, and that the ray enters the medium at the point $P_{i}=\left (0,0\right )$ (see Fig. 6). Then,

whileAs stated above, the GRIN only depends on the cartesian variable $y$, and our physical ray tracing algorithm can be applied using Eqs. (14), (22)–(24), as it is shown in Fig. 3. For instance, in Fig. 7, it is shown the ray obtained from the physical ray tracing algorithm (continuous line) in an axial GRIN described by

The physical ray tracing results are compared with the corresponding analytical solutions (dashed line in Fig. 7), which are given by [18]

where $\beta =n_{0}\cos {\theta _{0}}$, and $\theta _{0}$ is the incident angle of the ray at the point $x=0$ and $y=0$.From this comparison, it can be observed that our physical ray tracing ray (continuous line) is, for all practical purposes, spliced with the analytical solution (dashed line), in the same way as in the spherical GRIN lens (see Fig. 5).

## 5. Physical ray tracing in a cylindrical fiber with parabolic GRIN

The cylindrical case is specially interesting. In contrast to the spherically symmetric case where there is only one invariant, in the cylindrical geometry there are two invariants associated to rotational and translational symmetries respectively. The radial cylindrical GRIN possesses an axial symmetric index profile. Its corresponding isoindicial surfaces are concentric cylinders due to the axial symmetry. [1]. It can be represented by a function $n=n\left (r\right )$. Then, its optical Lagrangian in cylindrical coordinates $\left (r,\phi ,z\right )$, is given by [48]

where $\phi _{r}=\textrm {d}\phi /\textrm {d}r$, and $z_{r}=\textrm {d}z/\textrm {d}r$. The invariants corresponding to Eq. (5) for this medium areNow, consider a parabolic refractive index that its maximum value is located at the optical axis, from which it smoothly decreases along the transverse direction towards the periphery,

#### 5.1 Meridional rays

The meridional rays are those that their paths are confined into a plane, i.e., $\phi =\textrm {const}$ and $\tilde {l}=0$. This leads us to a single invariant, namely, $\tilde {\beta }$. Note that Eq. (29) can be rewritten as

where $\varphi _{C}$ and $\theta$ are complementary angles (see Fig. 8). By comparing Fig. 6 and Fig. 8, it can be observed that and consequently, where the invariant $\tilde {\beta }$ is the generalized Snell’s law for any plane containinig the propagation axis.The physical ray tracing is presented in Fig. 9, where the incident ray passes through the point $P_ {0}=\left (x_ {0}, y_{0}\right )$, and enters the medium at the point $P_{i}=\left (0,0\right )$. Our solution is compared with the corresponding analytical solution, which is given by [18]

#### 5.2 Helical rays

Throughout these sections, we have shown that our method is able to trace plane rays. However, our method is not limited to this kind of curves, it can be employed to trace rays that are not confined in a plane, that is, skew rays. Since the invariants $\tilde \beta$ and $\tilde {l}$ are manifestations of the translational and rotational invariance of the GRIN, they can be used to trace skew rays. The parameter $\tilde {l}$ is usually referred to as the skewness parameter [18].

Let us consider an incident ray that passes through the point $P_{0}=\left (a',y_{0},z_{0}\right )$, and enters the medium at the point $P_{i}=\left (a',0,0\right )$. Thus, the incident ray initially propagates on the $y$-$z$ plane, making an angle $\theta$ with the $z$-axis. The invariants $\tilde \beta$ and $\tilde {l}$ are given by

and where $\theta '=\frac {\pi }{2}-\sin ^{-1}\left \{\frac {n\left (r>a\right )}{n\left (a'\right )}\sin ^{-1}\left [\frac {\pi }{2}-\tan ^{-1}\left (\frac {y_{0}}{z_{0}}\right )\right ]\right \}$. Solving for $\theta$ and $\vartheta$ (see Fig. 8) it is obtained andUnlike the case of planar rays, in which Eq. (14) was useful, in the case of skew rays the angles $\theta$ and $\vartheta$ must be calculated through

By using Eqs. (36), (37), (38), (39) and (40), it is possible to perform the ray tracing using the algorithm described in Fig. 3, as can be observed in Fig. 10. From this Figure it is possible to observe that the path ray is a circular helix of radius $a'$, which is precisely the prediction given by the analytical solution [18], namely

where $\Gamma =\frac {n_{1}\sqrt {2\Delta }}{a\beta }$. Also, in Fig. 10 it is possible to observe that in the $x$-$y$ plane (Fig. 10(a)), in the $x$-$z$ plane (Fig. 10(b)), and in the three-dimensional space (Fig. 10(c)) the analytical solution (dashed line) for all practical purposes splices with the solution given by our physical ray tracing (continuous line).## 6. Conclusions

In this work we proposed a general ray tracing method for GRIN media with symmetrical properties, by means of Fermat’s ray invariants. The method was presented in detail and the corresponding algorithm benchmarked against several important GRIN media used for many applications. The algorithm is remarkably simple and efficient, and we showed that it works for any GRIN distribution possessing a symmetry in a given curvilinear orthogonal coordinate system. One of the main advantages of our algorithm is that it does not rely on approximations of any kind, that is, it is an exact numerical method. The algorithm gains importance since it significantly simplifies the representation of ray tracing in spherical symmetric systems, which are the least known among the ones presented in this work due to their complexity. A remarkable advantage over other methods is that it can be applied to cases where the ray bends backwards to the direction of the source or spirals around the center of the sphere. The setup of the method requires no explicit solution to the ray equation, only the GRIN distribution and the initial conditions of the incident ray as well. We provide evidence that for the cases discussed in this work, the physical ray tracing yields results with outstanding agreement with the analytical solution, as expected for an exact numerical method. Furthermore, the method is applicable to media where the refraction index has discontinuities of any order, including step discontinuities, which makes it reasonable to extend the applicability of the principles discussed here to the realm of acoustics for modeling acoustic devices such as those reported in [51–54]. Also, the physical principles and calculations allow for further generalization of Snell’s law to other variational principles, such as the Hamilton-Lagrange principle to describe the path of a particle subject to a potential, or the optical current streamlines of structured beams [55,56], or even to applications such as ray tracing in GRIN media describing from gravitational-like behavior [57], to human crystalline lens [58–60].

## Funding

Instituto Tecnológico y de Estudios Superiores de Monterrey.

## Acknowledgments

We respectfully dedicate this work to the memory of José Adrián Carbajal Domínguez, who was a researcher at Universidad Juárez Autónoma de Tabasco, Mexico. The authors acknowledge the support of the Tecnologico de Monterrey Vicepresidency of Research in this publication.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

No data were generated or analyzed in the presented research.

## References

**1. **C. Gómez-Reino, M. V. Pérez, and C. Bao, * Gradient-Index Optics: Fundamentals and Applications* (Springer, 2002).

**2. **E. W. Marchand, * Gradient Index Optics* (Academic Press, 1978).

**3. **M. Born and E. Wolf, * Principles of optics: Electromagnetic theory of propagation, interference and diffraction of light* (Cambridge University Press, 1999), 7th ed.

**4. **E. W. Marchand, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. **60**(1), 1–7 (1970). [CrossRef]

**5. **E. W. Marchand, “Rapid ray tracing in radial gradients,” Appl. Opt. **27**(3), 465–467 (1988). [CrossRef]

**6. **E. W. Marchand, “Ray tracing in cylindrical gradient-index media,” Appl. Opt. **11**(5), 1104–1106 (1972). [CrossRef]

**7. **D. T. Moore, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. **65**(4), 451–455 (1975). [CrossRef]

**8. **W. Streifer and K. B. Paxton, “Analytic solution of ray equations in cylindrically inhomogeneous guiding media. 1: Meridional rays,” Appl. Opt. **10**(4), 769–775 (1971). [CrossRef]

**9. **K. B. Paxton and W. Streifer, “Analytic solution of ray equations in cylindrically inhomogeneous guiding media. part 2: Skew rays,” Appl. Opt. **10**(5), 1164–1171 (1971). [CrossRef]

**10. **H. A. Buchdahl, “Rays in gradient-index media: separable systems,” J. Opt. Soc. Am. **63**(1), 46–49 (1973). [CrossRef]

**11. **A. Sharma, D. V. Kumar, and A. K. Ghatak, “Tracing rays through graded-index media: a new method,” Appl. Opt. **21**(6), 984–987 (1982). [CrossRef]

**12. **A. Sharma, “Computing optical path length in gradient-index media: a fast and accurate method,” Appl. Opt. **24**(24), 4367–4370 (1985). [CrossRef]

**13. **T. Sakamoto, “Ray trace algorithms for GRIN media,” Appl. Opt. **26**(15), 2943–2946 (1987). [CrossRef]

**14. **B. Richerzhagen, “Finite element ray tracing: a new method for ray tracing in gradient-index media,” Appl. Opt. **35**(31), 6186 (1996). [CrossRef]

**15. **D. A. Atchison and G. Smith, * Optics of the human eye* (Butterworth-Heinemann, 2002).

**16. **J. Aguilar-Gutiérrez, M. Arroyo Carrasco, and M. Iturbe-Castillo, “Slices method to describe ray propagation in inhomogeneous media,” Opt. Commun. **383**, 208–214 (2017). [CrossRef]

**17. **V. V. Krylov, “Transmission of a Rayleigh wave across smooth large-scale surface irregularities,” Soviet Physics Acoustics **34**(6), 613–618 (1988).

**18. **V. Lakshminarayanan, A. K. Ghatak, and K. Thyagarajan, * Lagrangian optics* (Springer, 2002).

**19. **A. J. Dragt, E. Forest, and K. B. Wolf, “Foundations of a lie algebraic theory of geometrical optics,” in * Lie Methods in Optics*, J. Sánchez Mondragón and K. B. Wolf, eds. (Springer Berlin Heidelberg, Berlin, Heidelberg, 1986), pp. 105–157.

**20. **S. V. Biryukov, Y. V. Gulyaev, V. V. Krylov, and V. P. Plessky, * Surface acoustic waves in inhomogeneous media* (Springer-Verlag Berlin Heidelberg, 1995).

**21. **J. W. Blaker and M. A. Tavel, “The Application of Noether’s Theorem to Optical Systems,” Am. J. Phys. **42**(10), 857–861 (1974). [CrossRef]

**22. **J. Evans and M. Rosenquist, ““F=ma” optics,” Am. J. Phys. **54**(10), 876–883 (1986). [CrossRef]

**23. **S. A. Khan, “Quantum methods in light-beam optics,” Opt. Photonics News **27**(12), 47 (2016). [CrossRef]

**24. **D. Ambrosini, A. Ponticiello, G. S. Spagnolo, R. Borghi, and F. Gori, “Bouncing light beams and the Hamiltonian analogy,” Eur. J. Phys. **18**(4), 284–289 (1997). [CrossRef]

**25. **S. A. Khan, “Hamilton’s optical–mechanical analogy in the wavelength-dependent regime,” Optik **130**, 714–722 (2017). [CrossRef]

**26. **H. Ohno, “Symplectic ray tracing based on Hamiltonian optics in gradient-index media,” J. Opt. Soc. Am. A **37**(3), 411–416 (2020). [CrossRef]

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

**28. **H. Ohno and T. Usui, “Points-connecting neural network ray tracing,” Opt. Lett. **46**(17), 4116–4119 (2021). [CrossRef]

**29. **R. K. Luneburg, * Mathematical Theory of Optics* (University of California Press, 1964).

**30. **G. B. Arfken and H. J. Weber, * Mathematical methods for physics*; 6th Ed. (Academic Press, 2005).

**31. **D. D. Holm, * Geometric Mechanics* (IMPERIAL COLLEGE PRESS, 2011), 2nd ed.

**32. **E. W. Weisstein, “Scale Factor, From MathWorld–A Wolfram Web Resource,” https://mathworld.wolfram.com/ScaleFactor.html.

**33. **D. Neuenschwander, * Emmy Noether’s Wonderful Theorem* (Johns Hopkins University Press, 2017).

**34. **J. R. Flores, “Estudio de elementos ópticos de gradiente de índice de simetría esférica,” Ph.D. thesis, Universidade de Santiago (1992).

**35. **A. S. Gutman, “Modified Luneberg lens,” J. Appl. Phys. **25**(7), 855–859 (1954). [CrossRef]

**36. **J. E. Eaton, “An extension of the Luneburg-type lenses,” Naval Research laboratory, Report No. 4110, 9–11 (1953).

**37. **E. E. Narimanov and A. V. Kildishev, “Optical black hole: Broadband omnidirectional light absorber,” Appl. Phys. Lett. **95**(4), 041106 (2009). [CrossRef]

**38. **Y. N. Demkov and V. N. Ostrovskiĭ, “Internal symmetry of the Maxwell “Fish-eye” problem and the fock group for the hydrogen atom,” Sov. Phys. JETP **6**, 1083–1087 (1971).

**39. **Q. Lei, R. Foster, P. S. Grant, and C. Grovenor, “Generalized Maxwell fish-eye lens as a beam splitter: A case study in realizing all-dielectric devices from transformation electromagnetics,” IEEE Trans. Microwave Theory Tech. **65**(12), 4823–4835 (2017). [CrossRef]

**40. **J. C. Maxwell, “Solutions of problems (prob. 3, vol. viii, p. 188),” Cambridge and Dublin Mathematical Journal **9**, 9–11 (1854).

**41. **J. E. Gómez-Correa, V. Coello, A. Garza-Rivera, N. P. Puente, and S. Chávez-Cerda, “Three-dimensional ray tracing in spherical and elliptical generalized Luneburg lenses for application in the human eye lens,” Appl. Opt. **55**(8), 2002–2010 (2016). [CrossRef]

**42. **J. M. Gordon, “Spherical gradient-index lenses as perfect imaging and maximum power transfer devices,” Appl. Opt. **39**(22), 3825–3832 (2000). [CrossRef]

**43. **S. P. Morgan, “General solution of the Luneberg lens problem,” J. Appl. Phys. **29**(9), 1358–1368 (1958). [CrossRef]

**44. **J. A. Lock, “Scattering of an electromagnetic plane wave by a Luneburg lens. i. ray theory,” J. Opt. Soc. Am. A **25**(12), 2971–2979 (2008). [CrossRef]

**45. **Y. G. Ma, C. K. Ong, T. Tyc, and U. Leonhardt, “An omnidirectional retroreflector based on the transmutation of dielectric singularities,” Nat. Mater. **8**(8), 639–642 (2009). [CrossRef]

**46. **J. A. Grzesik, “Focusing properties of a three-parameter class of oblate, Luneburg-like inhomogeneous lenses,” J. Electromagn. Waves Appl. **19**(8), 1005–1019 (2005). [CrossRef]

**47. **G. Beliakov, “Numerical evaluation of the Luneburg integral and ray tracing,” Appl. Opt. **35**(7), 1011–1014 (1996). [CrossRef]

**48. **A. Ghatak, E. Sharma, and J. Kompella, “Exact ray paths in bent waveguides,” Appl. Opt. **27**(15), 3180–3184 (1988). [CrossRef]

**49. **A. Ankiewicz and C. Pask, “Geometric optics approach to light acceptance and propagation in graded index fibres,” Opt. Quantum Electron. **9**(2), 87–109 (1977). [CrossRef]

**50. **A. Ankiewicz, “Geometric optics theory of graded index optical fibres,” Ph.D. thesis, The Australian National University (1978).

**51. **S.-H. Kim, “Sound focusing by acoustic Luneburg lens,” arXiv preprint arXiv:1409.5489 (2014).

**52. **C. M. Park, C. H. Kim, H. T. Park, and S. H. Lee, “Acoustic gradient-index lens using orifice-type metamaterial unit cells,” Appl. Phys. Lett. **108**(12), 124101 (2016). [CrossRef]

**53. **Y. Xie, Y. Fu, Z. Jia, J. Li, C. Shen, Y. Xu, H. Chen, and S. A. Cummer, “Acoustic imaging with metamaterial Luneburg lenses,” Sci. Rep. **8**(1), 16188 (2018). [CrossRef]

**54. **R. Fuentes-Domínguez, M. Yao, A. Colombi, P. Dryburgh, D. Pieris, A. Jackson-Crisp, D. Colquitt, A. Clare, R. J. Smith, and M. Clark, “Design of a resonant Luneburg lens for surface acoustic waves,” Ultrasonics **111**, 106306 (2021). [CrossRef]

**55. **M. V. Berry, “Optical currents,” J. Opt. A: Pure Appl. Opt. **11**(9), 094001 (2009). [CrossRef]

**56. **A. Jaimes-Nájera, J. E. Gómez-Correa, M. D. Iturbe-Castillo, J. Pu, and S. Chávez-Cerda, “Kepler’s law for optical beams,” Opt. Express **28**(21), 31979–31992 (2020). [CrossRef]

**57. **W. Xiao, S. Tao, and H. Chen, “Mimicking the gravitational effect with gradient index lenses in geometrical optics,” Photonics Res. **9**(7), 1197–1203 (2021). [CrossRef]

**58. **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**(17), 3990–3993 (2015). [CrossRef]

**59. **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**(7), 3699–3716 (2020). [CrossRef]

**60. **A. Jaimes-Nájera, J. E. Gómez-Correa, V. Coello, B. K. Pierscionek, and S. Chávez-Cerda, “A single-function model for the eye’s crystalline lens,” Opt. Photonics News **31**(12), 54 (2020). [CrossRef]