Black Hole Glimmer Signatures of Mass, Spin, and Inclination

Gravitational lensing near a black hole is strong enough that light rays can circle the event horizon multiple times. Photons emitted in multiple directions at a single event, perhaps because of localized, impulsive heating of accreting plasma, take multiple paths to a distant observer. In the Kerr geometry, each path is associated with a distinct light travel time and a distinct arrival location in the image plane, producing black hole glimmer. This sequence of arrival times and locations uniquely encodes the mass and spin of the black hole and can be understood in terms of properties of bound photon orbits. We provide a geometrically motivated treatment of Kerr glimmer and evaluate it numerically for simple hotspot models to show that glimmer can be measured in a finite-resolution observation. We discuss potential measurement methods and implications for tests of the Kerr hypothesis.


INTRODUCTION
Spinning supermassive black holes likely power astrophysical relativistic jets via the Blandford-Znajek mechanism (Blandford & Znajek 1977). To probe the spin-jet connection through observation, it is necessary to understand the properties of spacetime near the hole and to have an accurate model of the accreting plasma around the hole. If spacetime is described by the Kerr metric, then its properties are uniquely determined by the mass and angular momentum of the central black hole.
Many spin measurement methods propose measuring the size of the accretion disk, the qualities of disk oscillations, or the deviation in line profiles due to black hole spin (e.g., Hanawa 1989;Kojima 1991;Laor 1991;Kato 2001;Miller 2007). These measurement techniques rely on an accurate understanding of the accretion flow and are therefore subject to uncertainties in the plasma physics model. The 2017 Event Horizon Telescope observation of the black hole at the center of the galaxy M87 provided the first direct horizon-scale observation of a black hole (Event Horizon Telescope Collaboration et al. 2019); however, its ability to constrain the spin of the hole is also limited by the modeling uncertainties.
In the Kerr metric, photons can travel along bound (but unstable) orbits in space around black holes. The properties of these orbits are determined solely by the spacetime geometry, i.e., by the mass and angular momentum of the black hole. In the image plane, the bound orbits produce a characteristic critical curve whose size and shape are also set by the spacetime geometry. Since the critical curve is determined solely by the spacetime, it is independent of the accretion model, and thus it provides a consistent, measurable signature that can directly probe the hole's properties. Measurement strategies to infer spin from the shape of the curve have been proposed in the past (e.g., Falcke et al. 2000;Takahashi 2004;Bambi & Freese 2009;Hioki & Maeda 2009;Younsi et al. 2016;Johnson et al. 2020). Related proposals have suggested testing the hypothesis that spacetime is Kerr by looking for deviations in the curve's shape (e.g., Amarilla et al. 2010;Tsukamoto et al. 2014;Amarilla & Eiroa 2013;Medeiros et al. 2020).
The bound orbits allow light near a black hole to orbit it multiple times before escaping to an observer. Since emitters can radiate in multiple directions simultaneously, two rays produced by the same source may orbit the hole a different number of times. Light signals from subsequent orbitings will be separated by a time delay set by the length of a complete winding around the hole. These delays will cause the source to echo in the image plane. Since the echo period is a function of path length, it is intrinsically tied to the underlying bound orbits and can be measured to infer spin. Constraining spin by measuring dominant echoes has been considered in the past (e.g., Broderick & Loeb 2005;Moriyama & Mineshige 2015;Saida 2017;Thompson 2019;Moriyama et al. 2019;Gralla & Lupsasca 2020).
In detail, the Kerr geometry produces a rich spectrum of echo time delays associated with resonant bound orbits. These resonant echo delays are closely related to the black hole quasinormal-mode spectrum in the eikonal limit (see Yang et al. 2012), which has been studied in the context of gravitational-wave ringdown and measuring mass and spin (e.g., Berti et al. 2006;Buonanno et al. 2007;Berti et al. 2007).
Each echo maps to a distinct arrival location in the image plane; taken together, the set of echoes produces a characteristic black hole glimmer that encodes the mass and spin of the hole. Since the mechanism that produces these echoes is driven purely by the spacetime, the black hole glimmer signature is separable from the source emission model. If glimmer can be measured precisely, it is possible to test the Kerr hypothesis and infer the black hole mass and spin, even without a detailed understanding of the emission source. We provide a geometrically motivated treatment of Kerr glimmer and demonstrate that the glimmer signature of a hotspot can be measured even in a finite-resolution observation.
This article is structured as follows. In Section 2, we review the salient features of the Kerr geometry in the context of bound orbits, and in Section 3 we describe the observable properties of the bound orbits. We describe Kerr glimmer and provide an example measurement in Section 4, and we provide a brief discussion in Section 5.

