Adaptive Analytical Ray Tracing of Black Hole Photon Rings

Introduction

According to general relativity, black holes display unique lensing behavior: for instance, any two spatial points outside the event horizon are connected by infinitely many light rays, each of which executes a different number of orbits around the black hole under its extreme gravitational pull [1, 2, 3, 4, 5]. As a result, images of a black hole surrounded by a non-spherical, optically thin emission region decompose into a sequence of superimposed layers indexed by photon half-orbit number $n$, with each layer consisting of a full lensed image of the main emission [5, 6, 7, 8, 9, 10, 11]. The direct ($n = 0$) layer typically displays a central dark area—the black hole—encircled by the weakly lensed, primary image of the accretion flow onto the hole, whose details depend sensitively on astrophysical conditions. On the other hand, the higher-$n$ layers arise from photons on highly bent trajectories that are strongly lensed to form a series of narrow “photon rings.”¹ These rings are usually stacked on top of the broader $n = 0$ emission and their shape rapidly converges (exponentially fast in $n$) to that of the “Kerr critical curve” [14]: a theoretical curve in the observer sky that corresponds to the apparent image of asymptotically bound photon orbits. In contrast to the astrophysics-dependent $n = 0$ image, this “$n \to\infty$ photon ring” is completely determined by the Kerr black hole—depending only on its mass, spin and inclination—and delineates its cross-sectional area in the sky.

Recently, 1.3 mm interferometric observations by the Event Horizon Telescope have resolved the horizon-scale emission from sources in the immediate vicinity of two nearby supermassive black holes: M87* [15], the central compact object at the core of our neighboring galaxy Messier 87, and Sgr A* [16], our own black hole at the center of the Milky Way. Their reconstructed images display a central brightness depression within a thick ring consistent with theoretical expectations for the direct image of the surrounding accretion flow [17, 18, 19, 20, 21].

However, these observations have thus far not provided any evidence for the presence of a lensed photon ring [22]. They are instead dominated by $n = 0$ photons [6, 7], which form the image layer that is more sensitive to the astrophysics of the flow than to purely relativistic effects [23, 24, 25]. As a result, the bounds placed on possible deviations from general relativity are on the order of several percent [26, 27], and comparable to the constraints derived from gravitational-wave observations of stellar-mass, binary black holes [28], or x-ray spectroscopy measurements of low-mass binaries and active galactic nuclei [29]. By contrast, future detections of orbiting ($n \ge1$) photons will open a new window into strong gravity and enable higher-precision probes of the Kerr geometry, since it is this orbiting light which forms the part of the image—the photon ring—that belongs to the black hole itself, rather than to its plasma.

There are three major obstacles to measuring a photon ring. First, since the photon rings are exponentially narrow (in $n$) features, resolving them requires interferometric observations on exponentially long baselines [7]. At the current observing frequency of 230 GHz, even Earth-spanning baselines are too short to resolve the first ring, but it should become accessible to next-generation space-based interferometers. In particular, SALTUS (the Single Aperture Large Telescope for Universe Studies) is a bold proposal—currently a contender for NASA’s upcoming Probe mission—to launch within the next decade a spacecraft far enough to access the first two rings of M87*.

Optical depth poses a second hurdle: even though photons could in principle circumnavigate the black hole indefinitely (albeit unstably), in practice, those that traverse its emission region multiple times are eventually reabsorbed by the matter they intersect—an effect that cuts off image layers past some $n > 0$. Nevertheless, since absorptivity decreases with photon energy, the first few rings still ought to be present in images taken at sufficiently high frequencies. Simple models suggest that at 230 GHz, the $n = 1$ ring is always visible while the $n = 2$ ring may only be marginally observable, whereas both rings should be clearly visible at 345 GHz [11]. State-of-the-art simulations of general-relativistic magnetohydrodynamic (GRMHD) flows [30] also confirm these expectations [31].

