## Abstract

A universal method to design gradient-index (GRIN) optical elements is proposed here for a given desired light ray bundle. Fermat’s principle can be transformed into a spatial parametric ray equation where a spatial Cartesian coordinate is used as a parameter of the equation. The ray equation can thus be written in a time-independent form, which ensures that a refractive index distribution is in principle obtainable from a spatial light ray distribution. Based on the ray equation, an iterative GRIN mapping method using the neural network (NN) is then constructed to map a refractive index distribution that enables light rays to trace corresponding desired paths. Maxwell’s fisheye lens is used to demonstrate how well the GRIN mapping method works. The refractive index distribution is shown to be well reconstructed from only knowledge of the light ray paths.

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

## 1. Introduction

Gradient-index (GRIN) optics has widely been used for lenses with flat surfaces and for aberration-free lenses that are unlike traditional spherical lenses where aberration is inevitable. Several types of GRIN optical elements such as the optical black hole [1, 2], invisibility cloak [3], and silicon waveguides [4, 5] have been designed and have attracted attention. Although several ray-tracing methods [6–17] in GRIN media have been developed, there are relatively few methods for inversely finding the refractive index distribution from a given desired light ray bundle. Conformal mapping [18] is one method that provides a powerful tool for designing a refractive index distribution. However, it is often difficult to find an analytical function needed for the conformal mapping from virtual space to physical space.

Equation-driven unsupervised neural networks (NNs) have recently been developed to solve differential equations relevant to several physical systems [19–27]. The NN can be optimized through approximating a numerical loss, defined as a deviation of numerical calculation from the differential equation, to zero. Furthermore, the universal approximation theorem [24] of the NNs states that the NNs can approximate any function with arbitrary accuracy. This makes the NNs suitable approach for solving complicated problems governed by the differential equations.

Here, a universal method for designing GRIN optical elements is therefore proposed that is applicable to a given desired light ray bundle. The method is called a GRIN mapping method. Fermat’s principle can be transformed into a spatial parametric ray equation where a spatial Cartesian coordinate is used as a parameter of the equation. For the GRIN mapping method, a numerical loss can be defined as a deviation of numerical calculation from the spatial parametric ray equation. A refractive index distribution for the given desired light bundle can then be found by approximating the numerical loss to zero using the NNs. The GRIN mapping method can therefore provide the refractive index distribution, which enables light rays to trace corresponding paths of the given desired light ray bundle.

## 2. Concept of GRIN mapping

Figure 1 shows a conceptual overview of the GRIN mapping. Once a light ray bundle in a Cartesian coordinate of (*x*, *y*, *z*) is given, the GRIN mapping can provide a refractive index distribution that makes light rays trace the corresponding desired light ray paths. Fermat’s principle can be transformed into a spatial parametric ray equation where one spatial Cartesian coordinate, namely, *z*, is used as a parameter of the equation. A three-dimensional position vector ** Q** of a light ray can thus be described by a two-dimensional position vector

**of (**

*q**x*,

*y*) that is parametrized with

*z*. A three-dimensional unit directional vector

**of the light ray can also be described by a two-dimensional directional vector**

*U***that is also parametrized with**

*u**z*. Note that the magnitude of the unit directional vector

**is 1. The ray equation can thus be written in a time-independent form, which ensures that a refractive index distribution is in principle obtainable from the spatial light ray distribution.**

*U*## 3. GRIN mapping method

#### 3.1 Spatial parametric ray equation

A formulation of the GRIN mapping method is derived as follows. The traversal time, *T*, of a light ray in a GRIN medium between two given points A and B can be written with an infinitesimal time element, *dt*, as

*dt*can be transformed into an infinitesimal line element,

*ds*, with a refractive index

*n*and the speed of light

*c*in vacuum as Note that the speed of light

*c*in vacuum is constant and independent of space and time. The infinitesimal line element

*ds*can be further transformed into the Cartesian coordinate system of (

*x*,

*y*,

*z*) as where

*x*′ and

*y*′ are derivatives with respect to