KERR GEOMETRY AND BOUND ORBITS
Black holes (in vacuum) in general relativity are described by three quantities: mass, angular momentum, and charge, although it is unlikely that supermassive black holes with a dynamically important charge exist in nature. Charge-neutral holes with non-zero angular momentum are described by the Kerr metric. We use geometrized units with G = c = 1 and write angular momentum J in terms of the conventional dimensionless spin parameter a * ≡ J/M 2 where M is the mass of the hole. Hereafter, we set M = 1. In Boyer-Lindquist coordinates x µ = (t, r, θ, φ), the Kerr line element is (Bardeen et al. 1972) The Kerr geometry admits a set of spherical bound orbits-orbits with fixed radial coordinate-for both massive and massless (photon) particles near black holes. The set of all spherical photon orbits is known as the photon shell. The properties and consequences of the photon shell on black hole images have been studied both in the context of non-spinning holes (see especially Luminet 1979) and the more general spinning Kerr hole (e.g., Bardeen 1973). One of the first full treatments of spherical photon orbits was presented by Teo (2003), who also provided a convenient categorization of the types of such orbits. In this section, we review several key features of the bound photon orbits. Appendix A details the results described in this section. We do not consider black holes with extremal spin.
For a given spin, the set of all bound spherical photons orbits can be parameterized by the radii of the orbits. These radii lie continuously in the range r − ≤ r ≤ r + , where Only the two extremal orbits at r ± are confined to the midplane. The prograde orbit at r = r − rotates around the hole in the same direction as its spin, and the retrograde orbit rotates opposite the spin at r = r + . The other orbits at intermediate radii oscillate between symmetric minimum and maximum latitudes θ ± . For holes with nonzero spin, the latitudinal oscillation period is different from the azimuthal φ period, so after a full φ orbit around the hole (φ → φ + 2π), the geodesic will not return to the same θ coordinate. The magnitude of the precession can be written in terms of the deviation in φ from one latitudinal cycle to the next. In the case of no precession, ∆φ would be 2π. In general, it is given by where K and Π are complete elliptic integrals 1 of the first and third kinds, the roots of the latitudinal poten-tial are and the constants of motion Φ(r) and Q(r) are given by Figure 1 plots the trajectory of bound orbits in the θφ plane and illustrates this precession effect, where one complete latitudinal cycle does not correspond to an azimuthal displacement of 2π. The spread in ∆φ increases with a * . By studying ∆φ, we can identify which orbits are also closed-which geodesics return to their original positions and orientations. Closed (resonant) orbits correspond to the case that ∆φ/2π is a rational number = p/q in simplest form. Infinitely many closed orbits exist, but orbits with small p, q are the most interesting since they have the shortest path lengths. Figure 1 shows the trajectories of several closed orbits both in threedimensions and as projected on the latitude-azimuth plane.
We also compute the time delay for one complete latitudinal cycle of the orbit: where E is the complete elliptic integral of the second kind. In the case of closed orbits, the time delay to complete a full cycle is a function of the number of latitudinal oscillations per complete cycle ∆T = q ∆t 2π/∆φ.

