# Images and photon ring signatures of thick disks around black holes

Frederic H. Vincent, Samuel E. Gralla, Alexandru Lupsasca, Maciek Wielgus

► To cite this version:

Frederic H. Vincent, Samuel E. Gralla, Alexandru Lupsasca, Maciek Wielgus. Images and photon ring signatures of thick disks around black holes. Astronomy & Astrophysics - A&A, 2022, 667, pp.A170. 10.1051/0004-6361/202244339 . hal-03719961

HAL Id: hal-03719961

https://hal.science/hal-03719961v1

Submitted on 8 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Images and photon ring signatures of thick disks around black holes

F. H. Vincent<fnref n="1" />, S. E. Gralla<fnref n="2" />, A. Lupsasca<fnref n="3" />, and M. Wielgus<fnref n="4" />

Received 24 June 2022 / Accepted 26 July 2022

## ABSTRACT

**Context.** High-frequency very-long-baseline interferometry (VLBI) observations can now resolve the event-horizon-scale emission from sources in the immediate vicinity of nearby supermassive black holes. Future space-VLBI observations will access highly lensed features of black hole images – photon rings – that will provide particularly sharp probes of strong-field gravity.

**Aims.** Focusing on the particular case of the supermassive black hole M 87*, our goal is to explore a wide variety of accretion flows onto a Kerr black hole and to understand their corresponding images and visibilities. We are particularly interested in the visibility on baselines to space, which encodes the photon ring shape and whose measurement could provide a stringent test of the Kerr hypothesis.

**Methods.** We developed a fully analytical model of stationary, axisymmetric accretion flows with a variable disk thickness and a matter four-velocity that can smoothly interpolate between purely azimuthal rotation and purely radial infall. To determine the observational appearance of such flows, we numerically integrated the general-relativistic radiative transfer equation in the Kerr spacetime, taking care to include the effects of thermal synchrotron emission and absorption. We then Fourier transformed the resulting images and analyzed their visibility amplitudes along the directions parallel and orthogonal to the black hole spin projected on the observer sky.

**Results.** Our images generically display a wedding cake structure composed of discrete, narrow photon rings ($n = 1, 2, ...$) stacked on top of broader primary emission that surrounds a central brightness depression of model-dependent size. At 230 GHz, the $n = 1$ ring is always visible, but the $n = 2$ ring is sometimes suppressed due to absorption. At 345 GHz, the medium is optically thinner and the $n = 2$ ring displays clear signatures in both the image and visibility domains. We also examine the thermal synchrotron emissivity in the equatorial plane and show that it exhibits an exponential dependence on the radius for the preferred M 87* parameters.

**Conclusions.** The black hole shadow is a model-dependent phenomenon – even for diffuse, optically thin sources – and should not be regarded as a generic prediction of general relativity. Observations at 345 GHz are promising for future space-VLBI measurements of the photon ring shape, since at this frequency the signal of the $n = 2$ ring persists despite the disk thickness and nonzero absorption featured in our models. Future work is needed to investigate whether this conclusion holds in a larger variety of reasonable models.

**Key words.** gravitation – accretion, accretion disks – black hole physics – radiative transfer – relativistic processes

## 1. Introduction

In 2019, the Event Horizon Telescope collaboration released 1.3 mm interferometric observations of the supermassive black hole M 87* at the center of the galaxy Messier 87 (EHT; Event Horizon Telescope Collaboration 2019a, hereafter EHT L1), achieving an effective angular resolution comparable to the black hole size. These observations revealed the presence of a bright ring of approximately 40 µas in diameter that surrounds a much darker central region. While these basic image features are undisputed and have in fact been independently confirmed (Arras et al. 2022; Carilli & Thyagarajan 2022; Lockhart & Gralla 2022), several aspects of their theoretical interpretation remain open. One example is whether the dark area should be associated with the “black hole shadow” (Falcke et al. 2000), as originally proposed (EHT L1), or with the apparent position of the equatorial event horizon, as the data now seem to suggest (Chael et al. 2021).

A plausible range of observational appearances for M 87* is depicted in Fig. 1. Under the currently favored assumption of optically thin emission concentrated very near the event horizon

(e.g., Event Horizon Telescope Collaboration 2019d, hereafter EHT L5, 2021), the characteristic appearance ranges between two extremes: a narrow photon ring surrounding a dark region inside the critical curve (a “shadow”, see Falcke et al. 2000), and a series of discrete photons rings stacked on top of broader emission outside the apparent equatorial horizon (a “wedding cake”, see Gralla et al. 2019 Fig. 1, and Johnson et al. 2020 Fig. 3). Models of spherically symmetric, infalling matter lead to shadows, while models with equatorial, orbiting matter generate wedding cakes. General-relativistic magnetohydrodynamic (GRMHD) models generally favor the wedding cake over the shadow (Johnson et al. 2020; Chael et al. 2021), but the debate has yet to be settled (Bronzwaer & Falcke 2021).

The purpose of this paper is to more fully flesh out the model space between the shadow and wedding cake extremes. Having a broad range of models is crucial, not just for understanding the source, but also for accurately forecasting its potential for future observations. We are especially motivated by the promise of photon ring (orbiting light) measurements on long interferometric baselines (Johnson et al. 2020; Gralla 2020), which could provide a stringent test of the Kerr hypothesis (Gralla et al. 2020,

<footnote n="1">¹ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. de Paris, Sorbonne Paris-Cité, 5 Place Jules Janssen, 92195 Meudon, France
e-mail: frederic.vincent@obspm.fr</footnote>

<footnote n="2">² Department of Physics, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="3">³ Princeton Gravity Initiative, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="4">⁴ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>
hereafter GLM 20; Wielgus 2021). A recent proposal to measure the precise shape of the $n = 2$ photon ring (formed by light that executes a full orbit around the black hole before escaping to a detector) relied exclusively on a purely equatorial and perfectly absorption-free model of the emission (Gralla et al. 2020). It is imperative to check whether this exciting prospect survives the introduction of geometrical thickness, absorption, and other potential complications of the real astrophysical flow (see also Paugnat et al. 2022).

For this study, we constructed a set of semi-analytic models that bridge the gap between the shadow and wedding cake extremes, while also including a “realistic” level of absorption. For select density, temperature, and velocity profiles, we computed the synchrotron emission and absorption based on the assumption of a uniformly magnetized disk. Tuning parameters to match the total 230 GHz horizon-scale flux density reported by EHT, we performed radiative transport to determine the observational appearance at both 230 GHz and 345 GHz, which is the planned frequency of future observations (Doeleman et al. 2019; Event Horizon Telescope Collaboration 2019b, hereafter EHT L2). Finally, we Fourier transformed the images and checked whether the photon ring signatures are observable on moderate and long baselines. We considered both infalling and orbiting matter, as well as geometrically thin, thick, and fully spherical emission regions. We fixed the observer inclination to a value of 160°, which is deemed to be very likely for M 87* (Walker et al. 2018), while varying the black hole spin from near-zero to near-maximal spin.

Our results have both theoretical and practical implications. On the theory side, we confirm that the size of the central brightness depression is highly model-dependent, while the presence of discrete photon rings is generic (Gralla et al. 2019; Johnson et al. 2020). We also confirm that the classic black hole shadow (a sharp intensity drop inside the critical curve) is caused by extreme special-relativistic redshifting (Narayan et al. 2019; Gralla 2021), depending only indirectly on general-relativistic effects. Furthermore, we find that the intensity drop occurs only in fully spherical, infalling models, whereas the dark area is highly distinct from the critical curve even for very thick disks assumed to contain purely infalling matter. That is, far from being a generic prediction for optically thin flows (as claimed, e.g., in Psaltis 2019; Narayan et al. 2019), the appearance of a dark shadow filling the interior of the critical curve is in fact a remarkably fine-tuned phenomenon. This provides further evidence that the EHT observations should not be regarded as the black hole shadow in the strict sense given by Falcke et al. (2000).

On the practical side, our results indicate that the $n = 2$ photon ring is only marginally observable at 230 GHz. For strongly magnetized disks, we reproduced previous results (Johnson et al. 2020; Chael et al. 2021) showing a prominent $n = 2$ photon ring, but we found that in the weak-magnetization regime, the strength of the $n = 2$ signal is sensitive to the astrophysical details of the source, and especially to the black hole spin. In some cases, the $n = 2$ signal even vanishes entirely due to absorption. However, we found that at 345 GHz, the $n = 2$ signal returns to levels broadly consistent with previous estimates. This suggests that 345 GHz is an appropriate target for future space-VLBI observations of M 87*, since the signal at that frequency is robust to the astrophysical parameters we vary. We stress that these results only mark the beginning of a systematic study of the reasonable parameter space. In particular, our models exclude the possibilities of tilted disks, highly inclined observers, or jet-base emission.

The paper is organized as follows. We present our accretion disk model in Sect. 2 and display images of some accretion flows in Sect. 3, before discussing the resulting visibility amplitudes in Sect. 4 and summarizing our conclusions in Sect. 5.

## 2. Accretion flow model

In this paper, we considered accretion onto a black hole described by the Kerr metric in Boyer-Lindquist coordinates $(t, r, \theta, \phi)$. The dimensionless spin parameter is labeled $a$. We used natural units such that $G = c = 1$, and radii are thus expressed in units of the black hole mass $M$. We always use $\rho$ to denote the cylindrical radius $\rho = r \sin\theta$ (rather than a density) and $z = r \cos\theta$ to denote the cylindrical height above the equatorial plane $\theta = \pi/2$.

### 2.1. Disk geometry, physical quantities, and emission profile

Our goal is to develop a very general and fully analytical model of geometrically thick disks, extending the one in Vincent et al. (2021). We restricted our attention to axisymmetric disks, so there is no dependence on $\phi$ in any of our physical quantities.

We considered a population of thermal electrons that fills the spacetime outside the event horizon. The (fluid-frame) electron number density $n_e(\rho, z)$ is specified in the equatorial plane at the (cylindrical) horizon radius $\rho_H = r_H = 1 + \sqrt{1 - a^2}$,

$$n_{e;H} = n_e(\rho_H, z = 0), \qquad (1)$$

as is the electron temperature $T_e(\rho, z)$,

$$T_{e;H} = T_e(\rho_H, z = 0). \qquad (2)$$

We defined density and temperature profiles via a prescription similar to the one used in the analytical radiatively inefficient accretion flow (RIAF) models (see, e.g., Broderick et al. 2011):

$$n_e(r, z) = n_{e;H} \left(\frac{r}{r_H}\right)^{-2} \exp\left(-\frac{z^2}{2(\alpha\rho)^2}\right), \qquad (3a)$$

$$T_e(r) = T_{e;H} \left(\frac{r}{r_H}\right)^{-1}. \qquad (3b)$$

The dependence of the density profile on the cylindrical height $z$ is thus a simple Gaussian with standard deviation $s_z = \alpha\rho$, where $\alpha$ is a free parameter that sets the opening angle of the accretion disk. Indeed, if we define the surface of the disk to be a cone of height $s_z$ above the equatorial plane, then $\alpha$ corresponds to the tangent of its opening angle. As $\alpha \to 0$, we recovered equatorial models similar to those of GLM 20 and Paugnat et al. (2022), but here we kept $\alpha \neq 0$ to explore the effects of disk thickness, which is expected to be significant for realistic RIAF flows (Yuan & Narayan 2014).

The magnitude of the magnetic field, which is necessary to compute the synchrotron radiation, is prescribed by imposing a constant magnetization throughout the disk. That is, we fixed the ratio of the magnetic field and particle energy densities

$$\sigma = \frac{B^2/4\pi}{m_p c^2 n_e}, \qquad (4)$$

where $B$ is the magnetic field magnitude, $m_p$ the proton mass, and $c$ the speed of light (which we have restored here for clarity).