*z*, and are written respectively as Using Eqs. (2) and (3) with the two-dimensional position vector

**of (**

*q**x*,

*y*), the traversal time

*T*of Eq. (1) can be written in terms of the Lagrangian:

*z*, is used as a parameter in the integrand of Eq. (7) instead of the time

*t*. Fermat’s principle indicates that the light ray path remains stationary under variation within a family of nearby paths. The following equations can then be derived by applying the variational method to Eq. (7) as [28–30]

*z*is used as the spatial parameter of the equation.

The two-dimensional directional vector ** u** can be defined as

**corresponding to a partial derivative of the Lagrangian with respect to**

*p***′ can be defined with both the refractive index**

*q**n*and the directional vector

**as Equations (8) and (9) can then be written in a unified form with the momentum vector**

*u***as**

*p**u*denotes the magnitude of the directional vector

**.**

*u*A GRIN mapping operator $\hat{{\boldsymbol \varepsilon }}$ for a function *n* of a position (** q**,

*z*) is defined here with a vector function

**of the position (**

*u***,**

*q**z*) as

**and the refractive index**

*u**n*as Note that Eq. (15) depends only on the spatial distributions of the directional vector

**and the refractive index**

*u**n*, which indicates that the distribution of the refractive index is obtainable using Eq. (15) once the distribution of the light ray direction is given.

#### 3.2 Iterative calculation using neural network

A numerical loss *g* for a function *n* is defined using the GRIN mapping operator with the directional vector ** u** as

*n*is equal to an ideal refractive index for a given desired light ray bundle, which can be easily seen from Eq. (15). In other words, the refractive index that makes light rays trace the corresponding desired paths can be found by approximating the numerical loss

*g*to zero. An iterative GRIN mapping method using Eq. (16) with a neural network can then be constructed to map a refractive index distribution for a given desired light ray bundle.

Equation (16) can be discretized on a regular grid using the derivatives calculated with central finite differencing. The total loss *G* over the entire regular grid can be written using Eq. (16) in terms of a regular grid number *i* of 1, 2, 3, …, *M* as

*q**, Δ*

_{i}*x*, Δ

_{i}*y*, and Δ

_{i}*z*denote step sizes, and

_{i}

*u**and*

_{i}*n*are defined respectively as The total loss

_{i}*G*represented by Eq. (17) is here called a total GRIN mapping loss.

Figure 2 shows a schematic GRIN mapping process using the NN [23] with the total GRIN mapping loss *G* represented by Eq. (17). A set of two-dimensional directional vectors of a light ray bundle on the regular grid, namely, *u*_{1}, *u*_{2}, …, *u** _{M}*, is set as input. A set of positions on the regular grid, namely, (

*q*_{1},

*z*

_{1}), (

*q*_{2},

*z*

_{2}), …, (

*q**,*

_{M}*z*), is also set as input. A set of real number variables of

_{M}*γ*

_{1},

*γ*

_{2}, …,

*γ*on the regular grid is provided to a front-end process by the NN. The front-end process provides a discretized refractive index

_{M}*n*with the

_{i}*i*of 1, 2, 3, …,

*M*using an update function

*f*that is a function of both the discretized position (

*q**,*

_{i}*z*) and the variable

_{i}*γ*. Note that the update function

_{i}*f*can be set freely to satisfy a desired boundary condition. The NN can be optimized iteratively by approximating the total GRIN mapping loss to zero. The iteration is counted by an epoch number

*l*. The discretized refractive index can be taken as output.

## 4. Results

One well-known GRIN optical element is Maxwell’s fisheye lens [31]. Maxwell’s fisheye lens is then used to demonstrate how well the GRIN mapping method works. The refractive index of the fisheye lens in a cross-sectional plane at *y* = 0 can be written as

*n*with the

_{i}*i*of 1, 2, 3, …,

*M*will be iteratively updated using the update function

*f*in the front-end process as shown in Fig. 2. In this work, the update function

*f*is defined with both the discretized position (

*q**,*

_{i}*z*) and the real number variable of

