## Abstract

Light propagation obeys Fermat’s principle, and an important inference of Fermat’s principle is the optical Lagrange equation, from which the light trace can be determined with a given refractive index. Here, we consider the inverse problem of how to derive the refractive index distribution of a planar geometric optical system once the trace is predetermined. Based on the optical Lagrange equation, we propose a dynamic equation model which associates the refractive index with the light trace. With the consideration of a certain trace, we illustrate the process of solving the partial differential equation of refractive index through first integral method. By setting the distribution function of a gradient-refractive-index (GRIN) medium, one can control the light traveling along a desirable curve, adjust the incoming and outgoing rays, and also use the trace to paint geometrics. This method develops the Lagrangian optics in the application of ray dynamic system design, such as lens, beam splitter, metasurface and optical waveguide. It provides a theoretical guidance to manipulate the ray in a GRIN medium.

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

## 1. Introduction

The desire to control the propagation of light by using gradient-refractive-index (GRIN) medium has received considerable attention during recent years. Because the GRIN elements can remarkably improve the properties of optical elements while reducing the size, weight, and the number of optical components, they can greatly promote the integrated level of optical systems and are widely used in the manufacture of lens [1–6], metasurface [7–10], biological optics [11–14], and visible systems [15–18]. In addition, because of the similar principles between light and acoustic wave, its extraordinary performance is also introduced in the design of phononic crystals [19,20].

Due to its wide application, it is fundamentally important to study how to properly design a GRIN system, so that the light trace can exhibit some expected characters. Actually, this is a geometric optics problem. As known, in geometric optics, light trace can be determined by the Fermat’s principle. Its counterpart in mechanics is the famous Hamilton’s principle. Because of the similarity between these two principles, there is a very close analogy between analytical mechanics and geometric optics [21]. There are two direct inferences of Fermat’s principle, known as the Snell’s law and the optical Lagrange (or Hamilton) equations. To anticipate the light trace, the former is applicable for the refraction in a piecewise constant refractive index, while the latter is convenient to handle the continuous cases, that is, GRIN medium [22–26]. Aiming to the design of on-chip lens and bio-photonics applications, previous studies mostly focused on the study of optical system with a given refractive index function [2,6,7,11,14,27]. These systems are usually designed by experience. Since the lineshape of refractive index is already known, the only procedure is to modify the undetermined coefficients in a given function to match the technical requirement [2,14,27]. But usually, what we expect is a more purposeful design for optical systems. That is to say, we expect the ray can travel as imagined (its predetermined path may be an arbitrary curve), rather than to solve the light trace with a given refractive index function. Therefore, we propose such a problem: to design a certain distribution function (it is beforehand unknown and thus cannot be obtained by fitting parameters), so that the ray can travel along a predetermined path. However, the general method for this inverse process, to the best of our knowledge, has not been paid much attention.

In this paper, we will seek for a general method on how to design the GRIN optical system when the light trace is predetermined, so that we can precisely manipulate the light trace and let it travel along a desirable geometric curve. This may provide some reference models for the design of optical systems. The situation we consider here is a pure geometric optical system, that is, we ignore the dispersion, dissipation, interference, diffraction and other factors, only need to consider the beam’s refraction in a GRIN medium. At some point, this model is similar to the corresponding classical mechanics system [21]. Based on variation principle and optical Lagrange equation, we will illustrate the general procedure of solving the refractive index function with a given trace through several examples.

This paper is organized as follows: In Sec. 2, we briefly review the Fermat principle and optical Lagrange equation and derive its Cartesian-coordinate form. This allows us to perform the ray dynamic system design through a predetermined light trace. In Sec. 3, we choose the light trace as a given function, and seek for the solution of the background refractive index through first integral method. We will show the general procedures to design the background medium according to the light trace, and illustrate their applications. Sec. 4 summarizes the paper.

## 2. Basic equations

We begin our discussion by reviewing the basic Fermat principle: the optical path reaches the extreme value when light propagates through a background medium, which is mathematically presented as the variation of optical path functional equals to zero, that is

