## Abstract

A method based on Hamiltonian optics for ray tracing through gradient-index (GRIN) media is proposed. The ray equation that describes light-ray paths can be written in the form of the Hamiltonian equations. Although the Hamiltonian equations can be numerically calculated using a finite-difference explicit method, deviations from the exact equations are generally inevitable at subsequent time steps. An optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to such a separable optical Hamiltonian system and makes the optical Hamiltonian equations form invariant at each time step of numerical calculations. Accuracies of light-ray paths calculated using the first-order symplectic ray tracing in GRIN lenses approximate those calculated on the basis of the fourth-order Runge–Kutta algorithm, which shows the promising potential of the symplectic-ray-tracing method.

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

## 1. INTRODUCTION

Light-ray paths can be controlled by refractive index fields. Since novel design techniques of gradient-index (GRIN) media based on the coordinate transformation have rapidly been developed, new types of GRIN media have been created [1–6]. In designing such GRIN media, there is a growing need for a highly efficient and practical ray-tracing method. Physical quantities such as temperature, pressure, strain, and density in media can generate GRIN fields. Several methods to measure the distributions of such physical quantities in transparent media have been rapidly developed by means of measuring deflection angles of light rays passing through the generated GRIN fields [7,8]. These measuring methods also require a highly efficient and practical ray-tracing method.

The ray equation that is derived based on Fermat’s principle can describe light-ray paths in a GRIN medium. However, because the light-ray paths in the GRIN medium are arbitrarily curved in general, the ray equation should be solved numerically except for a few cases in which analytical solutions are possible. Several ray-tracing techniques based on finite-difference explicit methods have been developed for GRIN media [9–13]. In particular, numerical integration based on the Runge–Kutta algorithm uses a standard numerical technique, is practically treatable, and has proven to be a faster and more accurate numerical method than the other methods [14].

The ray equation can be transformed into the Hamiltonian equations. Several forms of the Hamiltonian equations have been extensively studied [15]. Although the Hamiltonian equations can be numerically calculated using finite-difference explicit methods, deviations from their exact equation are inevitable at subsequent time steps. This causes the accumulation of error. The symplectic integrator that is able to preserve the equation forms has therefore been developed for celestial mechanics, nonlinear Hamiltonian systems, and quantum physics [16–19]. The symplectic integrator has also revealed its worth in optical Hamiltonian systems to increase accuracies of numerical calculations of the finite-difference time-domain (FDTD) method and wavelet transformation [20–22]. A method based on the symplectic integrator algorithm for ray tracing in GRIN media is thus proposed here. The optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to such a separable optical Hamiltonian system and ensures that the optical Hamiltonian equations are almost form invariant at each time step of numerical calculations. This should increase accuracy of the ray tracing. The light-ray paths in the cylindrical and spherical GRIN lenses are calculated using the symplectic ray tracing method and are compared with those calculated on the basis of the fourth-order Runge–Kutta algorithm. It is shown that the deviations of light-ray paths from ideal paths are supremely reduced by the symplectic ray-tracing method.

## 2. SYMPLECTIC-RAY-TRACING METHOD BASED ON HAMILTONIAN OPTICS

The path of light rays in a GRIN medium can be described by the ray equation that can be derived from the Fermat’s principle as

A point ${\boldsymbol a}$ in the phase space can be written with ${\boldsymbol q}$ and ${\boldsymbol p}$ as

Another phase space point at a next time step, ${\boldsymbol A}$, can be written with an infinitesimal time $\Delta t$ as

Equation (11) is called the symplectic conditions that ensure the phase space volume conservation, and they can be written using Eqs. (5) and (8) as

where {} denotes the Poisson bracket, $i$ and $j$ denote components of each vector, and $\delta^{ij}$ indicates the Kronecker’s delta.The ${\boldsymbol P}$ and ${\boldsymbol Q}$ of the phase point ${\boldsymbol A}$ can be written with the ${\boldsymbol p}$ and ${\boldsymbol q}$ of the phase point ${\boldsymbol a}$ in a finite-difference explicit form using Eqs. (3) and (4) as

On the other hand, instead of the finite-difference explicit form of Eqs. (15) and (16), the following equations in the form of the first-order symplectic integrator (i.e., semi-implicit Euler form) are therefore adopted as