THE IMAGE OF THE PHOTON SHELL
We now consider the signature of the photon shell as seen by an observer far away from the hole. Since the orbits comprising the shell are unstable, photons whose paths deviate slightly from the precise trajectories described above will either fall onto the hole or escape to infinity where they can be captured by an observer. The exponential instability can be described in terms of how the deviation δr between a geodesic's radial position and the radial position of the corresponding bound orbit increases (or decreases) after n φ azimuthal cycles around the hole This equation defines the Lyapunov exponent γ φ , which is given by where γ θ is taken to be consistent with Johnson et al. (2020) and governs the deviation after one latitudinal θ half-cycle (from midplane to extremum back to midplane) Rather than consider the source-to-observer model, in which we start with all geodesics that are emitted from a source and select only those that make it to the observer, it is convenient to switch to the observerto-source model, where we begin with the set of all geodesics that intersect the image plane. For a camera at infinite distance, all photons incident on the camera arrive parallel.
Following Bardeen (1973), we parameterize geodesics that intersect the image plane according to their impact parameters x and y (α and β in Bardeen). By convention, the y axis is aligned with the projection of the black hole spin axis on the image. It is sometimes more convenient to write in terms of polar coordinates on the image ρ = x 2 + y 2 and ϕ = arctan y/x (see the top left panel of Figure 2).
Geodesics at large impact parameter (far from the image center) remain far from the hole and barely feel its influence. As the impact parameter decreases, however, the geodesics become increasingly bent, and eventually they wrap around the hole and undergo latitudinal oscillations. The set of impact parameters for the geodesics that intersect the midplane n times defines the nth subring on the image.
Nested subrings produce demagnified images of space. The magnification is given by the image area of the subring and scales according to the Lyapunov treatment above. If most of the emission is produced near the black hole, then it will be sampled n times in the nth subring, and so the ratio of flux (area-integrated intensity over the subring on the image) between the n and n + 1 subrings will be given by In the limit that n goes to infinity, geodesics wrap around the hole infinitely many times and are effectively trapped by the bound orbits described above. The set of all impact parameters with n → ∞ defines the critical curve on the image. The region within the critical curve is sometimes called the black hole shadow.
Points on the critical curve correspond to different bound orbits-different radii-within the photon shell, and so the path of the curve can be parametrized by r (see Bardeen 1973 and also Johannsen 2013): The shape of the critical curve is thus a function of two parameters: the spin of the hole (which defines the set of allowed bound orbits) and the viewing inclination i (which sets effective upper and lower limits on the accessible bound orbit radii according to their angular momenta-see Appendix B for more detail). Figure 2 shows the shape of the critical curve for (top) black holes with different spins and (bottom) the same hole observed at different inclinations. It also shows the mapping between points along the critical curve (parametrized by ϕ) and the bound orbits they correspond to (parametrized by r).

MEASURING KERR GLIMMER
Although messy gas dynamics determine the larger image features, the photon shell produces a separable, unique signature in the image domain. This signature is independent of the gas dynamics model, so it provides a direct way to measure the properties of the underlying black hole. The presence of an infinite number of subrings on the image means that an emission source will be imaged an infinite number of times, albeit with exponentially decreasing flux from one instance to the next. The image in the nth subring comes from light that has gone around the hole n times, and so the subring images will echo with a period equal to the light travel time around the hole.
The aggregate signature produced by an emission source depends on the characteristics of the source, but we can provide an initial analysis by making two re- marks in the context of a simplified model. First, since different positions along the critical curve correspond to different bound orbits and path lengths, the delay between subsequent imagings will be a function of position on the curve. Second, if the source emission is localized in space, then an orbit that probes the source must return to the same localized area in order for an echo to be excited along that geodesic. Since bound orbits precess, it may take many revolutions around the hole before a trajectory passes through the source again, and since flux decreases exponentially with n φ , the precession can render some orbits practically echoless. The latter criterion is the most restrictive and is relaxed in the second model we present. In practice, the echo response function to an individual emission event is complicated since it depends on the size and duration of the source. The details are complicated further by the initial transient response, which is set by source size, duration, and position. Rather than attempt a full treatment of the emission variability near the hole, we consider simplified hotspot models.
Each hotspot is taken to be a transient, isolated Gaussian blob with a time-and position-dependent emissivity where δr is the distance to a point in the midplane of the system on a Keplerian orbit at radius r 0 , t 0 ≡ t − δt sets the time when the hotspot is brightest, σ r and σ t describe the width of the hotspot in space and time, and where Ψ(t) is a bump function in time which forces a smooth decay to zero emissivity.  Figure 3. Echo delay times due to resonant orbits as a function of black hole spin. As spin increases, the number of accessible short-time-delay echoes increases. Color encodes the location of the echo on the critical curve. Only perfect resonances are shown in this plot. Bound orbits that undergo an even-to-odd ratio of latitude-to-azimuth oscillations will produce half-period echoes for emitters in the midplane.