For simplicity, we consider the two-dimensional case, that is, we assume the medium is*z*-direction homogeneous and the ray is launched in

*xOy*plane. Using Cartesian coordinates, Eq. (1) can be rewritten as

We substitute $\frac{\partial L}{\partial {y}^{\prime}}=n\frac{{y}^{\prime}}{\sqrt{1+{{y}^{\prime}}^{2}}},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial L}{\partial y}=\frac{\partial n}{\partial y}\sqrt{1+{{y}^{\prime}}^{2}}$ into Eq. (3) and further simplify it to

It should be emphasized that, for a predetermined trace $y=f(x)$, although the refractive index function $n=n(x,y)$ can be solved from Eq. (4), it does not ensure that the light must travel along the trace $y=f(x)$, because when we substitute $n=n(x,y)$ back into Eq. (4), $y=f(x)$ is not surely the unique solution. Physically, a specific $n(x,y)$ may support more than one propagation mode, and it depends on the position and the direction of light source (initial conditions). In addition, once we substitute solution $n=n(x,y)$ back to Eq. (4), the initial conditions of Cauchy problem $y({x}_{0})={y}_{0},{y}^{\prime}({x}_{0})={{y}^{\prime}}_{0}$ may conflict with the predetermined trace $y=f(x)$. In such a case, the variation problem has no solution, and the predetermined light trace just means the asymptotic line. The light will travel, to some extent, close to this curve, but not coincide. Therefore, in order to obtain predetermined trace $y=f(x)$, one need not only to find a certain $n(x,y)$ satisfies Eq. (4) (this may not be unique), but also to properly set the light source so that its launch condition (initial condition) satisfies $y=f(x)$.

Generally, to derive the PDE about $n(x,y)$, the light equation should be considered. This indicates that the solution of the PDE about $n(x,y)$ may be added an additional condition $y=f(x)$. In other words, $n(x,y)$ only needs to satisfy Eq. (4) on the predetermined light trace$y=f(x)$, and no need to be valid on the whole plane area. This condition is physically weak, but for Eq. (4), it is a constraint condition and thus mathematically strong. Nevertheless, for $n(x,y)$ the constraint condition $y=f(x)$is not necessary, since a given $n(x,y)$ which satisfies the equation in the whole area is certainly valid on a specific curve. That is to say, we can also physically strengthen the condition, so that $n(x,y)$ in Eq. (4) is valid in the whole plane area, and the conditions that $n(x,y)$ should satisfy are mathematically weakened. Making use of the trace equation may be, sometimes, beneficial to solve the problem.

Based on the first-order linear PDE theory, we can “play with the beam” and depict some common geometrics by properly setting the refractive index. Below are some simple examples.

## 3. Examples

Here we show some examples to solve Eq. (4) under a predetermined light trace, and simply introduce their applications. In order to conveniently illustrate this method, the chosen examples are all analytically solvable. For some more complicated light traces, we can also refer to numerical method to solve Eq. (4) and get the refractive index function.

#### 3.1 Straight line with certain slope (parallel lines)

As is familiar to us, light travels along a straight line in a homogeneous medium. This can be also seen in Eq. (4). As $n=n(x,y)\equiv const$ yields ${y}^{\u2033}=0$. However, can we add some restrictions so that only light propagates along a certain direction can go straight? This can be answered from Eq. (4). Assume the trace equation is $y=kx+b$ and substitute it into Eq. (4)

The character equation of Eq. (5) is $\frac{dx}{k}=\frac{dy}{-1}$, and we can derive its first integral$x+ky=C$_{.}The general solution of Eq. (5) isIn such a case, if the initial condition of $y=f(x)$ satisfies ${y}^{\prime}({x}_{0})=k$, the light will travel in straight line. In other words, if the light source’s inclination with respect to horizontal axis satisfies $\mathrm{tan}\alpha =k$, the light trace is straight. However, if the source does not satisfy this condition, the initial conditions will contradict with the light equation (whatever

*b*is), and the variation problem has no solution, which means $y=kx+b$ is not the light trace but an asymptotic line of the light trace. That is to say, if the refractive index is a non-constant function of

*x*+

*ky*, only light with slope