As one of the Hamiltonians that satisfy the condition of Eq. (21) for an optical system, the following optical Hamiltonian, ${H_{\rm opt}}$, that consists of two independent terms, i.e., one term dependent on the position, ${\boldsymbol q}$, and the other term dependent on the momentum, ${\boldsymbol p}$, is adopted as

Using the optical Hamiltonian, the ray equation of Eq. (1) can be transformed with the light-ray momentum represented by Eq. (2) as

In this paper, $t$ is called time, which is not actually physical time.

The derivative of the position ${\boldsymbol q}$ with respect to $t$ can be derived using the optical Hamiltonian ${H_{\rm opt}}$ as

Equations (24) and (26) are in the form of the Hamilton’s equations with canonical coordinates of (${\boldsymbol q}$, ${\boldsymbol p}$). The initial condition of the optical Hamiltonian can be derived from Eq. (2) as

Using the optical Hamiltonian of Eq. (22), the first-order symplectic integrator of Eqs. (18) and (19) can be written as

Note that the second term on the right-hand side of Eq. (29) is set to ${\boldsymbol P}({\boldsymbol q},\;{\boldsymbol p})$ instead of ${\boldsymbol p}$. Equations (28) and (29) represent a symplectic ray tracing on the basis of Hamiltonian optics, which ensures that the Hamiltonian equations corresponding to Eqs. (24) and (26) are almost form invariant. The canonical transformation that satisfies the symplectic conditions ensures that the phase space volume is conserved. In geometrical optics, the four-dimensional phase space volume conservation is called etendue conservation, which is widely used for lighting lens designs [24,25].

## 3. NUMERICAL CALCULATIONS OF LIGHT-RAY PATHS USING SYMPLECTIC RAY TRACING IN GRIN LENSES

#### A. GRIN Lens with Cylindrical Symmetry

GRIN lenses having cylindrical symmetry and flat surfaces are widely used for several applications, such as imaging optics, photocopiers, and scanners. A typical refractive index of a cylindrical GRIN lens [26] can be written in a cylindrical coordinate system ($\rho $, $\theta $, $z$) as

where $\rho $ is a radius, ${n_0}$ is a refractive index of center core of the GRIN lens, and $A$ is a positive constant. In the Cartesian coordinate system ($x$, $y$, $z$), the radius $\rho $ can be written as The analytical solution of the ray equation [i.e., Eq. (1)] for the refractive index of Eq. (30) leads to the following trajectory of a light ray initially parallel to the $z$ axis at the position of $x = {x_0}$ and $z = {0}$ in a cross-sectional $xz$-plane at $y = {0}$: where $\Lambda$ denotes a period and can be written asIn this work, the values of the parameters of $A$, ${n_0}$, and ${x_0}$ are set to ${0.5}\;{{\rm mm}^{ - 1}}$, 1.564, and 0.5 mm, respectively. These values of the parameters are taken from Ref. [26]. The analytical value of the period $\Lambda $ is 12.16733603 mm.

Figure 1 shows the light-ray paths numerically calculated using the several methods in the cross-sectional $xz$-plane. The horizontal axis indicates the $z$ axis normalized by the period $\Lambda $, and the vertical axis indicates the $x$ axis normalized by the initial $x$ coordinate of ${x_0}$. The light rays are incident on the cylindrical GRIN lens in a direction parallel to the $z$ axis. An ideal light-ray path arrives at $(x,z) = ({x_0}, \Lambda) $. Numerical calculations using the first-order explicit Euler method represented by Eqs. (15) and (16), the first-order symplectic method represented by Eqs. (28) and (29), and the fourth-order symplectic method [18,23] are, respectively, plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square. The analytical solution is also plotted by a dashed line. The time step $\Delta t$ for each method is set to 0.5 mm (i.e., $\Delta t/\Lambda = {0.04109363}$). Note that the time is not actual time and has the dimension of length in this work. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order explicit Euler method and approximates the analytical solution. The fourth-order symplectic method is found to be more accurate than the first-order symplectic method, as expected.

#### B. GRIN Lens with Spherical Symmetry

One of the widely known GRIN lenses having spherical symmetry is the Luneburg lens [27]. The refractive index of the Luneburg lens having a radius of ${R_0}$ can be written as