For this reason, SALTUS is slated to simultaneously observe at both frequencies. M87* makes for a particularly exciting prospective target because a measurement of its $n = 2$ ring diameter could deliver a stringent test of the Kerr hypothesis, which predicts a definite shape for its higher-$n$ rings: photons orbiting just outside the horizon of a black hole can probe its extreme gravity and carry away information about its spacetime geometry, encoded in the observable shape of the rings that these photons produce in their observer’s sky [32, 33, 34].

Time variability introduces a third significant complication. While time-averaged GRMHD-simulated movies have shown that the photon rings are persistent, sharp features that come to dominate observations with very-long-baseline interferometry (VLBI) after averaging over sufficiently long timescales [7], it remains to be understood how clearly visible the rings will be to a realistic, near-future, space-VLBI mission like SALTUS, which will be limited in the number of snapshots it can collect. In each snapshot, such an interferometer—with a single space leg—can only sample the radio visibility on one space-ground baseline, thereby only measuring the projected diameter $d_\varphi$ of the photon ring at one angle $\varphi$ in the image.² To compensate for its sparse baseline coverage, the instrument can observe at regular intervals along its orbit around the Earth, eventually filling in every angle $\varphi$ in the Fourier domain, with each $d_\varphi$ thus sampled twice per orbit. The baseline lengths over which the ring signature dominates the signal fix the orbital radius (about lunar distance for the $n = 2$ ring of M87*) and hence the orbital period ($\sim1$ month), which in turn sets the cadence of these snapshots: $\sim40M_{\mathrm{M87}^*}$, or roughly every two weeks. As it is evidently impractical to maintain coherence over such timescales, the snapshots must be incoherently time-averaged; moreover, only $\sim24$ snapshots of $d_\varphi$ can be sampled per year.

In a single snapshot, the “clean” interferometric signature of the ring—a periodic ringing in the visibility amplitude—is typically “polluted” by noise from both the instrument and from astrophysical fluctuations (plasma flares, emission ropes, or other ring mimickers), which can obscure the signal. The key question is then: Can the interferometric signature of a photon ring—and hence its projected diameter—be recovered from an incoherent time-average over $N \approx20$ snapshots of its visibility amplitude on very long space-ground baselines?

An affirmative answer to this question would open the door to a consistency test of the Kerr hypothesis—a cornerstone of general relativity (GR) in the strong-field regime—via space-VLBI measurements of the photon ring shape. The paper [34] (henceforth: GLM) took the first steps towards establishing the viability of such a test for M87*. GLM examined a range of models of stationary, axisymmetric, equatorial disks that reproduce the time-averaged observational appearance of GRMHD-simulated flows, and found that the observable shape of their $n = 2$ ring always follows a specific functional form, independent of the source model. They concluded that this ring shape is a robust prediction of strong-field GR.

In other words, observations of the $n = 2$ photon ring can in principle disentangle gravitational and astrophysical effects that are otherwise commingled in the direct image. Moreover, GLM simulated interferometric data of the kind that could be collected by a mission like SALTUS, and were able to extract this ring shape from the visibility amplitude on space-ground baselines. Their experimental forecast achieved a sub-percent level of precision for the resulting test of the Kerr hypothesis, suggesting that M87* holds exceptional promise as a target for such a test in practice. This analysis was recently reviewed in depth and extended to an even larger selection of equatorial disk models, supporting these conclusions [10].

While these early results are encouraging, demonstrating the feasibility of a ring shape measurement requires further theoretical work. Crucially, even though the GLM analysis did include realistic instrument noise, it only considered time-averaged images of equatorial disks. The latter limitation was recently tackled with a study of geometric thick-disk models [11], but to date a detailed investigation of source fluctuations and time variability has yet to be carried out.

