We propose a simple geometric optics analog of a gravitational lens with a refractive index equal to one at large distances and scaling like $n(r)2=1+C2/r2$, where C is a constant. We obtain the equation for ray trajectories from Fermat's principle of least time and the Euler equation. Our model yields a very simple ray trajectory equation. The optical rays bending, reflecting, and looping around the lens are all described by a single trigonometric function in polar coordinates. Optical rays experiencing fatal attraction are described by a hyperbolic function. We use our model to illustrate the formation of Einstein rings and multiple images.
I. INTRODUCTION
While knowledge of the general theory of relativity is required to fully explain curved spacetime, geometric optics is often used to illustrate the phenomena of gravitational lensing and Einstein rings. The logarithmic lens is a wellstudied example.^{1–7} One side of the logarithmic lens is flat, and its thickness is described by the expression $a\u2009ln\u2009(b/r)$ where a and b are constants and r is the distance from the axis of the lens. Experiments with a logarithmic lens illustrate the formation of Einstein rings from a point (star), making it a useful tool in the classroom. Interestingly, the base of a wine glass resembles a logarithmic lens and can be used for a quick and simple gravitational lens experiment.^{5,6} A more sophisticated experiment that involves multiple lenses can be realized by shaping a water surface using small rods to pull up a water–air interface.^{7} For small rods, the water surface is logarithmic close to the rods. Graded index optics offers more room for innovation and offers truly threedimensional spherical gradients.
The trajectory equations for the profile in Eq. (3) in polar coordinates (r, $\theta $) are remarkably simple and described by a single trigonometric function, $csc\u2009\theta =1/\u2009sin\u2009\theta $, or a hyperbolic function, $csch\u2009\theta =1/sinh\u2009\theta $. This simplicity is not entirely surprising. In mechanical systems in which the motion of particles mimics the motion of photons, the potential U(r) is directly related to the squared refractive index: $U(r)=E\u2212n(r)2/2$, where E > 0 is a constant.^{9,11} Hence, trajectory equations for particles in mechanical central force systems with $U(r)=E\u22121/r2$ are the same as for photons moving in an index described by Eq. (3).^{12,13} A more exotic example is that of a neutral atom and a charged wire where a $1/r2$ potential is present.^{14} For a Schwarzschild black hole, the ray trajectories are described in terms of Jacobi elliptic functions.^{15–17} The event horizon is present and some photons vanish within the black hole. While black hole mathematics is substantially more complex than the model presented here, similar behavior emerges: trajectories can be bent, looped back, or around a gravitational lens, or vanish within a gravitational lens at a singularity. Our model does not reach the depth of the general theory of relativity. However, just like the base of the wine glass, it provides a simple educational tool. In a simple way, it helps to explain the Eddington experiment, Einstein rings, Hubble and JWST images, and effects of dark matter.^{18–21} This simplicity makes this model accessible to a wide audience. Understanding of this paper does not require knowledge of the general theory of relativity.
The paper is organized as follows: Section II provides a derivation of the ray trajectory equation for a spherical refractive index using Fermat's principle and the calculus of variations. Section III offers a detailed analysis of ray trajectories and illustrates ray bending, reflection, looping, and fatal attraction. Section IV elaborates on delays encountered by optical rays looping around a gravitational lens. Section V illustrates the formation of Einstein rings. Sections VI and VII explore imaging in 2D and 3D, respectively.
II. SPHERICAL GRADED INDEX LENS
This section contains an analysis of the spherical graded refractive index profile. For completeness, we review the derivation of the ray trajectory equation.^{22} In addition to the ray trajectory equation, the concept of the flight invariant is introduced.
III. RAY TRAJECTORIES
The following observations can be made about rays described by Eq. (14):

A full trajectory is described for the argument of the cosecant ranging over one halfperiod, e.g., between 0 and π.

The argument of the cosecant is equal to zero for θ = 0. At that angle, the cosecant and the radius r approach infinity. Thus, $\theta 1=0$ can be interpreted as an entry angle.
 The radius r also approaches infinity when the cosecant argument is equal to π. Thus, the angle θ_{2}, such that $\theta 2B2\u2212C2/B=\pi $, can be interpreted as an exit angle, which yields$\theta 2=\pi BB2\u2212C2.$

The turning angle (the ray bend angle) is related to the exit angle and, depending on the sign of B, is equal to $\theta 2\u2212\pi $ or $\theta 2+\pi $.