Point source emission
We start by considering the response function produced by an anisotropic point source emitter as measured on the critical curve. Since echoes are produced when geodesics sample the same emission source multiple times, then in the limit as σ r → 0, the only localized (non-transient) echoes that are produced must come  As σr increases, more bound orbits pass through the emission region and thus more echoes are excited. The periodic signal observed at ϕ ≈ 7π/16 is produced by a half-period orbit. Bottom row: response produced by the same hotspot (with σr = 0.4 M) for different black hole spins. As above, the echoes near π/2 for the first two panels and near 3π/8 for the last one are due to half-period orbits.
from closed bound orbits, since they are the only ones that return to the exact same point. Other geodesics will either miss the emitter or not arrive at the correct observer location. The potential for a source to be imaged but not echo is what differentiates Kerr from Schwarzschild, since all spherical orbits in Schwarzschild are closed. Even though an infinite number of closed orbits exist, only orbits with short path lengths will produce observable echoes since the sensitivity required to detect echoes from multiple turnings increases exponentially with the number of turnings. Thus, at each spin, there is a limited number of fixed, position-dependent observable echo time delays. Figure 3 plots the allowed delays due to closed orbits as a function of spin. Because each delay is due to a particular bound orbit, it maps to a distinct position along the critical curve. The allowed delays in the figure are colored according to the angle ϕ along the critical curve where they appear. The time delays for the prograde and retrograde orbits are always accessible since they lie within the equatorial plane and thus always pass through the same points. The closed orbits can be computed by identifying which radii in Equation 4 correspond to (the reciprocals of) the first few levels of a Stern-Brocot tree.
The radius-dependent precession means that the time delay for a geodesic that circles the hole twice will not be twice that for a geodesic that only circles it once before closing. Since flux decreases quickly with time delay, it is challenging to observe the higher order resonances.
Finally, we note that emitters confined to the midplane are special since they admit a secondary set of half-delay echoes. Closed orbits with an even-to-odd ratio between the number of azimuthal and latitudinal cycles return to the same point (θ, φ) after having fol-lowed only half of the full closed path. Although the geodesic will have returned to the same position, it will follow a midplane-reflected trajectory since it has only completed half of its full orbit. We caution that light emitted in such different directions may not be coherent and thus may be harder to observe. retrograde-centered ϕ ∈ (π/4, −π/4) (blue), progradecentered (green) ϕ ∈ (3π/4, 5π/4), and two remaining (red) regions. The prograde orbit echoes are strongest because they have the smallest Lyapunov exponents. The initial transient to t ≈ 60 GM/c 3 is produced during the time that the hotspot is active and orbiting the hole.

Finite-width sources
Real world emission sources have non-zero width. In the context of glimmer, increasing the width of an emission sources softens the condition that a geodesic must return to the same point in (θ, φ) space for an echo to be produced, since geodesics will resample the fiducial source feature as long as the deviations in their positions are smaller than the size of the feature. These more permissive conditions broaden the set of accessible time delays, and therefore the set of ϕ that will exhibit clean echoes increases. Thus, as the source width is increased, the echo response becomes richer as the signal spreads across the bound orbits. One can understand the echo response function by plotting the intensity of light observed along the critical curve as a function of time to visualize how different echoes are excited by a source impulse.
The top row of Figure 4 shows how the echo response changes as the source size grows. If the source is smaller than the width of the photon shell, it cannot echo along the entire critical curve. As the source width increases, the ϕ width of each response blip increases since nearly closed orbits will begin to resample the source. This effect can be seen as vertical blip height increases from left to right across the panels in the figure. Since the precession of near-miss orbits is a function of their radii, the blip profiles skew with time.
The bottom row of Figure 4 shows the echo response of an emission source with fixed width σ r = 0.4 and fixed r = 2.4 radial coordinate but for black holes with different spin parameters ranging from 1/2 to 31/32. As spin increases, the range of angles on the critical curve that are excited changes from the prograde ϕ = π position to the retrograde ϕ = 0 position. The increased spread in ∆φ as a * increases can be seen as the echo features at large time delay become warped.
While the clearest glimmer signatures are encoded in the time-and position-dependent intensities along the critical curve, it may not be feasible to make a measure with sufficient resolution to measure them. Nevertheless, one can still probe the echoes by comparing light curves in different regions of a full, under-resolved image. Figure 5 shows an angular decomposition of the light curve produced by a hotspot (r 0 = 3 M, σ r = 0.8 M, and σ T = 20 M) orbiting around black holes with three different spins a * = 1/2, 3/4, and 31/32. After the initial transient from the hotspot orbiting the hole decays, the full light curve signal is dominated by the prograde echo period.
Nevertheless, as Figure 5 shows, it may be possible to resolve different echoes if the image is divided into quadrants and the light curves in each quadrant can be measured independently. Even if the precise time delay between peaks in a single light curve is not measurable, it may be possible to differentiate between the echo periods of the different curves. A rough comparison of the delay periods can be used to identify which side of the image corresponds to the prograde orbit and thus infer the orientation of the spin axis of the hole.