The present work is the first attempt to remedy this lacuna. The main obstruction is technical: as high-order photon rings are exponentially narrow compared to the overall structure of a black hole image, resolving them in the image plane incurs a large computational cost. More precisely, their presence in the image introduces a large separation of scales between the pixel grid size (which must be large enough to capture the entire field of view) and the pixel spacing (the grid must achieve a sufficiently fine resolution to see the narrow rings). While a brute-force approach—pumping millions of pixels in the grid—can overcome this scale separation for a handful of images, it becomes intractable when dealing with a black hole movie consisting of several hundred snapshots, in which case adaptive ray tracing is necessary [35, 36] (and less wasteful).

Here, we present a numerical tool³ designed to efficiently compute high-resolution “slow-light” movies of generic (non-stationary and non-axisymmetric) equatorial sources around a Kerr black hole, together with their associated radio visibility on very long baselines; we present example outputs in Figure 1.

The code was developed with the intent to maximize speed while still guaranteeing the accuracy of its output, particularly the radio visibility on very long baselines, which encodes the high-frequency components of a snapshot’s Fourier transform and is therefore extremely sensitive to its most minute image features. The code’s structure is highly modular and most of its individual components reproduce pre-existing capabilities; its main novelty is arguably to combine all these routines into one single and convenient-to-use (we hope!) package.

One new and important technique from which AART derives much of its power deserves special mention: it turns the very feature of photon rings that makes them so difficult to fully resolve—namely, their thinness—to its advantage. It does so by decomposing the full image into multiple layers labeled by half-orbit number $n$, and ray tracing in each layer the image of the $n^{\mathrm{th}}$ photon ring with exponentially high (in $n$) resolution.

More precisely, the lensing behavior of a black hole forces the $n^{\mathrm{th}}$ photon ring to lie within an exponentially small (in $n$) region of the image plane, dubbed the $n^{\mathrm{th}}$ lensing band, with each band completely fixed by the Kerr geometry [10, 34]. A key innovation of AART is to first compute (once and for all) the lensing bands associated with a given black hole spin and inclination, and then to ray trace, for the $n^{\mathrm{th}}$ layer, only pixels lying in the $n^{\mathrm{th}}$ lensing band. Since each lensing band contains an exponentially demagnified image of the main emission, by also increasing the resolution in each band exponentially,⁴ the code is guaranteed to resolve the source with roughly the same effective resolution in every layer. In other words, AART adapts its ray-tracing resolution (grid spacing) to each layer, but also adjusts the ray-tracing region (grid size), so as to resolve the increasingly fine features present in higher layers using only a fixed number of pixels per layer.

Other ray tracers use adaptive mesh refinement to increase pixel density in regions where they detect long geodesic path lengths [35] or large gradients [36, 37]. This results in a new, non-uniform grid for every new snapshot, which is recursively refined until a desired criterion is met, or else the number of recursions exceeds a hard-set limit. Therefore, small image features can sometimes be missed when this limit is hit.

By contrast, the adaptiveness provided by the lensing bands is determined by the Kerr geometry alone and results in the same, uniform grid for every image of a given black hole spin and inclination. These grids therefore need be computed only once, and provided that their resolution increases at the proper rate set by the demagnification factor, they cannot miss any feature that is already resolved in the direct $n=0$ image.

Finally, high-$n$ layers comprise photons that execute many orbits in the Kerr photon shell [38] where their radial potential almost develops a double root and geodesic integrals diverge logarithmically, leading to a growing risk of numerical error. To minimize error, AART performs analytical ray tracing using an exact solution of the Kerr null geodesic equation recently given in terms of Legendre elliptic integrals and Jacobi elliptic functions [39], similar in spirit to previous implementations based on Carlson symmetric forms [40, 41].

To summarize: AART uses lensing bands to prevent a large separation of scales (between the grid size and its spacing) and thereby maintain its speed as it increases the pixel density in higher-$n$ photon rings, which is necessary to ray trace the fine image features that (due to the nonlocal character of the Fourier transform) influence the visibility on long baselines; moreover, the rings are ray traced analytically to avoid errors.

