Ray-tracing simulation of gravitational lensing using a gradient-index model

Abstract. We present a method for modeling gravitational lensing as gradient-index lenses in the ray-tracing software CODE V. When applied to gravitational lensing by the Sun, these models yield results in agreement with theoretical and experimental results. The extension of this type of model to “lumpy” astronomical objects whose lensing behavior is difficult to predict with traditional methods is discussed.


Introduction
The techniques and insights particular to optical engineering have been applied to a variety of disciplines beyond optics. 1 In particular, the analysis of "ray" paths in analogs to gradient-index (GRIN) media 2 has yielded interesting results in fields as diverse as oceanography and coastal engineering, 3,4 acoustics, 5 seismic waves, 6 and chemical physics. 7 In this paper, we apply the tools of optical engineering to a problem that is, at once, familiar-since the rays here are indeed light rays-and unfamiliar-since the rays are bent not by materials but by the presence of massive gravitating bodies.
Gravitational lensing (i.e., the bending of light rays by a massive object) was one of the first predictions of Einstein's General Theory of Relativity to be experimentally confirmed. Eddington et al. 8 measured the gravitational deflection of light rays by the Sun during the total eclipse of the Sun in 1919. Results from these measurements (1.61 AE 0.30 arc sec and 1.98 AE 0.12 arc sec) were consistent with the value predicted by General Theory of Relativity of 1.75 arc sec, differing from the Newtonian prediction of 0.87 arc sec by a factor of nearly two. In the intervening years, agreement between theoretical and experimental values for ray deflection by the sun has been refined through the use of radio astronomy to better than three parts in ten thousand. 9 Knowledge of the lensing effects allows for estimation of the total mass of the galaxy acting as the lens. 10 Given that the amount of light (visible) matter in those galaxies is substantially less than the predicted mass of the galaxies, this method provides an estimate for the amount of dark matter present in them. Since the galaxies forming the lenses are typically several orders of magnitude further away from Earth and from the source of the light rays than they are thick, astronomers typically treat them as thin lenses (i.e., they assume that all the bending takes place in a single plane 11 ), using the models often taught in introductory physics courses to describe ray bending.
This method is, however, poorly suited to situations when a group of galaxies forms a "lumpy" lens, as is the case with visually striking examples such as the Cheshire Cat, 12 which consists of three galaxies acting as a lens for the light produced by several more distant galaxies. (The Cheshire Cat indeed takes its name from the grinning cat in Lewis Carroll's Alice in Wonderland. See it at NASA's Astronomy Picture of the Day, 2015 November 27 13 .) In addition to the simple deflection of light rays by gravitating bodies, astronomers have observed gravitational lensing of very distant galaxies by a somewhat less distant galaxy, including the production of multiple images. 14 In these cases, ray tracing cannot be accomplished by the simple application of the law of refraction since the index of refraction cannot be modeled as homogeneous. Instead, Fermat's principle governs ray propagation in such media: rays travel along paths of stationary optical path length, where the optical path length is the product of the refractive index and the path length. For example, the existence of four stationary paths is demonstrated in the four images formed by the lensing galaxy in the famous case of the Einstein Cross (see it at NASA's Astronomy Picture of the Day, 2017 December 17 15 ).
One solution to this difficulty is to model gravitational lenses as GRIN lenses. 16,17 In this paper, we present a method for modeling gravitational lenses as GRIN lenses in the ray tracing software CODE V. Although there are several methods to numerically integrate ray paths in a gradient index medium, CODE V uses a standard Runge-Kutta method. 18 Practical considerations, including the algorithmic efficiency and accuracy, of ray tracing in GRIN materials with this software are well documented. 19 We demonstrate that the method presented herein produces results consistent with theory and experiment for the case of gravitational deflection of grazing rays by the Sun. Then we discuss its potential for use for modeling a group of galaxies that form a lumpy lens.