DISCUSSION
We have described black hole glimmer, a position-and time-dependent effect in black hole images that is due to bound and closed photon orbits. The set of echo periods and arrival locations that comprise glimmer uniquely encodes the properties of the spacetime and provides a direct way to measure black hole mass and angular momentum. Since glimmer is determined only by the underlying geometry, its signature is separable from astrophysical and plasma uncertainties in the emission model. Black hole glimmer provides a robust, independent probe of black hole mass and angular momentum and makes precise predictions that can directly test the Kerr hypothesis. Evaluating the feasibility of a precise measurement is beyond the scope of this paper; however, we make several remarks.
Since it may be difficult to measure time-dependent intensity along the critical curve with high spatial precision, we have shown that the approximate magnitude and orientation of the spin axis can be inferred from glimmer by comparing echo periods and magnitudes in low-resolution components of a light curve decomposition. If astrophysical processes are not correlated on glimmer timescales, then such a measurement could be performed by autocorrelating resolution elements in an image or by correlating visibility amplitudes taken on long baselines. The energy received in subsequent echoes decreases exponentially, but glimmer is a robust, stable, achromatic feature, so measurements from multiple epochs, locations, and frequencies can be integrated to boost the measurement signal. An analytic analysis of the measurability of correlations along the critical curve as viewed by a nearly polar observer will be presented in Hadar et al. (2020, in preparation).
If the emission spectrum has a characteristic frequency, then it may be possible to measure glimmer by comparing the dominant echo periods of different components of the spectrum. Since each bound orbit follows a different latitudinal profile, different orbits will intersect the source at different angles relative to the source's net movement. Frequency shift is controlled by this angle, and so different orbits (and thus different echo periods) will peak along different characteristic frequencies.
In our treatment, we assumed that the emission was independent of frequency; in sources with non-trivial frequency dependence, redshift effects can change the intensity of light received along different segments of the critical curve and may influence the sensitivity required for a measurement. We also neglected the effect of optical depth, which decreases the strength of high order echoes compared to the analytic Lyapunov treatment. Since optical depth is a function of wavelength, it may be desirable to observe at frequencies where the optical depth of the plasma is minimal.
Our toy hotspot emission model may not be representative of true astrophysical scenarios, and the condition that the hotspot remain in the equatorial plane means it will excite echoes at frequencies that are inaccessible to non-midplane emitters. Moreover, our treatment neglected the initial transient (which is a strong function of the position, shape, and dynamics of the source). An analysis of general relativistic magnetohydrodynamic simulations may be required to study the detailed complexities of a true observation; an analytic treatment of correlations must faithfully reproduce the spatial structure of the emission source, since source position determines which echoes are excited.
It may also be possible to measure glimmer produced by an emitter far from the hole. If a coherent light source is located directly behind the hole relative to the observer, then the light paths between the source and the observer will undergo the same lensing effects as they pass the hole. This sort of scenario might be used to measure properties of systems without strong emission sources, like stellar mass black holes.

ACKNOWLEDGMENTS
I am grateful to Avery Broderick, Paul Chesler, Shahar Hadar, Alex Lupsasca, Kotaro Moriyama, Ramesh Narayan, Nico Yunes, and especially Charles Gammie for conversations and comments that greatly improved the text. This work was supported by NSF grant AST 17-16327 and by a Donald C. and F. Shirley Jones Fellowship.