These features of AART are specially tailored to the photon ring and its interferometric signature; with this tool in hand, it becomes feasible to investigate the effects of time variability on measurements of the $n=2$ ring shape, and we can begin to answer the key experimental question posed above.

To study the effects of source fluctuations on the observed visibility amplitude, we call upon both AART and inoisy [42]: a code that can rapidly generate realizations of a 2D Gaussian random field with Matérn covariance. Such a field can provide a simple statistical model for a generic (non-stationary and non-axisymmetric) stochastic source in the equatorial plane of a Kerr black hole, with prescribed two-point function. Since a “realistic” choice of autocorrelation structure is not yet known (and will require additional research into the plasma physics of the fluid), our goal will be to vary the statistics of the model within a wide range of “reasonable” possibilities (informed by GRMHD simulations), so as to parameterize our uncertainty in the expected variability of signals from sources like M87*.

Completing such a parameter survey is a large undertaking beyond the scope of this first paper. Here, we will be content with a proof-of-concept demonstration that AART is up to the task for such a study. To showcase its capabilities, we ray trace 100 snapshots of an inoisy source with statistics set by the “best-guess” parameters for M87* (listed in Table II). Sample snapshots and their associated visibility amplitude are shown in Figure 1, while Figure 2 includes all of the snapshots, with the left panel displaying their time-averaged image and the right panel all the individual visibility amplitudes together with their incoherent time average (solid lines). As expected, we find that the astrophysical fluctuations wash out from the time average, leaving an image that is visibly dominated by a prominent photon ring with clear $n = 1$ and $n = 2$ subrings.

Parameter	Default value	Description
$\psi$	1.07	BH mass-to-distance ratio
$a/M$	94%	BH spin
$\theta_{\mathrm{o}}$	$17^\circ$	Observer inclination
$\zeta$	1.5	Geometrical factor
$\sigma$	0.4	Fluctuation scale
$\xi$	0.95	Sub-Keplerian factor
$\beta_{r}$	0.95	Radial velocity factor
$\beta_{\phi}$	0.95	Angular velocity factor
$\theta_{\ell}$	$20^\circ$	Anisotropy direction
$\lambda_{0}$	$2\pi/\Omega$	Temporal correlation
$\lambda_{1}$	$5r_{\mathrm{s}}$	Spatial correlation in the $x_{\mathrm{s}}$ direction
$\lambda_{2}$	$0.1\lambda_{1}$	Spatial correlation in the $y_{\mathrm{s}}$ direction
$\mu$	$r_{-}$	Controls location of profile peak
$\vartheta$	$1/2$	Controls profile width
$\gamma$	$-3/2$	Controls profile asymmetry

Correspondingly, the incoherently time-averaged visibility amplitude is dominated on long baselines by the perfectly clean interferometric signature of the $n = 2$ ring. In particular, we find that the signal in the range $u \in[425, 440]$ G$\lambda$—which a satellite at lunar distance from the Earth could access with 345 GHz observations—exactly follows the periodic ringing pattern predicted for a thin ring (black overlay in panel inset).

As we show in Figure 3, the projected diameter $d_{\varphi}$ of the $n = 2$ ring can then be extracted from the frequency of this ringing in the visibility amplitude $|V(u,\varphi)|$. The top panels display the ring diameter $d_{\varphi}$ inferred from an average over $N$ snapshots, with $N = 5$, 10, 20, and finally 100, by which time the smooth shape of the ring has emerged. The bottom panels display the relative deviation from this fiducial shape that is induced by astrophysical fluctuations, whose noise clearly averages out of the image; it is encouraging to see this noise is also beat down in measurements of the ring shape using only a few snapshots.