*k*can travel in straight line, as the simulation shown in Fig. 1(a). These rays will be approximately parallel after a certain propagation distance. This is in contrast with homogeneous case.

By making use of this effect, we can transform the divergent rays into the directional ones, as shown in the schematics (Fig. 1(b)). This way is more convenient than traditional convex lens and concave mirror, since there is no need to place the source at the focus, only need that the GRIN medium is thick enough. Also, compared with the Wood lens [28], the design of the medium’s parameter has more choices, because we only need to design the refractive index as the form of Eq. (6), and specially, it can be a linear function of coordinates. Noticeably, when we take advantage of this effect, the medium should become denser along the propagation direction of light, so that it can collimate the beam. However, if the medium becomes rarer, the beam will be divergent, since the light path is invertible. This inverse process can be also used to focus the parallel beam. Compared with the traditional lens, if the thickness of GRIN medium is adjustable, it can realize a changeable focus.

Actually, this result is explainable from the perspective of refraction. From Eq. (6), the direction of refractive index gradient $\nabla n=\frac{\partial n}{\partial x}i+\frac{\partial n}{\partial y}j={\varphi}^{\prime}(x+ky)(i+kj)$ is in accordance with the propagation direction. This can be understood as the light always points toward the “normal incidence” direction from one medium to another. Refraction does not happen along the normal incidence direction and thus its direction does not change. Therefore, if we set the initial direction of source as $i+kj$, the light trace is apparently straight.

#### 3.2 Circle

The above example is, in a sense, to seek for the general condition under which the light can travel in straight line. From this one we begin to discuss how to control the light propagating along a certain geometric shape. Assume the light propagates as a trace of circle ${x}^{2}+{y}^{2}={R}^{2}$, and therefore ${y}^{\prime}=-\frac{x}{y}$ and ${y}^{\u2033}=-\frac{{R}^{2}}{{y}^{3}}$. Equation (4) becomes

Using polar coordinates, $n=n(r,\theta )$, $\frac{\partial n}{\partial x}=\mathrm{cos}\theta \frac{\partial n}{\partial r}-\frac{\mathrm{sin}\theta}{r}\frac{\partial n}{\partial \theta}$, $\frac{\partial n}{\partial y}=\mathrm{sin}\theta \frac{\partial n}{\partial r}+\frac{\mathrm{cos}\theta}{r}\frac{\partial n}{\partial \theta}$, and Eq. (7) becomesThe solution of Eq. (8) iswhere ${r}_{0}(\theta )$ is an undetermined function. If the distribution is axisymmetric, we can set ${r}_{0}(\theta )\equiv {r}_{0}$, and the solution becomesConsidering the fact that the refractive index is usually greater than 1 and in case it goes to infinite, we modify Eq. (10) asWe choose the Radii of circles as R = a, 2a, 3a. To ensure the circles are located in the GRIN medium, ${r}_{0}<a<\frac{5}{3}{r}_{0}$ is required. We show the result in Fig. 2. The initial conditions, according to the trace equation, should be ${y}^{\prime}|{}_{x={x}_{0}}=-\frac{{x}_{0}}{{y}_{0}}$ and ${x}_{0}^{2}+{y}_{0}^{2}={R}^{2}$. This means, to obtain a circle trace, the launch direction of the source should be perpendicular to the radial direction, and the radius of circle is determined by the launch point.

Undoubtedly, this result can be also explained from the perspective of refraction. The gradient of refractive index $\nabla n=-\frac{1}{{r}^{2}}{e}_{r}$ is along the radial direction, and the light propagating along the tangential direction ${e}_{\theta}$(perpendicular to the radial direction) can be understood as “grazing incidence”. Because the angular distribution of refractive index is uniform, the deflection of light should be equal everywhere when $r$ is fixed. This means, if the light initially travels towards the tangential direction ${e}_{\theta}$, it will point to this direction forever, indicating that the trace is a circle. This is similar to a particle’s circular motion under a centripetal force (or a gravitational potential), since the gradient field always points to the center. As the gradient gets smaller with the increase of distance, the refraction is weakened and the trace’s curvature is also getting smaller, resulting in a larger circle trace, as expected.