We included the effects of thermal synchrotron emission and absorption using the phenomenological expression derived by
<figure><img src="image_2.png" /><figcaption>Fig. 1. Intensity cuts of equatorial and spherical models. Under the assumption of optically thin emitting matter concentrated very near the horizon, the range of reasonable appearances for models of accretion onto a Kerr black hole can be bracketed by two extreme toy models: equatorial, orbiting matter (brown), which produces a “wedding cake” structure (Gralla et al. 2019; Johnson et al. 2020), and spherical, infalling matter (purple), which produces a “shadow” (Falcke et al. 2000). Here, we show a 230 GHz intensity cut parallel to the spin axis of M 87* (taken to be a rapidly spinning black hole with a = 0.94), including photons that orbit at most one full orbit around the black hole (up to n = 2 half-orbits), with each model normalized to its peak intensity. The brown and purple curves were ray traced from analytic models (see Sect. 3.4 for details), while the green curve is numerical data from Johnson et al. (2020) produced with the simulation pipeline described in Wong et al. (2022). The horizontal bars indicate the location of the Kerr critical curve and the apparent position of the equatorial event horizon.</figcaption></figure>

Leung et al. (2011). We presented this expression and discussed its accuracy in Appendix A. As shown in Appendix B, the radial profile of thermal synchrotron emission in our model is approximately

$$j_v(r) \propto \exp\left(-\zeta \frac{r}{r_H}\right), \qquad (5)$$

where $\zeta$ is a model-dependent number that takes a typical value near 3 for our models. However, we emphasize that Eq. (5) is used for interpretation only; our computations are performed with the more precise expression B.1 of Leung et al. (2011).

## 2.2. Dynamics

### 2.2.1. Orbiting motion

To compute images, we still need to specify the four-velocity of the emitting electrons in the accretion flow. We introduced a linear combination of an orbiting matter component and an infalling matter component that allows one to interpolate between the two extremes illustrated in Fig. 1.

First, we defined the orbiting component by an azimuthal four-velocity field with vanishing radial and polar components. As discussed in Appendix C.1, it is nontrivial to find a prescription that is simple, everywhere smooth and becomes Keplerian

at large distances from the black hole. We followed the prescription of Gold et al. (2020) and wrote the four-velocity 1-form as

$$\begin{gathered} u_{\mu}^{\text{circ}} \mathrm{d}x^{\mu} = -u_{t}^{\text{circ}} (-\mathrm{d}t + \ell \mathrm{d}\phi), \\ \ell = \frac{\rho^{3/2}}{1 + \rho}, \end{gathered} \qquad (6)$$

which is a special case of the general profile (Eq. (C.16)) described in Appendix C.1. This choice gives rise to a mild divergence $\Omega \sim \rho^{-1/2}$ at the poles, but is otherwise well-defined outside the horizon. The polar divergence is irrelevant for our models of thin and thick disks, as they have effectively no emission from the poles, and it has no noticeable effect even in the limit of fully spherical emission for infalling matter (Fig. 6 left).

Unit-normalization of the four-velocity (Eq. (6)) requires

$$-u_t^{\text{circ}} = \frac{1}{\sqrt{-(g^{tt} - 2g^{t\phi}\ell + g^{\phi\phi}\ell^2)}}. \qquad (7)$$

We note that this circular motion is not geodesic.

### 2.2.2. Infalling motion

Next, we defined the radially infalling matter component. Here, we simply assumed that the motion is geodesic, with no azimuthal angular momentum and an asymptotically vanishing
<table><caption>**Table 1.** Parameters of our model.</caption><thead><tr><th>Symbol</th><th>Value</th><th>Property</th></tr></thead><tbody><tr><td><i>M</i></td><td>6.2 × 10<sup>9</sup> <i>M</i><sub>☉</sub></td><td>Compact object mass</td></tr><tr><td><i>D</i></td><td>16.9 Mpc</td><td>Compact object distance</td></tr><tr><td><i>a</i></td><td>{0.01, 0.94}</td><td>BH spin parameter</td></tr><tr><td>α = tan θ<sub>op</sub></td><td>{0.1, 1}<br>circ/thin {1.5, 5}<br>rad/thin {5, 20}</td><td>Disk opening angle</td></tr><tr><td><i>n</i><sub>e;H</sub></td><td>circ/thick {0.7, 2}<br>rad/thick {2, 10}</td><td>Max density of electrons</td></tr><tr><td><i>T</i><sub>e;H</sub></td><td>10<sup>11</sup> K</td><td>Max electron temperature</td></tr><tr><td>σ</td><td>0.01</td><td>Magnetization</td></tr><tr><td><i>i</i></td><td>160°</td><td>Inclination angle</td></tr><tr><td>ν<sub>obs</sub></td><td>[230, 345] GHz</td><td>Observing frequency</td></tr><tr><td><i>f</i></td><td>100 µas</td><td>Field of view</td></tr><tr><td><i>N</i> × <i>N</i></td><td>20 000 × 20 000</td><td>Image resolution</td></tr></tbody></table>

**Notes.** The maximum electron number density varies across models in order to ensure an observed flux of ≈0.5 Jy at 230 GHz. It is expressed in units of 10<sup>6</sup> cm<sup>-3</sup>, with the first and second numbers referring to BH spins of *a* = 0.01 and *a* = 0.94, respectively. The models are labeled as follows: “circ” means “circular rotation” (Sect. 2.2.1), “rad” means “radial infall” (Sect. 2.2.2), “thin” refers to a disk with opening-angle parameter *α* = 0.1, and “thick” to a disk with *α* = 1.

velocity. As reviewed in Appendix C, the resulting four-velocity takes the form

$$u_{\text{rad}}^{\mu} \partial_{\mu} = u_{\text{rad}}^{t} \partial_{t} + u_{\text{rad}}^{r} \partial_{r} + u_{\text{rad}}^{\phi} \partial_{\phi}, \quad (8)$$

where

$$ \begin{aligned} u_{\text{rad}}^{t} &= -g^{tt}, \\ u_{\text{rad}}^{r} &= -\sqrt{(-1 - g^{tt}) g^{rr}}, \\ u_{\text{rad}}^{\phi} &= -g^{\phi t}. \end{aligned} \quad (9) $$

We note that the *φ* component is nonzero due to frame-dragging.

### 2.2.3. Combined motion

In this paper, we only considered purely circular motion or purely radial motion, but here we show that it is easy to combine them.

Indeed, one can linearly combine these motions to obtain the total four-velocity $u^{t} \partial_{\mu} = u^{t} \partial_{t} + u^{r} \partial_{r} + u^{\phi} \partial_{\phi}$. Following the notation of Pu et al. (2016), we introduced $\Omega = u^{\phi}/u^{t}$ and wrote

$$u^{r} = (1 - \beta_{r}) u_{\text{rad}}^{r}, \quad (10a)$$

$$\Omega = \Omega_{\text{circ}} + (1 - \beta_{\phi})(\Omega_{\text{rad}} - \Omega_{\text{circ}}), \quad (10b)$$

where $0 \le \beta_{r} \le 1$ and $0 \le \beta_{\phi} \le 1$ parametrize the superposition of the circular and radial components (Eqs. (6) and (8)). We note that only $u_{\text{rad}}^{r}$ is present in the first equation because our prescribed orbital four-velocity (Eq. (6)) has vanishing radial component. Unit-normalization fixes the time component of the four-velocity to

$$u^t = \sqrt{\frac{1 + g_{rr} (u^r)^2}{g_{tt} + 2g_{t\phi} \Omega + g_{\phi\phi} \Omega^2}}, \quad (11)$$

and its *φ* component is then given by $u^{\phi} = \Omega u^{t}$.

Again, we will not examine this combined motion here, so the parameters $\beta_{r}$ and $\beta_{\phi}$ will not be discussed in the remainder of this paper – this section is simply meant to highlight that our model allows for very general flow dynamics.

### 2.3. Parameter choices

In this paper, we are primarily interested in the effects of nonzero disk thickness, the motion of the emitting material, and the black hole spin. We thus considered extreme cases corresponding to a thin ($\alpha$ = 0.1) or thick disk ($\alpha$ = 1), purely azimuthal (Sect. 2.2.1) or purely radial (Sect. 2.2.2) motion, and low ($\alpha$ = 0.01) or high ($\alpha$ = 0.94) spin. The other source parameters are fixed to take likely – or at least reasonable – values for M 87*. We chose an inclination of $i$ = 160° and a mass-distance pair of $M$ = 6.2 × 10<sup>9</sup> M<sub>☉</sub> and $D$ = 16.9 Mpc (EHT L5; Event Horizon Telescope Collaboration 2019e). For the magnetization $\sigma$, we adopted a low value of 0.01, meaning that we considered a weakly magnetized disk<fnref n="1" />. The normalizations of the electron density and temperature are chosen such that the 230 GHz flux is on the order of 0.5 Jy for all configurations, in accord with the analysis (and assumptions) of the 2017 EHT data (Event Horizon Telescope Collaboration 2019c; EHT L5). The model is illustrated in Fig. 2, and all the values of its parameters are given in Table 1.

## 3. Image of the accretion flow

### 3.1. Ray tracing and image orders

We performed general-relativistic ray tracing in the Kerr spacetime with the model described in Sect. 2. We used the open-source code GYOTO (Vincent et al. 2011) to trace null geodesics backwards in time from a distant observer. The code integrates the radiative transfer equation along the null geodesics to evolve the specific intensity $I_{\nu}$ across the accretion flow. The output is an image, that is, a map of $I_{\nu}$. As discussed in Appendix D, the GYOTO integration parameters are chosen to ensure the highest possible numerical precision for the ray tracing and radiative transfer computations, while maintaining a reasonable computation time.

We will frequently differentiate between image orders. The 0th-order ($n$ = 0) primary image is defined as the image produced by selecting the part of each null geodesic that extends from the observer screen to its first angular turning point (following the ray backwards into the geometry). The 1st-order ($n$ = 1) image is produced by selecting the part of each null geodesic that extends between its first and second angular turning points, and likewise the 2nd-order ($n$ = 2) image arises by retaining the contributions from the portion of each geodesic between its second and third angular turning points; we did not consider higher image orders. These image orders are illustrated in Fig. 2, and Appendix E describes the technical details of their computation.

### 3.2. Resulting images

The images of our models for spins $a$ = 0.01 and $a$ = 0.94 are presented in Figs. 3 and 4, respectively. All of them share the same qualitative appearance, and they all display the three main features of black hole images:

– A central dark area, whose size depends on the astrophysical assumptions (e.g., Gralla et al. 2019; Chael et al. 2021).

<footnote n="1">¹ A commonly made distinction for black hole accretion flows is that between Standard and Normal Evolution (SANE) versus a Magnetically Arrested Disk (MAD). Prior investigations of high-order photon rings have focused on the MAD regime (Johnson et al. 2020; Chael et al. 2021), and we are able to qualitatively reproduce these results using strongly magnetized ($\sigma$ = 1) thin disks ($\alpha$ = 0.1). Throughout this paper, we adopted a weaker magnetization $\sigma$ = 0.01 corresponding to the SANE regime (see, e.g., Fig. 1 of Porth et al. 2019).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Density profiles of the “thick” ($\alpha$ = 1, left) and “thin” ($\alpha$ = 0.1, right) disk models we consider. The plots display the poloidal ($\rho$, $z$) plane, with all azimuthal angles $\phi$ projected to one single point in the plane. The filled black area is the black hole event horizon. The solid black lines delineate the region in which bound photon orbits exist (the “photon shell”), with the circular-equatorial orbits (prograde and retrograde) indicated by black dots. The red color scale encodes the log-scale profile of the electron number density, with a floor set at 1% of the maximum density (this floor is only applied to this figure for readability; it is not applied in our model). The two solid red contours correspond to a density of 1% and 10% of the maximum. The dashed red lines enclose the points located at a height less than $3s_c$ above the equatorial plane, where $s_c$ is the standard deviation of the Gaussian distribution controlling the electron number density above the equator (Eq. (3)). In the right panel, we also show a high-order null geodesic in green, with blue letters marking the distant observer’s screen ($S$) and the two $\theta$ turning points ($T_1$ and $T_2$) along the geodesic. Very sharp changes of direction appear at the $\theta$ turning points. These are due to the 2D projection in the ($\rho$, $z$) plane of the three spatial dimensions of the null geodesic. The geodesic represented in 3D space would look perfectly smooth. The $n$ = 0 part of the geodesic extends between $S$ and $T_1$, the $n$ = 1 between $T_1$ and $T_2$, and the $n$ = 2 between $T_2$ and the black hole. The yellow arrows highlight the portions of the null geodesic that are responsible for most of the $n$ = 0, $n$ = 1, and $n$ = 2 emission.</figcaption></figure>