Of course, whether these conclusions are likely to hold for M87* has yet to be established, and a systematic investigation of astrophysical fluctuations remains to be done. In a soon-to-be-released paper, we will initiate such a study by repeating this analysis for multiple models of M87*, varying both the parameters of the black hole (its spin and inclination) and of the source (the inoisy model). We also hope to report on the signature of the $n = 1$ ring and how it may encode the spin.

We now turn to a brief summary of the rest of the paper.

Summary

In Sec. II, we review the problem of light propagation in the Kerr spacetime. We write down the null geodesic equation and its exact analytical solution as it is implemented in AART, and describe the key concept of lensing bands (Figure 4). We then illustrate the lensing behavior of the black hole by plotting its “transfer functions”: mappings of directions in the observer sky to the spacetime points where the corresponding light rays intersect the equatorial plane. The transfer functions for polar coordinates $(r_s,\phi_s)$ in the plane are shown in Figs. 5, 6, and 9 for the $n=0$, $n=1$, and $n=2$ images, respectively. Likewise, the time lapse $\Delta t$ between source and observer is also shown in Figs. 7 and 8 for the $n=0$ and $n=1$ images, respectively. The $n^{\mathrm{th}}$ image always fills out the $n^{\mathrm{th}}$ lensing band, as expected.

We also give an approximate formula for light bending that was derived by Beloborodov [43] for nonrotating black holes in the weak-deflection regime. We express his result in a very simple form [Eq. (43)] that proves to be remarkably accurate for the computation of $n=0$ images of axisymmetric sources, for most black hole spins and inclinations. This observation, which is illustrated in Figure 10, highlights the fact that $n=0$ photons barely carry an imprint of the black hole spin, as they do not spend enough time near it to be strongly affected by its gravitational field. Beloborodov’s approximation can thus be regarded as a generalization to all inclinations of the “just add one” prescription for the transfer function of a spin-aligned observer, for whom the impact parameter $\rho$ is simply related to emission radius $r_s$ by adding one: $\rho\approx r_s + M$ [5, 44]. AART can also parallel transport linear polarization. We check that Beloborodov’s approximation is adequate for ray tracing $n=0$ polarimetric images (Figure 11), which can be done using simple algebraic equations that we write down explicitly.

In Sec. III, we describe our model of equatorial emission. First, we review stationary and axisymmetric equatorial disk profiles that can reproduce the time-averaged observational appearance of M87* in GRMHD simulations [9, 34]. Then, we add in astrophysical fluctuations with prescribed statistics.

In Sec. IV, we use noisy to simulate a variable source. With AART, we ray trace its instantaneous snapshots (Figure 12) and compute its light curve (Figure 13). Great care must be taken in the choice of resolution and field of view used in each layer: Sec. V discusses these issues in depth (Figure 15, Figure 16, and Figure 14). We can then produce movies of a stochastic, non-stationary, non-axisymmetric source. In Sec. VI, we compute movies of the associated visibility amplitude and use this synthetic data to reconstruct the projected diameter $d_{\varphi}$ of the $n=2$ ring, before concluding with a brief discussion of future prospects in Sec. VII. We relegate some details to Apps. A, B, and C. Throughout the paper, we work in $(-,+,+,+)$ metric signature with geometric units in which $G=c=1$. Our conventions for Legendre elliptic integrals are listed in App. A of Ref. [45].

II Theoretical framework

The special integrability properties of the Kerr spacetime reduce its geodesic equation to a problem of quadratures [46], resulting in elliptic integrals that are expressible in Legendre normal form [39, 47, 48]. Modern computers can evaluate the Legendre elliptic integrals very fast and to arbitrary precision, reducing the computational cost of ray tracing in Kerr. In this section, we review the exact solution of the Kerr null geodesic equation in Legendre form [5, 39] as well as the approximate solution derived in Schwarzschild by Beloborodov [43], and we use them to plot the transfer functions mapping Bardeen’s coordinates in the observer sky [14] to the equatorial plane.

A Null geodesics of the Kerr exterior