#### 3.3 Equilateral hyperbola

Assume the light travels as a trace of equilateral hyperbola ${x}^{2}-{y}^{2}=\pm {c}^{2}$, and we have ${y}^{\prime}=\frac{x}{y}$ and ${y}^{\u2033}=\mp \frac{{c}^{2}}{{y}^{3}}$, and Eq. (4) becomes

We can simply illustrate the application of this example. Actually, this optical system can be used as a beam splitter. Imagine a wide beam is incident along one asymptotic line from infinity, and it will be split into two beams. With the increase of the propagation distance, they will be totally separated, as shown in Fig. 4. Compared with Fig. 3, the optical system here is rotated by 45 degrees. If the GRIN medium is thick and wide enough, they are close to the other asymptotic line according to the above result. We can externally connect the homogenous medium to the GRIN medium, so that they will finally go to opposite directions.

Examples 3.2 and 3.3 are very similar to a particle moves in a central gravitational force and a central repulsive force. Principally, similar procedures can be performed for an ellipse and a general hyperbola if the first integrals of these traces are obtained.

#### 3.4 Parabola

We set the predetermined light trace as a parabola $y=a{x}^{2}(a>0)$, and we have ${y}^{\prime}=2ax,{y}^{\u2033}=2a$_{.} Equation (4) becomes

*y*, the refractive index increases while the gradient $\nabla n$ decreases. This indicates that, when the light direction is inclined to the + y direction, it is equivalent that the light continuously travels from the denser medium to the thinner one, and the slope of the trace will increase. Because $\nabla n\to 0$ when $y\to +\infty $, the medium is approximately homogeneous and the light approximately rises in a straight line when y is large. Definitely, if the light direction is inclined to the

*-y*direction, we can derive the opposite conclusion. This is in accordance with the characters of parabolic trace.

#### 3.5 Exponential (Logarithmic) curve

Through similar procedures, we can also manipulate the ray to travel as an exponential curve. Substituting $y={y}^{\prime}={y}^{\u2033}={e}^{x}$ into Eq. (4) ${y}^{\prime}\frac{\partial n}{\partial x}-\frac{\partial n}{\partial y}+\frac{{y}^{\u2033}}{1+{{y}^{\prime}}^{2}}n=0$ yields

The characteristic equation of Eq. (19) is $\frac{dx}{y}=-dy=-\frac{1+{y}^{2}}{ny}dn$, giving the first integral $2x+{y}^{2}={C}_{1},\frac{{n}^{2}}{1+{y}^{2}}={C}_{2}$. Therefore the implicit general solution is $F\left(2x+{y}^{2},\frac{{n}^{2}}{1+{y}^{2}}\right)=0$, which givesAs above, we set $\varphi \equiv 1$ and obtain the refractive index functionWhen the refractive index is set as Eq. (21), it is easy to see $y={e}^{-x}$,$y=-{e}^{x}$ and $y=-{e}^{-x}$ also satisfy Eq. (4) because of the symmetry. Therefore, depending on the initial conditions, we can obtain four exponential curves. If we exchange*x*and

*y*coordinates, the corresponding logarithmic curves can be derived. We depict the result in Fig. 6.

Similar to the Example 3.3, this system can be also used as a beam splitter. If the light is incident along the *x* axis from infinity, the original beam will split into two branches $y={e}^{x}$ and $y={e}^{-x}$. Compared with hyperbolic trace used in Fig. 4, these two branches will be separated in a shorter distance, since the exponential function varies more rapidly than hyperbola. In contrast with the traditional beam splitter, the beam-splitting angle between two branches can be adjusted by properly setting the thickness of GRIN medium. If we truncate the GRIN medium at $x={x}_{0}$, we can compute that the beam-splitting angle is $\theta =2\mathrm{arctan}{e}^{{x}_{0}}$.

From the above examples, we see that, the gradient of refractive index is very similar to the force imposed on a moving particle. It always deflects the light trace towards its direction, as a particle’s movement always bends to the direction of force. These examples indicate that the geometric optics and the classical mechanics are intrinsically related. As they can both be described by Lagrangian dynamics, they share the same essence. Therefore, we can understand the gradient of refractive index as a generalized force which can act on the light.