– A bright, narrow ring produced by strongly lensed pho- tons that execute multiple orbits around the black hole on their way to the observer. This thin ring can be decomposed into a series of $n$ = 1, 2, ... photon rings (often collectively referred to as “the photon ring”), each of which is a lensed image of increasingly higher order $n$ of the accretion disk. These subrings, which are observable, exponentially converge to the theoretical critical curve, which is not (Bardeen 1973). The geometry of the critical curve is a pure function of the black hole mass and spin and of the observer inclination. However, the geometry and flux of observable photon rings still depend on the astrophysical assumptions (e.g., Paugnat et al. 2022).

– A thick annular region of primary emission (produced by $n$ = 0 photons travelling “straight” from the source to the observer, without orbiting around the black hole) that strongly depends on the astrophysical assumptions. In our images, it lies primarily within the thin photon ring.

For each of the models in Figs. 3 and 4, we present the full image (left column) and images of the $n$ = 1 and $n$ = 2 rings only (second and third columns). The right column displays intensity along horizontal and vertical (i.e., perpendicular and parallel to the spin axis) cuts of the full image, zoomed in around the $n$ = 2 contribution, which is also displayed by itself. Clearly, the details of the $n$ = 2 contribution are quite sensitive to the astrophysical assumptions. We now analyze the origin of these differences.

### 3.3. Intensity contribution from the $n$ = 2 ring

Let us consider first the low-spin images, whose $n$ = 2 contributions are shown in the rightmost column of Fig. 3. Two features

are noteworthy: (i) the left and right peaks are of approximately equal intensity; (ii) the peaks in the models of radially infalling matter have lower intensity overall. Feature (i) is a consequence of the weak frame-dragging at low spin, which implies that the left and right geodesics follow a very similar path through the flow (Fig. 5 top panel). The suppressed intensity of the radial case relative to the circular case arises from a difference in their redshift factors that can be attributed to a complicated interplay of various competing effects, which we now explore in detail.

To this end, we define a local increment of specific intensity loaded onto any given geodesic via radiative transfer. First, we introduce a local emissivity $j_v$ and absorptivity $\alpha_v$, as well as the Planck function $B_v$, which depends only on the local temperature and which equals (by virtue of Kirchhoff’s law) the ratio of emission to absorption. We also need the optical depth

$$ \tau_v = \int_{\text{optical path}} \alpha_v ds, \qquad (12) $$

where the integral is taken over the portion of the ray connecting the source to the observer, so that $\exp(-\tau_v)$ is the transmission. The local intensity increment then reads

$$ \delta I_v = \frac{j_v}{\alpha_v} [1 - \exp(-\alpha_v \delta s)] \times \exp(-\tau_v) \times g^3 \qquad (13a) $$

$$ = \underbrace{B_v [1 - \exp(-\alpha_v \delta s)]}_{\text{self-absorbed local emission}} \times \underbrace{\exp(-\tau_v)}_{\text{transmission}} \times \underbrace{g^3}_{\text{redshift}} \qquad (13b) $$

$$ \approx B_v \exp(-\tau_v) g^3, \qquad (13c) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Images and photon rings at 230 GHz. Three leftmost columns: inner 50 µas of the low-spin ($a = 0.01$) brightness temperature maps. The two left columns share the same linear color scale, which goes up to a brightness temperature of $3.7 \times 10^{10}$ K. The third column is in logarithmic scale, with the overall scale varying across panels to enable better visualization. The total specific flux of each image, as well as the maximum brightness temperature of each $n = 2$ image (in units of $10^9$ K), are indicated in yellow font. The white dashed curve in the left column shows the primary image of the equatorial event horizon. Rightmost column: Horizontal (red) and vertical (blue) cuts of the full and $n = 2$ brightness temperature profiles, centered around the $n = 2$ peaks regions. The temperature ratios $r$ between the two $n = 2$ peaks are provided.</figcaption></figure>

where $\delta s$ is the increment of length along the ray as measured by the emitter, and $g$ denotes the redshift factor

$$g = \frac{u_{\text{obs}} \cdot p_{\text{obs}}}{u_{\text{em}} \cdot p_{\text{em}}}. \qquad (14)$$

Here, $u_{\text{obs/em}}$ denotes the four-velocity of the observer or emitter and $p_{\text{obs/em}}$ the photon momentum at the observer or emitter. All the quantities appearing in Eq. (13) depend on the local value of the emission frequency, which itself depends on the observed redshift, which in turn depends on the flow dynamics.

To illustrate this, we show in the top row of Fig. 5 some of the geodesics that most contribute to the $n = 2$ peaks visible in the rightmost column of Fig. 3. For the circular-thin model, these

light rays collect much more intensity than the analogous light rays for the radial-thin model, whose contributions are greatly suppressed by the redshift factor. This explains point (ii) above.

To summarize, the heights of the $n = 2$ peaks appearing in the rightmost columns of Figs. 3 and 4 are controlled by three factors that contribute to the RHS of Eq. (13), namely:

1. The emission, which in our model depends mostly on radius, as it is highly concentrated near the equatorial plane (near
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Same as Fig. 3 for high spin ($a = 0.94$). The radial-thin and radial-thick horizontal cuts do not show the left peak of the $n = 2$ ring because it is too small to be visible (the associated temperature ratios are $r = 10^4$ and $r = 10^8$). The right column has a different scale compared to Fig. 3.</figcaption></figure>

the $n = 2$ equatorial crossing of each light ray); however, it also depends on frequency and therefore redshift (see 3. below);

2. The absorption and transmission, which behave similarly to the emission, with the important difference that they are non-local, as they depend on the full history of the geodesic as it travels from the $n = 2$ emission site to the observer;

3. The redshift, which has a complicated dependence on radius (the gravitational redshift is linked to the metric and blows up at the horizon), as well as on the relative orientation of the emitted photon and emitting particle velocity.

The interplay of these three ingredients strongly depends on the observing frequency, on the dynamics of the flow, on the local physical conditions at the emission site (density, magnetic field, temperature), as well as on the physical conditions along the complete path followed by the photons. It is therefore not a simple task to predict the properties of the $n = 2$ contribution, and simulations are certainly the only way to investigate them. Figure 5 provides the typical values of these three

quantities in the region of $n = 2$ emission for each of our models.

At high spin, the $n = 2$ peaks (shown in the rightmost column of Fig. 4) are very different from their nonrotating counterparts. They also display two striking features: (i) in stark contrast with the nonrotating case, the left and right peaks in the horizontal cuts (perpendicular to the spin vector, shown in red) have very different heights, indicating a large brightness asymmetry in the image; (ii) this brightness asymmetry ratio (the intensity ratio of the two horizontal peaks) is extremely high for the radial models.

These features can be understood by examining the bottom row of Fig. 5, which shows that the left and right geodesics in the high-spin case probe very different parts of the inner flow, resulting in a different collected intensity. This is closely related to the fact that as black hole spin increases, the critical curve is both distorted (from a circle) and displaced (relative to the primary radiation, see, e.g., Chan et al. 2013 Fig. 5 left panel): since the $n = 2$ photon ring shares these properties, the left and
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Intensity cuts for our models. We display cuts along the directions perpendicular to spin (red) and parallel to spin (blue) for spin $a = 0.01$ for the various thick-disk models used in this article (four central panels), as well as for a spherical model obtained from the thick-disk model in the limit of large disk thickness (top left panel). The size of the critical curve (thin line) and of the primary ($n = 0$) image of the equatorial horizon (thick line) are shown in green.</figcaption></figure>

### 3.4.1. Extreme models: spherical and equatorial

Figure 1 of the introduction presents the intensity cuts of two analytical models that are not based on the thick-disk model that we develop here. These two models are defined as follows:

- The spherical infall model consists of a spherically symmetric distribution of emitting matter located everywhere outside the event horizon. This matter is falling into the black hole with four-velocity defined in Sect. 2.2.2. The number density and temperature of the emitting fluid were chosen at the event horizon and evolve with radius following the power laws defined in Eq. (3) (i.e., $n_e(r) \propto r^{-2}$, $T_e(r) \propto r^{-1}$). The magnetic field was prescribed by demanding that the magnetization $\sigma$ defined in Eq. (4) remain constant, at the same value as in our thick-disk model, $\sigma = 0.01$. We considered thermal synchrotron emission and ignored self-absorption.

- The equatorial orbiting model is that of GLM 20 with an ad-hoc emission profile (their emission model 1). There is also no absorption considered in this model.

### 3.4.2. Canonical thick-disk models

We now discuss the intensity cuts of our four canonical thick-disk models, which are displayed in Fig. 6. We focus on the case of low spin, as the high-spin case is qualitatively similar.

The “spherical limit” in the top-left panel is the $\alpha \to \infty$ limit of our disk model with radially infalling matter. This large-thickness limit of our thick-disk model exactly coincides with the spherical infall model discussed in the previous subsection.

The only difference is that in the upper-left panel of Fig. 6, we take synchrotron self-absorption into account. This spherical limit is similar to the models used by Falcke et al. (2000) and Narayan et al. (2019), except that we self-consistently include absorption in our model of synchrotron emission.

The $\alpha \to 0$ limit of our thick-disk is numerically tricky, and the emission profile of the analytical equatorial model depicted in Fig. 1 differs from that obtained with our thick-disk models. For this reason, we do not show an “equatorial limit” intensity cut in Fig. 6.

Figure 6 shows a continuum of models spanning a large parameter space, from near-equatorial to spherical emission, and from circularly orbiting matter to radially infalling matter. We highlight a few important observations:

- The radial spherical model (top left) shows a prominent photon ring surrounding a large shadow that essentially occupies the full interior of the critical curve. The existence of the shadow is due to the strong redshift effect on the radiation emitted by matter that falls toward the black hole<fnref n="2" />. The spherical model is the only case where the intensity profile does not decompose into discrete layers ($n = 0, 1, 2$).

- Our thick-disk models produce a wide range of profiles that interpolate between the spherical limit and the equatorial model depicted in Fig. 1. Like the equatorial model, these models all produce intensity profiles that decompose into distinct $n = 0$ and