Entry and exit angles are interchangeable (reciprocity).

A minimum positive value of cosecant is equal to 1 when its argument is equal to $\pi /2$. At that point, the ray trajectory is tangent to the lines of constant refractive index, and the distance from the origin is at minimum (periapsis) and is equal to $d=B2\u2212C2$.
The following observations can be made about rays described by Eq. (15):

A full trajectory is described for the argument of the hyperbolic cosecant ranging between 0 and infinity.

The argument of the hyperbolic cosecant is equal to zero for θ = 0. At that angle, the hyperbolic cosecant and the radius r approach infinity. Thus, $\theta 1=0$ can be interpreted as an entry angle.

As the argument of the hyperbolic cosecant approaches infinity, its value approaches zero. Thus, the trajectory terminates at the origin. This is the case of fatal attraction.
Two families of rays described by Eqs. (14) and (15) are shown in Fig. 2 for the constant C = 1. All rays enter at a zero angle $\theta 1=0$ from the right forming a bundle of parallel rays. Their initial distance from the x axis (the impact parameter) is equal to the flight invariant B whose values are stepped from 0.05 to 1.95 in increments of 0.1.
Those rays whose flight invariant B is smaller than C = 1 belong to the family described by Eq. (15) and follow a spiral trajectory that terminates in the origin. The decay rate of the spirals is slower as the flight invariant B (the initial distance from the x axis) approaches C = 1. The rays from this family experience fatal attraction and vanish at the singularity after tracing an infinite number of loops.
Those rays whose flight invariant B is greater than C = 1 are bent or looped and escape at the exit angle θ_{2} described by Eq. (18). They belong to the family described by Eq. (14). The bend is smaller for the rays that enter further away from the x axis. As the distance B from the x axis approaches C = 1, the bend can become arbitrarily large and can take a form of a loop or even multiple loops. Interestingly, the rays from this family can pass arbitrarily close to the origin and still escape.
The figure contains two highlighted regions for this family: yellow and green (shown in color online and as different shades in print). As parallel rays enter from the right, both highlighted regions have the same width $\Delta B=0.1$. However, the range of exit angles θ_{2} is different. The green region, bounded by impact parameters B = 1.25 and B = 1.35, covers a substantially wider range of exit angles θ_{2} than the yellow region bounded by impact parameters B = 1.85 and B = 1.95. This has interesting consequences when propagation in the opposite direction is considered. This time light originates in the plane of the figure and exits on the right. An observer located in the far distance along the positive x axis sees all point sources/stars that are located within the green region at the green exit on the right, and all the sources/stars located in the yellow region at the yellow exit on the right as marked in Fig. 2. In general, exits (viewing regions) on the right that are closer to C = 1 contain more star images. This property is explored further in a numerical example in Sec. VI.
IV. OPTICAL PATH
V. EINSTEIN RINGS
In Fig. 3, the observer is on the x axis in the far distance at the lower right corner of the figure. Two stars (red and blue) are on the x axis on the other side of the lens. The optical rays from the red star, which is closer to the origin than the blue star, form the red inner ring seen by the observer. The optical rays from the blue star form the outer blue ring. The radii of the rings are equal to the flight invariants B. The bend of the rays forming a ring is related to the exit angle of Eq. (18) and equal to $\theta 2\u2212\pi $. The positions of the stars are described by Eq. (28).
A gravitational lens does not focus a parallel bundle of rays into a single focal point. All rays of a bundle aligned with the x axis are contained in a plane that also contains the x axis, and all of them cross the x axis. However, the locations of the crossings depend on the impact parameter (flight invariant B). The rays whose impact parameter $B$ is smaller than C all terminate at the origin. The rays whose impact parameter $B$ is larger than C cross the x axis at the point P_{s} described by Eq. (28). Thus, a bundle of rays parallel with an axis and with the impact parameter $B>C$ is focused onto this axis^{6} predominantly on the opposite side of the lens. (Predominantly, because some rays may loop around the gravitational lens and cross the axis multiple times.) By reciprocity, the light source (star) positioned on the x axis (or any other axis) is observed as a ring and not a as full bundle of rays. Multiple stars on the x axis are observed as multiple rings of different diameter. A ring of diameter 2B originates from a star positioned at the point P_{s} described by Eq. (28). When a point source (a star) or an object of some diameter (a galaxy) is not on the x axis, the solutions become more complicated. They are described in the next two sections.
VI. 2D IMAGING
Four lowest order solutions for the point $P(\u22122,2)$ are illustrated in Fig. 4. In numerical calculations, the constant C is assumed to be equal to one. The flight invariant B can be positive or negative. Positive values describe cw trajectories, negative ccw trajectories. As illustrated in the figure, the observer positioned in the far distance on the x axis can see multiple images below and above the x axis. The images are located at distance $B$ from the axis. For the observer, this maps to angles B/R. The images that appear farthest away from the x axis correspond to the smallest bend of the optical rays. As the bend curvature increases, the images get closer to $\xb1C$, in our case ±1. This especially applies to the looped trajectories since their bend exceeds $2\pi $. For clarity of the figure, multiple loop rays are not shown, however, they produce images that are very close to the single loop images and close to $\xb1C$. There are no images that appear between –C and C.
In Fig. 5, a set of 100 random point sources/stars together with ray trajectories and the histogram of the exit points are shown. A uniform distribution in a twodimensional space was used to generate this set of stars. Of the multiple trajectories that connect each star to the exit point, only the trajectories with simple cw and ccw bends and no loops are shown. Each star can be observed above and below the x axis (blue or magenta trajectories). The displacement of the image depends on the curvature of the trajectory. The images of the stars that are located in the upper halfplane are also in the upper halfplane and their displacement is relatively small. However, the images of the stars that are in the lower halfplane appear in the upper halfplane as well after their ray trajectories are bent by the gravitational lens. These images are densely clustered above y = 1. Similar effect takes place for the images below the x axis. The images of the stars from the lower halfplane remain in the lower halfplane while their displacement is relatively small. However, the images of the stars from the upper halfplane that are imaged by the gravitational lens into the lower halfplane are strongly clustered below y = – 1. As a result, the distribution of images is highly nonuniform with strong peaks above y = 1 and below y = – 1 as illustrated in the histogram of exit points.
This 2D illustration can be generalized to 3D just by rotating the image in Fig. 5 about the x axis. Thus, in the 3D space very few rays emerge within the circle of the radius C centered on the x axis. These are the images of stars that are within that region between the gravitational lens and the observer. Just beyond the circle of radius C, the region is most crowded with star images gathered from a wide angle on the opposite side of the x axis. Further beyond, where the effect of a gravitational lens diminishes, the density of images also decreases and more closely approximates the real position of the sources. Mathematically, the density of images is related to Eq. (19): the acceptance angle $\Delta \theta 2$ is larger per $\Delta B$ as B gets closer to C. This once again indicates that the images (and their light) tend to concentrate around a gravitational lens. This is true for any observer looking toward the gravitational lens.
VII. 3D IMAGING
Since each trajectory is in a plane that contains the x axis, the position of each point of the image is described in polar coordinates by two variables: a distance from the x axis that is equal to the flight invariant B found from a numerical solution of Eq. (14), and a rotation angle ψ; in polar coordinates $(B,\psi )$. Figure 7 illustrates a circular object (disk) and its lowest order crescentshaped images. If the blue circle represents the outline of a galaxy positioned behind a gravitational lens, the red crescent shapes represent the observed warped images. The similarity to the recorded astronomical images is evident.^{21}
VIII. CONCLUSIONS
We analyzed a refractive index profile that qualitatively mimics the effects of a Schwarzschild gravitational lens. The profile provides two very simple ray trajectory equations using a single trigonometric or a hyperbolic function, simple enough to make it accessible to a wide audience. One solution describes rays that are bent but escape the lens; the other, rays that are captured by the gravitational lens. The escaping rays bend increasingly more as the lens center is approached. The deflection is unbounded, ranging from a miniscule bend to a reflection, a full loop around the lens, and multiple loops. There is a critical impact parameter $B=C$ below which the rays experience the fatal attraction and are captured within the lens. Analysis of the gravitational lens imaging in 2D reveals that images may be multiple and that the gravitational lens tends to gather images of stars from all directions and from every star in the sky. This is true for all observers. Analysis in 3D reveals that the images may be multiple and warped, as observed in images from Hubble and the JWST.
ACKNOWLEDGMENTS
We thank the reviewers for their valuable suggestions and our colleagues Carsten Langrock and Gerry Owen.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts of interest.