© H. Paugnat et al. 2022

1. Introduction

The existence of black holes (BH) is a key prediction of general relativity (GR) in the strong-field regime. More specifically, the theory posits that the spacetime geometry around these compact objects is described by the Kerr metric. This “Kerr hypothesis” underlies a considerable amount of astrophysics and also plays a driving role in theoretical physics, but it has yet to be directly tested. A number of experiments are now poised to change this. The recent successes obtained in gravitational-wave and radio interferometry by LIGO-Virgo, GRAVITY, and EHT (LIGO Scientific Collaboration & Virgo Collaboration 2016; GRAVITY Collaboration 2018; Event Horizon Telescope Collaboration 2019a) herald the advent of a new era of black hole astronomy, and their future extensions promise to deliver observations that will enable tests of the Kerr hypothesis with unprecedented precision.

In April 2019, the EHT collaboration published the very first images resolving the nucleus of the galaxy M87 (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f). Reconstructed from 1.3 mm very-long-baseline interferometry (VLBI) observations, these images achieved an angular resolution comparable to the expected size of M87* – the supermassive compact object at the center of M87 – and revealed a thick, non-uniformly bright ring with a diameter of ∼40 μas encircling a central brightness deficit. These observed features are roughly compatible with those appearing in simulated black hole images (e.g., Luminet 1979; James et al. 2015; Event Horizon Telescope Collaboration 2019e), and indeed M87* may be consistently modeled as a Kerr BH surrounded by an accretion flow whose emission originates within a few Schwarzschild radii of the event horizon. Direct (weakly lensed) images of such a flow present an asymmetric bright ring that is consistent with the 2017 EHT observations (e.g., Gralla et al. 2019; Johnson et al. 2020; Chael et al. 2021).

Due to the limited resolution of these first EHT images of M87*, this consistency test of the Kerr hypothesis remains weak; instead, the Event Horizon Telescope Collaboration (2019f) has thus far assumed the Kerr nature of M87*, so as to constrain its mass from observations. However, as technology continues to improve, increasingly better images will be obtained, and these should enable correspondingly sharper probes of the spacetime geometry around M87*. Nonetheless, it is not immediately clear how to best extract relativistic signatures from such probes, which raises the question: how could sharper images of M87* in principle lead to more stringent tests of the Kerr hypothesis?

Addressing this question is made particularly difficult by the inherent asymmetry of the imaging problem. Thanks to curved-spacetime ray tracing codes such as Gyoto (Vincent et al. 2011) and others (Gold et al. 2020), solving the forward problem has become a fairly routine task, and determining the observational appearance of a compact object with specified properties (spacetime metric, emission model) is now relatively straightforward. The inverse problem, on the other hand, remains nigh intractable: inferring the nature of an astrophysical source from its visual appearance requires one to disentangle effects of the spacetime geometry from properties of the emission, a problem subject to enormous degeneracy. A veritable zoo of theoretically proposed objects (besides the Kerr BH) have been shown to also produce simulated images with features compatible with the recent EHT observations¹, and many of these alternative models cannot yet be ruled out by observations with current EHT resolution (e.g., Vincent et al. 2021).

This present blurriness is not the only limitation: even with higher resolution, images that resolve only the direct emission will likely never separate astrophysical and geometrical effects. (This partly explains the present lack of robust constraints on the spin of M87*.) As such, these images may always be compatible with models of a Kerr BH surrounded by some (possibly exotic) emission, and hence they could never lead one to definitively conclude a violation of the Kerr hypothesis, as this would require achieving greater confidence in the astrophysics than in GR – this could of course break the degeneracy in principle, but seems very unlikely in practice (e.g., Gralla 2020; Bauer et al. 2022).

Recent work suggests a promising path forward. The strong gravity of a Kerr BH creates around it a “photon shell” – a region outside its event horizon in which its gravity is so strong that light rays can become trapped on unstably bound orbits (Bardeen 1973; Teo 2003). As a result, sources in the vicinity of the black hole can produce multiple relativistic images in the observer sky: provided that the emission is optically thin, the primary image (consisting of weakly lensed photons that travel directly to the observer after being emitted) is generically superimposed with a series of mirror images arising from highly bent photons that circumnavigated the photon shell multiple times on their way from source to observer (Johnson et al. 2020; Gralla & Lupsasca 2020a). If the emission is not spherical (so that it does not fully surround the hole), then these successive images form a discrete sequence of “photon rings” that are labeled by photon half-orbit number n and converge to a “critical curve” in the observer sky corresponding to the light rays that asymptote to unstably bound orbits in the photon shell. Crucially, GR predicts that these rings must display a self-similar structure governed by Kerr critical exponents γ, δ, and τ that respectively control the successive demagnification, rotation, and time delay of successive images (Johnson et al. 2020; Gralla & Lupsasca 2020a). This suggests that observing this characteristic lensing pattern may provide a way to extract information about the black hole geometry, and to discriminate between different spacetimes (e.g., Wielgus 2021).

The photon ring has not yet been conclusively observed, but its GR-predicted substructure has already been numerically confirmed with state-of-the-art (general relativistic magneto-hydrodynamics) GRMHD-simulated images of M87*, in which at least n = 3 subrings can be distinguished (Johnson et al. 2020). Although these rings (which are observable) are distinct from the theoretical critical curve (which is not observable), they closely track its shape, which is insensitive to the astrophysical details of the surrounding emission and entirely fixed by the Kerr metric (Bardeen 1973; Gralla & Lupsasca 2020c). Critical curves have been extensively discussed in the literature and spacetimes other than Kerr are known to produce different shapes (e.g., Johannsen & Psaltis 2010; Johannsen 2013; Cunha & Herdeiro 2018; Okounkova et al. 2019; Medeiros et al. 2020). Thus, measuring the shape of a high-n photon ring of M87* may be a good proxy for probing the shape of its critical curve, and hence its geometry.

This idea was recently made precise by Gralla et al. (2020), who established that, in a certain class of models, the shape of the n = 2 subring is extremely close to that of the critical curve. Moreover, they showed that on very long baselines u ∼ 300 Gλ, which could be accessible with space-based² extensions of the EHT, this subring presents a clean interferometric signature, as expected from Fourier theory (Johnson et al. 2020; Gralla 2020). They argued that a space mission targeting M87* could resolve the shape of its n = 2 photon ring with sufficient precision to enable a sharp comparison with the GR prediction: simulations of their experiment achieved sub-percent precision, forecasting a stringent consistency test of the Kerr hypothesis. We will refer to their proposed procedure as the GLM method.

Of course, this first study is subject to multiple limitations, as it relies on a restricted class of emission models consisting of stationary, axisymmetric disks composed of circular-equatorial orbiters. More precisely: 1) it only investigated a choice of 3 toy emission profiles, viewed from observer inclinations of 10° − 30° (the range of expected relevance to M87*); 2) it considered only geometrically thin disks; 3) it ignored astrophysical fluctuations, which amounts to studying (coherently) time-averaged images.

It is crucial to assess the robustness of the GLM method in greater generality. The aim of this paper is to extend the GLM analysis to many more configurations; that is, to broaden its scope by performing a parameter survey of emission profiles, as well as BH spins and inclinations, in order to remove the first limitation outlined above. This work remains confined to the framework of thin, equatorial accretion disks; a parallel study of the impact of disk thickness is ongoing and will soon address the second limitation (Vincent et al. 2022). In the meantime, it is worth noting that recent work by Chael et al. (2021) suggests that geometrically thin models yield effective approximations to numerically modeled GRMHD accretion flows.

Finally, another objective of this paper is to initiate the study of black-hole-spin estimates based on astrophysics-independent shape measurements of the n = 2 photon ring. A first step is taken in this direction by establishing that photon rings arising from astrophysically plausible configurations constitute only a small subset of the full range of theoretically allowed shapes.

The paper is organized as follows. First, Sect. 2 reviews the theoretical foundations of the GLM method: the critical curve is defined and contrasted with the concept of “black hole shadow,” the notion of “lensing bands” is introduced, and the expected interferometric signature of the photon ring is described. Next, Sect. 3 presents a numerical implementation of the GLM method that improves upon the original technique. Then, Sect. 4 describes how the choice of baseline may affect a shape measurement of the photon ring, especially if its thickness varies substantially with angle. After that, Sect. 5 presents the results of a survey over emission profiles while Sect. 6 examines the influence of BH spin and inclination. We find that a measurement is always possible at low inclinations, while new complications that may arise at high inclination necessitate further refinements to the method. We provide further discussion in Sect. 7, before ending with our conclusions and perspectives in Sect. 8.

2. Theoretical overview

2.1. Shadow, critical curve, and photon rings

In the introduction, we asserted that the photon ring generically decomposes into a stack of discrete subrings. For completeness, here we list the (mild) astrophysical assumptions required for the presence of this substructure; this also provides us with an opportunity to distinguish the photon ring from the “black hole shadow” and critical curve, terms deserving of disambiguation.

First, the emission region must be optically thin, so that light may cross it multiple times before being reabsorbed; otherwise, highly lensed (n > 0) images may be obscured (e.g., Beckwith & Done 2005). The validity of this assumption depends on both the source and the observing frequency. As indicated in Fig. 3 of Johnson et al. (2020), GRMHD predicts it to likely hold for 230 GHz photons collected from M87* by present-day EHT, and it is even more likely true at the higher 345 GHz frequency that ngEHT will target; by contrast, the core of Centaurus A is only expected to be optically thin past 1–5 THz (Janssen et al. 2021).

Second, the emission region must have a gap rather than fully surround the hole; otherwise, the discrete structure visible in the intensity cross-sections displayed in Johnson et al. (2020), Chael et al. (2021), which resemble the layers of a wedding cake, will not be present. Instead, one will see a continuous cross-section with a mild (logarithmic) divergence near the the critical curve: the image-plane curve, first derived by Bardeen (1973; who used the term “apparent boundary”), corresponding to light rays that asymptote to unstably bound spherical photon orbits in the photon shell. The resulting image will then display a distinctive feature known as the “black hole shadow”: a central brightness deficit precisely bounded by the critical curve. This effect was first highlighted by Falcke et al. (2000) in the context of a BH surrounded by a spherical, radially infalling accretion flow, and later revisited by Johannsen & Psaltis (2010).

However, as recently pointed out by Gralla et al. (2019) and Narayan et al. (2019), this spherical-infall scenario is highly fine-tuned, and in fact no longer viewed as relevant for M87*: the latest EHT constraints derived from the polarimetric image of M87* (Event Horizon Telescope Collaboration 2021a,b) strongly favor “magnetically arrested disk” (MAD) models (Narayan et al. 2003; Igumenshchev et al. 2003) with emission concentrated near the midplane. Such models do not display the traditional “shadow” feature (i.e., a central brightness depression filling the critical curve), but rather a different “inner shadow” feature associated with the inner edge of the emission, which produces a smaller but even darker central brightness depression contained well within the critical curve (Chael et al. 2021). For equatorial disks extending to the horizon, the inner edge of the emission in principle coincides with the lensed equatorial event horizon (though it may appear larger due to redshift effects).

To summarize, the boundary of the black hole shadow is a mathematical curve which is not in itself observable, and only happens to coincide with an image feature – the central brightness depression – in special configurations currently disfavored by EHT data on M87*. On the other hand, the photon ring is a visible feature that always arises in images of a BH surrounded by an optically thin emission region. Moreover, if this emission is gapped, then this photon ring decomposes into subrings that converge (exponentially in n) to the theoretical critical curve, which can be thought of as the n → ∞ subring. The most easily accessible n = 1 and n = 2 subrings are close to, but nonetheless still distinct from, the exact critical curve. Therefore their shape is constrained to be similar, but not exactly equal, to that of the shadow. This observation is key for the GLM method.

2.2. Bound geodesics form the photon shell

Before turning to the study of the observable photon ring shape, we first elucidate the shape of the theoretical image-plane curve comprising the null geodesics that asymptote to unstably bound orbits around a Kerr BH (Bardeen 1973). The region of spacetime containing these orbits – the photon shell – is extensively reviewed in Teo (2003), Johnson et al. (2020). Here we describe it with Boyer-Lindquist coordinates (t, r, θ, ϕ) on the spacetime of a Kerr BH with mass M and angular momentum J = a_*M:

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2& =\frac{\mathrm{\Delta }}{\mathrm{\Sigma }}{(\mathrm{d}t-a_{\ast }{sin}^2\theta \mathrm{d}\varphi )}^2+\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}\mathrm{d}{r}^2+\mathrm{\Sigma }\mathrm{d}{\theta }^2\hfill \\ \hfill & +\frac{{sin}^2\theta }{\mathrm{\Sigma }}{[(r^2+a_{\ast }^2)\mathrm{d}\varphi -a_{\ast }\mathrm{d}t]}^2,\hfill \end{array}$$(1)

with $\mathrm{\Delta }(r)=r^2-2Mr+a_{*}^2$ and $\mathrm{\Sigma }(r,\theta )=r^2+a_{*}^2{cos}^2\theta$. We will often use a dimensionless spin parameter a ≡ a_*/M ∈ [ − 1, 1].

A Kerr photon with four-momentum p^μ has a conserved energy-rescaled angular momentum and Carter constant

$$\begin{array}{c}\hfill \lambda =\frac{p_{\varphi }}{-p_t},\eta =\frac{p_{\theta }^2}{p_t^2}-a_{\ast }^2{cos}^2\theta +{\lambda }^2{cot}^2\theta .\end{array}$$(2)

Bound photon orbits in Kerr are “spherical”: they evolve at a constant Boyer-Lindquist radius in the range ${\stackrel{\sim }{r}}_{-}\le r\le {\stackrel{\sim }{r}}_{+}$, where

$$\begin{array}{c}\hfill {\stackrel{\sim }{r}}_{\pm }=2M[1+cos(\frac{2}{3}arccos(\pm a))]\end{array}$$(3)

is the radius of the retrograde (upper sign) or prograde (lower sign) circular-equatorial orbit. As they rotate around the BH, these bound orbits undergo polar librations between angles $[{\stackrel{\sim }{\theta }}_{-},{\stackrel{\sim }{\theta }}_{+}]$, where

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\theta }}_{\pm }& =arccos(\mp \sqrt{{\stackrel{\sim }{u}}_{+}})\gtrless \frac{\pi }{2},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{u}}_{\pm }& =\frac{r}{a_{\ast }^2{(r-M)}^2}[-r^3+3M^{2}r-2a_{\ast }^{2}M\hfill \\ \hfill & \pm 2\sqrt{M\mathrm{\Delta }(2r^3-3M{r}^2+a_{\ast }^{2}M)}],\hfill \end{array}$$(5)

and their characteristic conserved quantities are given by

$$\begin{array}{c}\hfill \stackrel{\sim }{\lambda }=\frac{M(r^2-a_{\ast }^2)-r\mathrm{\Delta }}{a_{\ast }(r-M)},\stackrel{\sim }{\eta }=-a_{\ast }^2{\stackrel{\sim }{u}}_{+}{\stackrel{\sim }{u}}_{-}.\end{array}$$(6)

Bound orbits are either prograde or retrograde, depending on the sign of $\stackrel{\sim }{\lambda }$. Orbits at generic radii $r\in [{\stackrel{\sim }{r}}_{-},{\stackrel{\sim }{r}}_{+}]$ are not closed, as they do not return to the same location after a full rotation by Δϕ = 2π, but rather explore the entire allowed region within their shell of fixed r; resonant orbits which are closed form an exceptional set of measure zero in the photon shell (Wong 2021). The photon shell occupies a 3D region of space given at any time t by ${\stackrel{\sim }{r}}_{-}\le r\le {\stackrel{\sim }{r}}_{-}$, θ₋ ≤ θ ≤ θ₊, and 0 ≤ ϕ < 2π. In the Schwarzschild case (a = 0), it reduces to the sphere at r = 3M.

2.3. Asymptotically bound geodesics form the critical curve

A photon that is initially outside the photon shell but has exactly the same conserved quantities $(\lambda ,\eta )=(\stackrel{\sim }{\lambda },\stackrel{\sim }{\eta })$ as a photon bound at orbital radius r asymptotes to that orbit in its far future or past. The critical curve is defined as the theoretical curve in the sky of an observer corresponding to these asymptotically bound photons. It delineates the region of photon capture by the BH (inside the curve) from that of photon escape (outside the curve).

For the case of an observer at large distance D ≫ M and inclination i ≠ 0 from the spin axis of a BH, Bardeen (1973) showed that a photon with conserved quantities (λ, η) appears in the sky at a position given by³

$$\begin{array}{cc}\hfill \rho & =\frac{1}{D}\sqrt{\eta +{\lambda }^2+a_{\ast }^2{cos}^{2}i},cos{\phi }_{\rho }=-\frac{\lambda }{D\rho sini},\hfill \end{array}$$(7)