where $r$ denotes a radius from the center of the Luneburg lens. A light-ray path can be solved analytically in a cross-sectional $xz$-plane at $y = {0}$ with an initial $x$ position of ${x_0}$ and an initial direction parallel to the $z$ axis asFigure 2 shows light-ray paths numerically calculated using several methods in the cross-sectional $xz$-plane at $y = {0}$. The horizontal axis indicates the $z$ axis normalized by the lens radius ${R_0}$, and the vertical axis indicates the $x$ axis normalized by the radius ${R_0}$. The light rays are incident on the Luneburg lens in a direction parallel to the $z$ axis. An ideal focus point is $(x,z) = ({0},{R_0})$. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively.

The analytical solution is also plotted by dashed line. The time step $\Delta t$ for each method is set to 0.01 ${R_0}$. The $x$ positions of ${x_0}$ of incident light rays (i.e., $x$ positions of light rays at the periphery of the Luneburg lens) are set to 0.15 ${R_0}$ and 0.75 ${R_0}$, respectively. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order explicit Euler method and approximates the analytical solution. The first-order symplectic method also approximates the fourth-order symplectic method.

## 4. RESULTS

Errors of numerically calculated light-ray paths using several methods for the cylindrical and spherical GRIN lenses are calculated and are compared with each other.

First, errors of light-ray paths for the cylindrical GRIN lens as mentioned in Section 3.A are calculated. The error is defined as a ratio of an absolute value of deviation of $x$ from analytical solution ${x_{\rm ana}}$ to ${x_{\rm ana}}$ at $z = \Lambda $. Figure 3 shows errors with respect to the time step $\Delta t$ normalized by $\Lambda $. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method [14] is also plotted by dashed line. The time steps, $\Delta t$, for both the symplectic method and Runge–Kutta method are set to the same value where the definitions of the time in both methods are the same. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order Euler method. The first-order symplectic method also approximates the fourth-order Runge–Kutta method. The fourth-order symplectic method agrees well the fourth-order Runge–Kutta method, whereas there is a slight difference in the larger $\Delta t$ of more than 0.09 $\Lambda $.

Second, errors of light-ray paths for the spherical GRIN lens (i.e., the Luneburg lens) as mentioned in Section 3.B are calculated. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution to the lens radius ${R_0}$ at the focal point of $z = {R_0}$. Figure 4 shows errors with respect to the time step $\Delta t$ normalized by ${R_0}$. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order Euler method. The first-order symplectic method also approximates the fourth-order symplectic method. The first-order symplectic method, the fourth-order symplectic method, and the fourth-order Runge–Kutta method are the same order of magnitude. The error of the first-order symplectic method can be less than that of the fourth-order Runge–Kutta method.

The first-order symplectic method represented by Eqs. (28) and (29) is simple enough to reduce the computational cost and approximates higher-order methods. Thus, it should be a promising candidate for ray tracing in GRIN media.

## 5. CONCLUSIONS

A method for ray tracing through GRIN media with a simple formulation and high accuracy is proposed on the basis of Hamiltonian optics. The optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to the separable Hamiltonian system. It can be shown that the Hamiltonian equations are ensured to be almost form invariant at each time step of numerical calculations using the symplectic integrator, which can increase accuracy of the ray tracing. The light-ray paths in GRIN lenses are calculated using the symplectic-ray-tracing method of the first and fourth order, which approximate those calculated on the basis of the fourth-order Runge–Kutta algorithm. It is shown that the deviations of light-ray paths from the ideal paths are supremely reduced by the symplectic-ray-tracing method in comparison with the explicit Euler method. The first-order symplectic ray tracing represented by Eqs. (28) and (29) is simple enough to reduce the computational cost and approximates higher-order methods. Thus, symplectic ray tracing should be a promising candidate for ray tracing in GRIN media and should enlarge the range of the method choice.

## Disclosures

The author declares no conflicts of interest.

## REFERENCES

**1. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic
field,” Science **312**,
1780–1782 (2006). [CrossRef]

**2. **U. Leonhardt, “Optical conformal
mapping,” Science **312**, 1777–1780
(2006). [CrossRef]