A. KERR EQUATIONS OF MOTION
The Kerr metric is cyclic in the t and φ coordinates, and so it has two killing fields ∂ t and ∂ φ . These corre-spond to the conserved conjugate momenta p t and p φ that are canonically associated with (negative) energy at infinity and angular momentum about the spin axis of the hole. Carter (1968) identified a third conserved quantity Q that is associated with a second-order killing tensor field and is physically related to the motion of a particle as it passes through the midplane.
For a test particle with four-momentum p µ , the three conserved quantities (A1)-(A3) along with the particle's mass −µ 2 = p µ p µ uniquely determine geodesics in the Kerr spacetime.
Photons follow null geodesics x µ (λ), where λ is an affine parameter such that u µ =ẋ µ ≡ dx µ /dλ and u µ u µ = 0. The equations of motion for null geodesics around a black hole are scale invariant with respect to the mass of the hole are can be written as a oneparameter set of ordinary different equations (Carter 1968;Bardeen et al. 1972) Since photons are massless and their paths are independent of their energies E, it is convenient to normalize both L z and Q by E to define new constants of motion Φ and Q. These two constants can be written in terms of the orbit radius r and are often used to parameterize the bound orbits around a hole of a given spin (see the equivalent Equations 6 and 7).
By definition, spherical orbits lie at fixed, unchanging radii,ṙ =r = 0. By solving Equation A5 for these two conditions and rejecting non-physical solutions (see Teo 2003 for more detail), we find that spherical orbits must lie between r ± = 2 1 + cos 2 3 cos −1 ±a * .
It is convenient to introduce a new variable for the latitudinal coordinate θ. Substituting (A8) and (A9) and writing u ≡ cos θ, Equation A6 can be written as a fourth order polynomial in u: The roots of the above expression are given by (see also Equation 5) In terms of u, the two other equations of motion arė To compute the ∆φ azimuthal precession identified in §2, we integrate dφ/du over four quarter cycles Here, we have written the answer in terms of complete elliptic integrals by expressing du in terms of the product written in Equation A13. The ∆t integral reported in Equation 8 is solved in the same way. For a * = 0, the ratio ∆φ/2π cannot be one, and thus a complete azimuthal cycle will not correspond to a complete latitude cycle. The sign of ∆φ corresponds to the net displacement of the orbit (as either prograde or retrograde) and mirrors the sign of Φ. The radius of the polar orbit, which has Φ = 0, is given by (e.g., Teo 2003) Varying the inclination angle between the spin axis of the black hole and line of sight to the observer changes the shape and size of the critical curve by restricting the set of bound orbits that are accessible to the observer. The radial range of accessible bound orbits decreases Left panel: Range of allowed Boyer-Lindquist radii probed by the critical curve for different black hole spins (vertical axis) and different inclinations (shades of blue in 15 degrees steps). The set of all bound orbits for a given spin lie within the horizontal range delineated by the two black lines. Center/right panels: The relation between inclination and the map between echo periods and position on the critical curve. The full range of ϕ on the curve, given by the colored lines in the center panel, decreases as inclination goes from 90 to 0 degrees. from the edge-on case at i = π/2 to the single Φ = 0 orbit (Equation A18) at i = 0.
For an arbitrary inclination, the minimum and maximum visible bound radii r are given by roots of Equation B19. Figure 6 shows the allowed radii as a function of spin and inclination and also plots which portion of the ring response that inclination will sample. As inclination is decreased, the arrival time of the primary signal along the critical curve may be increasingly warped due to the azimuthal location of the source. 0 = η + a 2 * cos 2 (i) − ξ 2 cot 2 (i) = a 2 * cos 2 (i) − a 2 * (r + 1) + r 2 (r − 3) 2 cot 2 (i) + r 3 r(r − 3) 2 − 4a 2 * a −2 * (r − 1) −2 (B19)

C. NUMERICAL IMAGE GENERATION
The ray traced results presented in this paper were produced using a custom version of the ipole code (Mościbrodzka & Gammie 2018). Most observer-toemitter codes like ipole solve the radiation transport equations along a single geodesic per resolution element (pixel). This approach is reasonable when the difference between neighboring geodesic trajectories is small (such as when the pixels are small or when the geodesics do not pass close to the photon sphere). Since we study the neighborhood of the critical curve, we must resolve differences between geodesics as they begin to wind around the hole. Thus in our case, the pixel-centered method with a fixed grid fails since the widths of the subrings decrease exponentially. We deal with this issue by adaptively concentrating resolution elements near the critical curve.
The modified code first identifies the geodesics that have the longest path lengths relative to their neighbors and then constructs a connected set of pixels (the set containing the identified geodesics) to refine. Each pixel is refined into a 3x3 set of pixels centered around the original geodesic. The process is repeated until a stopping criterion is met. Figure 7 shows: the image produced by ray tracing on a grid with eight refinement levels, a schematic of refinement levels near the critical curve, and the multiple self-similar subrings produced by each of the imaged turnings.