where, following Johnson et al. (2020), we replaced Bardeen’s original Cartesian coordinates (α, β) with (dimensionless) polar coordinates (ρ, φ_ρ) that are centered about the “line of sight to the black hole” and such that φ_ρ = 90° corresponds to the BH spin axis projected onto the plane perpendicular to this line of sight. Hence, for such an observer, the critical curve is the parametric curve $C_{\gamma }=(\stackrel{\sim }{\rho },{\stackrel{\sim }{\phi }}_{\rho })$ obtained by tracing

$$\begin{array}{cc}\hfill \stackrel{\sim }{\rho }& =\frac{1}{D}\sqrt{{\stackrel{\sim }{\lambda }}^2+a_{\ast }^2({cos}^{2}i-{\stackrel{\sim }{u}}_{+}{\stackrel{\sim }{u}}_{-})},\hfill \end{array}$$(8a)

$$\begin{array}{cc}\hfill cos{\stackrel{\sim }{\phi }}_{\rho }& =-\frac{\stackrel{\sim }{\lambda }}{D\stackrel{\sim }{\rho }sini},\hfill \end{array}$$(8b)

over the radial extent $r\in [{\stackrel{\sim }{r}}_{-},{\stackrel{\sim }{r}}_{+}]$ of the photon shell. Counter-intuitively, the angle ${\stackrel{\sim }{\phi }}_{\rho }$ is parameterized not by the angle ϕ around the BH, but rather by the radius r within the photon shell. This effect highlights the warped nature of the Kerr spacetime.

Since Eq. (8b) admits two solutions for every choice of r, each orbit in the photon shell corresponds to two angles ${\stackrel{\sim }{\phi }}_{\rho }$ around the critical curve: this corresponds to the fact that an asymptotically bound photon can reach the observer on either an upwards or downwards libration. Only the equatorial observer with i = 90° can see the full set of orbits within the photon shell; other observers see only the subset of orbits that have sufficient inclination to reach them (that is, such that ${\stackrel{\sim }{\rho }}^2\ge 0$). Examples of critical curves for various BH spins and inclinations are plotted in many papers (e.g., Fig. 1 in Farah et al. 2020).

Finally, we emphasize that the notion of a critical curve is not limited to observers far from the BH. For instance, the critical curve for observers on timelike circular-equatorial orbits around a Kerr BH is analyzed in Gates et al. (2021).

2.4. Nearly bound photons form the photon ring

A photon with conserved quantities $(\lambda ,\eta )\approx (\stackrel{\sim }{\lambda },\stackrel{\sim }{\eta })$ very close to those of a photon bound at orbital radius $\stackrel{\sim }{r}$ can describe multiple oscillations in $\theta \in [{\stackrel{\sim }{\theta }}_{-},{\stackrel{\sim }{\theta }}_{+}]$ near the photon shell before leaving it. Since Kerr bound orbits are unstable, their nearby, nearly bound geodesics display an exponential rate of radial deviation with respect to their (fractional) number N of polar half-orbits,

$$\begin{array}{c}\hfill |r-\stackrel{\sim }{r}|\propto e^{\gamma N},\end{array}$$(9)

where one full orbit is defined as a complete oscillation from a turning point θ_± back to itself, while fractional N is precisely defined in Eq. (36) of Gralla & Lupsasca (2020a). The Lyapunov exponents governing the orbital instability depend on both the BH spin and photon shell radius (Johnson et al. 2020):

$$\begin{array}{c}\hfill \gamma =\frac{4}{a_{\ast }}\sqrt{r^2-\frac{Mr\mathrm{\Delta }}{{(r-M)}^2}}{\displaystyle {\int }_0^1}\frac{\mathrm{d}t}{\sqrt{(1-t^2)({\stackrel{\sim }{u}}_{+}{t}^2-{\stackrel{\sim }{u}}_{-})}}.\end{array}$$(10)

By definition, when a nearly bound photon reaches a distant observer, it appears at a position close to the critical curve (Eq. (8)): $(\rho ,{\phi }_{\rho })\approx (\stackrel{\sim }{\rho }(\stackrel{\sim }{r}),{\stackrel{\sim }{\phi }}_{\rho }(\stackrel{\sim }{r}))$ where $\stackrel{\sim }{r}$ is the radius of the nearby bound orbit. More precisely, its arrival position can be labeled by the coordinates $(\stackrel{\sim }{r},d)$, where d denotes the perpendicular distance from the nearest point $(\stackrel{\sim }{\rho }(\stackrel{\sim }{r}),{\stackrel{\sim }{\phi }}_{\rho }(\stackrel{\sim }{r}))$ on the critical curve. Johnson et al. (2020) heuristically showed, and then Gralla & Lupsasca (2020a) rigorously proved, that a light ray shot backwards into the geometry from position $(\stackrel{\sim }{r},d)$ in the sky executes

$$\begin{array}{c}\hfill N(\stackrel{\sim }{r},d)\approx -\frac{1}{\gamma (\stackrel{\sim }{r})}log\left|d\right|\end{array}$$(11)

half-orbits around the BH before either returning to asymptotic infinity, if it was outside the critical curve (d > 0), or crossing the event horizon, if it was inside (d < 0). Gralla & Lupsasca (2020a) also compute subleading 𝒪(d⁰) corrections to $N(\stackrel{\sim }{r},d)$.

While Eq. (11) is a purely geometric statement about the lensing behavior of a Kerr BH, it also constrains the structure of BH images. To describe these constraints, we will from now on restrict our attention to emission regions that are optically thin, stationary and axisymmetric. The assumption of stationarity and axisymmetry amounts to considering only time-averaged images (averaged over sufficiently long timescales), and the assumption of infinite optical depth allows one to neglect absorption effects; as discussed in Sect. 2.1, this approximation becomes exact in the limit of large observation frequency, and could already provide a good approximation at the frequencies of current or near-future observations of M87*.

Since the emission region around the BH is optically thin, a light ray that passes through it 2N times (each time at essentially the same inclination) collects roughly twice as many photons as a neighboring light ray that passes through it N times. That is,

$$\begin{array}{c}\hfill I(\stackrel{\sim }{r},d)\propto N(\stackrel{\sim }{r},d)\sim -\frac{1}{\gamma (\stackrel{\sim }{r})}log\left|d\right|.\end{array}$$(12)

As such, the lensing formula (Eq. (11)) implies a mild (logarithmic) divergence in the observed intensity I(ρ, φ_ρ) near the critical curve. The photon ring is defined as the intensity bump caused by this divergence; in practice, it is cut off by absorption effects (finite optical depth) as N grows large. In summary: GR predicts that embedded within every image of a BH surrounded by an optically thin emission region, there lies a narrow photon ring.

Due to Eq. (11), any axisymmetric ring emitting isotropically produces an infinite sequence of lensed images within the photon ring (Gralla & Lupsasca 2020a). Successive images within this sequence are produced by photons describing increasingly many half-orbits around the black hole: if the first (direct) image arises from photons with fractional half-orbit number 0 ≤ N₀ < 1, then the nth lensed image will arise from photons with fractional half-orbit number N ≈ N₀ + n. Since the entirety of the emission surrounding the black hole can be decomposed into such source rings, it follows that the photon ring consists of a series of lensed images of the main emission labeled by half-orbit number n ∈ ℕ.

Moreover, Eq. (11) implies that these successive images are increasingly lensed towards the critical curve by an exponential demagnification factor: a source ring whose nth image appears at perpendicular distance $d_n(\stackrel{\sim }{r})$ from the critical curve will produce its (n + 1)th image at an exponentially smaller distance

$$\begin{array}{c}\hfill d_{n+1}(\stackrel{\sim }{r})=e^{-\gamma (\stackrel{\sim }{r})}{d}_n(\stackrel{\sim }{r}).\end{array}$$(13)

Depending on whether the emission region is gapped or not, this lensing pattern can result in one of two very distinctive image substructures within the photon ring.

If the BH is entirely immersed in its surrounding emission (as in the case of a spherically symmetric accretion flow onto the hole, for instance), then successive images of the emission region continuously blend together and the resulting intensity profile logarithmically rises to a smooth peak centered about the critical curve. Such an intensity profile is displayed in Fig. 1 of Narayan et al. (2019), with the corresponding image shown in Fig. 2 therein: a bright photon ring encircling a central brightness deficit (the interior of the critical curve). Moreover, as shown in that paper, if the emitting matter is infalling, then beaming effects strongly attenuate the intensity inside the critical curve, creating the distinctive shadow effect reviewed in Sect. 2.1 above. The reason is that light rays appearing inside the critical curve cannot encounter a radial turning point, so any photons collected along such rays had to have been emitted opposite to the infalling matter’s direction of motion, thereby incurring a strong redshift.

On the other hand, if the emission is gapped (as in the case of emission localized around the equator with gaps near the poles, for instance), then each successive image of the emission region displays a corresponding gap, and the resulting intensity profile no longer consists of a smooth logarithmic peak (Eq. (12)), but rather a sequence of distinct peaks accumulating near the critical curve. Such an intensity profile is displayed in Fig. 1 of Gralla et al. (2019), with the corresponding image shown in Fig. 5 therein. In such images, the photon ring decomposes into a stack of distinct photon subrings, each of which is labeled by n and corresponds to a single lensed image of the main emission, demagnified by e^−γn. As reviewed in Sect. 2.1, the latest EHT observations favor M87* models with near-equatorial emission, and thus we expect this subring structure to be present in actual images of M87*.

High-resolution, time-averaged, 230 GHz images of M87* produced using state-of-the-art GRMHD simulations display at least n = 3 subrings and numerically confirm this GR-predicted substructure (Johnson et al. 2020), which is also present in all the images that we obtained from our toy models of M87*.

2.5. Photon subrings lie in lensing bands

To study the subring structure of the photon ring, it is helpful to introduce the (purely geometric) notion of a lensing band (Gralla et al. 2020; Chael et al. 2021): a Kerr observer’s nth lensing band is defined as the region in the sky corresponding to geodesics that cross the equatorial plane at least n + 1 times before escaping to asymptotic infinity (if d > 0) or crossing the event horizon r = r₊ of the BH (if d < 0), where the inner and outer horizon radii are

$$\begin{array}{c}\hfill r_{\pm }=M\pm \sqrt{M^2-a^2}.\end{array}$$(14)

We present a general method for computing Kerr lensing bands in Appendix A. Figure 1 illustrates the n = 1 and n = 2 lensing bands for an observer at inclination i = 45° from a BH of spin a = 0.99.

Each n > 0 lensing band consists of a region bounded by two concentric, closed curves, and therefore takes an annular shape. The inner and outer edges of the nth lensing band correspond to the nth lensed images of the equatorial circles of radius r = r₊ and r → ∞, respectively. The inner edge is always contained within the critical curve, while the outer edge always encircles it on the outside. According to Eq. (13), these edges converge to the critical curve exponentially fast in n, with widths scaling as

$$\begin{array}{c}\hfill w_{n+1}(\stackrel{\sim }{r})\approx e^{-\gamma (\stackrel{\sim }{r})}{w}_n(\stackrel{\sim }{r})\approx e^{-n\gamma (\stackrel{\sim }{r})}{w}_1(\stackrel{\sim }{r}),\end{array}$$(15)

where ${\mathit{w}}_1(\stackrel{\sim }{r})$ is the angle-dependent width of the first (n = 1) lensing band (here parameterized by photon shell radius) and the first relation becomes exact in the limit of large n → ∞. In short, the lensing bands form a stacked sequence of annular shapes, each straddling the critical curve, with the (n + 1)th lensing band strictly contained within the nth one.

From now on, we further restrict our attention to a thin disk of isotropic emitters filling the equatorial plane. In this setting, the lensing bands are particularly useful because they coincide with the photon subrings: the nth lensed image of the disk, which produces the nth photon subring, exactly fills the nth lensing band. It follows that the photon subrings also display a stacked structure with widths scaling as in Eq. (15). Hence, measuring the angle-dependent ratio ${\mathit{w}}_{n+1}(\stackrel{\sim }{r})/{\mathit{w}}_n(\stackrel{\sim }{r})\approx e^{-\gamma (\stackrel{\sim }{r})}$ of successive subrings could (in this setting) allow for the Lyapunov exponents $\gamma (\stackrel{\sim }{r})$ to be experimentally determined (Johnson et al. 2020).

2.6. Interferometric signature of a narrow curve in the sky

Having described the subring structure of BH images, we next turn to a description of the corresponding subring structure in the visibility (Fourier) domain. We begin by reviewing some facts about the interferometric signature of a narrow feature in the sky.

We consider an infinitely thin, closed light curve in the sky. Assuming this plane curve is convex, we may parameterize it by normal angle φ ∈ [0, 2π) as r(φ) = (α(φ),β(φ)), where we temporarily revert to Cartesian coordinates (α, β). (Non-convex curves decompose into multiple such segments r_i(φ), but we will not need this generalization here.) The inward-pointing unit normal is $\widehat{\mathit{n}}(\phi )=-(cos\phi ,sin\phi )$. Following Gralla & Lupsasca (2020c), we define the “projected position” f(φ) by

$$\begin{array}{c}\hfill f(\phi )\equiv -\mathit{r}(\phi )·\widehat{\mathit{n}}(\phi )=\alpha (\phi )cos\phi +\beta (\phi )sin\phi .\end{array}$$(16)

This function encodes the whole curve r(φ) = (α(φ),β(φ)), which may be reconstructed as

$$\begin{array}{cc}& \alpha (\phi )=f(\phi )cos\phi -f^{\prime }(\phi )sin\phi ,\hfill \end{array}$$(17a)

$$\begin{array}{cc}& \beta (\phi )=f(\phi )sin\phi +f^{\prime }(\phi )cos\phi .\hfill \end{array}$$(17b)

The projected position f(φ) thus encodes all other properties of the curve, such as its radius of curvature ℛ(φ) = f(φ)+f″(φ). The real motivation for introducing the quantity f(φ) is that it naturally connects to the interferometric signature of the curve. To see this, we first note that the projected position can be uniquely decomposed into parity-even and parity-odd parts

$$\begin{array}{c}\hfill f(\phi )=\frac{d_{\phi }}{2}+C_{\phi },\end{array}$$(18)

where

$$\begin{array}{cc}& d_{\phi }=f(\phi )+f(\phi +\pi ),\hfill \end{array}$$(19a)

$$\begin{array}{cc}& C_{\phi }=\frac{1}{2}[f(\phi )-f(\phi +\pi )].\hfill \end{array}$$(19b)

The even part d_φ = d_φ + π is the “projected diameter” of the curve at angle φ, while the odd part C_φ = −C_φ + π is its centroid motion. Gralla & Lupsasca (2020c) illustrate several examples of this shape decomposition.

We now arrive at the key point: if the curve is not infinitely thin, but still very narrow relative to its diameter (w ≪ d_φ), then its Fourier transform V(u, φ) (where (u, φ) are polar coordinates in the Fourier plane) exhibits a simple behavior in the regime

$$\begin{array}{c}\hfill \frac{1}{d}\ll u\ll \frac{1}{w},\end{array}$$(20)

where it adopts the “universal” form, valid to leading order in an expansion in 0 < w/d ≪ 1 (Gralla 2020),

$$\begin{array}{c}\hfill V(u,\phi )\approx \frac{e^{-2\pi i{C}_{\phi }u}}{\sqrt{u}}[{\alpha }_L(\phi )e^{-\frac{i\pi }{4}}{e}^{i\pi d_{\phi }u}+{\alpha }_R(\phi )e^{\frac{i\pi }{4}}{e}^{-i\pi d_{\phi }u}],\end{array}$$(21)

with α_L, R(φ) = α_R, L(φ + π) > 0 encoding the intensity profile of the curve and ensuring that V(u, φ + π) = V^*(u, φ), as required by definition (Eq. (29) below). In the context of interferometry, the complex function V(u, φ) is known as the radio visibility, with the visibility plane (u, φ) usually referred to as the baseline plane. Intuitively, short baselines u ≪ 1/d are not sufficient to resolve the shape of a bright curve, which then just appears as a blob, while on very long baselines u ≫ 1/w, the smooth profile of the curve has been resolved and its signal therefore dies off exponentially; on the other hand, in the universal regime (Eq. (20)), the curve appears infinitely thin and as a result, its signal follows a very weak u^−1/2 power-law fall-off (Eq. (21)).

In the context of realistic high-frequency VLBI observations, it is in practice significantly easier to measure the amplitude |V| of the complex visibility – the “visamp” – than its absolute phase, as the latter is more susceptible to contamination by a variety of atmospheric and instrumental effects. While the full complex visibility (Eq. (21)) encodes both d_φ and C_φ, the visibility amplitude depends on d_φ only, and obeys |V(u, φ + π)| = |V(u, φ)|:

$$\begin{array}{cc}\hfill |V(u,\phi )|& \approx \sqrt{\frac{{\alpha }_L^2(\phi )+{\alpha }_R^2(\phi )+2{\alpha }_L(\phi ){\alpha }_R(\phi )sin\left(2\pi d_{\phi }u\right)}{u}}.\hfill \end{array}$$(22)

To summarize, the clearest interferometric signature of a bright curve in the sky is its visibility amplitude (Eq. (22)), which is strong in the universal regime (Eq. (20)) and only encodes partial information about the shape, namely its projected diameter d_φ (but not the full shape information, which also requires the centroid motion C_φ that is only available in the harder to measure visibility phase).

2.7. Interferometric signature of the photon ring

Combining the results of the previous sections, we are now in a position to describe the interferometric signature produced by the photon ring. We expect the photon subrings to be the only narrow features present in time-averaged BH images, since other fine features (such as emission ropes and other flares) should be transient and wash out after sufficient time-averaging. Thus, we expect that the nth subring will dominate the visibility amplitude in the regime (Johnson et al. 2020)

$$\begin{array}{c}\hfill \frac{1}{w_{n-1}}\ll u\ll \frac{1}{w_n},\end{array}$$(23)

where, according to Eq. (15), ${\mathit{w}}_n(\stackrel{\sim }{r})=e^{-\gamma (\stackrel{\sim }{r})}{\mathit{w}}_{n-1}(\stackrel{\sim }{r})$. Beyond, its signal starts to die off exponentially and the (n + 1)th subring begins to dominate. This results in a distinctive cascade structure illustrated in Fig. 2. In each range of Eq. (23), the nth subring produces a damped oscillation (Eq. (22)) governed by its projected diameter $d_{\phi }^{(n)}$. All that remains for us to do is to describe the diameters $d_{\phi }^{(n)}$, beginning with the projected diameter ${\stackrel{\sim }{d}}_{\phi }\equiv d_{\phi }^{(\infty )}$ of the critical curve, which they exponentially converge to in the limit n → ∞.

2.8. Theoretical shape of the critical curve

A “phoval” (for “photon ring oval”) is a geometric shape with

$$\begin{array}{cc}& \frac{d_{\phi }}{2}=R_0+\sqrt{R_1^2{sin}^2\phi +R_2^2{cos}^2\phi },\hfill \end{array}$$(24a)

$$\begin{array}{cc}& C_{\phi }=(X-\chi )cos\phi +arcsin(\chi cos\phi ).\hfill \end{array}$$(24b)

More precisely, this is a family of shapes parameterized by five parameters R₀, R₁, R₂, X, and χ, with the interesting property that the Kerr critical curve is extremely well-approximated by some member of the family for any choice of observer inclination and BH spin (Gralla & Lupsasca 2020c).

We briefly summarize the geometric interpretation of these parameters. First, the projected position f_circ(φ) = R₀ describes a circle centered at the origin with radius R₀. Next, the projected position $f_{\text{ellipse}}^2(\phi )=R_1^2{sin}^2\phi +R_2^2{cos}^2\phi$ describes an ellipse centered at the origin with radii R₁ and R₂. Adding these two functions results in the projected position f_circlipse(φ) = d_φ/2 given in Eq. (14), which defines a shape known as a “circlipse.” It follows from the last paragraph that the projected diameter of the critical curve is always exquisitely close to that of a circlipse,

$$\begin{array}{c}\hfill {\stackrel{\sim }{d}}_{\phi }\equiv d_{\phi }^{(\infty )}\approx 2R_0+2\sqrt{R_1^2{sin}^2\phi +R_2^2{cos}^2\phi },\end{array}$$(25)

with the best-fit parameters R₀, R₁, and R₂ depending in a rather intricate way on the mass-to-distance ratio M/D, BH spin a, and observer inclination i. Finally, adding the centroid motion (Eq. (24b)) to the projected function produces the full phoval shape, which is translated along the α axis (orthogonal to the spin axis projected onto the plane perpendicular to the line of sight) by a horizontal offset X ∈ ℝ, and carries a left-right asymmetry induced by a parameter χ ∈ [ − 1, 1] that matches the warping of the critical curve at high inclinations and high spins. With these additional parameters, the best-fit phoval reproduces the Kerr critical curve to a part in 10⁵ over the vast majority of parameter space, and to a part in 10³ in the extremal limit a → 1 (Gralla & Lupsasca 2020c). Away from extremality, the (parity-odd) centroid motion gives a minor correction to the shape of the critical curve, which is almost entirely encoded in its (parity-even) circlipse diameter.

A critical curve with best-fit circlipse parameters R₀, R₁, and R₂ has diameters d_⊥ ≡ d₀ and d_∥ ≡ d_π/2 along the directions perpendicular (α axis) and parallel (β axis) to the projected spin

$$\begin{array}{c}\hfill d_{\perp }=2(R_0+R_2),d_{\Vert }=2(R_0+R_1),\end{array}$$(26)

which define a fractional asymmetry (Johnson et al. 2020)

$$\begin{array}{c}\hfill f_A\equiv 1-\frac{d_{\perp }}{d_{\Vert }}=\frac{R_1-R_2}{R_0+R_1}\in [0,1].\end{array}$$(27)

The range 0 ≤ f_A ≤ 1 follows from the observation that d_⊥ ≤ d_∥ over all parameter space, which also implies that R₂ ≤ R₁. Moreover, we can see from Eq. (24a) that d_⊥ and d_∥ are the minimal and maximal diameters of the circlipse, respectively.

Finally, we point out that the two quantities d_∥ and f_A encode the BH spin a and observer inclination i, as illustrated in Fig. 7 of Johnson et al. (2020). The map (a, i)→(d_∥, f_A) is not exactly bijective, however, but rather two-to-one, as (a, i) and ( − a, i + π) both map to the same (d_∥, f_A). Putting aside this discrete twofold degeneracy (which can be broken by other means), this suggests that one could in principle infer a BH’s spin and inclination from the projected diameter ${\stackrel{\sim }{d}}_{\phi }$ of its critical curve (Eq. (25)), or even just its minimum and maximum d_⊥ and d_∥ (without any need for the parameters X and χ encoding its centroid).

However, this suggestion is misleading because the critical curve is purely theoretical and not in itself observable. Nonetheless, the photon subrings – which are observable – converge to it exponentially fast and moreover, while the first subring is still noticeably different to the naked eye, the second subring tracks it quite closely. It is therefore tempting to extract the BH spin and inclination from the projected diameter $d_{\phi }^{(2)}$. Unfortunately, we have found this approach to be only moderately successful, as the residual astrophysics-dependence of the n = 2 ring shape introduces significant uncertainties (see Sect. 5.3 below).

2.9. Testing GR with the photon ring shape

As we have explained, the projected diameter of the Kerr critical curve depends only on the spacetime geometry: ${\stackrel{\sim }{d}}_{\phi }$ encodes the BH spin and inclination, and is fully determined by them. However, this theoretical curve, and hence its diameter ${\stackrel{\sim }{d}}_{\phi }$, are not directly observable. Instead, what we may (in principle) observe is a sequence of photon subrings labeled by half-orbit number n and approaching the critical curve as n → ∞, with projected diameters $d_{\phi }^{(n)}\to d_{\phi }^{(\infty )}={\stackrel{\sim }{d}}_{\phi }$. Thus, in the large-n limit, the rings enter a “universal” regime in which their dependence on astrophysical conditions drops out and they converge to the GR-predicted, astrophysics-independent shape ${\stackrel{\sim }{d}}_{\phi }$. This suggests that a shape measurement of a “large-n” subring could provide a test of GR and a means to infer the BH parameters.

The first (n = 1) subring is certainly not yet in the universal regime, as the photons comprising it are not sufficiently bent and its shape still carries a substantial imprint of the astrophysical details of the surrounding source. On the other hand, the shape of the higher n ≥ 3 subrings is much less sensitive on astrophysics, but these rings are also much harder to detect experimentally, for at least two reasons: first, their signature only comes to dominate on exponentially long baselines (Eq. (23)), which (given present observation frequencies) could only be accessed using a space element at unrealistic separation from the Earth (e.g., Fig. 5 of Johnson et al. 2020); second, this problem is compounded by the need for exceptional sensitivity, as the flux in each subring is also exponentially suppressed in n (Eq. (15)). As such, the second subring seems to be the most promising target for a shape test of GR, as it may sit in the “sweet spot” where it remains accessible with current or near-future technology (Gralla et al. 2020), while n = 2 is “sufficiently large” that $d_{\phi }^{(2)}\approx {\stackrel{\sim }{d}}_{\phi }$.

Gralla et al. (2020) initiated a quantitative investigation of this idea and found that in a wide class of toy models of M87* (subject to the limitations listed in the introduction), the n  =  2 subring indeed produced the universal visibility amplitude (Eq. (22)). Moreover, they showed that it is possible to recover $d_{\phi }^{(2)}$ from this signature, even in the presence of simulated instrument noise. They found that, while the n = 2 ring is still noticeably different from its n → ∞ limit (the critical curve), it nonetheless adopts the same shape: more precisely, even though the n = 2 ring and the critical curve have different projected diameters $d_{\phi }^{(2)}\ne {\stackrel{\sim }{d}}_{\phi }$, they both take the shape of a circlipse. That is, $d_{\phi }^{(2)}$ must also follow the functional form of Eq. (25).

We emphasize that, while the best-fit radii R₀, R₁, and R₂ for ${\stackrel{\sim }{d}}_{\phi }$ are determined only by the BH parameters M/D and a, as well as the observer inclination i, the best-fit radii R₀, R₁, and R₂ for $d_{\phi }^{(2)}$ additionally depend on astrophysical conditions as well. This is illustrated in Fig. 7 of Gralla et al. (2020), which shows how these best-fit parameters vary with the astrophysical model, keeping the BH parameters fixed. Thus, the n = 2 ring carries some astrophysics-dependence (which makes it harder to extract the BH parameters from its diameter), but it is weak enough that $d_{\phi }^{(2)}$ does not break out of the functional form

$$\begin{array}{c}\hfill \frac{d_{\phi }^{(2)}}{2}\approx R_0+\sqrt{R_1^2{sin}^2(\phi -{\phi }_0)+R_2^2{cos}^2(\phi -{\phi }_0)},\end{array}$$(28)

which includes an additional parameter φ₀ to account for the fact that the orientation of the image is generically not aligned with the Bardeen coordinate system. This observation led Gralla et al. (2020) to propose a novel GR test based on the n = 2 ring shape: a space-based VLBI mission targeting M87* could extract $d_{\phi }^{(2)}$ and check to what extent it follows the form of Eq. (28).

If the deviation of $d_{\phi }^{(2)}$ from its best-fit circlipse is large, then the Kerr hypothesis fails the test; if it is small, then one can report a test of the Kerr hypothesis and GR with a precision given by the root-mean-square deviation of $d_{\phi }^{(2)}$ from its best fit. With their simulated experimental forecast, Gralla et al. (2020) attained a stringent, sub-percent (0.04%) level of precision.

Finally, we wish to clarify in what sense we regard the GLM method as a test of the Kerr hypothesis. Such a test can either be:

• A consistency test aiming to confirm that some observable is consistent with the prediction of the Kerr spacetime, in which case there is no notion of model comparison; or,

• A model-comparison test aiming to compare predictions of the Kerr hypothesis with alternative predictions derived from other spacetime geometries, so as to determine which one best explains the data (in the Bayesian sense).

In this article, we discuss the GLM method as a consistency test of the Kerr hypothesis and thus do not consider any other spacetime but that of Kerr. Performing an unambiguous model-comparison test of the Kerr hypothesis is a difficult challenge, as one must be careful to properly take into account all complex astrophysical degeneracies (e.g., Bauer et al. 2022).

Additional work is needed to establish the viability of this GLM method for a wider range of configurations. In this paper, we explore whether the method remains viable when some of the limitations of the original GLM analysis are removed; this will be the object of our parameter survey in Sects. 5 and 6 below.

3. Implementation

To check the robustness and precision of the Kerr hypothesis test reviewed in Sect. 2, we simulated high-resolution images of Kerr BHs using various types of emission and confirmed that they all display a photon ring with n = 1 and n = 2 subrings. Then, we inferred the n = 2 ring diameter $d_{\phi }^{(2)}$ from its characteristic damped oscillation (Eq. (22)) in the visibility amplitude on long baselines, and verified that it follows the predicted functional form (Eq. (28)), thereby establishing the viability of the GLM method as a test of the Kerr hypothesis. We also studied the discrepancy between the n = 2 subring and the critical curve, and checked that a phoval with the minimal and maximal diameters measured from the n = 2 ring could fit within the n = 2 lensing band. Our simulations were performed with the relativistic ray tracing code Gyoto (General relativitY Orbit Tracer of the Observatoire de Paris), whose details are presented in Vincent et al. (2011).

3.1. Image simulation

Gyoto simulates an image, that is, a map I(ρ, φ_ρ) of the specific intensity in the observer sky. To reproduce the complex visibility V(u, φ) that is directly sampled via VLBI observations, we had to additionally compute its 2D Fourier transform,

$$\begin{array}{c}\hfill V(u,\phi )={\displaystyle {\int }_0^{\infty }{\int }_0^{2\pi }I(\rho ,{\phi }_{\rho }+\phi )e^{-2\pi iu\rho cos{\phi }_{\rho }}\rho \mathrm{d}\rho \mathrm{d}{\phi }_{\rho },}\end{array}$$(29)

which by definition satisfies V(u, φ + π) = V^*(u, φ). Instead of computing this 2D FT, we followed Gralla et al. (2020) and made use of the projection-slice theorem to directly compute V(u, φ) along slices of fixed angle φ in the Fourier plane: for each φ, we computed the Radon projection (i.e., the integrals of I(ρ, φ) along lines perpendicular to the slice of constant φ across the image), and then applied a 1D fast Fourier transform (FFT) to obtain the visibility V(u, φ) evaluated at that angle φ.

We made the following assumptions in our computations:

• BH mass: M = 6.2 × 10⁹ M_⊙, in the range of values for M87* favored by stellar dynamical measurements (e.g., Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019f);

• Distance between the BH and observer: D = 16.9 Mpc, in the range of values for M87* (e.g., Event Horizon Telescope Collaboration 2019f);

• Observation wavelength: λ = 1.3 mm, matching the 230 GHz frequency of current observations Event Horizon Telescope Collaboration 2019b;

• We neglected absorption effects, a suitable approximation for optically thin accretion flows, as discussed in Sect. 2 above;

• We used Bardeen’s coordinates (Eq. (7)) in the sky (such that the projected spin axis points in the direction φ_ρ = 90°);

• We terminated the ray tracing of geodesics after their second equatorial crossing, so as to avoid parasitic pixels produced by an under-resolved, partially imaged n = 3 ring.

From these values of M and D, one obtains a conversion ratio from microarcseconds to units of M consistent with EHT priors (Event Horizon Telescope Collaboration 2019f):

$$\begin{array}{c}\hfill {(M/D)}_{\mathrm{M}{87}^{\ast }}=3.6212\mathrm{\mu }\mathrm{as}.\end{array}$$(30)

We took our field of view to be of 180 μas ≈ 50M and ray traced most images with a resolution of 10 000 × 10 000 pixels, though to resolve the n = 2 ring, we found it necessary in some cases (when the ring is highly peaked) to double the linear resolution.

3.2. Image treatment and analysis

3.2.1. Apodization

We find that the visibility amplitude profile |V(u, φ)| obtained from the FFT of a Radon slice indeed presents (on sufficiently long baselines) the expected periodicity linked to the n = 2 ring. Without further treatment, however, this ringing signature does not quite take the “clean” form of Eq. (22) because it is “polluted” by some other oscillation with higher period and lesser amplitude.

This effect arises from an artificial discontinuity in the Radon projection, which is set to zero outside our field of view but does not exactly vanish at its edges, where it instead decays to a tiny (but still nonzero) intensity that is ≲10⁻⁴ times the maximum intensity. Though small, this step function nonetheless produces its own ringing on long baselines, with periodicity set by the size of the field of view, and this interferes with our signal.

To eliminate this spurious artifact, one solution would be to increase the field of view until the intensity at the edges is so small that this effect vanishes. However, in order to still resolve the n = 2 ring, these larger images would have to be ray traced at a correspondingly higher resolution. This would be too costly, especially since the bulk of the computation would be dedicated to exterior pixels of very little interest to us. An ingenious way to overcome this difficulty would be to use adaptative ray tracing, either by employing algorithms designed for this purpose (e.g., Gelles et al. 2021), or else by decomposing the total image into a sum of multiple layers, each with its own resolution. The latter approach was adopted by Gralla et al. (2020), who divided the image into layers labeled by half-orbit number n: since the nth layer only has nonzero pixels in the nth lensing band, which is exponentially small, they were able to exponentially increase the resolution in each layer while keeping fixed the number of pixels computed in each layer.