**3. **Y. A. Urzhumov, N. B. Kundtz, D. R. Smith, and J. B. Pendry, “Cross-section comparisons of
cloaks designed by transformation optical and optical conformal
mapping approaches,” J. Opt. **13**, 024002 (2011). [CrossRef]

**4. **U. Leonhardt, “Notes on conformal
invisibility devices,” New J. Phys. **8**, 118 (2006). [CrossRef]

**5. **C. Lu and Z. L. Mei, “Multi-functional lens based
on conformal mapping,” Opt. Express **23**, 19901–19910
(2015). [CrossRef]

**6. **H. Ohno and T. Usui, “Gradient-index dark hole
based on conformal mapping with etendue conservation,”
Opt. Express **27**,
18493–18507 (2019). [CrossRef]

**7. **H. Ohno and K. Toya, “Reconstruction method of
axisymmetric refractive index fields with background-oriented
schlieren,” Appl. Opt. **57**, 9062–9069
(2018). [CrossRef]

**8. **H. Ohno and K. Toya, “Scalar potential
reconstruction method of axisymmetric 3D refractive index fields with
background-oriented schlieren,” Opt.
Express **27**,
5990–6002 (2019). [CrossRef]

**9. **E. W. Marchand, “Fifth-order analysis of GRIN
lenses,” Appl. Opt. **24**, 4371–4374
(1985). [CrossRef]

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

**11. **M. Sluijter, D. K. G. de Boer, and H. P. Urbach, “Ray-optics analysis of
inhomogeneous biaxially anisotropic media,” J.
Opt. Soc. Am. A **26**,
317–329 (2009). [CrossRef]

**12. **Y. Nishidate, “Closed-form analytical
solutions for ray tracing in optically anisotropic inhomogeneous
media,” J. Opt. Soc. Am. A **30**, 1373–1379
(2013). [CrossRef]

**13. **J. E. Gomez-Correa, V. Coello, A. Garaza-Rivera, N. P. Puente, and S. Chavez-Cerda, “Three-dimensional ray tracing
in spherical and elliptical generalized Luneburg lenses for
application in the human eye lens,” Appl.
Opt. **55**,
2002–2010 (2016). [CrossRef]

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

**15. **G. W. Forbes, “On variational problems in
parametric form,” Am. J. Phys. **59**, 1130–1140
(1991). [CrossRef]

**16. **E. Hairer, C. Lubich, and G. Wanner, “Geometric numerical
integration illustrated by the Störmer/Verlet method,”
Acta Numer. **12**,
399–450
(2003). [CrossRef]

**17. **R. D. Ruth, “A canonical integration
technique,” IEEE Trans. Nucl. Sci. **30**, 2669–2671
(1983). [CrossRef]

**18. **E. Forest and R. D. Ruth, “Fourth-order symplectic
integration,” Physica D **43**, 105–117
(1990). [CrossRef]

**19. **H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by
simulating adiabatic bifurcations in nonlinear Hamiltonian
systems,” Sci. Adv. **5**, eaav2372
(2019). [CrossRef]

**20. **T. Hirono, W. W. Lui, K. Yokoyama, and S. Seki, “Stability and numerical
dispersion of symplectic fourth-order time-domain schemes for optical
field simulation,” J. Lightwave
Technol. **16**,
1915–1920 (1998). [CrossRef]

**21. **T. Hirono and Y. Yoshikuni, “Accurate modeling of
dielectric interfaces by the effective permittivities for the
fourth-order symplectic finite-difference time-domain
method,” Appl. Opt. **46**, 1514–1524
(2007). [CrossRef]

**22. **H. Fan and H. Lu, “Symplectic wavelet
transformation,” Opt. Lett. **31**, 3432–3434
(2006). [CrossRef]

**23. **H. Yoshida, “Construction of higher order
symplectic integrators,” Phys. Lett.
A **150**, 262–268
(1990). [CrossRef]

**24. **R. Winston, J. Minano, and P. Benitez, *Nonimaging Optics*
(Academic,
2005).

**25. **H. Ohno, “Design of a coaxial light
guide producing a wide-angle light distribution,”
Appl. Opt. **56**,
3977–3983 (2017). [CrossRef]

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

**27. **R. K. Luneburg and Brown
University, *Mathematical Theory of Optics*
(Providence,
1944).