_{i}*γ*as

_{i}*r*of 1 where the radius is a distance from the coordinate origin. The refractive index distribution of the fisheye lens is therefore ensured to have a spherical boundary where the refractive index is 1. Note that the update function of Eq. (21) can be replaced by other functions to satisfy desired conditions.

A light ray path in the fisheye lens can be analytically derived using the ray equation with a constant *α* as

*x*,

*z*) = (

*x*

_{0},

*z*

_{0}), the

*α*can be written using Eq. (22) as The derivative of

*x*with respect to

*z*at the point (

*x*

_{0},

*z*

_{0}) can then be written using Eqs. (22) and (23) as

*g*for a function

*n*in the cross-sectional plane at

*y*= 0 can be written using Eq. (16) as

*u*can be written using Eq. (10) as The derivative of

_{x}*x*with respect to

*z*in Eq. (26) for Maxwell’s fisheye lens is represented by Eq. (24). Note that the derivative is always calculable if a light ray bundle is given.

The refractive index distribution can then be reconstructed from the light ray direction distribution of Eq. (26) with the gradient of Eq. (24) by approximating the total GRIN mapping loss to zero. Figure 3 shows reconstructed results by applying the GRIN mapping method to the light ray direction distribution. PyTorch [32] is used to build the NN with the Python language. In the NN, the number of the hidden layers having 50 neurons is set to 2. The total grid number *M* is set to 10000 where grid numbers in *x*- and *z*-directions are both set to 100. Perspective views of refractive index distributions reconstructed using the GRIN mapping method for various epoch numbers of *l* = 0, 4, 16, 1024, and 8192 are shown in (a). A perspective view of an ideal refractive index distribution is shown in (b). Each refractive index distribution *n* is color-contoured in the *z*-*x* plane (i.e., the cross-sectional plane at *y* = 0). A white solid line denotes the radius *r* of 1. The reconstructed refractive index at epoch number *l* of 8192 agrees well with the ideal refractive index.

Figure 4 shows a cross-sectional view of the refractive index distribution. The refractive index distributions with respect to the *z*-axis for various cross-sectional lines parallel to the *z*-axis, namely, *x* = 0, 0.2, 0.4, 0.6, and 0.8, are plotted. The horizontal axis denotes the *z*-axis, and the vertical axis denotes the refractive index. A solid black line denotes the refractive index reconstructed using the GRIN mapping method. An orange dashed line denotes the ideal refractive index. Note that the refractive index distributions are reconstructed using the GRIN mapping method from only the distribution of the light ray direction represented by Eq. (26). The reconstructed distributions agree well with the ideal distributions.

Here, a numerical deviation *σ* is set to a standard deviation of the reconstructed refractive index distribution from the ideal distribution within a radius of 1. Figure 5 shows the numerical deviations of the *σ* and *G* (the total GRIN mapping loss) with respect to the epoch number *l*. The horizontal axis denotes the epoch number, and the vertical axis denotes the numerical deviation. A solid black line denotes *σ*, and a dotted black line denotes *G*. The *σ* becomes less than 10^{−2} for epoch number *l* of more than 10^{3}. The total GRIN mapping loss *G* becomes less than 10^{−3} for epoch number *l* of more than 10^{3}.

## 5. Conclusions