Instead of computing images with an unnecessarily large field of view or using adaptive ray tracing, we dealt with this effect by simply multiplying the Radon projection by a window function. This method is known as apodization.

More specifically, we used the window function (also known as an apodization or tapering function)

$$\begin{array}{c}\hfill W(x)=p(r(1-x))p(r(1+x)),\end{array}$$(31)

where x ∈ [ − 1, 1] ranges over the Radon projection, ϵ = 1/r is the range of the cutoff, and p is a “smooth plateau” function,

$$\begin{array}{c}\hfill p(x)=\frac{f(x)}{f(x)+f(1-x)},f(x)=\{\begin{array}{cc}{e}^{-1/x^2}\hfill & \mathrm{if}x>0,\hfill \\ 0\hfill & \mathrm{if}x\le 0.\hfill \end{array}\end{array}$$(32)

This function is infinitely smooth (C^∞), exactly equal to unity for x ∈ [ − 1 + ϵ, 1 − ϵ], and exactly vanishes at the edges x = ±1. As a result, multiplying the Radon transform with it smoothens the FFT while retaining the periodicity that we are interested in.

3.2.2. Visibility amplitude profile fit

Once the visibility amplitude of a given image has been obtained, the next step in the GLM test of the Kerr hypothesis is to extract the diameter $d_{\phi }^{(2)}$ of the n = 2 ring at every angle φ_ρ around the image. This diameter can be inferred from the ringing signature displayed by |V(u, φ)| at the corresponding angle φ = φ_ρ in the Fourier plane, which is predicted to take the universal form (Eq. (22)) in the baseline regime appropriate to the n = 2 ring (Eq. (23)). Naïvely, one would simply determine the best-fit parameters α_L, α_R, and d_φ, and take the latter to provide the desired measurement of $d_{\phi }^{(2)}$.

In practice, however, a simple fitting method – such as the FindFit routine implemented in Mathematica – works only when the signal is extremely “clean,” in the sense that |V(u, φ)| displays no artifacts and follows the precise functional form of Eq. (22). This happens in ideal cases where the image perfectly resolves the photon ring and the baselines chosen for the fit correspond exactly to the domain over which the contributions from the n = 2 subring are dominant (Eq. (23)). Unfortunately, these cases are quite rare and the signal is typically less clean. As a result, the naïve fitting methods fail, even if – as in Fig. 3 – the periodicity is still manifestly present. A more robust approach is clearly called for.

Markov chain Monte-Carlo (McMC) methods offer such an approach and prove effective even when the signal is buried in noise. For instance, the data simulated by Gralla et al. (2020) for their experimental forecast contained so much noise so as to render the periodicity invisible to the eye (see panels (a) & (c) in their Fig. 8), and yet an MCMC analysis enabled the true ring diameter to be extracted with great precision (see panels (b) & (d) in their Fig. 8). While extremely robust, MCMC methods are also computationally intensive and thus a “middle-ground” approach is more desirable, especially for models without noise.

We devised such a procedure, which proved highly effective. We now describe it using Fig. 3 as an illustrative example:

1. For each φ, choose a baseline window in which to perform the fit (in Fig. 3, we chose the range u ∈ [2500, 2600] Gλ).

2. Identify the local extrema (both maxima and minima) of the visibility amplitude profile within this window.

3. Interpolate between maxima to obtain the superior envelope e_max(u) of the amplitude (green dashed curve in Fig. 3).

4. Interpolate between minima to obtain the inferior envelope e_min(u) of the amplitude (blue dashed curve in Fig. 3).

5. Using a simple fitting routine, obtain the ring diameter d_φ as the parameter d for which the refined model

   $$\begin{array}{c}\hfill V_{\mathrm{fit}}(u;d)=\sqrt{{\alpha }_L^2(u)+{\alpha }_R^2(u)+2{\alpha }_L(u){\alpha }_R(u)sin\left(2\pi du\right)}\end{array}$$(33)

   best fits the data, where the functions α_L/R(u) are defined as

   $$\begin{array}{c}\hfill {\alpha }_{L/R}(u)=\frac{e_{\max }(u)\pm e_{\min }(u)}{2}.\end{array}$$(34)

The precise interpolation method used to obtain the envelopes e_max/min(u) is not important; we resorted to cubic interpolation with the interp1d routine provided with SciPy. Likewise, the particular fitting method used in the last step is not important either, as the model is by construction very close to the data; we made use of the curve_fit routine, also provided with SciPy.

Some words of explanation are in order. First, the use of the refined model (Eq. (33)) allows us to consider signals that display the same periodic behavior as predicted by Eq. (22), but with variable extrema e_min and e_max. This relaxation significantly improves the method’s robustness and enables it to determine the periodicity of even moderately resolved rings, which can display non-monotonic modulations in their envelopes. We thus expect the method to remain effective even at lower resolution, thereby enabling a noticeable speed-up of parameter surveys. A trade-off of the method is that, in order to attain sufficient precision on the determination of the periodicity, we are required to examine the signal over a relatively larger interval of baselines; in practice, we used a window of 100 Gλ, corresponding to ∼15 periods.

This method also derives significant power from the flexible fall-off rate built into the visamp model (Eq. (33)), which generalizes the u^−1/2 power-law fall-off exhibited by the analytic formula (Eq. (22)) to arbitrary damping rates. Such a step was already taken by Johnson et al. (2020), who noticed that their signal presented damped oscillations with envelope e_max(u) = u^−1/2e^−(uw)ζ, with two additional parameters w and ζ needed to obtain a good fit (see their Fig. 4). The present approach is more general still. The use of an envelope tailored to the oscillation damping provides more robustness and allows one to ignore other ringing patterns that could arise from pixel effects, or even beating produced by interactions with other signals of larger periodicity: as shown in Fig. 3, this method works even when such beats produce a local increase of the visamp. Its drawback is that astrophysical noise (ignored here) may result in ill-defined envelope functions for Eq. (34).

We also wish to emphasize that Eq. (21) is only the leading term in an expansion in powers of 0 < w/d ≪ 1. Interestingly, preliminary analysis suggests that subleading corrections set by the angle-dependent, nonzero thickness of the ring w(φ) could also predict a beating pattern; analytic study of these corrections will be the subject of future work.

Finally, we note that while Eq. (33) is symmetric under L ↔ R interchange, this symmetry is broken by the definitions of Eq. (34), according to which α_L measures the oscillations’ mean and α_R their amplitude, with α_R < α_L by convention.

3.2.3. Multi-peaked distribution for the diameter

The above fitting procedure returns a “best-fit” diameter d_φ that minimizes the normalized root-mean-square deviation (RMSD)

$$\begin{array}{c}\hfill {\mathrm{RMSD}}_u(d)=\frac{\sqrt{{\langle {[V_{\mathrm{fit}}(u;d)-|V(u)|]}^2\rangle }_u}}{{\langle V_{\mathrm{fit}}(u;d)\rangle }_u},\end{array}$$(35)

where ⟨⋅⟩_u denotes an average over the chosen baseline window.

Often, this d_φ is only a local (rather than global) minimum of RMSD_u(d), in which case it is not the absolute best-fit diameter; even so, it still gives a good fit so long as 0 < RMSD_u(d_φ)≪1. This is illustrated with the example in Fig. 4, for which the above procedure yields a diameter d_φ = 37.312 μas that is a local (but not global) minimum of Eq. (35). Nevertheless, the resulting model V_fit(u; d_φ) (dashed yellow curve) closely tracks the visamp |V(u)| (solid red curve), as indicated by their small deviation (shown in the left panel inset) and as measured by the correspondingly small numerical value of RMSD_u(d_φ).

Still, the global minimum of RMSD_u(d) is d_φ = 37.115 μas, and with this diameter the model V_fit(u; d_φ) (dashed black curve) provides an even better fit to |V(u)|: their deviation (left inset) is even narrower, as measured by the smaller value of RMSD_u(d_φ).

By definition, the global minimum of RMSD_u(d) is always the absolute best-fit diameter. Equivalently, it is also determined as the global maximum of the “goodness-of-fit” measure

$$\begin{array}{c}\hfill g(d_{\phi })\equiv e^{-{\mathrm{RMSD}}_u(d_{\phi })}.\end{array}$$(36)

We show g(d_φ) for the above example in the right panel of Fig. 4, where we recognize the tallest peak to be the global maximum d_φ = 37.115 μas, with the next local maximum corresponding to the naïve best-fit parameter d_φ = 37.312 μas. Surprisingly, we also observe an entire periodic sequence of local maxima (with several providing a good fit) separated by a gap of Δd_φ ≈ 0.2 μas.

This multi-peak structure was already observed by Gralla et al. (2020) in their experimental forecast. Its physical origin can be intuitively understood as follows.

Over a sufficiently narrow baseline window, one may view u ≈ u_w ≫ 1/d as fixed. Within such a window, the universal form (Eq. (22)) of the visibility amplitude is periodic in the diameter, as |V(u)| is approximately invariant under shifts d_φ → d_φ + k/u_w for integer k ∈ ℤ. As a result, the probability distribution for d_φ is multi-peaked, with peaks separated roughly by Δd_φ ≈ 1/u_w. The example in Fig. 4 confirms this argument: the observed gap of Δd_φ ≈ 0.2 μas between peaks corresponds to a baseline length

$$\begin{array}{c}\hfill u_w\approx \frac{1}{2\times {10}^{-7}\text{arcsec}}\approx \frac{1.03\times {10}^{12}}{\text{rad}}\approx 1030\text{G}\lambda ,\end{array}$$(37)

which, as expected, lies well within the chosen baseline window of u ∈ [1000, 1045] Gλ.

Another intuitive way to understand this degeneracy comes from the observation that the number of damped oscillations (or “hops”) from the origin u = 0 to the baseline window u ≈ u_w is

$$\begin{array}{c}\hfill N_w\approx u_w{d}_{\phi }.\end{array}$$(38)

Hence, shifts d_φ → d_φ + k/u_w are equivalent to N_w → N_w + k, and this degeneracy in d_φ implies that the number N_w of hops from the origin cannot be measured with perfect precision using only the visibility amplitude on very long baselines u ≈ u_w: adding or substracting a few units to this period number N_w would barely shift the extrema within our fitting window, and as a result the fit V_fit(u, d_φ + k/u_w) remains acceptable for several values of k ∈ ℤ.

In practice, we accounted for this degeneracy by adjuting our fitting procedure to keep track of several of the most prominent peaks of the multi-peak distribution, as follows:

1. For each φ, obtain a first estimate of d_φ by computing the mean distance between local maxima of the visamp |V(u)|.

2. Pick an interval of some width Δd_φ ≳ 10/u_w centered around this estimate and plot the goodness of fit g(d_φ) defined by Eqs. (39) and (40), as in the right panel of Fig. 4.

3. Find the peaks of g(d_φ) in this interval and keep the values of d_φ corresponding to local maxima above a chosen peak height, together with the associated goodness of fit g(d_φ).

The resulting list should contain several diameters d_φ that locally maximize g(d_φ), including the global maximum. It is important to note, however, that the latter is not always the “true” diameter d_φ of the ring as measured from its image, which can sometimes correspond to one of the smaller peaks in g(d_φ). There are also cases where the multi-peaked distribution g(d_φ) is relatively flat, so we do not obtain a good measurement at that angle φ. Despite such failures, one can usually still infer the true diameter d_φ by carrying out this analysis at multiple angles, as we next explain.

3.2.4. Circlipse fit and test of the Kerr hypothesis

Having determined the absolute best-fit ring diameter d_φ at every angle φ around the image, the final step of the GLM test is to fit it to the GR-predicted functional form of the ring (Eq. (28)),

$$\begin{array}{c}\hfill d_{\mathrm{fit}}(\phi ;\mathit{R})=2R_0+2\sqrt{R_1^2{sin}^2(\phi -{\phi }_0)+R_2^2{cos}^2(\phi -{\phi }_0)},\end{array}$$(39)

with parameters R = {R₀,R₁,R₂,φ₀}. In practice, we carried out this fitting by examining 36 angles in the range [0° ,180° ) spaced at regular intervals of 5°. That is, we used φ ∈ {0°,5°,…,175°}. The best-fit parameters R_fit are then obtained by minimizing

$$\begin{array}{c}\hfill {\mathrm{RMSD}}_{\phi }(\mathit{R})=\frac{\sqrt{{\langle [d_{\mathrm{fit}}(\phi ;\mathit{R})-d_{\phi }{]}^2\rangle }_{\phi }}}{{\langle d_{\mathrm{fit}}(\phi ;\mathit{R})\rangle }_{\phi }},\end{array}$$(40)

where ⟨⋅⟩_φ denotes an average over the chosen baseline angles.

If 0 < RMSD_φ(R_fit)≪1, then the circlipse fit is good and the Kerr hypothesis passes the GLM test at a level of precision given by the normalized root-mean-square deviation RMSD_φ(R_fit)⁴. Otherwise, the fit is poor and the Kerr hypothesis fails the test.

As explained in the last section, this naïve approach is not robust because of the multi-peaked distribution for d_φ: keeping only the absolute best-fit diameter d_φ may yield a suboptimal fit, since it may force us to use a wrong diameter at some angles φ.

However, we found a simple way to remedy this by using the data collected at all angles simultaneously, as follows:

1. For each angle φ, compute g(d_φ), identify several peak values of d_φ, and display them all together as in Fig. 5.

2. Separate these peaks into multiple circlipse-shaped subsets 𝒞_i, each with “joint goodness-of-fit” measure

   $$\begin{array}{c}\hfill g({\mathcal{C}}_i)\equiv {\displaystyle \prod _{\phi }g(d_{\phi }).}\end{array}$$(41)

3. Finally, fit each 𝒞_i to a circlipse shape (Eq. (39)) with normalized RMSD(𝒞_i)≡RMSD_φ(R_fit) given by Eq. (40).

This “multi-fit method” results in several circlipses 𝒞_i that are not equally likely, as measured by their different values of g(𝒞_i). In favorable examples such as the one in Fig. 5, the likeliest 𝒞_i (i.e., the one with maximal g(𝒞_i), shown in dark blue in the right panel) does turn out to be the “true” circlipse shape of the actual n = 2 ring in the image⁵.

That said, there is no reason for this to be the case in general, and in some cases the “true” circlipse in the underlying image may well differ from the likeliest circlipse 𝒞^max as determined by this analysis. Of course, there is no other way to infer the true circlipse in a real experiment, so in principle we must take it to be the likeliest one 𝒞^max, even if it does not provide the best fit to the circlipse shape, that is, even if RMSD(𝒞^max)≠min_iRMSD(𝒞_i).

In practice, the least favorable examples present a handful of circlipses ${\mathcal{C}}_i^{\text{max}}$ sharing similar near-maximal likelihoods g(𝒞_i). In such cases, we may only infer the true ring diameter up to a degeneracy of a few periods Δd_φ ≈ 1/u_w, but we may still report a test of the Kerr hypothesis at a level of precision given by the maximal RMSD(${\mathcal{C}}_i^{\text{max}}$) within the set of best-fit circlipses ${\mathcal{C}}_i^{\text{max}}$.

4. Importance of the choice of baseline window

Before reporting the results of our parameter surveys, we first describe some new and important features of the GLM method that we noticed in the course of our investigation. In order to determine the diameter $d_{\phi }^{(2)}$ of the n = 2 ring, it is crucial to sample the visamp |V(u, φ)| in the appropriate regime (Eq. (23)), that is, on baselines long enough to resolve the width of the n = 1 ring (but not that of the n = 2 ring) at image angle φ_ρ = φ:

$$\begin{array}{c}\hfill \frac{1}{w_1({\phi }_{\rho })}\ll u\ll \frac{1}{w_2({\phi }_{\rho })}.\end{array}$$(42)

It is only in this regime that the signature of the n = 2 ring can dominate the visamp and produce damped periodic oscillations (Eq. (22)) that encode the ring diameter $d_{\phi }^{(2)}$. In the example of Fig. 2, this requires us to choose a baseline window u ≳ 100 Gλ.

The angle-dependence of Eq. (42) is important. At low inclinations, this dependence is relatively weak, and the baseline threshold b₂(φ)≈1/w₁(φ) past which the n = 2 ring dominates |V(u, φ)| is approximately constant in φ. As a result, there exist baseline windows for which a ring diameter measurement works for all angles φ. On the other hand, at high inclinations, b₂(φ) may vary so much that no choice of baseline can satisfy the condition (Eq. (42)) for all angles simultaneously; worse, the visamp may even cease to encode the n = 2 ring diameter altogether.

These difficulties can be quantitatively fleshed out as follows. The GLM method can only succeed for a fixed choice of baseline window u ≈ u_w that satisfies Eq. (42) for all angles, so that

$$\begin{array}{c}\hfill \underset{\phi }{max}{b}_2(\phi )\ll u_w\ll \underset{\phi }{min}{b}_3(\phi ).\end{array}$$(43)

Here, b_n(φ)≈1/w_n − 1(φ) denotes the threshold past which the visamp contribution from the nth ring overtakes its predecessor’s. These thresholds are well-defined when there is a clear transition between the baseline regimes in which each subring’s signature dominates; however, these are not always cleanly delineated. To understand why, we recall that |V(u, φ + π)| = |V(u, φ)|, which implies that the baseline thresholds must obey (when they exist)

$$\begin{array}{c}\hfill b_n(\phi +\pi )=b_n(\phi ).\end{array}$$(44)

On the other hand, the subring widths w_n(φ_ρ) are only subject to the periodicity constraint w_n(φ_ρ + 2π) = w_n(φ_ρ) and may display significant variation with image angle φ_ρ ∈ [0, 2π). In particular, w₁(φ_ρ) and w₁(φ_ρ + π) may be very different at high inclinations, in which case b₂(φ_ρ) fails to be sharply defined. Worse, while it must always be the case that w₂(φ_ρ) < w₁(φ_ρ), it may happen that w₂(φ_ρ)∼w₁(φ_ρ + π) for some angles φ_ρ, in which case the periodic ringing of the visamp |V(u, φ)| at φ = φ_ρ is affected by both the n = 1 and n = 2 rings but encodes neither’s diameter. We now describe these inclination-dependent effects in detail.

4.1. Low to moderate inclinations

At low-to-moderate inclinations, each subring produces a distinct signature in its own baseline regime, and so the threshold b₂(φ) exists. This is the case illustrated in Figs. 2 and 6. According to Eq. (15), the next thresholds are related to it by

$$\begin{array}{c}\hfill b_{n+1}(\phi )\approx e^{\gamma (\phi )}{b}_n(\phi ),\end{array}$$(45)

an approximate Kerr-lens equation that becomes exact as n → ∞ and is already an excellent approximation for n ≳ 2.

On the other hand, one cannot determine b₂(φ) from first principles, as the first subring’s width w₁(φ_ρ) must be computed for each model separately. We may thus rewrite Eq. (43) as

$$\begin{array}{c}\hfill \underset{\phi }{max}{b}_2(\phi )\ll u_w\ll \underset{\phi }{min}{e}^{\gamma (\phi )}{b}_2(\phi ),\end{array}$$(46)

where b₂(φ) is model-determined (astrophysics-dependent) and the demagnification factor e^γ(φ) is universal (GR-predicted).

Letting γ₀ = min_φγ(φ_ρ) denote the minimal demagnification factor, we see that a sufficient (though not necessary) condition for the existence of a baseline window u_w satisfying Eq. (46) is

$$\begin{array}{c}\hfill \underset{\phi }{max}{b}_2(\phi )\ll e^{{\gamma }_0}\underset{\phi }{min}{b}_2(\phi ).\end{array}$$(47)

At low inclinations, where b₂(φ)≈b₀ is relatively flat and e^γ₀ takes values of approximately 10 ∼ 20, this sufficient condition is manifestly satisfied, guaranteeing the existence of a suitable choice of baseline window over which to carry the measurement. At higher inclinations, however, e^γ₀ ∼ 1 (it is in fact exactly unity for equatorial inclination and maximal spin a = 1). Moreover, b₂(φ) varies significantly, so this condition is clearly violated. This does not necessarily imply that Eq. (43) cannot be satisfied, but it does mean that some choices of baseline window u_w that are suitable at a given angle φ may not be suitable at other angles.

This effect is illustrated in Fig. 6, which displays a model from our parameter survey with BH spin a = 0.5 and moderate inclination i = 45°. Panels B–D illustrate the angle-dependence of the threshold b₂(φ), which varies from b₂(60° ) ≈ 200 Gλ to b₂(140° ) ≈ 700 Gλ. If one were to extract a diameter d_φ from the periodicity of the visamp in the baseline window u_w ≈ 600 Gλ, then one would be measuring the n = 2 ring diameter $d_{\phi }^{(2)}$ at φ = 60° (panel B) and the n = 1 ring diameter $d_{\phi }^{(1)}$ at φ = 140° (panel D). At an intermediate angle of φ = 102°, one would land in a transition region between the n = 1 and n = 2 regimes in which the visamp periodicity encodes neither ring’s diameter, but rather some “superposition” of the two.

Of course, at moderate inclinations, this issue can always be avoided by selecting a larger baseline window. For the example of Fig. 6, any choice u_w ≳ 1000 Gλ is sufficiently large: over such baseline windows, the visamp periodicity always encodes the n = 2 ring diameter (until the n = 3 signal takes over on even larger baselines u ≳ 2000 Gλ). Meanwhile, measurements of d_φ on insufficiently long baselines lead to tell-tale transitions between $d_{\phi }^{(1)}$ and $d_{\phi }^{(2)}$ that inevitably result in the characteristic, non-periodic pattern of data points shown in panel A of Fig. 6.

The reason is the following. For the angles in the blue zone, d_φ = $d_{\phi }^{(2)}$ is the “true” diameter of the n = 2 ring’s image, while for the angles in the gray zone, d_φ = $d_{\phi }^{(1)}$ is the “true” diameter of the n = 1 ring’s image. These regions are separated by green zones in which neither ring dominates the signal, so that d_φ is equal to neither $d_{\phi }^{(2)}$ nor $d_{\phi }^{(1)}$: in such regions, d_φ does not track the diameter of any actual image feature, and in particular it need no longer be a continuous curve. As a result, the absolute best-fit diameter d_φ – defined to maximize Eq. (36) – is allowed to “jump” across circlipse bands (and indeed, one can see two such jumps in the second green zone 150° < φ < 175° of panel A). In fact, such jumps are necessary to ensure the periodicity of the physical ring diameter $d_{\phi }^{(2)}$ in the blue zone, which must necessarily remain invariant under shifts φ → φ + π.

On the other hand, the continuous curves 𝒞_i that the human eye naturally joins points into need not be periodic as φ → φ + π; it is only the totality {C_i} of such curves that needs to respect that periodicity. It is therefore perfectly acceptable for these curves to continuously wrap into their successor as one looks up across panel A of Fig. 6. While these curves C_i are still separated by a gap Δd_φ ≈ 1/u_w, they are manifestly not circlipses, and their lack of periodicity is an indicator that the baseline window is inadequate for a full measurement of $d_{\phi }^{(2)}$.

To summarize, at low inclinations i ≲ 20°, the baseline regimes corresponding to each subring are sharply delineated; in particular, the baseline threshold b₂(φ) exists and is relatively flat, so that if a baseline window u_w is suitable for a measurement of the n = 2 ring diameter $d_{\phi }^{(2)}$ at some angle φ, then it is almost surely suitable at every other angle around the baseline plane. At higher but still moderate inclinations i ≲ 50°, the baseline threshold b₂(φ) still exists but now exhibits substantial variation. Because of this, some baseline windows u_w may be suitable for a measurement of the n = 2 ring diameter $d_{\phi }^{(2)}$ at certain angles φ, but not others, where they instead pick up the n = 1 ring diameter $d_{\phi }^{(1)}$ or some unphysical “superposition” of the two. Nonetheless, it should be possible to tell if this is happening from the shape of d_φ, which would display the general pattern shown in panel A of Fig. 6. Moreover, at these inclinations, this difficulty can always be avoided by selecting a sufficiently large baseline window.

4.2. High inclinations

At higher inclinations i ≳ 60°, subrings do not always produce clearly delineated signatures in individual baseline regimes, and as a result the baseline threshold b₂(φ) does not exist at all angles φ. Physically, this can occur because the baseline threshold must be π-periodic (Eq. (44)), while the subring images need only be 2π-periodic. In particular, while w₂(φ_ρ) < w₁(φ_ρ) for every angle φ_ρ around the image, at high inclinations it is usually the case that w₂(φ_ρ)∼w₁(φ_ρ + π) for angles φ_ρ opposite the spin axis. In that case, the visamp V(u, φ) at baseline angle φ = φ_ρ (or equivalently, φ = φ_ρ + π) is sensitive to the diameters of both the n = 1 and n = 2 subrings.

This phenomenon is clearly visible in the example shown in Fig. 7, which corresponds to a BH of spin a = 0.5 observed from a near-equatorial inclination i = 89°. Since w₂(φ_ρ)≈w₁(φ_ρ + π) for φ_ρ = 90°, it follows that the visamp profile |V(u, 90° )|, which is plotted in red in Fig. 8, cannot display a sharply delineation between regimes dominated by the n = 1 and n = 2 subrings. Nonetheless, this profile still displays oscillations. According to Eq. (15), their periodicity ought to be set by the diameter of a ring of fixed width in the image, which in this case must be the composite ring consisting of the upper half of the n = 1 ring joined together with the lower half of the n = 2 ring. In other words, we would expect the visamp for this model to still adopt the universal form (Eq. (22)) on long baselines, but with a projected diameter d_φ corresponding to a hybrid n = 1&2 ring of uniform width. In practice, this expectation is not quite realized in this model because of the direct (n = 0) image of the near-side of the equatorial disk, which cuts across the photon ring. Since this horizontal line is itself thin and bright like the hybrid n = 1&2 ring, the vertical intensity cross-section consists of three spikes; as a result, the visamp profile at angles φ ≈ 90° is sensitive to multiple projected diameters, including that between the upper half of the n = 1 ring and the n = 0 emission line, and that between the n = 0 emission line and the lower half of the n = 2 ring. These diameters are roughly half that of the hybrid ring, and indeed the visamp consists of a beating pattern with roughly twice the period expected for the signature of the photon ring.

5. Parameter survey for the emission profile

5.1. Description of the models

In this paper, we restrict our attention to the “thin-disk” emission models presented in Sect. III.A of Gralla et al. (2020) and Sect. 3.2 of Chael et al. (2021), which consist of emission profiles that are stationary, axisymmetric, and confined to the equatorial plane. The emitters forming the disk are assumed to describe stable, circular-equatorial orbits all the way down to the radius of the innermost stable circular orbit (ISCO), past which they follow the prescription of Cunningham (1975) for infall to the horizon. This specifies the redshift factor g of light rays emitted from the disk, whose observed (bolometric) intensity is then computed as

$$\begin{array}{c}\hfill I(\alpha ,\beta )={\displaystyle \sum _{n=0}^{N(\alpha ,\beta )}}{f}_{n}J\left(r_s^{(n)}\right)g^4(r_s^{(n)},\alpha ),\end{array}$$(48)

where J(r) is an arbitrary radial emission profile, $r_s^{(n)}=r_s^{(n)}(\alpha ,\beta )$ denotes the radius at which a light ray crosses the equatorial plane for the nth time after being shot back into the geometry from position (α, β) in the observer sky, and N + 1 is the total number of times the light ray intersects the equatorial plane⁶. Here, we also introduced a “fudge factor”

$$\begin{array}{c}\hfill f_n=\{\begin{array}{cc}1\hfill & n=0,\hfill \\ 1.5\hfill & n\in \{1,2\},\hfill \\ 0\hfill & n\ge 3,\hfill \end{array}\end{array}$$(49)

whose effect is to remove n ≥ 3 subrings from the image (which would be underresolved), while enhancing the intensity of the n = 1 and n = 2 subrings relative to the direct n = 0 image. This enhancement can account for the effects of geometric thickness and results in images whose intensity cross-sections match those observed in time-averaged GRMHD-simulated images (Johnson et al. 2020; Chael et al. 2021). In this sense, the equatorial model (Eq. (48)) provides a good approximation for realistic accretion flows, and its simplicity and relatively low computational cost make it ideal for the survey of emission profiles undertaken here.

While the equatorial source profile J(r) can in principle take any form, in this survey we considered only profiles of the form

$$\begin{array}{c}\hfill J(r;\mu ,\gamma ,\sigma )=\frac{e^{-\frac{1}{2}{[\gamma +\mathrm{arcsinh}\left(\frac{r-\mu }{\sigma }\right)]}^2}}{\sqrt{{(r-\mu )}^2+{\sigma }^2}}.\end{array}$$(50)

This is a three-parameter subset of Johnson’s S_U-distribution (the fourth parameter δ is set to unity), already adopted by Gralla et al. (2020) for its versatility: this functional form provides a simple way to generate smooth profiles with emission typically concentrated near the event horizon of the BH, and displaying substantial variation with changes in μ, γ, and σ. In short, these three parameters respectively determine a location (the peak of the distribution), a shape (its asymmetry), and a scale (its width). We now describe in detail how their variation affects ring width.

First, μ influences the location of the emission profile’s peak. We considered values in the range [r₋, 2r₊], where r_± denote the horizon radii given in Eq. (14). At the lower end of this range (left column of Fig. 9), the peak is very close to or behind the event horizon, but it moves farther out near the upper end of the range (right column of Fig. 9). The observed thickness of the photon subrings is not significantly affected by this variation.

Second, γ controls the profile’s asymmetry. We considered values in the range [ − 2, 2]. When γ < 0 (left column of Fig. 9), the steep end of the distribution is to the left and (mostly) cut off by the horizon, leaving only a gentle-sloping tail in the emission profile. When γ > 0 (right column of Fig. 9), the reverse is true: the profile’s steep end lies to the right and remains visible outside the horizon, while its gentle-sloping side is cut off. Hence, the rings tend to become thinner as γ increases.

Third, σ encodes the profile’s width. We considered values in the range [0.25M, 1.5M]. At the upper end of this range (left column of Fig. 9), the profile is broad enough to produce thick rings, which makes the GLM test easy to carry out on a baseline window u ∈ [450,550] Gλ for all our considered values of μ and γ. At the lower end (right column of Fig. 9), the profile is a steep spike (provided that the distribution’s peak lies outside the horizon) that produces narrower rings: as a result, the n = 1 ring dominates on longer baselines, and a measurement of the n = 2 ring diameter may require (in the most extreme cases) a baseline window u_w ≳ 2000 Gλ.

Images with thinner rings have a larger baseline threshold b₂(φ) for carrying out a measurement of the n = 2 ring diameter $d_{\phi }^{(2)}$. This trend is illustrated in Fig. 9 with two extreme examples that display maximally thick rings (left column) and maximally thin rings (right column) within our considered parameter range.

5.2. Choice of parameters for the survey

In order to assess the robustness of the GLM method described in Sects. 2 and 3, we conducted a survey over the emission-profile parameter space {μ,γ,σ}. We used “intermediate” values for the BH spin and inclination of a = 0.5 and i = 45°, and examined a total of 100 emission models corresponding to all the possible combinations of parameters with values presented in Table 1.

As a cross-check, we repeated the survey using the currently favored values of BH spin and inclination for M87*. The power in the relativistic jet produced by M87* suggests that this BH is highly spinning, since most low-spin GRMHD models do not give rise to a sufficiently powerful jet (Event Horizon Telescope Collaboration 2019a). We therefore took a = 0.94, the highest value considered in EHT simulations (Event Horizon Telescope Collaboration 2019e). We also assumed that the spin axis of the BH and of its accretion disk are aligned with its forward jet, which has a measured inclination of i ≈ 17° (Walker et al. 2018). Since r₋ ≈ r₊/2 for this value of BH spin, we examined the same combinations of parameters except for those with μ = r₊/2; this amounted to a total of 80 emission profiles in our survey.

5.3. Results

We investigated these models to address the following questions:

• Q1.

  Is the GLM test possible in principle? In particular, does $d_{\phi }^{(2)}$ follow the GR-predicted functional form of a circlipse (Eq. (28))?

• Q2.

  If so, how do the best-fitting circlipse parameters R_fit vary with the astrophysical profile? In particular, what can we learn about the BH parameters from the photon ring shape?

Our survey led us to the following answers:

• A1.

  Yes. In all the models considered here, we identified a regime in which the n = 2 ring dominates the visamp and extracted the projected ring diameter $d_{\phi }^{(2)}$ from the signal’s periodicity. For reasons discussed in Sect. 4, this required using different baseline lengths, sometimes as low as u ≳ 150 Gλ (Fig. 2) or as high as u ≳ 2000 Gλ (bottom right panel of Fig. 9). We always found an excellent fit to the circlipse shape (Eq. (39)) as measured by the normalized RMSD (Eq. (40)): in all cases,

  $$\begin{array}{c}\hfill {\mathrm{RMSD}}_{\phi }({\mathit{R}}_{\mathrm{fit}})<{10}^{-4}.\end{array}$$(51)

• A2.

  By analogy with Eq. (26), it is helpful to repackage the best-fit parameters R_fit = {R₀,R₁,R₂,φ₀} describing the projected diameter of the n = 2 ring into maximal/minimal diameters⁷

  $$\begin{array}{c}\hfill d_{+}=d_{{\phi }_0+\pi /2}^{(2)}=2(R_0+R_1),d_{-}=d_{{\phi }_0}^{(2)}=2(R_0+R_2)\end{array}$$(52)

  of the associated circlipse shape (d₋ ≤ d₊ since R₁ ≥ R₂ by definition). We plot the set of (d₊, d₋) obtained from all our emission models in Figs. 10 and 11, and provide mean values and standard deviations for (d₊, d₋) in Tables 2 and 3. These results give a good overview of the variabilty in the circlipse shape obtained from different astrophysical source profiles. We return to the question of spin inference in Table 4 below.

As discussed in Sect. 2.9, one expects the n = 2 photon ring to closely track the critical curve. This is verified in practice, though we find that the choice of astrophysical source profile still has a significant impact on the shape of the n = 2 ring. In most of our models, the n = 2 ring is slightly larger in diameter than the critical curve, so that d₊ > d_∥ and d₋ > d_⊥. In particular, we find that broader emission profiles produce thicker n = 2 rings and therefore a larger difference in these diameters. Conversely, narrow profiles that produce thinner rings lead to diameters that are closer to (and sometimes, even smaller than) those of the critical curve. Either way, conflating the shape of the n = 2 ring with that of the critical curve can lead to erroneous conclusions.

Nevertheless, it is possible to approximate the critical curve diameters (d_∥, d_⊥) by (d₊, d₋), and to then use these parameters to estimate the BH spin and inclination (assuming a mass prior) as described in Sect. 2.8. However, the precision of this spin and inclination determination is significantly limited by the spread of results for the n = 2 ring diameters at a given (a, i) value. To illustrate this, Figs. 10 and 11 also display diameters of critical curves that are compatible with the measured n = 2 ring shape, but which correspond to (a, i) values that are different from the “true” values of the BH spin and inclination used in the model (green and pink stars). The constraints on these BH parameters derived from the inferred n = 2 ring diameters have substantial error bars, for which we provide lower bounds in Table 4.

5.4. Range of allowed phovals within the n = 2 lensing band

These results highlight that the critical curve is not necessarily the most relevant concept for the analysis of the n = 2 photon ring shape. For our thin-disk models, a more fruitful approach is to consider the n = 2 lensing band described in Sect. 2.5, since it contains all the possible geometric shapes of the n = 2 photon ring for a given BH spin and inclination. The critical curve lies within this band as well, but so do the critical curves of BHs with neighboring spins and inclinations, as Figs. 10 and 11 indicate.

Since Sect. 2.9, we have used “photon ring shape” to refer to the simplest observational target for a measurement of the n = 2 ring: its projected diameter $d_{\phi }^{(2)}$, which matches that of a circlipse (Eq. (24a)). In this section, we study the full geometric shape of the n = 2 ring, which typically includes centroid motion effects (Eq. (24b)) and is therefore not a circlipse in general. Instead, we expect this geometric shape to be close to a phoval (Eq. (24)); that is, the n = 2 ring’s projected position should closely match

$$\begin{array}{cc}& f_{\mathrm{ph}}(\phi ;\mathit{R})=R_0+\sqrt{R_1^2{sin}^2(\phi -{\phi }_0)+R_2^2{cos}^2(\phi -{\phi }_0)}\hfill \\ & +(X-\chi )cos(\phi -{\phi }_0)+arcsin[\chi cos(\phi -{\phi }_0)],\hfill \end{array}$$(53)

for some choice of parameters R = {R₀,R₁,R₂,X,χ,φ₀}. (Here, as in Eq. (28), we have included a rotation angle φ₀ to account for the arbitrary orientation of the BH spin axis in the image.) This expectation is based on the observation (reviewed in Sect. 2.8) that the Kerr critical curve is always a phoval to high accuracy, together with the fact that the photon rings converge very rapidly (exponentially fast in n) to the critical curve (Eq. (13)).

Since we expect the geometric shape of the n = 2 photon ring to be that of a phoval, it is interesting to consider the range of all phoval shapes that can possibly fit inside the n = 2 lensing band. Since phovals are described by six parameters, it is convenient to project this 6D region of allowed phoval shapes down onto the 2D plane of diameters (d₊, d₋) via Eq. (52). The resulting region in the (d₊, d₋) plane gives a good sense of the allowed range of possible n = 2 ring shapes.

For the two choices of BH spin and inclination in our survey of emission profiles, this projection results in the blue regions displayed in Figs. 10 and 11. Their numerical computation is nontrivial and we relegate the details to Appendix B. We note that since d₋ ≤ d₊, these regions must lie below the diagonal⁸.

As expected, all the emission profiles in our survey produce n = 2 photon rings that lie within the n = 2 lensing band, with diameters (d₊, d₋) lying within the blue region of allowed phoval shapes that fit inside the band. The band’s outer and inner edges can themselves be fitted to the phoval shape (Eq. (53)) with excellent agreement, and their diameters (d₊, d₋) – shown in Figs. 10 and 11 as blue and orange dots, respectively – lie on opposite corners of this blue region. Remarkably, all of our astrophysical source profiles produce phovals whose diameters (d₊, d₋) lie within a surprisingly narrow subregion of the full allowed region.

More precisely, these diameters (the gold dots in Figs. 10 and 11) form a relatively straight line connecting the corners of the blue region that correspond to the outer and inner edges of the lensing band. In addition, they are not uniformly scattered along this line; rather, they are tightly grouped near its center. Both of these facts may be intuitively understood as follows. We recall from Sect. 2.5 that the inner and outer edges of the nth lensing band are the nth lensed images of the equatorial circles of radius r = r₊ (the event horizon) and r → ∞, respectively. Images of intermediate radii r₊ < r < ∞ form curves that fill the lensing band and interpolate between its inner and outer edges. From Fig. 6 of Gralla & Lupsasca (2020a), these curves roughly appear to be dilated versions of each other (with γ(φ) controlling the angular-dependent dilation) and indeed, we find that they map to a series of points in the (d₊, d₋) plane tracing a diagonal line that connects the blue and orange dots (shown in gray in Figs. 10 and 11): the same line that the gold dots are concentrated around (see Appendix B for details of its computation).

Since an equatorial disk is a superposition of equatorial rings, the n = 2 photon ring is also a superposition of the n = 2 images of said rings. That is, its position in the (d₊, d₋) plane should be a weighted average of points along the diagonal line, and must therefore also lie on this line. Moreover, this position may be associated to an “effective” radius of equatorial emission lying somewhere in the line’s middle, away from its extremities. We intend to further explore this connection in future work.

6. Parameter survey for BH spin and inclination

6.1. Choice of parameters for the survey

To better assess the robustness of the GLM method described in Sects. 2 and 3, we also conducted a parameter survey over BH spins and inclinations (a, i). Simulating 100 different emission models for each (a, i), as we did in Sect. 5.2, would have been prohibitively expensive here; instead, we focused on a handful of “representative” emission profiles of the form of Eq. (50) that spanned the full range of outcomes observed in that survey:

• A “typical 1” profile for which the diameters (d₊, d₋) of the n = 2 ring tended to be near their average value (that is, near their mean value averaged over all emission models);

• An “overestimated” profile for which the diameters (d₊, d₋) of the n = 2 ring were among the largest;

• An “underestimated” profile for which the diameters (d₊, d₋) of the n = 2 ring were among the smallest;

• A “typical 2” profile for which they were slightly subaverage;

• A “narrow” profile with a sharp peak (small σ) for which the circlipse fit was more difficult (required longer baselines) and led to the smallest values for (d₊, d₋) when successful.

Table 5 lists the parameters {μ,γ,σ} of these emission profiles.

In this survey, we considered all the combinations (a, i) with a ∈ {0, 0.25, 0.5, 0.75, 0.99} and i ∈ {1° ,22° ,45° ,67° ,89° } (avoiding the fine-tuned cases 0° and 90°). In the Schwarzschild case a = 0, the inclination is defined relative to the equatorial disk, whose spin is always aligned with the BH spin.

6.2. Results

We investigated these models to address the two questions posed in Sect. 5.3. Our survey led us to the following answers:

• A1.

  Yes, up to i ≲ 45°. In all the models considered here, we identified a regime in which the n = 2 ring dominates the visamp and extracted the projected ring diameter $d_{\phi }^{(2)}$ from the signal’s periodicity. For reasons discussed in Sect. 4, the diameters inferred from images at inclinations of i = 67° and i = 89° sometimes presented features that did not allow for a full circlipse fit as described in Sects. 2 and 3. However, for low-to-moderate inclinations i ≲ 45°, we always found an excellent fit to the circlipse shape (Eq. (39)) as measured by the normalized RMSD (Eq. (40)): in all such cases,

  $$\begin{array}{c}\hfill {\mathrm{RMSD}}_{\phi }({\mathit{R}}_{\mathrm{fit}})<5\times {10}^{-4}.\end{array}$$(54)

  For most configurations, these best fits were obtained in a window u ∈ [1000,1100] Gλ. However, the “narrow” profile produced n = 2 rings that were so thin that they only became dominant on much longer baselines, requiring the fit to be carried out in a baseline window u_w ≈ 2000 Gλ. This also became necessary at i = 45° for the “underestimated” profile.

• A2.

  As in Sect. 5.3, we extracted a pair of diameters (d₊, d₋) for each emission model. This enabled us to produce, for each BH spin and inclination, a plot like the one shown in Fig. 12 for the case of spin a = 0.5 and inclination i = 45°. Since our models were handpicked to form a “representative” sample spanning the range of possible diameters (d₊, d₋) that one can obtain for each spin and inclination, it is informative to draw a bounding box containing them (like the green box in Fig. 12) and then measure its size. The dimensions of these bounding boxes are listed for each (a, i) in Table 6 below. As we explain in Sect. 6.3, this analysis bears relevance for the “inverse problem” of how to infer BH spin and inclination from the n = 2 photon ring shape.

We first note that this survey over BH spins and inclinations confirms one of our crucial observations from Sect. 5.4: for each choice of (a, i), all our emission models produce n = 2 rings with diameters (d₊, d₋) that are very tightly clustered near the line connecting the inner and outer edges of the n = 2 lensing band. In other words, the five observed ring shapes corresponding to the models in Table 5 occupy a very small subregion (green box in Fig. 12) of the full region of allowed phoval shapes (blue box). As explained in Sect. 5.4, diameters along this line correspond to n = 2 images of equatorial circles of different radii r_s ∈ [r₊, ∞), so we may equivalently say that the n = 2 rings produced by our emission models have the same diameters (d₊, d₋) as images of equatorial rings of emission concentrated near a fixed radius r_s.

Rather than pursue this association with “effective” emission radii, we employed a more agnostic parameterization of this line in the (d₊, d₋) plane and used an affine parameter s such that s = 0 corresponds to the inner edge of the lensing band and s = 1 to its outer edge. We also introduced δ, a perpendicular distance from the line, thereby defining a map (d₊, d₋)↔(s, δ). Then, the diameters produced by our models are all contained in a bounding box extending from some s_min to some s_max in the direction parallel to the line, with a perpendicular width of |δ|_max on each side on the line. For each (a, i), the resulting triplet (s_min, s_max, |δ|_max) characterizes the spread of the ring diameters obtained from our models. These triplets are listed in Table 6.

We note that, aside from the extreme case with a = 0.99, i = 45°, we typically have s_min ≳ 0.18 and s_max ≲ 0.75. Hence, our models produce rings that are concentrated away from the edges of the lensing band, as already observed in Sect. 5.4. This is consistent with an “effective” emission radius 2M ≲ r_s ≲ 7M, in line with the choice of emission peak locations in Table 5.

6.3. Inferring BH spin and inclination from photon ring shape

In this survey, we have solved the “forward problem” of how to compute the diameters (d₊, d₋) of the n = 2 photon ring given a choice of emission profile, and a BH spin and inclination (a, i). Purely geometric considerations constrain the observed photon ring to lie within the n = 2 lensing band, which corresponds to the blue regions in Figs. 10–12. Perhaps suprisingly, we have found that the ring shapes produced by our emission models occupy a much smaller subregion (the green bounding box in Fig. 12) of this full region of geometrically allowed shapes.

This observation is highly relevant to the “inverse problem” of how to infer BH spin and inclination from a measurement of the n = 2 ring diameters (d₊, d₋). Given such a data point, one could try to determine which blue regions (n = 2 lensing bands) it is contained within. The answer would provide the set of (a, i) for which lensing would allow the observed photon ring shape.

Evidently, given the considerable size of these blue regions, these purely geometric constraints on (a, i) (which are derived from lensing alone) are relatively weak. However, if one were to make additional assumptions about the astrophysics of the source (such as the class of emission models being considered), then one would obtain narrower subregions for each (a, i) (the smaller green boxes), and a given data point would therefore be compatible with a smaller range of BH spins and inclinations.

Since our five models were picked to span the full range of s observed in our survey over emission profiles in Sect. 5, we expect s_min and s_max to be close to the values they would take in more general situations, at least so long as the emission peaks near an equatorial source ring at small radius outside the horizon. On the other hand, our values of |δ|_max are probably less reliable, and would likely increase if we were to consider geometrically thick models in which the emission is not well-approximated by an equatorial source ring.

Thus, we do not claim here to have definitively computed the bounding boxes for “realistic” choices of astrophysical profiles. Rather, we hope to have provided useful organizing principles for tackling the problem of spin inference: the choice between blue and green boxes clarifies the interplay between geometric and astrophysical assumptions, with the additional constraining power provided by the latter made manifest in the much smaller size of the green regions relative to the blue ones.

An additional benefit of this framework is that it enables us to quantify the uncertainty in a given spin estimate by computing the range of (a, i) that produce green or blue boxes containing a given measurement (d₊, d₋) of the n = 2 photon ring shape. We leave a more detailed examination of this problem to the future.

7. Discussion

7.1. Limitations of this work

We have revisited the proposal made by Gralla et al. (2020) for a test of the Kerr hypothesis based on the photon ring shape. We have established the robustness of their proposed test in a large class of models, provided that the BH inclination is not too high. This confirms that M87* holds great promise as an observational target for this test. However, as pointed out in the introduction, this conclusion remains subject to some caveats.

First, we have only considered simple equatorial emission models. It is important to test whether the GLM method remains viable in the presence of more realistic accretion flows (such as fluctuating, geometrically thick structures), for which effects such as absorption can become very significant (for instance, the interferometric signature of the photon ring is suppressed as the emission region grows thicker and more opaque). Encouraging preliminary results in this direction will be reported elsewhere (Vincent et al. 2022).

Second, we have not presented a complete, quantitative study of the range of baseline lengths for which the method works. In particular, we have not identified an “ideal” baseline window in which the GLM test succeeds for the largest set of configurations (astrophysical profiles, BH spins, and inclinations), and we have largely ignored practical concerns in our theoretical study of the baseline domain. With the current wavelength of observations (1.3 mm), a space antenna would have to be almost 1.5 million km away from Earth – close to L₂! – to sample the visibility at u_w ∼ 1000 Gλ (Johnson et al. 2020).

Third, we have not included any noise in our simulations. Generating mock data with both instrumental and astrophysical noise is crucial, not only for demonstrating the feasibility of this proposed test, but also for developing a more robust, Bayesian approach to the multi-fit method described in Sect. 3. In their original analysis, Gralla et al. (2020) included instrument noise but not did not model astrophysical fluctuations, which can have a considerable impact on the photon ring signatures present in instantaneous snapshots. A new code tailored to include such fluctations – AART, for “Analytic Adaptive Ray Tracing” – has recently been developed (Cárdenas-Avendaño et al. in prep.) to tackle this specific problem and will inform the requirements for (incoherent) time-averaging.

While amplitude gain errors can be problematic for VLBI data analyses (e.g., Event Horizon Telescope Collaboration 2019d,f), it is worth noting that they may not pose a significant problem for a measurement of the photon ring diameter, which would only require measuring the periodicity of oscillations in the visibility amplitude (rather than their absolute scale).

7.2. Perspectives: Inferring black hole spin and inclination

If a future space mission targeting M87* successfully measured the shape of its n = 2 photon ring and confirmed that it matched the Kerr prediction (Eq. (39)), then the natural next step would be to assume the Kerr nature of the source and try to infer the BH spin and inclination (a, i) from the measured photon ring shape.

Our survey over emission profiles in Sect. 5 shows that this shape displays a residual dependence on the astrophysical source profile (even for the n = 2 ring). Though small, this dependence nonetheless introduces significant uncertainty in estimates of the spin based on the n = 2 ring shape alone, as evidenced by the lower bounds in Table 4. (This uncertainty is undoubtedly much greater for spin estimates based on the n = 1 ring shape.)

Spin estimates can be tightened with additional assumptions about the nature of the emission. As discussed in Sects. 5.4 and 6.2, the photon ring shapes produced by our emission models inhabit only a small subset of the full range of geometric shapes allowed by BH lensing. Thus, assuming that an observed photon ring shape is produced by a model within this class considerably narrows down the set of BH spins and inclinations compatible with the measurement, excluding many more values of (a, i) than could be ruled out based on the Kerr geometry alone.

It would be extremely interesting to quantify the reduction in uncertainty enabled by different astrophysical assumptions. This requires one to tackle the inverse problem of how to determine the values of (a, i) that are compatible with a given photon ring shape; we hope that our framing of this problem in Sect. 6.3 will prove useful for future discussions of spin inference. We note that other approaches based on the identification of the photon ring with the critical curve either fail to quantify the significant systematic error in their spin estimates (see Figs. 10 and 11), or else they effectively incorporate some version of our analysis.

Finally, spin estimates can also be substantially improved by combining the constraints derived from multiple image features. Early work suggests that measuring both the photon ring shape and that of the central brightness depression could significantly constrain BH spin and inclination (Chael et al. 2021). Of course, simultaneously resolving more than one ring (such as the n = 1 and n = 2 subrings, for instance) should also enable sharper spin estimates, though this question has so far only been explored in the simplest cases with circular subrings (Broderick et al. 2022).

7.3. Perspectives: Other photon ring tests

The prospect of measuring multiple subrings naturally opens the door to even more stringent tests of the Kerr hypothesis (or more precise techniques for determining the BH parameters).

For instance, comparing the widths of successive subrings at different angles around the image could in principle provide a measurement of the angle-dependent demagnification factors e^−γ(φ), and therefore of the Lyapunov exponents $\gamma (\stackrel{\sim }{r})$ that offer a direct probe of the Kerr geometry (Johnson et al. 2020; Gralla & Lupsasca 2020a). If the demagnification factors were not e^−γ(φ), then this might suggest that the spacetime is not Kerr (Wielgus 2021). Conversely, if γ(φ) did match the form predicted from Kerr, then these exponents would provide a function’s worth of constraints on the BH spin and inclination (a, i), which could allow these parameters to be inferred with great precision.

Likewise, polarimetric imaging of successive subrings could also be used to test the Kerr hypothesis, or conversely, to yield an estimate of the BH parameters (a, i) (Himwich et al. 2020), while also providing a probe of the plasma and magnetic field around the BH (Event Horizon Telescope Collaboration 2021b). Lastly, the photon ring substructure predicts a characteristic pattern of autocorrelations whose detection could also provide a window into the Lyapunov exponents γ (and other Kerr critical exponents τ and δ), even if the rings themselves remain unresolved (Hadar et al. 2021). A related promising idea involves timing signatures such as quasi-periodic oscillations (Johannsen & Psaltis 2010).

Finally, a wide variety of other techniques for testing the Kerr hypothesis and strong-field GR via VLBI observations have been proposed (e.g., Psaltis 2019), many of which do not rely upon the photon ring. Similarly, some constraints on BH inclination can be derived from other image features, such as its brightness asymmetry (Medeiros et al. 2022).

7.4. Open questions and future directions

We have thoroughly investigated the shape of the n = 2 photon ring and its corresponding interferometric signature in a certain class of equatorial models. We can now confirm the prediction that the ring should follow a circlipse shape at low-to-moderate inclinations, but have yet to obtain a precise characterization of its shape at high inclinations, where its image can display kinks or discontinuities (Fig. 7) that produce characteristic jumps in its signature in the visibility domain (Fig. 6). Understanding these patterns is crucial for devising a test of the Kerr hypothesis based on the photon ring shape that remains valid at high inclinations.

Reversing the logic, one can assume the Kerr nature of the source and try to infer its parameters from the photon ring shape. We hope to further explore the question of spin inference by solving the inverse problem described in Sect. 6.3. This approach could most likely also be used to infer the BH mass with much greater precision than the EHT has so far achieved without being able to probe higher spatial frequencies.

As discussed in Sect. 7.1, it is important to determine whether our conclusions generalize to models with geometrically thick emission, with noise created by inevitable source fluctuations, with different astrophysical conditions (e.g., a radially infalling flow rather than a circularized one), and also with tilted disks that are misaligned with the BH spin. For instance, our observation that the measured diameters (d₊, d₋) of the n = 2 ring lie on a straight line connecting the inner and outer edges of the n = 2 lensing band could well depend sensitively on any one of these assumptions, with direct consequences for estimates of BH spin.

Finally, the most pressing task is to revisit all these questions in the context of the first (n = 1) subring, which is likely to be observed first. What is the predicted shape of the n = 1 ring? How well can the BH parameters be recovered from its shape?

From a theoretical standpoint, these questions are harder to tackle because the n = 1 ring is likely too thick to display a universal regime (Eq. (20)); that is, its width-to-diameter ratio w/d is likely not small enough for the separation of scales (Eq. (20)) to arise. If so, then there is no baseline regime in which its interferometric signature follows the universal form (Eq. (22)), and it may well be necessary to compute the subleading correction in w/d to this visibility amplitude. We take a first stab at these issues in Appendix C.

8. Summary

Images of a BH are typically dominated by a bright, narrow ring encircling a central brightness depression. This “photon ring” decomposes into a sequence of self-similar subrings indexed by half-orbit number n, with each subring consisting of a full, lensed image of the main emission surrounding the BH (Sect. 2). Here, we analyzed a method proposed by Gralla et al. (2020) to test the Kerr hypothesis via an interferometric measurement of the n = 2 ring shape. To check the viability of this shape test, we applied it to a wide range of simulated images of a Kerr BH surrounded by a thin, equatorial disk (Sect. 3), varying both the emission profile (Sect. 5) and the BH spin and inclination (Sect. 6). These two parameter surveys confirm the robustness of the test at low-to-moderate inclinations i ≲ 45°, indicating that M87* (which likely lies at an inclination i ≈ 17° from Earth) is a promising target for such a test via space-VLBI observation.

At higher inclinations i ≳ 45°, our study uncovers the emergence of qualitatively new phenomena that complicate ring shape measurements (Sect. 4). These effects arise because the ring width develops significant angular variation, which impacts the range of baselines dominated by the ring’s interferometric signature. In particular, we find that observations made at fixed baseline length in the visibility domain cannot measure the shape of a single ring at every angle: such observations are sensitive to a fixed ring width and must therefore either jump between the n = 1 and n = 2 rings (at moderately high inclinations; Fig. 6), or else pick up the signature of a “hybrid n = 1&2 ring” of roughly uniform width (at very high inclinations; Fig. 7). A more refined, sophisticated shape test has yet to be devised for these situations.

Our investigation also reveals that for a given choice of BH spin and inclination, the observable photon ring shape depends noticeably on the astrophysical properties of the emitting source. In particular, the most directly accessible (n = 1 and n = 2) rings can differ significantly from the theoretical critical curve, and in many models they are in fact closer to the critical curves of neighboring BH spins and inclinations (Figs. 10 and 11). As such (and as discussed throughout Sects. 5 and 6), the lensing bands introduced in Sect. 2.5 are more relevant than the critical curve to the problem of BH spin inference from the photon ring shape. We have outlined a promising approach for tackling this problem in Sect. 6.3, with some encouraging preliminary results (Fig. 12).

---

Acknowledgments

All of the images in this paper were ray traced on the tycho computer cluster of the Paris Observatory. We thank Samuel Gralla, Alejandro Càrdenas-Avendaño, Andrei Lobanov, and Hengrui Zhu for helpful comments and discussions. A.L. also gratefully acknowledges support from Will and Kacie Snellings, and M.W. thanks Alexandra Elbakyan for her contributions to the open science initiative.

References

Appendix A: Computing Kerr lensing bands

We consider a null geodesic connecting a source at (t_s, r_s, θ_s, φ_s) to an observer at (t_o, r_o, θ_o, φ_o). Its behavior can be characterized by its conserved quantities (λ, η) [Eq. (2) above] via the radial and angular potentials (e.g., Gralla & Lupsasca 2020b)

$$\begin{array}{cc}\hfill \mathcal{R}(r)& ={(r^2+a_{\ast }^2-a_{\ast }\lambda )}^2-\mathrm{\Delta }(r)[\eta +{(\lambda -a_{\ast })}^2],\hfill \end{array}$$(A.1a)

$$\begin{array}{cc}\hfill \mathrm{\Theta }(\theta )& =\eta +a_{\ast }^2{cos}^2\theta -{\lambda }^2{cot}^2\theta ,\hfill \end{array}$$(A.1b)

whose zeros define the radial and angular turning points of the ray. Integrating the geodesic equation along the trajectory gives

$$\begin{array}{c}\hfill I_r{\displaystyle {\equiv }_{r_s}^{r_o}\frac{\mathrm{d}r}{{\pm }_r\sqrt{\mathcal{R}(r)}}{=}_{{\theta }_s}^{{\theta }_o}\frac{\mathrm{d}\theta }{{\pm }_{\theta }\sqrt{\mathrm{\Theta }(\theta )}}\equiv G_{\theta },}\end{array}$$(A.2)

where the signs ±_r and ±_θ switch at radial and angular turning points, respectively. These expressions keep track of the total Mino time τ = I_r = G_θ elapsed along the trajectory.

We consider the case of a distant observer (r_o = D ≫ M) at an inclination θ_o = i from the BH spin axis who shoots photons backwards into the geometry from a polar position (ρ, φ_ρ) in the sky. Such a photon has conserved quantities

$$\begin{array}{c}\hfill \lambda =-D\rho cos{\phi }_{\rho }sini,\eta ={(D\rho )}^2-{\lambda }^2-a_{\ast }^2{cos}^{2}i,\end{array}$$(A.3)

in terms of which the four roots r_j of the quartic radial potential ℛ(r) are given by [Eqs. (79)–(95) of Gralla & Lupsasca (2020b)]:

$$\begin{array}{cc}\hfill r_j& =-z-{ϵ}_j\sqrt{-\frac{\mathcal{A}}{2}-z^2+{\stackrel{\sim }{ϵ}}_j\frac{\mathcal{B}}{4z}},\hfill \\ \hfill {ϵ}_j& =\{\begin{array}{cc}1\hfill & \mathrm{if}j\in \{1,3\},\hfill \\ -1\hfill & \mathrm{if}j\in \{2,4\},\hfill \end{array}{\stackrel{\sim }{ϵ}}_j=\{\begin{array}{cc}1\hfill & \mathrm{if}j\in \{1,2\},\hfill \\ -1\hfill & \mathrm{if}j\in \{3,4\},\hfill \end{array}\hfill \\ \hfill \mathcal{A}& =a_{\ast }^2-\eta -{\lambda }^2,\mathcal{B}=2M[\eta +{(\lambda -a_{\ast })}^2],\mathcal{C}=-a_{\ast }^2\eta ,\hfill \\ \hfill \mathcal{P}& =-\frac{{\mathcal{A}}^2}{12}-\mathcal{C},\mathcal{Q}=-\frac{\mathcal{A}}{3}[{\left(\frac{\mathcal{A}}{6}\right)}^2-\mathcal{C}]-\frac{{\mathcal{B}}^2}{8},\hfill \\ \hfill {\omega }_{\pm }& ={(-\frac{\mathcal{Q}}{2}\pm \sqrt{\frac{{\mathcal{P}}^3}{27}+\frac{{\mathcal{Q}}^2}{4}})}^{1/3},z=\sqrt{\frac{{\omega }_{+}+{\omega }_{-}-\mathcal{A}/3}{2}}.\hfill \end{array}$$

Letting r_ij = r_i − r_j, we define $A=\sqrt{r_{32}{r}_{42}}$, $B=\sqrt{r_{31}{r}_{41}}$, and introduce the antiderivatives of I_r (Gralla & Lupsasca 2020a,b):

$$\begin{array}{cc}\hfill {\mathcal{I}}_r^{+}(r)& =\frac{2}{\sqrt{r_{31}{r}_{42}}}F\left(arcsin\sqrt{\frac{r-r_4}{r-r_3}\frac{r_{31}}{r_{41}}}|\frac{r_{32}{r}_{41}}{r_{31}{r}_{42}}\right),\hfill \end{array}$$(A.4)

$$\begin{array}{cc}\hfill {\mathcal{I}}_r^{-}(r)& =\frac{1}{\sqrt{\mathit{AB}}}F\left(arccos\left(\frac{1-\frac{r-r_2}{r-r_1}\frac{B}{A}}{1+\frac{r-r_2}{r-r_1}\frac{B}{A}}\right)|\frac{{(A+B)}^2-r_{21}^2}{4AB}\right).\hfill \end{array}$$(A.5)

• If $\rho >\stackrel{\sim }{\rho }({\phi }_{\rho })$, so that the photon is shot from a position outside the critical curve, then all its roots are real and ordered as r₄ > r₃ > r₊ > r₋ > r₂ > r₁. As such, it must encounter a radial turning point at r = r₄ and eventually be deflected back to infinity. The Mino time τ(r) elapsed along its trajectory is

  $$\begin{array}{c}\hfill I_r^{{\mathcal{C}}^{+}}(r)={\mathcal{I}}_r^{+}(r_o)\mp {\mathcal{I}}_r^{+}(r),\end{array}$$(A.6)

  with upper/lower sign before/after reaching the turning point. The total Mino time elapsed along the full light ray is

  $$\begin{array}{c}\hfill I_r^{\mathrm{total}}=\frac{4}{\sqrt{r_{31}{r}_{42}}}F\left(arcsin\sqrt{\frac{r_{31}}{r_{41}}}|\frac{r_{32}{r}_{41}}{r_{31}{r}_{42}}\right).\end{array}$$(A.7)

• If $\rho <\stackrel{\sim }{\rho }({\phi }_{\rho })$, so that the photon is shot from a position inside the critical curve, then it cannot encounter a radial turning point and must therefore asymptotically reach the horizon at r = r₊. The Mino time τ(r) elapsed along its trajectory is

  $$\begin{array}{c}\hfill I_r^{{\mathcal{C}}^{-}}(r)={\mathcal{I}}_r^{\pm }(r_o)-{\mathcal{I}}_r^{\pm }(r),\end{array}$$(A.8)

  with the upper sign when all roots are real (in which case r₊ > r₋ > r₄ > r₃ > r₂ > r₁) and the lower sign when $r_3={\overline{r}}_4$ are complex conjugate roots with r₊ > r₋ > r₂ > r₁.⁹ The total Mino time elapsed along the full light ray is

  $$\begin{array}{c}\hfill I_r^{\mathrm{total}}=I_r^{{\mathcal{C}}^{-}}(r_{+}).\end{array}$$(A.9)

Regardless of where the photon is shot back from, it crosses the equatorial plane a maximum number N + 1 of times given by

$$\begin{array}{cc}\hfill N& =\lfloor \frac{I_r^{\mathrm{total}}\sqrt{-u_{-}{a}_{\ast }^2}+\mathrm{sign}(sin{\phi }_{\rho })F_o}{2K}\rfloor -H(sin{\phi }_{\rho }),\hfill \\ \hfill u_{\pm }& ={\mathrm{\Delta }}_{\theta }\pm \sqrt{{\mathrm{\Delta }}_{\theta }^2+\frac{\eta }{a_{\ast }^2}},{\mathrm{\Delta }}_{\theta }=\frac{1}{2}(1-\frac{\eta +{\lambda }^2}{a_{\ast }^2}),\hfill \\ \hfill K& =K\left(\frac{u_{+}}{u_{-}}\right),F_o=F\left(arcsin\left(\frac{cos(i)}{\sqrt{u_{+}}}\right)|\frac{u_{+}}{u_{-}}\right),\hfill \end{array}$$(A.10)

where H, F, and K respectively denote the Heaviside function, the incomplete elliptic integral of the first kind F(u|k), and its completion K(k)≡F(π/2|k). u_± are zeros u = cos²θ of the angular potential Θ(θ); when evaluated on critical geodesics with $(\lambda ,\eta )=(\stackrel{\sim }{\lambda },\stackrel{\sim }{\eta })$ [Eq. (6)], they reproduce ${\stackrel{\sim }{u}}_{\pm }$ in Eq. (5).

The photon intersects the equatorial plane for the (n + 1)^th time (for n ≤ N) when the Mino time τ(θ) elapsed along its trajectory is [Eq. (81) of Gralla & Lupsasca (2020a)]

$$\begin{array}{c}\hfill G_{\theta }(n)=\frac{2mK-\mathrm{sign}(sin{\phi }_{\rho })F_o}{\sqrt{-u_{-}{a}_{\ast }^2}},\end{array}$$(A.11)

where m = n + H(sin φ_ρ) is the number of angular turning points encountered along the trajectory. We note that N is by definition the non-negative integer such that $G_{\theta }(N)\le I_r^{\text{total}}\le G_{\theta }(N+1)$.

The n^th image of an equatorial source ring of constant radius r_s is the set of points (ρ, φ_ρ) in the observer sky that satisfy

$$\begin{array}{c}\hfill \boxed{{\mathrm{I}}_{\mathrm{r}}^{{\mathcal{C}}^{\pm }}({\mathrm{r}}_{\mathrm{s}})={\mathrm{G}}_{\theta }(\mathrm{n})}\end{array}$$(A.12)

with G_θ(n) given in Eq. (A.11) and $I_r^{{\mathcal{C}}^{\pm }}(r_s)$ given in Eqs. (A.7) and (A.9) for points outside and inside the critical curve. These points always trace a closed, continuous curve 𝒞_n(r_s).

The contours 𝒞₀(r_s) and 𝒞₁(r_s) corresponding to a selection of radii r_s are shown for various BH spins and inclinations in Fig. 6 of Gralla & Lupsasca (2020a). While all these curves are continuous, some of them may cross the critical curve, in which case they are stitched together from segments solving Eq. (A.12) with opposite signs. Aside from this minor inconvenience, it is straightforward to numerically solve Eq. (A.12) for any 𝒞_n(r_s).

In particular, as explained in Sec. 2.5, the n^th lensing band is the annular region in the sky bounded within the contours 𝒞_n(r₊) and 𝒞_n(∞), which respectively trace its inner and outer edges. These contours are relatively simple to compute because they can never cross the critical curve: trajectories that are deflected back to infinity must obey Eq. (A.6), so the outer edge 𝒞_n(∞) is obtained by solving Eq. (A.12) with a + sign only; likewise, trajectories that reach the BH must obey Eq. (A.8), so the inner edge 𝒞_n(r₊) is obtained by solving Eq. (A.12) with a − sign only.

We note that when crossing the inner and outer edges 𝒞_n(r₊) and 𝒞_n(∞) of the n^th lensing band, the maximal number N + 1 of equatorial crossings (Eq. (A.10)) must by definition change by 1. This provides a means of checking that 𝒞_n(r₊) and 𝒞_n(∞) have been correctly computed (or an alternative way to compute them).

Finally, we note that Eq. (A.2) may also be directly solved for r_s(I_r) using a single “unified inversion formula” given in Gralla & Lupsasca (2020b). Then $r_{\mathrm{eq}}^{(n)}(\rho ,{\phi }_{\rho })\equiv r_s(G_{\theta }(n))$ are transfer functions that give the radii at which light rays shot back from position (ρ, φ_ρ) cross the equatorial plane for the (n + 1)^th time. Contours of these functions are precisely the 𝒞_n(r_s) defined above. Hence, lensing bands may also be computed using these transfer functions, as described in App. A of Chael et al. (2021).

Appendix B: Computing the range of allowed phovals within a lensing band

Here, we provide a detailed solution to the following problem. Given a lensing band ℬ ⊊ ℝ², how do we compute the subregion 𝒜 of the (d₊, d₋) plane spanned by the phovals that fit within ℬ?

We glossed over this question after introducing it in Sec. 5.4, but its answer is crucial for obtaining the blue regions in Figs. 10, 11, and 12. At first glance, this problem seems straightforward: we could scan the 6D space of phovals (Eq. (53)) and ask, for each R = {R₀,R₁,R₂,X,χ,φ₀} in this parameter space 𝒫₆, whether the associated phoval fits inside the lensing band ℬ. The phovals that fit inside the lensing band would form a 6D region 𝒜₆ ⊊ 𝒫₆ whose 2D projection Π : 𝒫₆ → ℝ² via Eq. (52),

$$\begin{array}{c}\hfill \mathrm{\Pi }:\mathit{R}\mapsto (d_{+},d_{-})=2(R_0+R_1,R_0+R_2),\end{array}$$(B.1)

would give precisely the sought-after region in the (d₊, d₋) plane:

$$\begin{array}{c}\hfill \mathcal{A}=\mathrm{\Pi }({\mathcal{A}}_6).\end{array}$$(B.2)

In practice, however, the computational cost of this “brute force” approach is too high to achieve sufficient precision in reasonable time: the high-dimensionality of the parameter space 𝒫₆ makes any detailed scan too slow. Hence, another approach is needed.

The tractable method we came up with is quite general and consists of two parts. First, we recast the question of whether a given point (d₊, d₋) belongs in 𝒜 into a minimization problem, which is easier to solve. Second, we devised an efficient way to scan across regions of the (d₊, d₋) plane while solving this minimization problem. We now describe each part separately.

B.1. Part 1: Minimization method

Let ℬ ⊊ ℝ² denote an annular subregion of the plane with inner edge ℐ and outer edge 𝒪. We can define a separation function characterizing the distance from a closed curve 𝒞 to ℬ as

$$\begin{array}{c}\hfill {\delta }_{\mathcal{B}}(\mathcal{C})=\underset{z\in \mathcal{C}}{max}d(z,\mathcal{B}),\end{array}$$(B.3)

where d is the usual distance between a point and a set, namely:

$$\begin{array}{c}\hfill d(z,\mathcal{B})=\{\begin{array}{cc}{min}_{a\in \mathcal{I}}|z-a|\hfill & \mathrm{if}z\mathrm{lies}\mathrm{inside}\mathrm{of}\mathcal{I},\hfill \\ {min}_{a\in \mathcal{O}}|z-a|\hfill & \mathrm{if}z\mathrm{lies}\mathrm{outside}\mathrm{of}\mathcal{O},\hfill \\ 0\hfill & \mathrm{if}z\in \mathcal{B}.\hfill \end{array}\end{array}$$(B.4)

For us, 𝒞 will be a phoval shape [Eq. (53)] and ℬ a Kerr lensing band with inner edge ℐ = 𝒞_n(r₊) and outer edge 𝒪 = 𝒞_n(∞), as in App. A.

With this definition, δ_ℬ(𝒞) = 0 if and only if 𝒞 lies within the lensing band ℬ. Thus, to check if a point (d₊, d₋) corresponds to a phoval 𝒞 in the lensing band ℬ, we can minimize the function

$$\begin{array}{c}\hfill \mathrm{\Delta }:(d_{+},d_{-},\mathit{r})\mapsto {\delta }_{\mathcal{B}}\left(\mathcal{C}(R_0,\frac{d_{+}}{2}-R_0,\frac{d_{-}}{2}-R_0,X,\chi ,{\phi }_0)\right),\end{array}$$

where 𝒞(R) is the phoval with the projected position given by Eq. (53), over the 4D space 𝒫₄ of parameters r = (R₀,X,χ,φ₀). In other words, we compute m : ℝ² → ℝ defined as

$$\begin{array}{c}\hfill m:(d_{+},d_{-})\mapsto \underset{\mathit{r}\in {\mathcal{P}}_4}{min}\mathrm{\Delta }(d_{+},d_{-},\mathit{r}).\end{array}$$(B.5)

If m(d₊, d₋) vanishes (or stays below a tolerance threshold that accounts for numerical artifacts), then there exists a phoval with diameters (d₊, d₋) that lies within the lensing band, and thus,

$$\begin{array}{c}\hfill \mathcal{A}=\{(d_{+},d_{-})\in {\mathbb{R}}^2:m(d_{+},d_{-})=0\}.\end{array}$$(B.6)

We have thus recast the problem of determining whether a given (d₊, d₋) belongs to 𝒜 to a problem of minimization over 𝒫₄, which is only a 4D space. This completes part 1.

B.2. Part 2: Random walk

To compute 𝒜, it now suffices to define a grid of values (d₊, d₋) (typically, between the values found for ℐ and 𝒪) and compute m(d₊, d₋) for each point on the grid. Since the computational cost of the minimization method in part 1 is still quite high, we must still devise an efficient method to scan across a large grid, ideally without having to evaluate m for each point. This can be done using a random walk with a reflecting barrier:

1. We start from the parameters R ∈ 𝒜₆ ⊂ 𝒫₆ of some phoval that we know to lie in the lensing band: for instance, some intermediate phoval between those approximating ℐ and 𝒪.

2. We randomly pick one of the radii {R₀,R₁,R₂} and add to this parameter a predefined increment ±b, with the sign picked at random (we recall our convention that R₁ ≥ R₂; if this changes, then we must interchange R₁ ↔ R₂ to maintain it).

3. After the increment, we test if the resulting phoval (with the same values of X, χ, and φ₀) still lies within the lensing band. If that is the case, then we store (d₊, d₋) = 2(R₀ + R₁, R₀ + R₂) as an accepted value in the region 𝒜 and then proceed with the random walk by repeating step 2. If not, then we compute m(d₊, d₋), with two possible outcomes:

• m(d₊, d₋) = 0 (or the minimum is under threshold), in which case we store (d₊, d₋) = 2(R₀ + R₁, R₀ + R₂) as an accepted value in the region 𝒜 and return to step 2 to continue the random walk; otherwise,

• m(d₊, d₋) > 0 (or the minimum is above threshold), in which case we store (d₊, d₋) = 2(R₀ + R₁, R₀ + R₂) as a rejected value, then restore the values of {R₀,R₁,R₂} to their previous state and return to step 2 to take another step of the random walk in a new direction.

Steps 2 and 3 can be repeated any number of times, resulting in a set of rejected values {(d₊,d₋)}_rejected and a set of accepted values

$$\begin{array}{c}\hfill \mathcal{A}={\{(d_{+},d_{-})\}}_{\mathrm{accepted}}.\end{array}$$(B.7)

One advantage is that the minimization routine m is only evaluated in the rejected region of the (d₊, d₋) plane, or inside the allowed region but very near its boundary. This boundary acts as a reflecting barrier: each time it is crossed, the random walk returns to a previous position in 𝒜. Thus, it cannot escape far in the rejected region, where the computational cost is higher.

The boundaries of the region 𝒜 (which corresponds to the blue boxes in Figs. 10, 11, and 12) can be computed to arbitrary precision by tuning the value of the increment b that controls the step size in the random walk, and also by running several instances of the random walk starting from different positions.

Appendix C: Shape of the n = 1 photon ring

The n = 1 ring of M87* will almost certainly be the first photon subring to be measured (since its signal dominates the visibility on the shortest baselines to space), and it is therefore interesting to predict its geometric shape and interferometric signature. Does the n = 1 ring also follow the circlipse shape (Eq. (39))? If so, to what precision? Can its projected diameter $d_{\phi }^{(1)}$ be recovered from the periodicity of the visibility amplitude |V(u, φ)|?

Answering these questions in detail would be the subject of a whole other paper. Here, we merely record some observations that arose from our study and should inform future analyses:

1. Only narrow rings with a very small width-to-diameter ratio w/d ≪ 1 admit a universal regime (Eq. (20)). For thicker rings, there is not a large enough separation between the scales 1/d and 1/w for universality to set in. The n = 1 ring tends to be rather thick in most models (though not all, since there can be substantial variation in its width with the choice of astrophysical profile, as discussed below Eq. (46) and in Sec. 5), with a typical width of w₁ ≈ 1M and diameter d ≈ 10M. Hence w₁/d ≈ 10%, implying that subleading corrections to the universal visamp profile (Eq. (22)) may be substantial. Indeed, in most of our n = 1 visamps, we observed a much faster decay than the predicted u^−1/2 universal power-law fall-off.

2. Geometrically, this also means that the periodicity $d_{\phi }^{(1)}$ of the n = 1 visamp no longer admits a simple interpretation as the projected diameter of the n = 1 ring in the image, as its diameter is not sharply defined to better than w₁/d ≈ 10%.

3. Even in the absence of a universal regime 1/d ≪ u ≪ 1/w with simple interferometric signature (Eq. (22)), the n = 1 ring can still dominate the visamp on some baselines b₁ ≪ u ≪ b₂ where we may hope to measure the periodicity of |V(u, φ)|, and thus $d_{\phi }^{(1)}$, but only so long as the visamp displays enough periods within this baseline domain of width Δu = b₂ − b₁.

While these issues may complicate a measurement of the n = 1 ring and its interpretation, none of them seems insurmountable (see Fig. C.1), and we briefly address them in reverse order.

Regarding point 3: a baseline region of width Δu contains a number of visamp hops N = Δu/T, where T is the period of a hop, which for a ring is T = 1/d_φ. Hence, the number of visamp hops in the n = 1-dominated regime b₁ ≪ u ≪ b₂ is

$$\begin{array}{c}\hfill N_1(\phi )\approx [b_2(\phi )-b_1(\phi )]d_{\phi }.\end{array}$$(C.1)

Likewise, for higher-order rings n ≳ 2, using Eq. (45), we have

$$\begin{array}{cc}\hfill N_n(\phi )& \approx [b_{n+1}(\phi )-b_n(\phi )]d_{\phi }\approx e^{\gamma (\phi )}{N}_{n-1}(\phi ),\hfill \end{array}$$(C.2)

so that the number of visamp hops within the regime dominated by the n^th ring grows exponentially with n at a rate controlled by the image demagnification factor e^γ. Thus, the numbers of hops N_n(φ) for higher-order rings are predicted by Kerr lensing once the number of hops N₁(φ) for the first ring is known, but this latter number depends on astrophysical details. For instance, for the model in Fig. C.1, we have Δu = b₂ − b₁ ≈ 80 Gλ at φ = 135°. Since the photon ring of M87* has an expected diameter of

$$\begin{array}{c}\hfill d_{\phi }\approx 40\mu \mathrm{as}\approx \frac{1}{5\mathrm{G}\lambda },\end{array}$$(C.3)

where we manipulated units as in Eq. (37), it follows that

$$\begin{array}{c}\hfill N_1({135}^{°})\approx 16,\end{array}$$(C.4)

and indeed, we count 16 visamp nulls within the gray region in Fig. C.1. This is in principle a sufficient number of hops to obtain a good estimate of the periodicity, and hence of $d_{\phi }^{(1)}$, which in the example of Fig. C.1 follows a circlipse shape, albeit with larger error bars than for the n = 2 ring.

Regarding point 2: in the limit w/d → 0, a ring’s projected diameter d_φ becomes sharply defined in its image. Otherwise, it is only defined up to ambiguities of order w/d. If the n = 1 ring displays a visamp periodicity 1/$d_{\phi }^{(1)}$, which feature of the ring image does that correspond to? A natural guess would be the distance $d_{\phi }^{\text{peak}}$ between peaks of the ring’s intensity profile, but in simple models (such as the Gaussian ring below) one can prove that $d_{\phi }^{(1)}<d_{\phi }^{\text{peak}}$, with the difference vanishing as w/d → 0. This effect can be attributed to the everywhere-positive curvature of a ring, which implies that its Radon transform has its center of mass shifted inward relative to its peaks, with the corresponding 1D Fourier transform |V(u, φ)| displaying a slightly smaller $d_{\phi }^{(1)}$.

Regarding point 1: it is helpful to consider an axisymmetric ring with purely radial intensity profile I(ρ, φ) = I(ρ), such as

$$\begin{array}{c}\hfill I(\rho )=\frac{1}{2\pi w^2}{e}^{-\frac{d^2}{8w^2}}{I}_0\left(\frac{d\rho }{2w^2}\right)e^{-\frac{{\rho }^2}{2w^2}},\end{array}$$(C.5)

with I₀ denoting the 0^th modified Bessel function of the first kind. This is known as the Gaussian ring because $I_0(x)\stackrel{x\to \infty }{\approx }{(2\pi x)}^{-\frac{1}{2}}{e}^x$, so that as d → ∞ the radial profile approaches a Gaussian of width w centered as ρ = d/2, with normalization chosen such that the total flux is unity, V(0) = 1. The visibility is exactly

$$\begin{array}{c}\hfill V(u,\phi )=V(u)=J_0(\pi du)e^{-2{\pi }^2{w}^2{u}^2}.\end{array}$$(C.6)

The key point is the following. If the ring is very thin, then there is a wide range of baselines u ≪ 1/w in which the exponential is negligible and the visibility is simply V(u)≈J₀(πdu), which matches the universal prediction (Eq. (22)) as soon as du ≳ 1 ≫ wu. If the ring is thick, then we may never ignore that exponential term, yet we still expect to see the universal visamp (Eq. (22)), only multiplied by that exponential; in general, we expect |V(u, φ)| to take the form of Eq. (22) times e^−c(φ)up (with p = 2 for a Gaussian ring).

All Tables

All Figures