## 4. Conclusions

In conclusion, we have shown the general method to manipulate the light trace in a GRIN system. Based on optical Lagrange equation, we propose a PDE model which associates the refractive index function with the light trace. By performing the first integral method, the refractive index function can be conveniently solved once a predetermined trace is expected. Through properly designing the refractive index function and setting the light source, we have shown some examples to manipulate the light travelling along specific geometric curves, such as straight line, circle, equilateral hyperbola, parabola, and exponential curve. The light trace in a GRIN medium is quite similar to a particle’s motion trace under a coordinate-dependent force. The solution of this PDE model provides us an efficient way for the planar geometric optical system design. Our work presents a versatile method to control the light propagation, incidence and emergence through designing the refractive index function of a GRIN medium. It provides a theoretical framework for the potential application of GRIN systems, such as beam splitter, lens antenna, metasurface, optical waveguide and optical fiber.

## Funding

National Natural Science Foundation of China (11864032), Youth Talent Scientific Project of Education Department of Guizhou Province (KY[2017]337, KY[2016]319), Talent Project of Qiannan Normal University for Nationalities (QNSYRC201720, QNSYRC201618, QNSY2018BS017), and the Foundation of Scientific Innovative Research Team of Education Department of Guizhou Province (201329).

## References

**1. **R. B. Greegor, C. Parazzoli, J. Nielsen, M. A. Thompson, M. H. Tanielian, and D. R. Smith, “Simulation and testing of a graded negative index of refraction lens,” Appl. Phys. Lett. **87**(9), 91114 (2005). [CrossRef]

**2. **X. Mao, S. C. Lin, M. I. Lapsley, J. Shi, B. K. Juluri, and T. J. Huang, “Tunable Liquid Gradient Refractive Index (L-GRIN) lens with two degrees of freedom,” Lab Chip **9**(14), 2050–2058 (2009). [CrossRef] [PubMed]

**3. **A. I. Hernandez-Serrano, M. Weidenbach, S. F. Busch, M. Koch, and E. Castro-camus, “Fabrication of gradient-refractive-index lenses for terahertz applications by three-dimensional printing,” J. Opt. Soc. Am. B **33**(5), 928 (2016). [CrossRef]

**4. **S. D. Campbell, J. Nagar, and D. E. Brocker, “On the use of surrogate models in the analytical decompositions of refractive index gradients obtained through quasiconformal transformation optics,” J. Opt. **18**(4), 044019 (2016). [CrossRef]

**5. **J. D. Ellis, D. R. Brooks, K. T. Wozniak, G. A. Gandara-Montano, E. G. Fox, K. J. Tinkham, S. C. Butler, L. A. Zheleznyak, M. R. Buckley, P. D. Funkenbusch, and W. H. Knox, “Manufacturing of Gradient Index Lenses for Ophthalmic Applications.” Design and Fabrication Congress 2017 (IODC, Freeform, OFT) **OW1B–3** (2017). [CrossRef]

**6. **H. F. Ma, B. G. Cai, T. X. Zhang, Y. Yang, W. X. Jiang, and T. J. Cui, “Three-dimensional gradient-index materials and their applications in microwave lens antennas,” IEEE Trans. Antenn. Propag. **61**(5), 2561–2569 (2013). [CrossRef]

**7. **E. Erfani, M. Niroo-Jazi, and S. Tatu, “A high-gain broadband gradient refractive index metasurface lens antenna,” IEEE Trans. Antenn. Propag. **64**(5), 1968–1973 (2016). [CrossRef]

**8. **N. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: generalized laws of reflection and refraction,” Science **334**(6054), 333–337 (2011). [CrossRef] [PubMed]

**9. **S. Sun, Q. He, S. Xiao, Q. Xu, X. Li, and L. Zhou, “Gradient-index meta-surfaces as a bridge linking propagating waves and surface waves,” Nat. Mater. **11**(5), 426–431 (2012). [CrossRef] [PubMed]