Implementation of Models in CODE V
In this work, the Sun is simulated as a Schwarzschild gravitating body (i.e., a static, spherical mass distribution). For such a body, the index of refraction in spherical coordinates is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 4 7 7 where r is the radial distance from the center of mass of the object causing the gravitational field and r g is the Schwarzschild radius (i.e., the radius of the largest black hole with the same mass as the the object). 17 This index distribution is plotted in Fig. 1. For the purpose of numerical evaluation, the unit for simulation of solar gravitational lensing is selected as the solar radius (SR), which is ∼6.957 × 10 8 m. The Schwarzschild radius of the Sun is 2.95 × 10 3 m, which is equivalent to 4.24 × 10 −6 SR. The change of index of refraction from unity decreases as the point of interest moves away from the center of the Sun. At an infinite distance, the index of refraction reduces to unity. It can be concluded that, in the case of solar gravitational lensing, the majority of the index of refraction variation from the Sun's gravity occurs within 100 SR from the center of the Sun. For a numerical ray tracing simulation with finite computational resources, this conclusion provides an important guideline for setting up the simulation in CODE V. However, the gravitational field GRIN does not follow any of the existing GRIN conventions in CODE V. In order to define this kind of arbitrary gradient in a material, we employ the most rapid and robust method to allow ray tracing through arbitrary GRIN materials-the user-defined gradient (UDG) subroutine in CODE V. The UDG subroutine must be programmed to compute the index of refraction at any point in space as well as the explicit derivatives of the index of refraction. For a gravitational field that is represented by a GRIN material, the GRIN coefficients are related to the astronomical constants and coefficients of the body that is under simulation, such as its Schwarzschild radius.
To simulate the light deflection due to the gravitational lensing of the Sun, a user-defined GRIN material representing as a Schwarzschild gravitating body with the Schwarzschild radius equal to that of the Sun is programmed. Additional details regarding this GRIN material are provided in the Appendix. The simulation region is programmed to be a cylindrical region, 2 SR in diameter and 1000 SR in length, symmetrically positioned around the center of the Sun. A schematic of this simulation is presented in Fig. 2. When a geometrical ray, aimed to graze the surface of the Sun, is traced through the simulation region, the ray deflection angle between the input and output rays is recorded as 1.749 arc sec, in agreement with theoretical (1.75 arc sec) and experimental results, which, as noted in the introduction, agree with theoretical predictions to three parts in a thousand. (The advent of radio astronomy has obviated the need to look only for grazing rays during solar eclipses, and hence these measurements have been made for rays traveling at various distances from the Sun's center. For a grazing ray, this would be equivalent to 1.750 AE 0.005 arc sec:) Sen's computational model using the same expression for the index of refraction yields a ray deflection angle of 1.78 arc sec, so we see that the robust and well-tested codes of the UDG subroutine in CODE V yield better agreement with theory and experiment.
In addition to the simulation of the Sun, the ray deflection angles of Schwarzschild gravitating bodies with differing ratios between Schwarzschild radius and physical radius are evaluated. Note that this ratio is ≥1 for black holes. The relationship between the deflection angle and the Schwarzschild-physical-radius-ratio is linear when the Schwarzschild radius is very small compared to the physical radius of the gravitational body (as is the case for main sequence stars and even galaxies), as shown in Fig. 3. However, as the Schwarzschild radius becomes a noticeable fraction of the physical radius (or equivalently, as the density becomes very high, as it does  for objects like neutron stars), the deflection angle becomes nonlinear with the Schwarzschild radius, as shown in Fig. 4.
It is also possible to model a Kerr (i.e., rotating spherical) body in a similar fashion, 20 although several simplifying assumptions are necessary to reduce the complexity of that problem, as discussed in the Appendix. Note that treating the Sun as a Kerr body does not make a meaningful difference in the ray deflection angle.

Limitations of the Model
The programmed GRIN model for both Schwarzschild and Kerr bodies can be applied to different gravitational bodies by tuning the parameters. Additionally, more complex models can be programmed by combining multiple displaced bodies with appropriate parameters. However, some limitations must be addressed before more complex systems can be simulated usefully using the aforementioned methodology.
When modeling the gravitational lensing effect in imaging optical-design software such as CODE V or OpticStudio, there is an intrinsic trade-off between the range and the precision, because the software is primarily designed for the purpose of modeling imaging lenses for which an extraordinarily large dynamic range is not required. Although optical values can be evaluated with 12 significant digits in CODE V in double-precision floating-point format, ray tracing can be problematic when simulating gravitational lensing because the object is often very far from the lens compared with its spatial extent. It is, therefore, difficult to simulate the propagation of light on a large scale while simultaneously modeling the lensing effects that occur on a small spatial scale. When modeling imaging optics, it is common to consider objects that are far away as objects at infinity, in which case the accurate representation of very large object distance is not critical to the ray tracing results. This is true because the object distance is usually much greater than the effective focal length of the lens. However, in the scenario of modeling gravitational lensing, the lensing effect from an astronomical body is usually very weak, creating an exceptionally large effective focal length for the lens. Thus although the object is far away on an absolute scale, it cannot be treated as object at infinity. The accurate representation of these extraordinarily large numbers becomes more important, and in some cases, this can become a challenge when the software has limited dynamic range in its data structure.

Concluding Remarks
Our current models using existing imaging optical design software are capable and accurate in forward calculations of the lensing effect. The simulation method described is dependent on accurate knowledge of the astronomical configurations and parameters of the gravitational system. This information is often unknown or needs to be derived based on observations of the gravitational lensing effects. Performing backward calculations to solve for astronomical facts Fig. 4 Log-log plot of deflection angle versus ratio of Schwarzchild radius to physical radius of a gravitating body. This relationship is no longer linear as the Schwarzchild radius becomes an appreciable fraction of the physical radius. The ratios for common astronomical objects, including the Sun, the Milky Way, a white dwarf, and a neutron star are indicated for reference. from the observed data is not practical using the GRIN model within existing optical design software.
Ideally, custom ray tracing software must be developed to fully exploit the potential of the GRIN model for simulation of gravitational lensing, allowing for modeling of "lumpy" lenses and other complex phenomena. Recent work on the modeling of freeform GRIN in two dimensions to achieve a predetermined and complex irradiance distribution 21 is an interesting inverse problem, whose solution in three dimensions may prove useful to the inverse problem of determining the dark matter composition of these lumpy astronomical lenses. Additionally, biological media have many similarities to complex astronomical lenses, consisting as they do of many elements of different compositions and shapes, despite the large differences in scale. 22,23 Similarly, propagation of light through the atmosphere presents many related challenges. 24,25 Continued work to improve and expand optical engineering codes will make problems from across science and engineering disciplines more tractable.

Appendix: Gradient-Index Models for Spherical Gravitating Objects
In this section, we consider both Schwarzschild (static spherical) and Kerr (rotating spherical) gravitating objects.

Schwarzschild Gravitating Body
For a static, spherically symmetric gravitational field, the index of refraction in spherical coordinates is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 5 2 where r is the radial distance from the center of mass of the object causing the gravitational field and r g is the Schwarzschild radius. 17 This can be rewritten in terms of the Cartesian coordinates for ease of use in CODE V E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 7 2 nðx; y; zÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi For CODE V to treat this index distribution as a GRIN, n∇n is also necessary; due to the symmetry of this index distribution the three terms will be fundamentally similar. A bit of algebra after taking the partial derivatives yields

Kerr Gravitating Body
For a Kerr gravitating body (i.e., a rotating spherically symmetric object), the index of refraction in spherical coordinates is given in Refs. 16