In this work, a universal method to design GRIN optical elements is proposed for a given desired light ray bundle. First, Fermat’s principle is transformed into a parametric ray equation where one of the spatial Cartesian coordinates is used as a parameter of the equation as represented by Eqs. (8) and (9). The ray equation can be further transformed using a GRIN mapping operator of Eq. (14) into Eq. (15), which depends only on spatial distributions of a light ray directional vector and a refractive index. The ray equation can thus be written in a time-independent form, which ensures that the refractive index distribution is, in principle, obtainable using Eq. (15) from the spatial distribution of the light ray direction. Second, an iterative GRIN mapping method using the ray equation with the neural network (NN) is constructed to map a refractive index distribution that enables light rays to trace the corresponding paths of the given desired light ray bundle. The total GRIN mapping loss of Eq. (17) is defined as the deviation of the numerical calculation from the ray equation, which can be approximated to zero by the GRIN mapping method. Third, Maxwell’s fisheye lens is used to demonstrate how well the GRIN mapping method works. As shown in Figs. 3 and 4, the refractive index distribution is well reconstructed from only knowledge of the distribution of the light ray direction of Eq. (26) with the gradient of Eq. (24). Back scattered light rays and magnetic optical materials are not taken into account in this work. The GRIN mapping method, however, should be applicable to wide variety of light ray bundles. Once a desired light ray bundle is given, the GRIN mapping method can be used to find a corresponding refractive index distribution.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the corresponding author upon reasonable request.

## References

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

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

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

**4. **S. H. Badri and M. M. Gilarlue, “Silicon nitride waveguide devices based on gradient-index lenses implemented by subwavelength silicon grating metamaterials,” Appl. Opt. **59**(17), 5269–5275 (2020). [CrossRef]

**5. **S. H. Badri, M. M. Gilarluea, and S. G. Gavganib, “Ultra-thin silicon-on-insulator waveguide bend based on truncated Eaton lens implemented by varying the guiding layer thickness,” Photonics and Nanostructures - Fundamentals and Applications **39**, 100766 (2020). [CrossRef]

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

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

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

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

**10. **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**(2), 317–329 (2009). [CrossRef]

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

**12. **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**(8), 1559–128X (2016). [CrossRef]

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

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

**15. **R. D. Ruth, ““A canonical integration technique,” Nuclear Science,” IEEE Trans. on. NS **30**(4), 2669–2671 (1983). [CrossRef]

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

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

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

**19. **I. E. Lagaris, A. Likas, and D. I. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” IEEE transactions on neural networks **9**(5), 987–1000 (1998). [CrossRef]

**20. **J. A. Sirignano and K. Spiliopoulos, “Dgm: A deep learning algorithm for solving partial differential equations,” J. Comput. Phys. **375**, 1339–1364 (2018). [CrossRef]

**21. **M. Magill, F. Qureshi, and H. W. Haan, “Neural networks trained to solve differential equations learn general representations,” 32nd Conference on Neural Information Processing Systems (NeurIPS) (2018).

**22. **R.T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural Ordinary Differential Equations,” 32nd Conference on Neural Information Processing Systems (NeurIPS) (2018).

**23. **M. Mattheakis, D. Sondak, A. S. Dogra, and P. Protopapas, “Hamiltonian Neural Networks for Solving Differential Equations,” arXiv:2001.11107v2 [physics.comp-ph] (2020).

**24. **K. Hornik, “Approximation capabilities of multilayer feedforward networks,” Neural Networks **4**(2), 251–257 (1991). [CrossRef]

**25. **N. Akashi, M. Toma, and K. Kajikawa, “Design of metamaterials using neural networks,” Proc. SPIE **11194**, 111940U (2019). [CrossRef]

**26. **M. Qiao, X. Liu, and X. Yuan, “Snapshot temporal compressive microscopy using an iterative algorithm with untrained neural networks,” Opt. Lett. **46**(8), 1888–1891 (2021). [CrossRef]

**27. **X. Liu, Z. Wei, M. Li, L. Wang, Z. Liu, C. Yu, L. Wang, Y. Luo, and H. Y. Fu, “Experimental investigation of 16.6 Gbps SDM-WDM visible light communication based on a neural network receiver and tricolor mini-LEDs,” Opt. Lett. **46**(12), 2888–2891 (2021). [CrossRef]

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

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

**30. **H. Ohno and K. Toya, “Localized gradient-index field reconstruction using background-oriented schlieren,” Appl. Opt. **58**(28), 7795–7804 (2019). [CrossRef]

**31. **C. T. Tai, “Maxwell Fish-eye treated by Maxwell Equations,” Nature **182**(4649), 1600–1601 (1958). [CrossRef]

**32. **PyTorch, https://pytorch.org/