**10. **D. Lin, P. Fan, E. Hasman, and M. L. Brongersma, “Dielectric gradient metasurface optical elements,” Science **345**(6194), 298–302 (2014). [CrossRef] [PubMed]

**11. **R. H. H. Kröger and A. Gislén, “Compensation for longitudinal chromatic aberration in the eye of the firefly squid, Watasenia scintillans,” Vision Res. **44**(18), 2129–2134 (2004). [CrossRef] [PubMed]

**12. **D. E. Nilsson, L. Gislén, M. M. Coates, C. Skogh, and A. Garm, “Advanced optics in a jellyfish eye,” Nature **435**(7039), 201–205 (2005). [CrossRef] [PubMed]

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

**14. **S. Ji, M. Ponting, R. S. Lepkowicz, A. Rosenberg, R. Flynn, G. Beadie, and E. Baer, “A bio-inspired polymeric gradient refractive index (GRIN) human eye lens,” Opt. Express **20**(24), 26746–26754 (2012). [CrossRef] [PubMed]

**15. **O. M. Efimov, L. B. Glebov, L. N. Glebova, K. C. Richardson, and V. I. Smirnov, “High-efficiency bragg gratings in photothermorefractive glass,” Appl. Opt. **38**(4), 619–627 (1999). [CrossRef] [PubMed]

**16. **O. M. Efimov, L. B. Glebov, V. I. Smirnov, and L. Glebova, “Process for production of high efficiency volume diffractive elements in photo-thermo-refractive glass,” US 09/648,293, (2003).

**17. **L. G. Atkinson, D. S. Kindred, D. T. Moore, and J. R. Zinter, “Negative abbe number radial gradient index relay and use of same,” US 08/017,034 (1994).

**18. **T. H. Tomkinson, J. L. Bentley, M. K. Crawford, C. J. Harkrider, D. T. Moore, and J. L. Rouke, “Rigid endoscopic relay systems: a comparative study,” Appl. Opt. **35**(34), 6674–6683 (1996). [CrossRef] [PubMed]

**19. **S. S. Lin, T. J. Huang, J.-H. Sun, and T.-T. Wu, “Gradient-index phononic crystals,” Phys. Rev. B Condens. Matter Mater. Phys. **79**(9), 094302 (2009). [CrossRef]

**20. **S. Tol, F. L. Degertekin, and A. Erturk, “Gradient-index phononic crystal lens-based enhancement of elastic wave energy harvesting,” Appl. Phys. Lett. **109**(6), 063902 (2016). [CrossRef]

**21. **H. Goldstein, *Classical Mechanics* (Addison-Wesley, Cambridge, Mass. 1950)

**22. **A. Romano, *Geometric Optics: Theory and Design of Astronomical Optical Systems Using Mathematica* (Birkhäuser, 2010).

**23. **A. L. Rivera, S. M. Chumakov, and K. B. Wolf, “Hamiltonian foundation of geometrical anisotropic optics,” J. Opt. Soc. Am. A **12**(6), 1380–1389 (1995). [CrossRef]

**24. **J. C. Kimball and H. Story, “Fermat’s principle, Huygens’ principle, Hamilton’s optics and sailing strategy,” Eur. J. Phys. **19**(1), 15–24 (1998). [CrossRef]

**25. **D. S. Lemons, “Gaussian thin lens and mirror formulas from Fermat’s principle,” Am. J. Phys. **62**(4), 376–378 (1994). [CrossRef]

**26. **V. Lakshminarayanan, A. Ghatak, and K. Thyagarajan, *Lagrangian optics*. (Springer Science & Business Media, 2013).

**27. **J. N. Mait, G. Beadie, R. A. Flynn, and P. Milojkovic, “Dispersion design in gradient index elements using ternary blends,” Opt. Express **24**(25), 29295–29301 (2016). [CrossRef] [PubMed]

**28. **A. Hussain, Q. A. Naqvi, and K. Hongo, “Radiation characteristics of the Wood lens using Maslov’s method,” Prog. Electromagnetics Res. **73**, 107–129 (2007). [CrossRef]