<footnote n="2">² Null geodesics arriving inside the critical curve have no radial turning points so the emitted radiation is directed toward increasing radii, while the emitter is heading toward decreasing radii, giving rise to strong redshift suppression of the observed intensity via the “headlight” effect.</footnote>
<table><caption>**Table 2.** Key features of the visibility amplitudes of the various models considered at $v_{\text{obs}} = 230$ GHz.</caption><thead><tr><th>Model</th><th>40 G&lambda;</th><th>100 G&lambda;</th><th>n = 1</th><th>n = 2</th></tr></thead><tbody><tr><td rowspan="2">a = 0.01 circular thin</td><td>9.2 mJy</td><td>2.4 mJy</td><td>(150 G&lambda;, 1 mJy)</td><td>(1400 G&lambda;, 25 &mu;Jy)</td></tr><tr><td>10 mJy</td><td>2.7 mJy</td><td>(140 G&lambda;, 2 mJy)</td><td>(1300 G&lambda;, 20 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.01 radial thin</td><td>14.4 mJy</td><td>2.2 mJy</td><td>(54 G&lambda;, 10 mJy)</td><td>(800 G&lambda;, 10 &mu;Jy)</td></tr><tr><td>12 mJy</td><td>3.3 mJy</td><td>(65 G&lambda;, 7.5 mJy)</td><td>(900 G&lambda;, 5.6 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.01 circular thick</td><td>8.6 mJy</td><td>2.4 mJy</td><td>(54 G&lambda;, 6.8 mJy)</td><td>(650 G&lambda;, 55 &mu;Jy)</td></tr><tr><td>7.7 mJy</td><td>2.8 mJy</td><td>(44 G&lambda;, 7.7 mJy)</td><td>(650 G&lambda;, 43.2 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.01 radial thick</td><td>12.2 mJy</td><td>2.2 mJy</td><td>(35 G&lambda;, 14.4 mJy)</td><td>(660 G&lambda;, 23 &mu;Jy)</td></tr><tr><td>11.2 mJy</td><td>2.1 mJy</td><td>(34 G&lambda;, 12.8 mJy)</td><td>(660 G&lambda;, 14.9 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.94 circular thin</td><td>12.0 mJy</td><td>2.0 mJy</td><td>(112 G&lambda;, 1.8 mJy)</td><td>(700 G&lambda;, 57 &mu;Jy)</td></tr><tr><td>11.2 mJy</td><td>1.8 mJy</td><td>(N/A, N/A)</td><td>(N/A, N/A)</td></tr><tr><td rowspan="2">a = 0.94 radial thin</td><td>10.4 mJy</td><td>1.3 mJy</td><td>(44.2 G&lambda;, 9.7 mJy)</td><td>(896 G&lambda;, 0.25 &mu;Jy)</td></tr><tr><td>7.1 mJy</td><td>1.4 mJy</td><td>(54 G&lambda;, 5.5 mJy)</td><td>(N/A, N/A)</td></tr><tr><td rowspan="2">a = 0.94 circular thick</td><td>10.7 mJy</td><td>1.8 mJy</td><td>(54 G&lambda;, 6 mJy)</td><td>(328 G&lambda;, 180 &mu;Jy)</td></tr><tr><td>7.3 mJy</td><td>0.7 mJy</td><td>(63.8 G&lambda;, 3.5 mJy)</td><td>(N/A, N/A)</td></tr><tr><td rowspan="2">a = 0.94 radial thick</td><td>7.4 mJy</td><td>0.6 mJy</td><td>(34.4 G&lambda;, 10.0 mJy)</td><td>(847 G&lambda;, 0.01 &mu;Jy)</td></tr><tr><td>6.2 mJy</td><td>0.5 mJy</td><td>(34.4 G&lambda;, 7.3 mJy)</td><td>(N/A, N/A)</td></tr></tbody></table>

**Notes.** For each model, we list the visibility amplitude $|V(u, \varphi)|$ at $u = 40$ G&lambda; and 100 G&lambda;, as well as the baseline thresholds $b_1(\varphi)$ and $b_2(\varphi)$ (and the corresponding visibility amplitudes) past which the $n = 1$ and $n = 2$ contributions start to dominate the signal (Eq. (15)). For each model, two lines are provided, corresponding to a baseline polar angle $\varphi = 90^\circ$ (aligned with spin, upper line) or $\varphi = 0^\circ$ (perpendicular to spin, lower line). The percent level $p$ of Eq. (15) is 5% everywhere, except for the numbers in red, where it is set to 10% (in these cases, the ring never dominates according to the 5% criterion). The “N/A” in the $n = 2$ column for the high-spin cases at $\varphi = 0^\circ$ mean that the $n = 2$ signal never dominates even for $p = 10\%$.

$n = 1$ contributions. The $n = 2$ photon ring is also clearly present for all the circular models; in the radial models, it is also present but (as discussed above) weaker and less easily discerned. The contrast between the photon rings and the $n = 0$ emission is stronger for the radial models because their $n = 0$ contribution is strongly suppressed (by the same redshift effect that produces a shadow in the spherical model). As a consequence, the central brightness depression is larger in the radial models than in the circular models; for the latter, as for the equatorial model, the central dark region is restricted to the interior of the apparent equatorial horizon.

## 4. Visibility signatures

One of the motivations for this paper is the prospect of an experimental detection of the $n = 2$ photon ring via its long-baseline “ringing” in the Fourier plane, which might be observable with future space-based VLBI missions (Johnson et al. 2020; Pesce et al. 2019; Gralla 2020; Gralla & Lupsasca 2020c; GLM 20). To check whether such a feature would be observable in our models, we computed the visibility along cuts in the Fourier plane in the directions parallel and perpendicular to the spin axis. Following GLM 20, we leveraged the projection-slice theorem to simplify the calculation. For a given orientation, we determined the line integrals on all lines perpendicular to the orientation (a 1D cut of the Radon transform of the image), and we Fourier transformed this 1D function to find the desired cut of the 2D Fourier transform, which is the radio visibility. The magnitude of this 1D complex visibility is the visibility amplitude, which we will simply refer to as the visibility. In the same way that we extracted from the full image its $n = 1$ and $n = 2$ components, we will compute the full Fourier transform and its $n = 1$

and $n = 2$ components. Our Fourier transform conventions agree with GLM 20.

Table 2 shows the main features of the visibilities associated to the various models considered in this study. In particular, it provides the visibility amplitude level at 40 and 100 G&lambda;, as well as the baseline thresholds (and corresponding visibility levels) past which the $n = 1$ and $n = 2$ rings start to dominate the signal. These thresholds were called $b_1$ and $b_2$ in Paugnat et al. (2022).

The following criterion is used to define when the $n$th ring starts to dominate the visibility. The full image decomposes into layers indexed by $n$, $I_v(x) = I_0(x) + I_1(x) + I_2(x) + \dots$, and the full visibility $V(u) = V_0(u) + V_1(u) + V_2(u) + \dots$ does too, by linearity of the Fourier transform. We considered a sliding baseline window $u \sim u_w$ of fixed width 10 G&lambda;, which is approximately twice the photon ring diameter. At each angle $\varphi$ in the baseline plane, we computed over this window $u \sim u_w$ the mean visibility $\langle|V(u, \varphi)|\rangle_w$ of the full image $I_v(x)$ and the mean visibility $\langle|V_n(u, \varphi)|\rangle_w$ of the $n$th ring image $I_n(x)$ only. We defined the baseline threshold $b_n(\varphi)$ past which the $n$th ring dominates the visibility to be the shortest baseline window $u_w \equiv b_n(\varphi)$ such that the percentage difference between these two means dips below a given threshold $p$:

$$ \frac{\langle|V_n(u, \varphi)|\rangle_w - \langle|V(u, \varphi)|\rangle_w}{\langle|V(u, \varphi)|\rangle_w} < p \Rightarrow \begin{array}{l} \text{n}^{th} \text{ ring dominates} \\ \text{on } u \ge b_n(\varphi) \equiv u_w. \end{array} \quad (15) $$

We will generally take $p = 5\%$, unless otherwise stated.

We again focus our comments on the $n = 2$ ring signature. In all low-spin models, the $n = 2$ ring clearly dominates the signal on sufficiently long baselines $u \ge b_2$, with the threshold $b_2$ ranging from ~600 to 1400 G&lambda;. The corresponding visibility amplitude displays clear oscillations at levels varying between ~10 and 50 &mu;Jy (see Table 2 and Fig. 7). Thus, these low-spin
Spin = 0.01

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Visibility amplitude profiles for all low-spin ($a = 0.01$) models. The full profile $|V(u, \varphi)|$ is shown in black, the $n = 1$ profile $|V_1(u, \varphi)|$ in green, and the $n = 2$ profile $|V_2(u, \varphi)|$ in red. The four upper panels correspond to an orientation $\varphi = 0^\circ$ (perpendicular to the spin axis, $\perp$), and the four lower panels to an orientation $\varphi = 90^\circ$ (parallel to the spin axis, $\parallel$). The accretion flow model is specified in the top-right corner of each panel.</figcaption></figure>
Spin = 0.94

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Same as Fig. 7 for the high-spin ($a = 0.94$) models. In all models, the $n = 2$ signal is very weak at $\varphi = 0^\circ$ (orientation perpendicular to the spin axis) due to the strong frame-dragging effect at high spin (see text for details). The vertical scale is different for the *two lower right panels*.</figcaption></figure>
<table><caption>**Table 3.** Same as Table 2, but at $v_{\text{obs}} = 345$ GHz.</caption><thead><tr><th>Model</th><th>40 G&lambda;</th><th>100 G&lambda;</th><th>n = 1</th><th>n = 2</th></tr></thead><tbody><tr><td rowspan="2">a = 0.94 circular thin</td><td>11.6 mJy</td><td>2.5 mJy</td><td>(112 G&lambda;, 1.8 mJy)</td><td>(455 G&lambda;, 176 &mu;Jy)</td></tr><tr><td>30.1 mJy</td><td>6.4 mJy</td><td>(83 G&lambda;, 8.6 mJy)</td><td>(593 G&lambda;, 150 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.94 radial thin</td><td>15.7 mJy</td><td>4.1 mJy</td><td>(34.4 G&lambda;, 18.5 mJy)</td><td>(377 G&lambda;, 209 &mu;Jy)</td></tr><tr><td>24.2 mJy</td><td>3.6 mJy</td><td>(34.4 G&lambda;, 26.3 mJy)</td><td>(397 G&lambda;, 110 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.94 circular thick</td><td>7 mJy</td><td>1.5 mJy</td><td>(45 G&lambda;, 5.4 mJy)</td><td>(384 G&lambda;, 140 &mu;Jy)</td></tr><tr><td>18.6 mJy</td><td>4.9 mJy</td><td>(45 G&lambda;, 18 mJy)</td><td>(254 G&lambda;, 372 &mu;Jy)</td></tr><tr><td rowspan="2">a = 0.94 radial thick</td><td>15.3 mJy</td><td>3.5 mJy</td><td>(35 G&lambda;, 17.4 mJy)</td><td>(284 G&lambda;, 227 &mu;Jy)</td></tr><tr><td>19.4 mJy</td><td>1.6 mJy</td><td>(35 G&lambda;, 23.2 mJy)</td><td>(374 G&lambda;, 100 &mu;Jy)</td></tr></tbody></table>

**Notes.** Only high-spin models are considered.

models match the idealized models considered in GLM 20 and Paugnat et al. (2022) reasonably well.

At high spin, by contrast, it is only on the baseline oriented parallel to the spin axis ($\varphi = 90^{\circ}$) that the $n = 2$ ring clearly dominates the signal (see the four lower panels of Fig. 8). However, even at this orientation, the models of radially infalling matter produce a very weak visibility amplitude in the $n = 2$-dominated regime $\approx 0.01-0.1$ $\mu$Jy. This is directly related to the low height of the $n = 2$ peaks in the rightmost column of Fig. 4: the blue profiles for the two models of radial infall display a contrast with the $n = 1$ radiation of $\approx 100-1000$.

On baselines perpendicular to the spin axis ($\varphi = 0^{\circ}$), all the high-spin models produce a very weak $n = 2$ signal, with essentially no oscillation at long baselines (see four upper panels of Fig. 8). This is the reason why some numbers are missing in Table 2: these entries correspond to cases in which the $n = 2$ oscillation is vanishingly small. To understand this effect, recall that the oscillations in visibility amplitude are caused by “interference” between opposite peaks in the image, with the amplitude of the oscillation set by the smaller of the two (Eq. (7) in Gralla 2020). Given that all high-spin, spin-perpendicular cuts display a large brightness asymmetry between the two $n = 2$ peaks (see the red profiles of the rightmost column in Fig. 4), the resulting visibility cuts exhibit essentially no oscillations associated with the $n = 2$ ring. This asymmetry between the $n = 2$ peaks is itself due to the interplay between absorption and redshift effects on the various parts of the $n = 2$ ring (see Sect. 3.3), as illustrated by the geodesics plotted in the bottom row of Fig. 5.

To summarize what we have discussed so far, the $n = 2$ peaks in the intensity cuts displayed in the rightmost column of Figs. 3 and 4 have two features of fundamental importance: (i) the ratio of flux in the $n = 2$ image $I_2(x)$ relative to the total flux in the full image $I_v(x)$ controls the average strength of the visibility in the $n = 2$-dominated regime; (ii) the brightness asymmetry (intensity ratio) between the left and right (or top and bottom) peaks in the $n = 2$ image cuts controls the amplitude of the visibility oscillations in the $n = 2$-dominated regime. These empirical observations are in perfect agreement with the analytical study of Gralla (2020). They allow one to guess whether a model will produce a strong $n = 2$ signal by simply computing a high-resolution intensity cut along a given orientation in the image plane, which is very cheap in terms of computational resources.

The fact that the $n = 2$ signal is absent from the spin-perpendicular visibility in the high-spin models is problematic from an experimental perspective. However, since absorption

is the main culprit in suppressing the $n = 2$ contribution, one might expect the signal to reappear at higher observation frequencies, where absorption is lower. While we have so far restricted our attention to the present-EHT observation frequency $v_{\text{obs}} = 230$ GHz, we can also consider a higher frequency $v_{\text{obs}} = 345$ GHz of planned future VLBI observations (Doeleman et al. 2019; EHT L2). Table 3 provides the same information as Table 2, for the high-spin models observed at $v_{\text{obs}} = 345$ GHz. This higher frequency contains a clear $n = 2$ signal in all cases, and the associated visibility level is in the regime $\sim 100$ $\mu$Jy that is considered potentially observable by GLM 20. The stronger $n = 2$ signal at higher frequency is also readily seen in the image cross-sections. For example, Fig. 9 displays intensity cuts for the circular-thin model at these two frequencies, and the higher-order peaks are clearly much sharper at 345 GHz.

## 5. Conclusion

In this article, we studied the images and visibility amplitude profiles of a variety of models of geometrically thin and thick disks, including the effects of absorption. We found that the size of the central brightness depression depends strongly on astrophysical assumptions, and that the photon ring is generically discrete rather than smooth, supporting the “wedding cake” heuristic over the “black hole shadow” paradigm. We also studied the long-baseline visibility amplitude, focusing on the signature of the $n = 2$ photon ring of M 87*. Although the $n = 2$ signal is always clear at low spin at 230 GHz, it can disappear at high spin due to a combination of frame-dragging and absorption along the line of sight. However, at 345 GHz, we find a clear $n = 2$ signal even at high spin, because the absorption is less prominent at this higher frequency.

We thus find that 345 GHz is a promising target for future space-VLBI observations of the photon ring of M 87*. However, we emphasize that this conclusion applies only for the particular astrophysical conditions that we have considered herein. Alternative reasonable choices of inclination, total compact flux, magnetization, density profile, etc., may very well lead to a different conclusion. We also cannot discount the possibility that the source geometry is qualitatively different, consisting for instance of a tilted disk (e.g., White et al. 2020), a jet-dominated profile (e.g., Kawashima et al. 2021), or emission from current sheets (e.g., Crinquand et al. 2022). Finally, we should keep in mind that the source itself may be variable in ways that we do not understand and cannot predict. For example, the total compact source flux density appears to have varied by a factor
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Intensity cuts dependence with observing frequency. We display here cuts along directions perpendicular (red) and parallel (blue) to spin for the circular-thin model at high spin, for two observing frequencies: 230 GHz and 345 GHz. The emission region is more optically thin at 345 GHz, so the higher-order peaks are more pronounced. In particular, the *n* = 2 left peak of the spin-perpendicular case is very sharp at 345 GHz, and fully absorbed at 230 GHz (see green ellipses).</figcaption></figure>

of 2 across a decade of VLBI observations, which would influ-
ence the density scale and hence the system optical depth
(Wielgus et al. 2020). It is therefore of paramount importance
to extend our analysis to a larger astrophysical parameter space
covering all reasonable possibilities for the source.

Acknowledgements. F.H.V. thanks F. Roy from Paris Observatory-LUTH for his help in compiling the GYOTO code on the MesoPSL cluster. We thank

Alejandro Cárdenas-Avendaño, Hadrien Paugnat, and George Wong for help-
ful comments and discussions. George Wong also provided the data for the green
curve in Fig. 1, and Alejandro Cárdenas-Avendaño computed the equatorial hori-
zon images in Figs. 3 and 4 using AART (Cárdenas-Avendaño et al., in prep.). This
work was supported in part by NSF grant PHY-1752809 to the University of
Arizona. It was also granted access to the HPC resources of MesoPSL financed
by the Région Île-de-France and the project Equip@Meso (reference ANR-10-
EQPX-29-01) of the Programme d'Investissements d'Avenir supervised by the
Agence Nationale pour la Recherche. Lastly, A.L. also gratefully acknowledges
support from Will and Kacie Snellings, while M.W. thanks Alexandra Elbakyan
for her contributions to the open science initiative.

References

Arras, P., Frank, P., Haim, P., et al. 2022, Nat. Astron., 6, 259

Bardeen, J. M. 1973, in Black Holes (Les Astres Occlus), eds. C. Dewitt, & B. S. Dewitt (Gordon and Breach Science Publishers), 215

Broderick, A. E., Fish, V. L., Doeleman, S. S., & Loeb, A. 2011, ApJ, 735, 110

Bronzwaer, T., & Falcke, H. 2021, ApJ, 920, 155

Carilli, C. L., & Thyagarajan, N. 2022, ApJ, 924, 125

Chael, A., Johnson, M. D., & Lupsasca, A. 2021, ApJ, 918, 6

Chan, C.-K., Psaltis, D., & Özel, F. 2013, ApJ, 777, 13

Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. A. 2022, Phys. Rev. Lett., accepted [arXiv:2202.04472]

Cunningham, C. T. 1975, ApJ, 202, 788

Doeleman, S., Blackburn, L., Dexter, J., et al. 2019, Bull. Am. Astron. Soc., 51, 256

Event Horizon Telescope Collaboration (Akiyama, K. et al.) 2019a, ApJ, 875, L1

Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L2

Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019c, ApJ, 875, L4

Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019d, ApJ, 875, L5

Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019e, ApJ, 875, L6

Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021, ApJ, 910, L13

Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13

Gold, R., Broderick, A. E., Younsi, Z., et al. 2020, ApJ, 897, 148

Gralla, S. E. 2020, Phys. Rev. D, 102, 044017

Gralla, S. E. 2021, Phys. Rev. D, 103, 024023

Gralla, S. E., & Lupsasca, A. 2020a, Phys. Rev. D, 101, 044031

Gralla, S. E., & Lupsasca, A. 2020b, Phys. Rev. D, 101, 044032

Gralla, S. E., & Lupsasca, A. 2020c, Phys. Rev. D, 102, 124003

Gralla, S. E., Holz, D. E., & Wald, R. M. 2019, Phys. Rev. D, 100, 024018

Gralla, S. E., Lupsasca, A., & Marrone, D. P. 2020, Phys. Rev. D, 102, 124004

Johnson, M. D., Lupsasca, A., Strominger, A., et al. 2020, Sci. Adv., 6, eaaz1310

Kawashima, T., Toma, K., Kino, M., et al. 2021, ApJ, 909, 168

Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21

Lockhart, W., & Gralla, S. E. 2022, MNRAS, 509, 3643

Narayan, R., Johnson, M. D., & Gammie, C. F. 2019, ApJ, 885, L33

Paugnat, H., Lupsasca, A., Vincent, F., & Wielgus, M. 2022, A&A, in press
[arXiv:2206.02781]

Pesce, D., Haworth, K., Melnick, G. J., et al. 2019, Bull. Am. Astron. Soc., 51, 176

Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26

Psaltis, D. 2019, Gen. Relativ. Gravit., 51, 137

Pu, H.-Y., Akiyama, K., & Asada, K. 2016, ApJ, 831, 4

Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics
(New York: Wiley)

Vincent, F. H., Paumard, T., Gourgoulhon, E., & Perrin, G. 2011, Class. Quant. Grav., 28, 225011

Vincent, F. H., Wielgus, M., Abramowicz, M. A., et al. 2021, A&A, 646, A37

Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128

White, C. J., Dexter, J., Blaes, O., & Quataert, E. 2020, ApJ, 894, 14

Wielgus, M. 2021, Phys. Rev. D, 104, 124058

Wielgus, M., Akiyama, K., Blackburn, L., et al. 2020, ApJ, 901, 67

Wong, G. N., Prather, B. S., Dhruv, V., et al. 2022, ApJS, 259, 64

Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529
# Appendix A: Relativistic thermal synchrotron emission and absorption

The synchrotron power per unit frequency emitted by a single ultrarelativistic electron is well-known (see Rybicki & Lightman 1979, Sec. 6.2) to be

$$ \frac{dE}{dt d\nu} = \sqrt{3} \frac{e^3 B \sin \theta}{mc^2} F\left(\frac{\nu}{\nu_{\text{crit}}}\right), \qquad (A.1) $$

where

$$ \nu_{\text{crit}} = \frac{3}{2} \gamma^2 \nu_{\text{cyclo}} \sin \theta, \qquad \nu_{\text{cyclo}} = \frac{eB}{2\pi mc}. \qquad (A.2) $$

Here, $\gamma$ denotes the Lorentz factor of the electron, $\theta$ the angle between the direction of emission and the magnetic field, and

$$ F(x) \equiv x \int_{x}^{\infty} K_{5/3}(u) du = \begin{cases} 2^{2/3} \Gamma(\frac{2}{3}) x^{1/3} & \text{for } x \ll 1, \\ \sqrt{\frac{\pi x}{2}} e^{-x} & \text{for } x \gg 1, \end{cases} \qquad (A.3) $$

with $K_{5/3}(u)$ a modified Bessel function of the second kind. We note that in these formulas, we have identified the pitch angle between the magnetic field and the direction of the electron's velocity with the emission angle $\theta$, which is a valid approximation for ultrarelativistic electrons (Rybicki & Lightman 1979). We also note that the derivation of Eq. (A.1) is essentially built on the strong beaming effect of the gyrating ultrarelativistic electrons, which is thus taken into account from the very start in this treatment.

The associated thermal synchrotron emissivity is

$$ j_v \equiv \frac{dE}{dt d\nu dV d\Omega} = \frac{\sqrt{3} e^3 B \sin \theta}{4\pi mc^2} \int_0^\infty F\left(\frac{\nu}{\nu_{\text{crit}}}\right) \frac{dn_e}{d\gamma} d\gamma, \qquad (A.4) $$

where

$$ \frac{dn_e}{d\gamma} = \frac{n_e \gamma (\gamma^2 - 1)^{1/2}}{\Theta_e K_2(\Theta_e^{-1})} \exp\left(-\frac{\gamma}{\Theta_e}\right) \qquad (A.5) $$

is the relativistic thermal (Maxwell-Jüttner) distribution, $K_2(u)$ is a modified Bessel function of the second kind, $\Theta_e = k_B T_e / mc^2$ is the dimensionless electron temperature, and the factor of $(4\pi)^{-1}$ implicitly assumes an isotropic distribution of the electrons' momenta. We note that the integral over $\gamma$ in Eq. (A.4) is evaluated over the range $[0, \infty)$, even though physically, $\gamma \ge 1$. Extending the range of integration to include the interval $[0, 1]$ simplifies the computation of the integral, while maintaining a good accuracy (as we will see below). Using the asymptotic expansions of $F(x)$ for $x \ll 1$ and $x \gg 1$, the emissivity (A.4) may be analytically expressed as (see, e.g., Leung et al. 2011)

$$ j_v \approx \begin{cases} \frac{2^{4/3} \pi n_e e^2 v_s}{3 c \Theta_e^2} X^{1/3} & \text{for } X \ll 1, \\ n_e \frac{\sqrt{2\pi} e^2 v_s}{6c\Theta_e^2} X e^{-X^{1/3}} & \text{for } X \gg 1, \end{cases} \qquad (A.6) $$

where

$$ X = \frac{\nu}{\nu_s}, \qquad \nu_s = \frac{2}{9} \nu_{\text{cyclo}} \Theta_e^2 \sin \theta. \qquad (A.7) $$

Leung et al. (2011) also provide a fitting function that bridges these two asymptotic regimes:

$$ j_v^{\text{approx}} = n_e \frac{\sqrt{2\pi} e^2 v_s}{3cK_2(\Theta_e^{-1})} (X^{1/2} + 2^{11/12} X^{1/6})^2 e^{-X^{1/3}}. \qquad (A.8) $$

To summarize, Eq. (A.8) is derived under the following list of approximations: we assume the electrons are ultrarelativistic and their population isotropic, we slightly extend the integration bounds on $\gamma$, and we match the asymptotic expansions of $F(x)$.

In our ray tracing code, we use the formula (A.8) for the emissivity, averaged over the emission direction $\theta$:

$$ \langle j_v^{\text{approx}} \rangle = \frac{1}{4\pi} \int j_v^{\text{approx}} d\Omega = \frac{1}{2} \int_0^{\pi} j_v^{\text{approx}} \sin \theta d\theta, \qquad (A.9) $$

where $d\Omega = \sin \theta d\theta d\phi$ and the emissivity is independent of the azimuthal angle $\phi$.

The exact expression for the emissivity, free from the above approximations, is (see Leung et al. 2011 for details)

$$ J_v^{\text{exact}} = \frac{2\pi e^2 \nu}{c |\cos \theta|} \sum_{n=n_{-}}^{\infty} \int_{\gamma_{-}(n)}^{\gamma_{+}(n)} \frac{1}{\beta} J_n(\gamma) \frac{dn_e}{d\gamma} d\gamma, \qquad (A.10) $$

where $n_{-} = \frac{\nu}{\nu_{\text{cyclo}}} |\sin \theta|$ and

$$ \gamma_{\pm}(n) = \frac{n_{\nu_{\text{cyclo}}} \pm |\cos \theta| \sqrt{(n_{\nu_{\text{cyclo}}})^2 - \sin^2 \theta}}{\sin^2 \theta}, \qquad (A.11) $$

where $n$ is the integer that is being summed over in Eq. (A.10). The dimensionless function $J_n(\gamma)$ is given in terms of the Bessel function of the first kind $J_n(z)$ and its derivative $J'_n(z)$ by

$$ J_n(\gamma) = [M J_n(z)]^2 + [N J'_n(z)]^2, \qquad (A.12a) $$

$$ M = \frac{\cos \theta - \beta \cos \alpha_0}{\sin \theta}, \qquad z = \frac{\nu \gamma \beta}{\nu_{\text{cyclo}}} \sin \theta \sin \alpha_0, \qquad (A.12b) $$

$$ N = \beta \sin \alpha, \qquad \cos \alpha_0 = \frac{1}{\beta \cos \theta} \left(1 - \frac{n \nu_{\text{cyclo}}}{\gamma} \right). \qquad (A.12c) $$

The sum over $n$ in Eq. (A.10) accounts for the helical motion of electrons along magnetic field lines, which makes each electron emit at every multiple of its Doppler-shifted gyrofrequency.

The accuracy of the approximation (A.8) can be assessed by comparing it with the exact expression (A.10). This is done in Figs. 10 and 11 of Leung et al. (2011), where the relative error between these two quantities at $\theta = 30^\circ$ and $80^\circ$ is displayed as a function of $\Theta_e$ and $\nu/\nu_{\text{cyclo}}$. This error is typically on the order of a few percent for $\Theta_e \in [1, 10]$ and $\nu/\nu_{\text{cyclo}} \in [10^3, 10^5]$, which includes the typical ranges for these parameters in the inner region of the accretion flow surrounding M87*. Indeed, the temperature of the source is $T_e \approx 10^{11}$ K, so $\Theta_e \approx 15$ (and the accuracy gets better with increasing $\Theta_e$), while the magnetic field is $B \approx 10$ G, which implies a cyclotron frequency $\nu_{\text{cyclo}} \approx 10^7$ Hz, so that $\nu/\nu_{\text{cyclo}} \approx 10^4$. On the other hand, the error typically exceeds 10% when $\Theta_e \le 1$, which means that the approximation (A.8) does not accurately model the emission from the outer region of the accretion flow. Nevertheless, this is not an issue for us because these contributions are negligible compared to the emission from the inner region of the flow.

The thermal absorption coefficient is simply related to the emissivity by Kirchhoff's law $\alpha_v = j_v / B_v$, where $B_v$ is the Planck function, so the approximations for the absorption are exactly the same as for the emission.
## Appendix B: Synchrotron emission profile for M87*

<figure><img src="image_10.png" /><figcaption>**Fig. B.1.** Validity of the approximated emissivity given by Eq. (B.3). This figure shows the ratio of $(X^{1/2} + 2^{11/12}X^{1/6})^2$ to X [or equivalently, of the emissivity (B.1) to its approximation (B.3)] as a function of magnetic field strength B and dimensionless temperature $Θ_e$, fixing $sin θ = 1$ and $ν = 230$ GHz (and neglecting redshift effects). The solid white curves are level sets of ratios of 1.5, 2 and 3. The dashed red curves are level sets of the X parameter corresponding to X = 10 and X = 100.</figcaption></figure>

### B.1. Thermal synchrotron

Since $Θ_e ≈ 15$ for M87*, we may assume that $Θ_e ≫ 1$ for this source. In that case, the approximation (A.8) for the emissivity can be further simplified using $K_2(x) \stackrel{x→0}{≈} 2x^{-2}$. This substitution results in the further approximation

$$j_v^{Θ_e \gg 1} \approx n_e \frac{\sqrt{2\pi e^2 v_{\text{cyclo}} \sin \theta}}{27c} (x^{1/2} + 2^{11/12}x^{1/6})^2 e^{-x^{1/3}}, \quad (B.1)$$

which is accurate to better than 10% provided that $Θ_e > 1.5$, and to better than 1% provided $Θ_e > 5$.

Plugging typical values for M87* yields an estimate for X of

$$X \approx \frac{1600}{\sin \theta} \left( \frac{\nu[\text{Hz}]}{10^{12}} \right) \left( \frac{10}{B[G]} \right) \left( \frac{10}{\bar{\Theta}_e} \right)^2, \quad (B.2)$$

where the frequency $ν$ and magnetic field strength B should be expressed in Hertz and Gauss, respectively. The factor of 1600 incorporates all the constant terms in the definition (A.7) of X.

Since this expression satisfies $X \gg 1$ for all typical values of M87* parameters, we can further simplify the approximation (B.1) for the emissivity by taking its $X \to \infty$ limit. In the regime $X \gg 1$, we may replace $(X^{1/2} + 2^{11/12}X^{1/6})^2$ by X to find

$$j_v^{Θ_e, X \gg 1} \approx n_e \frac{\sqrt{2\pi e^2 v_{\text{cyclo}} \sin \theta}}{27c} X e^{-X^{1/3}}, \quad (B.3)$$

which is a slight underestimation because $(X^{1/2} + 2^{11/12}X^{1/6})^2$ is always bigger than X. Their precise ratio is plotted as a function of B and $Θ_e$ in Fig. B.1 and is always of order ~ 1–10 for typical parameter values for M87*.

Next, since we are only interested in the overall profile of $j_v$, we may drop all constant prefactors and write

$$j_v \propto \frac{n_e \nu}{\Theta_e^2} e^{-x^{1/3}} = \frac{n_e \nu}{\Theta_e^2} \exp \left[ -C \left( \frac{\nu}{B \Theta_e^2} \right)^{1/3} \right], \quad (B.4)$$

where $C = (9\pi mc/e)^{1/3}$ is a constant.

Since we are assuming that the density and temperature fall off as $r^{-2}$ and $r^{-1}$, respectively, we may immediately conclude that the radial dependence of the emissivity is simply given by

$$j_v(r) \propto \exp[-X(r)^{1/3}], \quad (B.5)$$

provided that we neglect the radial dependence of the frequency, which is affected by a position-dependent redshift; this is only a reasonable assumption provided the emission is not too close to the horizon $r = r_+$, where the redshift diverges.

Returning to the estimate (B.2) for X, and using the assumed magnetic field fall-off $B \propto r^{-1}$ [Eq. (4)], we may now write

$$X \approx \frac{3.7 \times 10^5}{B_{\text{inner}} \Theta_{e;\text{inner}}^2 \sin \theta} \left( \frac{r}{r_{\text{inner}}} \right)^3, \quad (B.6)$$

where we have set $ν = 230$ GHz and once again neglected the frequency redshift. Here, $r_{inner}$ is the innermost radius of the disk, while $B_{inner}$ and $Θ_{e;inner}$ respectively denote the innermost values of the magnetic field and dimensionless temperature. Hence,

$$j_v(r) \propto \exp \left[ -\zeta \frac{r}{r_{\text{inner}}} \right], \quad \zeta = \left( \frac{3.7 \times 10^5}{B_{\text{inner}} \Theta_{e;\text{inner}}^2 \sin \theta} \right)^{1/3}. \quad (B.7)$$

To summarize, Eq. (B.7) is derived under the approximations listed below Eq. (A.8), together with the following additional assumptions: we assume that $Θ_e \gg 1$ and $X \gg 1$ (which puts us between the red contours of Fig. B.2), and we prescribe power-law fall-offs for the density, temperature, and magnetic field as described in our model ($n_e \propto r^{-2}$, $Θ_e \propto r^{-1}$, $B \propto r^{-1}$). Should the various power-law indices differ from this choice, one could derive a corresponding correction to Eq. (B.7) starting from the more general expression (B.4). We reiterate that this approximation underestimates the emissivity by a factor of a few, and that its quality increases with X (Fig. B.1).

The thermal-synchrotron emission profile for typical M87* parameters thus decays exponentially in the radius, with a slope $ζ$ depending on the exact conditions. Table B.1 lists the values of $ζ$ for our geometrically thin models, and Fig. B.2 displays $ζ$ as a function of the magnetic field strength and temperature.

For thermal synchrotron emission, the absorptivity may be obtained from Kirchhoff's law as

$$\alpha_v = \frac{j_v}{B_v(T_e)}, \quad (B.8)$$

where $B_v(T_e)$ is the Planck function at the local electron temperature $T_e$. At 230 GHz, photons are deep in the Rayleigh-Jeans regime where $hν \ll k_B T_e$, so that this relation simplifies to

$$\alpha_v \approx \frac{c^2}{2k_B T_e v^2} j_v. \quad (B.9)$$

Assuming as usual a temperature profile $T \propto r^{-1}$, we obtain

$$\alpha_v(r) \propto \frac{r}{r_{\text{inner}}} \exp \left[ -\zeta \frac{r}{r_{\text{inner}}} \right]. \quad (B.10)$$

The highly lensed emission originates from small radii $r \gtrsim r_{inner}$, for which the exponential profiles (B.7) and (B.10) provide good approximations to the emissivity and absorptivity, respectively.
<figure><img src="image_11.png" /><figcaption>**Fig. B.2.** Values of the parameter $\zeta$ [Eq. (B.7)] as a function of the magnetic field strength $B$ and dimensionless temperature $\Theta_e$. Here, we fix $\sin\theta = 1$. Our models always stay in the part of the plane located above the dashed horizontal red line ($\Theta_e > 5$), where the approximation (B.1) is good) and below the dashed red contour ($X > 10$, where the approximation (B.3) is excellent). The solid white contours are level sets of $\zeta$.</figcaption></figure>

<table><caption>**Table B.1.** Values of the parameter $\zeta$ in a few of our models. The emissivity profile is approximately $\exp(-\zeta r/r_{\text{inner}})$, provided that $r$ is far enough from the event horizon and redshift effects are neglected.</caption><thead><tr><td>Model</td><td>ζ</td></tr></thead><tbody><tr><td>a = 0.01 circ-thin</td><td>4.0</td></tr><tr><td>a = 0.01 rad-thin</td><td>3.4</td></tr><tr><td>a = 0.94 circ-thin</td><td>3.4</td></tr><tr><td>a = 0.94 rad-thin</td><td>2.7</td></tr></tbody></table>

## B.2. Power-law synchrotron

Although we only study thermal-synchrotron emission in this paper, it is interesting to also consider a nonthermal population of electrons (accelerated by magnetic reconnection, for instance) with synchrotron emission following a power law with index $p$. Still following Leung et al. (2011), the emissivity is then

$$j_v^{\text{PL}} = n_e \frac{e^2 v_{\text{cyclo}}}{c} \frac{3^{p/2}(p-1)\sin\theta}{2(p+1)(\gamma_{\text{min}}^{1-p} - \gamma_{\text{max}}^{1-p})} \times \Gamma\left(\frac{3p-1}{12}\right) \Gamma\left(\frac{3p+19}{12}\right) \left(\frac{\nu}{v_{\text{cyclo}}}\right)^{-(p-1)/2}, \quad (\text{B.11})$$

where $\gamma_{\text{max/min}}$ are the maximum or minimum values of the Lorentz factor for electrons in the distribution and $\Gamma(x)$ the gamma function. Following the same steps as in the previous section leads to

$$j_v^{\text{PL}} \propto \left( \frac{r}{r_{\text{inner}}} \right)^{-(p+5)/2}, \quad (\text{B.12})$$

that is, a power-law emissivity profile depending only on the power-law index of the electron distribution, and not on the physical conditions of the flow. For the particular value $p = 3$, $j_v^{\text{PL}} \propto r^{-2}$.

The power-law absorptivity $\alpha_v^{\text{PL}}$ is also provided by Leung et al. (2011) and also follows a power law that depends only on the power-law index of the electron distribution. For $p = 3$, it behaves as $\alpha_v^{\text{PL}} \propto r^{-9/2}$. Thus, it appears that the radial profile of the emissivity and absorptivity can vary significantly with the choice of electron distribution, and this could in turn have significant impact on the $n = 2$ ring signature (see Sec. 3.2). It would therefore be useful to also investigate this signature in models with nonthermal emission.

# Appendix C: Geometrically thick accretion flows

## C.1. Orbiting motion

In this section, we restricted ourselves to the Schwarzschild metric for simplicity. We want to prescribe a four-velocity for circularly orbiting matter in a thick disk. A circular four-velocity $u^\mu$ and its associated 1-form $u_\mu$ take the general form

$$u_{\text{circ}}^{\mu} \partial_{\mu} = u^{\prime} (\partial_{t} + \Omega \partial_{\phi}), \quad (\text{C.1a})$$

$$u_{\mu}^{\text{circ}} dx^{\mu} = -u_{t} (-dt + \ell d\phi), \quad (\text{C.1b})$$

where $\Omega = u^\phi/u^\prime$ denotes the angular velocity and $\ell = -u_\phi/u_t$ the specific angular momentum of the flow.

### C.1.1. Equatorial Keplerian motion

Timelike circular geodesic motion in the equatorial plane (i.e., Keplerian motion) is well-known (Bardeen 1973). The four-velocity takes the form (C.1) with

$$-u_t = \frac{r-2}{\sqrt{r(r-3)}}, \quad \ell = \frac{r^{3/2}}{r-2}, \quad \Omega = r^{-3/2}. \quad (\text{C.2})$$

It is obvious from these relations that Keplerian motion is not allowed inside the photon shell at $r = 3$. Moreover, Keplerian motion is only stable outside the innermost stable circular orbit (ISCO) located at $r_{\text{ISCO}} = 6$.

Below the ISCO, Cunningham (1975) prescribed that the geodesic constants of motion $-u_t$ and $u_\phi$ (the energy and spin angular momentum) keep their ISCO value, which ensures that the flow remains continuous across the ISCO as it spirals into the horizon. This prescription results in the four-velocity

$$u_{\mu}^{\text{C}^{75}} dx^{\mu} = u_{t}^{\text{ISCO}} dt + u_{r} dr + u_{\phi}^{\text{ISCO}} d\phi, \quad (\text{C.3})$$

where the $(t, \phi)$ components $u_{t,\phi}^{\text{ISCO}} \equiv u_{t,\phi}(r = r_{\text{ISCO}})$ take their Keplerian values from Eq. (C.2) evaluated at the ISCO radius, while the radial component is fixed by unit-normalization:

$$u_r^2 = g_{rr} \left[ -1 - g^{tt} (u_t^{\text{ISCO}})^2 - g^{\phi\phi} (u_\phi^{\text{ISCO}})^2 \right]. \quad (\text{C.4})$$

This defines a circular equatorial flow at every radius outside the event horizon $r_H$: its four-velocity is given by Eqs. (C.1)–(C.2) for $r \ge r_{\text{ISCO}}$, and by Eq. (C.3) for $r_H < r < r_{\text{ISCO}}$.

### C.1.2. Keplerian thick disk

We now wish to find a natural way of extending the equatorial Keplerian disk to a geometrically thick configuration.

For cylindrical radii outside the ISCO ($\rho > 6$), we simply defined the specific angular momentum $\ell$ to take the same value at every height $z$ above the equator as it does in the equatorial Keplerian disk model.
That is, for $\rho > 6$, we imposed the axially symmetric profile

$$ \ell(\rho, z) = \frac{\rho^{3/2}}{\rho - 2}, \qquad (C.5) $$

resulting in a circular four-velocity of the type (C.1),

$$ u_{\mu}^{\text{thick;}\text{ISCO}} dx^{\mu} = -u_t(\rho, z) \left( -dt + \frac{\rho^{3/2}}{\rho - 2} d\phi \right), \qquad (C.6) $$

whose unit-normalization fixes the time component to be

$$ -u_t(\rho, z) = \left( -g^{tt} - \frac{\rho}{(\rho - 2)^2} \right)^{-1/2}, \qquad (C.7) $$

where we used the fact that $g^{\phi\phi} = (r \sin\theta)^{-2} = \rho^{-2}$. We note that $g^{tt}$ depends on $z$ via the spherical radius $r = \sqrt{\rho^2 + z^2}$. It is easy to check that this quantity is well-defined for $\rho > 6$ (i.e., that the expression in parentheses remains positive).

For cylindrical radii inside the ISCO ($\rho < 6$), a natural extension of the Cunningham (1975) prescription is to require that, at any given height $z$, the components $u_t$ and $u_\phi$ keep their ISCO values for all radii $\rho < 6$. That is,

$$ u_{\mu}^{\text{thick;}\text{ISCO}} dx^{\mu} = u_{t}^{\text{ISCO}}(z) dt + u_{r}(\rho, z) dr + u_{\phi}^{\text{ISCO}}(z) d\phi, \qquad (C.8) $$

where the $(t, \phi)$ components $u_{t,\phi}^{\text{ISCO}}(z) \equiv u_{t,\phi}(\rho = \rho_{\text{ISCO}}, z)$ are obtained by evaluating Eqs. (C.6)–(C.7) at the ISCO radius $\rho_{\text{ISCO}} = 6$, while the radial component is again fixed by unit-normalization:

$$ [u_r(\rho, z)]^2 = g_{rr} \left( -1 - g^{tt} [u_t^{\text{ISCO}}(z)]^2 - g^{\phi\phi} [u_\phi^{\text{ISCO}}(z)]^2 \right). \qquad (C.9) $$

A simple numerical investigation of this relation reveals that it quickly turns negative as the height $z$ above the disk increases, so that this four-velocity is not well-defined. We must therefore find a new prescription that produces a well-defined flow everywhere.

### C.1.3. General thick disk

We will start from a general circular four-velocity of type (C.1),

$$ u_{\mu}^{\text{circ}}(\rho, z) dx^{\mu} = -u_t(\rho, z) (-dt + \ell(\rho) d\phi), \qquad (C.10) $$

which assumes that the specific angular momentum depends only on cylindrical radius. The function $\ell(\rho)$ is undetermined at this stage, but it is constrained by the requirement that this four-velocity remain well-defined everywhere outside the horizon.

Its angular velocity is

$$ \Omega \equiv \frac{u^{\phi}}{u^{t}} = \frac{g^{\phi\phi} u_{\phi}}{g^{tt} u_{t}} = \frac{\ell}{\rho^{2}} \left(1 - \frac{2}{r}\right), \qquad (C.11) $$

and unit-normalization fixes its time component to be

$$ -u_t = (-g^{tt} - g^{\phi\phi} \ell^2)^{-1/2} = \sqrt{\frac{1 - 2/r}{1 - \Omega \ell}}, \qquad (C.12) $$

which is well-defined outside the horizon ($r > 2$) provided that

$$ \Omega \ell < 1, \qquad (C.13) $$

or equivalently [by Eq. (C.11)],

$$ \frac{\ell^2}{\rho^2} \left(1 - \frac{2}{r}\right) < 1. \qquad (C.14) $$

Since $\rho = r \sin\theta$, we have $\rho - 2 \sin\theta > 0$ outside the horizon, where the condition (C.14) becomes

$$ |\ell(\rho)| < \frac{\rho^{3/2}}{\sqrt{\rho - 2 \sin \theta}}. \qquad (C.15) $$

We have thus far imposed no conditions on the function $\ell(\rho)$. At this point, we will assume that it converges to the Newtonian profile $\ell_{\text{Newton}} = \rho^{1/2}$ at spatial infinity. We also want to consider a family of profiles that (i) includes the Keplerian profile (C.2) as a special case and (ii) such that the constraint (C.15) leads to a simple condition. We are thus naturally led to consider the ansatz

$$ \ell_{\alpha}(\rho) = \frac{\rho^{3/2}}{\rho + \alpha}, \qquad (C.16) $$

where $\alpha$ is a real constant such that $\alpha = 2$ recovers the Keplerian profile (C.2). However, we want our profile to be defined at every spacetime point and to remain free of singularities. We will thus require that $\alpha > 0$, so that the denominator may never vanish. This also ensures that $\ell > 0$, and reduces condition (C.15) to

$$ (\rho + \alpha)^2 - \rho + 2 \sin \theta > 0. \qquad (C.17) $$

Demanding that the discriminant of this quadratic equation be strictly negative results us in a stronger condition

$$ \alpha > \frac{1}{4}, \qquad (C.18) $$

which is the final condition for our angular momentum profile to be well-defined everywhere outside the horizon. Following Gold et al. (2020), we consider $\alpha = 1$ in this article. We have numerically checked that this choice leads to a well-defined four-velocity even when the spin is nonzero.

Finally, we note that there is still a locus of spacetime where our chosen profile (C.16) produces singular behavior: the axis $\rho = 0$. Taking the limit $\rho \to 0$ keeping $z > 2$ fixed results in

$$ u^{t} \approx \left(1 - \frac{2}{r}\right)^{-1/2}, \qquad u^{\phi} \approx \rho^{-1/2} \left(1 - \frac{2}{r}\right)^{1/2}, \qquad (C.19) $$

so the $\phi$ component of the four-velocity, as well as the angular velocity $\Omega$, diverge in the limit $\rho \to 0$. Indeed, Eq. (C.11) shows that the angular velocity grows like $\Omega \approx \ell\rho^{-2} \approx \rho^{-1/2}$ near the axis for the Keplerian $\ell$, which behaves as $\ell \approx \rho^{3/2}$. Hence, any family $\ell(\rho)$ that includes the Keplerian profile as a member must have a divergent $\Omega$ near $\rho = 0$. However, since the orbiting models that we consider in this paper always have vanishingly small emission close to the axis, this pathology never causes problems.

### C.2. Infalling geodesic motion

In this section, we work in the Kerr metric and consider

$$ u_{\text{rad}}^{\mu} \partial_{\mu} = u^{t} \partial_{t} + u^{r} \partial_{r} + u^{\phi} \partial_{\phi}, \qquad (C.20a) $$

$$ u_t = -1, \qquad u_\phi = 0. \qquad (C.20b) $$

This four-velocity describes particles that fall into the black hole from spatial infinity, where they start out with vanishing velocity and angular momentum. The $(t, \phi)$ contravariant components are

$$ u^t = g^{tt} u_t = -g^{tt}, \qquad u^\phi = g^{t\phi} u_t = -g^{t\phi}. \qquad (C.21) $$
We note that $u^\phi$ must be nonzero when the black hole is rotating (frame-dragging prevents purely radial geodesic infall). Unit-normalization fixes the $r$ covariant component to be

$$u_r = -\sqrt{\frac{-1 - g^{tt}}{g^{rr}}}, \qquad (C.22)$$

where we picked the negative sign of the square root to describe infall, and then the $r$ contravariant component is

$$u^r = g^{rr} u_r = -\sqrt{(-1 - g^{tt}) g^{rr}}. \qquad (C.23)$$

These formulas are easier to interpret in the Schwarzschild case, where they simplify to

$$u^t = \left(1 - \frac{2}{r}\right)^{-1}, \qquad u^r = -\sqrt{\frac{2}{r}}, \qquad u^\phi = 0. \qquad (C.24)$$

This results in a purely radial velocity (since $u_\phi = 0$ at zero spin)

$$\mathbf{v} = v^r \mathbf{e}_r = \frac{dr}{dt} \mathbf{e}_r = \frac{u^r}{u^t} \mathbf{e}_r, \qquad (C.25)$$

where $\mathbf{e}_r$ the unit spacelike vector $\mathbf{e}_r = \partial_r / \sqrt{g_{rr}}$ and

$$v^r = -\sqrt{\frac{2}{r}} \left(1 - \frac{2}{r}\right) \stackrel{r \to \infty}{\approx} \sqrt{\frac{2}{r}}, \qquad (C.26)$$

which at spatial infinity reduces to the Newtonian formula for velocity derived from classical conservation of energy:

$$\frac{1}{2} (v^r)^2 - \frac{1}{r} = 0. \qquad (C.27)$$

## Appendix D: Precision of ray tracing computations

In this section, we examine the strongly lensed null geodesics in Kerr, that is, those that explore the shell of spherical photon orbits before escaping to produce high-order images. Their trajectories are very highly bent, as illustrated in Fig. D.1.

Kerr null geodesics obey (e.g., Gralla & Lupsasca 2020a,b)

$$I_r \equiv \int_{r_s}^{r_o} \frac{dr}{\pm r \sqrt{\mathcal{R}(r)}} = \int_{\theta_s}^{\theta_o} \frac{d\theta}{\pm \theta \sqrt{\Theta(\theta)}} \equiv G_{\theta}, \qquad (D.1)$$

where $\mathcal{R}(r)$ and $\Theta(\theta)$ are radial and angular geodesic potentials, and the integrals are evaluated along the entire trajectory with the signs $\pm_{r,\theta}$ flipping at turning points of the radial or angular motion, which correspond to zeros of their respective potentials.

We considered a geodesic traveling from a source at $(r_s, \theta_s)$ to an observer at $(r_o, \theta_o)$. To assess the precision of our geodesic integrator, we plot in Fig. D.2 the numerical evolution along the geodesic of the quantity $|I_r - G_{\theta}|$, as well as of the conserved energy $E = -p_t$ and azimuthal angular momentum $L = p_\phi$.

These quantities should remain constant along the geodesic, and the figure shows that their error is consistently below $10^{-13}$, except for a peak of $10^{-9}$ in $|I_r - G_{\theta}|$ close to an angular turning point. We have also checked that the radii of the three equatorial crossings match the analytical formula derived by Gralla & Lupsasca (2020b), which gives the radius of the equatorial crossing after $m$ angular turns in terms of the Jacobi elliptic sine $\operatorname{sn}(x|k)$:

$$r_{\text{eq}}^{(m)} = \frac{r_4 r_{31} - r_3 r_{41} \operatorname{sn}^2\left(\frac{1}{2} \sqrt{r_{31} r_{42}} \tau_{\text{eq}}^{(m)} - \mathcal{F}_o\right) k}{r_{31} - r_{41} \operatorname{sn}^2\left(\frac{1}{2} \sqrt{r_{31} r_{42}} \tau_{\text{eq}}^{(m)} - \mathcal{F}_o\right) k}. \qquad (D.2)$$

<figure><img src="image_12.png" /><figcaption>Fig. D.1. A null geodesic of order $n = 2$ (two $\theta$ turning points), in blue. The black disk is the event horizon and the yellow lines indicate the angular turning points $\theta_-$ and $\theta_+ = \pi - \theta_-$, which are computed from the analytical expression for $\theta_-$ as a function of black hole spin and the geodesic constants of motion (see, e.g., Gralla & Lupsasca 2020a). The three equatorial crossings of the geodesic are labeled by the black dots.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. D.2. Precision of the geodesic integration. The figure shows the evolution of the difference between the initial and current values of the geodesic constants $E = -p_t$ (blue) and $L = p_\phi$ (orange), as well as the quantity $|I_r - G_{\theta}|$ (green), as a function of the integration step. The peak of $|I_r - G_{\theta}|$ happens close to an angular turning point where the geodesic evolves in the Kerr spherical null geodesics region.</figcaption></figure>

Here, $r_{ij} = r_i - r_j$ with $\{r_1, r_2, r_3, r_4\}$ the four roots of the quartic potential $\mathcal{R}(r)$ (with exact expressions given in Eq. (A8) of Gralla & Lupsasca 2020a), and $\mathcal{F}_o$ is the elliptic integral of the first kind

$$\mathcal{F}_o = F \left( \arcsin \left| \frac{r_{31}}{r_{41}} \right| k \right), \qquad k = \frac{r_{32} r_{41}}{r_{31} r_{42}}. \qquad (D.3)$$

The Mino time elapsed up to the $m^{th}$ equatorial crossing is

$$\tau_{\text{eq}}^{(m)} = \frac{2mK - \operatorname{sign}(p_o^\theta) F_o}{a \sqrt{-u_{-}}}, \qquad (D.4)$$
where $p_o^\theta$ is the polar photon momentum at the observer,

$$u_{\pm} = \Delta_{\theta} \pm \sqrt{\Delta_{\theta}^{2} + \frac{\eta}{a^{2}}}, \qquad \Delta_{\theta} = \frac{1}{2} \left(1 - \frac{\eta + \lambda^{2}}{a^{2}}\right), \qquad (D.5)$$

are the zeros $u = \cos^2 \theta$ of $\Theta(u)$ in terms of the energy-rescaled angular momentum $\lambda = L/E$ and Carter constant $\eta = Q/E^2$, and

$$F_o = F\left(\arcsin\left(\frac{\cos\theta_o}{\sqrt{u_+}}\right)\bigg|\frac{u_+}{u_-}\right), \qquad K = F\left(\frac{\pi}{2}\bigg|\frac{u_+}{u_-}\right) \qquad (D.6)$$

are an elliptic integral of the first kind and its completion $K$.

We found that the three numerical and analytical values of $r_{eq}^{(m)}$ agree to within $10^{-4} - 10^{-5}M$. We note that the equatorial crossing labeled $\tau^{eq,2}$ in Fig. D.1 occurs after the backward ray traced geodesic has encountered two angular turning points; this shows that a very high level of numerical precision is maintained despite the spike in the error $|I_r - G_\theta|$ at that crossing. Finally, we have compared the value of the radial turning point $r_4$ (i.e., the minimum value of the radial trajectory) to the analytical expression given in Eq. (A8d) of Gralla & Lupsasca (2020a), and found an error $\lesssim 10^{-8}M$. All these results lead us to conclude that our null geodesic integrator is highly accurate, even within the photon shell of bound spherical photon orbits.

Besides the geodesic integration, $GYOTO$ also integrates the radiative transfer equation by discretizing the part of the geodesic that intersects the disk, using a constant step of size $\delta = 0.1M$. We have checked that reducing this step size by a factor of 10 changes the observed specific intensity by less than 0.1%, which is sufficient given that our synchrotron emissivity prescription is only at the percent level of precision (see App. A).

## Appendix E: Image orders

Null geodesics in Kerr can be conveniently parametrized by the so-called Mino time (Gralla & Lupsasca 2020a,b)

$$\tau = I_r = G_\theta, \qquad (E.1)$$

with $I_r$ and $G_\theta$ are defined in Eq. (D.1). Before the first angular turning point (i.e., between the points $S$ and $T_1$ in Fig. 2), the Mino time elapsed along the geodesic as it travels backwards from the observer at $\theta_o$ to $\theta = \theta_s$ is (Gralla & Lupsasca 2020b)

$$\tau(\theta) = \operatorname{sign}(p_s^\theta) (\mathcal{G}_\theta^o - \mathcal{G}_\theta^s), \qquad (E.2)$$

where $p_s^\theta$ is the polar photon momentum at $\theta = \theta_s$, while

$$\mathcal{G}_{\theta} = -\frac{1}{a\sqrt{-u_{-}}} F\left(\arcsin\left(\frac{\cos\theta}{\sqrt{u_{+}}}\right)\bigg|\frac{u_{+}}{u_{-}}\right), \qquad (E.3)$$

and $\mathcal{G}_\theta^o$ and $\mathcal{G}_\theta^s$ are evaluated at the $\mathcal{G}_\theta$ evaluated at $\theta = \theta_o$ and $\theta = \theta_s$, respectively. We can determine the Mino time $\tau_1$ elapsed at the first angular turning point by setting

$$\theta_s = \arccos(\pm\sqrt{u_+}), \qquad \pm = \operatorname{sign}(p_o^\theta), \qquad (E.4)$$

while the Mino time elapsed between two successive angular is

$$\Delta\tau = \frac{2K}{a\sqrt{-u_{-}}}. \qquad (E.5)$$

The order of any null geodesic is then defined as

* $n = 0$: between the observer and $\tau = \tau_1$,

* $n = 1$: between $\tau = \tau_1$ and $\tau = \tau_1 + \Delta\tau$,

* $n = 2$: between $\tau = \tau_1 + \Delta\tau$ and $\tau = \tau_1 + 2\Delta\tau$.

The integration is cut off past $n = 2$, so as not to pollute the image by unresolved $n > 2$ features.

We are thus able to provide not only the full image (containing all layers $n \le 2$), but also individual images of each layer $n \in \{0, 1, 2\}$ consisting of photons loaded onto each null geodesic at the corresponding order.