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

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

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

 

Fig. 1. The path of rays in a medium with gradient index.

Download Full Size | PDF

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)}$.

 

Fig. 2. Optical Path Length in a medium with constant refractive index.

Download Full Size | PDF

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.

 

Fig. 3. Algorithm of the physical ray tracing.

Download Full Size | PDF

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

Table 1. Refractive Gradient Index of spherical lenses

View Table

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]:

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

 

Fig. 4. Ray trajectories throught of Gutman lenses with (a) $f=0.75$, (b) $f=1$ (Luneburg lens), (c) $f=1.25$, (d) Eaton lens, (e) Concentrator Lens, Generalized Maxwell’s Fish-eye lens with (f) $m=1$ (Maxwell’s Fish-Eye lens), (g) $m=0.95$, (h) $m=0.85$, and (i) $m=0.75$.

Download Full Size | PDF

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

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}$.

 

Fig. 5. (a) Analytical and Physical ray tracing in Luneburg Lens. (b) Relative error of the physical ray tracing for several step values ($\Delta \mathbf {s}$). (c) Relative error of the four order Runge-Kutta algorithm for several step values. The horizontal axis is the propagation coordinate and the vertical axis indicates the numerical error. This error is calculated on the outermost ray of Fig. 5(a) because it is the ray that must have the highest error.

Download Full Size | PDF

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

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

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

 

Fig. 6. The ray tangent makes an angle $\varphi _{L}$ with the vertical.

Download Full Size | PDF

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,

As 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

 

Fig. 7. Physical ray tracing and analytical solution (with $n_{0}=1.5$, and $\alpha =0.1$) in the axial refractive index profile. The error magnitude is of $10^{-3}$ with $\Delta \mathbf {s}=0.0048$.

Download Full Size | PDF

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

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]

 

Fig. 8. The cylindrical parameters [49,50].

Download Full Size | PDF

Now, 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

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]

 

Fig. 9. Physical ray tracing and analytical solution (with $n_{1}=1.38$ and $\Delta =0.2$) in parabolic index fiber. The error magnitude is of $10^{-4}$ with $\Delta \mathbf {s}=0.027$.

Download Full Size | PDF

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

Unlike 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

 

Fig. 10. Physical ray tracing and analytical solution (with $n_{1}=1.38$ and $\Delta =0.2$) in parabolic index fiber with $P_{0}=\left (a'=4,0,0\right )$. Helical ray propagation observed from the $x$-$y$ ((a)), and $x$-$z$ ((b)) planes. (c) Three-dimensional helical ray propagation.

Download Full Size | PDF

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