© F. H. Vincent et al. 2022

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).

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.

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, 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, θ, ϕ). 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 ρ to denote the cylindrical radius ρ = r sin θ (rather than a density) and z = r cos θ to denote the cylindrical height above the equatorial plane θ = π/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 ϕ 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(ρ, z) is specified in the equatorial plane at the (cylindrical) horizon radius ${\rho }_{\mathrm{H}}=r_{\mathrm{H}}=1+\sqrt{1-a^2}$,

$$\begin{array}{c}\hfill n_{\mathrm{e};\mathrm{H}}=n_{\mathrm{e}}({\rho }_{\mathrm{H}},z=0),\end{array}$$(1)

as is the electron temperature T_e(ρ, z),

$$\begin{array}{c}\hfill T_{\mathrm{e};\mathrm{H}}=T_{\mathrm{e}}({\rho }_{\mathrm{H}},z=0).\end{array}$$(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):

$$\begin{array}{cc}\hfill n_{\mathrm{e}}(r,z)& =n_{\mathrm{e};\mathrm{H}}{\left(\frac{r}{r_{\mathrm{H}}}\right)}^{-2}exp(-\frac{z^2}{2{(\alpha \rho )}^2}),\hfill \end{array}$$(3a)

$$\begin{array}{cc}\hfill T_{\mathrm{e}}(r)& =T_{\mathrm{e};\mathrm{H}}{\left(\frac{r}{r_{\mathrm{H}}}\right)}^{-1}.\hfill \end{array}$$(3b)

The dependence of the density profile on the cylindrical height z is thus a simple Gaussian with standard deviation s_z = αρ, where α 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 α corresponds to the tangent of its opening angle. As α → 0, we recovered equatorial models similar to those of GLM 20 and Paugnat et al. (2022), but here we kept α ≠ 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

$$\begin{array}{c}\hfill \sigma =\frac{B^2/4\pi }{m_{\mathrm{p}}{c}^2{n}_{\mathrm{e}}},\end{array}$$(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 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

$$\begin{array}{c}\hfill j_{\nu }(r)\propto \exp (-\zeta \frac{r}{r_{\mathrm{H}}}),\end{array}$$(5)

where ζ 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 azimutal 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{array}{cc}\hfill u_{\mu }^{\mathrm{circ}}\mathrm{d}{x}^{\mu }& =-u_t^{\mathrm{circ}}(-\mathrm{d}t+\ell \mathrm{d}\varphi ),\hfill \\ \hfill \ell & =\frac{{\rho }^{3/2}}{1+\rho },\hfill \end{array}$$(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 Ω ∼ ρ^−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

$$\begin{array}{c}\hfill -u_t^{\mathrm{circ}}=\frac{1}{\sqrt{-(g^{\mathit{tt}}-2g^{t\varphi }\ell +g^{\varphi \varphi }{\ell }^2)}}.\end{array}$$(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 velocity. As reviewed in Appendix C, the resulting four-velocity takes the form

$$\begin{array}{c}\hfill u_{\mathrm{rad}}^{\mu }{\mathit{\partial }}_{\mu }=u_{\mathrm{rad}}^t{\mathit{\partial }}_t+u_{\mathrm{rad}}^r{\mathit{\partial }}_r+u_{\mathrm{rad}}^{\varphi }{\mathit{\partial }}_{\varphi },\end{array}$$(8)

where

$$\begin{array}{cc}\hfill u_{\mathrm{rad}}^t& =-g^{\mathit{tt}},\hfill \\ \hfill u_{\mathrm{rad}}^r& =-\sqrt{(-1-g^{\mathit{tt}})g^{\mathit{rr}}},\hfill \\ \hfill u_{\mathrm{rad}}^{\varphi }& =-g^{t\varphi }.\hfill \end{array}$$(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^μ∂_μ = u^t∂_t + u^r∂_r + u^ϕ∂_ϕ. Following the notation of Pu et al. (2016), we introduced Ω = u^ϕ/u^t and wrote

$$\begin{array}{cc}& u^r=(1-{\beta }_r)u_{\mathrm{rad}}^r,\hfill \end{array}$$(10a)

$$\begin{array}{cc}& \mathrm{\Omega }={\mathrm{\Omega }}_{\mathrm{circ}}+(1-{\beta }_{\varphi })({\mathrm{\Omega }}_{\mathrm{rad}}-{\mathrm{\Omega }}_{\mathrm{circ}}),\hfill \end{array}$$(10b)

where 0 ≤ β_r ≤ 1 and 0 ≤ β_ϕ ≤ 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

$$\begin{array}{c}\hfill u^t=\sqrt{-\frac{1+g_{\mathit{rr}}{(u^r)}^2}{g_{\mathit{tt}}+2g_{t\varphi }\mathrm{\Omega }+g_{\varphi \varphi }{\mathrm{\Omega }}^2}},\end{array}$$(11)

and its ϕ component is then given by u^ϕ = Ωu^t.

Again, we will not examine this combined motion here, so the parameters β_r and β_ϕ 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 (α = 0.1) or thick disk (α = 1), purely azimuthal (Sect. 2.2.1) or purely radial (Sect. 2.2.2) motion, and low (a = 0.01) or high (a = 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⁹ M_⊙ and D = 16.9 Mpc (EHT L5; Event Horizon Telescope Collaboration 2019e). For the magnetization σ, we adopted a low value of 0.01, meaning that we considered a weakly magnetized disk¹. 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.

Fig. 2. Density profiles of the “thick” (α = 1, left) and “thin” (α = 0.1, right) disk models we consider. The plots display the poloidal (ρ, z) plane, with all azimuthal angles ϕ 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 3sz above the equatorial plane, where sz 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 θ turning points (T1 and T2) along the geodesic. Very sharp changes of direction appear at the θ turning points. These are due to the 2D projection in the (ρ, 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 T1, the n = 1 between T1 and T2, and the n = 2 between T2 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.

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_ν across the accretion flow. The output is an image, that is, a map of I_ν. 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:

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 × 1010 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 109 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.

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 = 104 and r = 108). The right column has a different scale compared to Fig. 3.

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

– A bright, narrow ring produced by strongly lensed photons 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.

Fig. 5. Origin of n = 2 emission. We show geodesics corresponding to the maximum intensity of the n = 2 ring of each model along a horizontal cut, on the left (blue) and right (red) sides of the image (see inset of the upper-left panel). The accretion flow and spin parameters are specified in the upper-left corner of each panel. The blue and red ellipses encircle the region emitting most of the n = 2 photons loaded onto each geodesic. The blue geodesic is abruptly cut in the lower-right panel because the medium becomes optically too thick (defined in the code as a transmission smaller than 10−6). For each ellipse, the local values of the emitted frequency (νem), self-absorbed emission (em, in cgs units), transmission (tr), and redshift factor (g3) are provided, as well as the resulting increment of specific intensity δIν (see Eq. (13) for the definition of these quantities).

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_ν and absorptivity α_ν, as well as the Planck function B_ν, 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

$$\begin{array}{c}\hfill {\tau }_{\nu }={\displaystyle {\int }_{\text{optical path}}}{\alpha }_{\nu }\mathrm{d}s,\end{array}$$(12)

where the integral is taken over the portion of the ray connecting the source to the observer, so that exp(−τ_ν) is the transmission. The local intensity increment then reads

$$\begin{array}{cc}\hfill \delta I_{\nu }& =\frac{j_{\nu }}{{\alpha }_{\nu }}[1-exp(-{\alpha }_{\nu }\delta s)]\times exp(-{\tau }_{\nu })\times g^3\hfill \end{array}$$(13a)

$$\begin{array}{cc}& =\underset{\text{self-absorbed local emission}}{\underbrace{B_{\nu }[1-exp(-{\alpha }_{\nu }\delta s)]}}\times \underset{\text{transmission}}{\underbrace{exp(-{\tau }_{\nu })}}\times \underset{\text{redshift}}{\underbrace{g^3}}\hfill \end{array}$$(13b)

$$\begin{array}{cc}& \approx B_{\nu }exp(-{\tau }_{\nu })g^3,\hfill \end{array}$$(13c)

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

$$\begin{array}{c}\hfill g=\frac{u_{\mathrm{obs}}·p_{\mathrm{obs}}}{u_{\mathrm{em}}·p_{\mathrm{em}}}.\end{array}$$(14)

Here, u_obs/em denotes the four-velocity of the observer or emitter and p_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 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 nonlocal, 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 right geodesics visit different parts of the flow, which in turn leads to the large horizontal (perpendicular-to-spin) brightness asymmetry of the n = 2 ring. This explains point (i) above.

We note that for the vertical cuts (parallel to the spin vector, shown in blue in the rightmost column of Fig. 4), the brightness asymmetry is much milder. This makes sense because at high spin, the distortion and displacement of the critical curve and n = 2 ring primarily act in the direction orthogonal to the spin.

The extreme horizontal brightness asymmetry obtained in models of radial infall can be attributed to a strong absorption of the n = 2 radiation along one of the geodesics, as well as to an extreme redshift effect due to the fact that the n = 2 emission is loaded onto the left geodesic (in blue in Fig. 5 lower-right panel) very close to the event horizon. This explains point (ii) above.

3.4. Intensity cuts

We conclude this section by revisiting the framework proposed in the introduction, in which realistic models are regarded as lying on a continuum between the shadow and wedding cake extremes. To this end, here we discuss intensity cuts of all the models analyzed in this paper, starting with the extreme models (spherical and equatorial) and then moving on to our thick-disk models.

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)∝r⁻², T_e(r)∝r⁻¹). The magnetic field was prescribed by demanding that the magnetization σ defined in Eq. (4) remain constant, at the same value as in our thick-disk model, σ = 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.

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.

The “spherical limit” in the top-left panel is the α → ∞ 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 α → 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². 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 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λ, 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₁ and b₂ in Paugnat et al. (2022).

The following criterion is used to define when the nth ring starts to dominate the visibility. The full image decomposes into layers indexed by n, I_ν(x) = I₀(x)+I₁(x)+I₂(x)+…, and the full visibility V(u) = V₀(u)+V₁(u)+V₂(u)+… does too, by linearity of the Fourier transform. We considered a sliding baseline window u ∼ u_w of fixed width 10 Gλ, which is approximately twice the photon ring diameter. At each angle φ in the baseline plane, we computed over this window u ∼ u_w the mean visibility ⟨|V(u, φ)|⟩_w of the full image I_ν(x) and the mean visibility ⟨|V_n(u, φ)|⟩_w of the nth ring image I_n(x) only. We defined the baseline threshold b_n(φ) past which the nth ring dominates the visibility to be the shortest baseline window u_w ≡ b_n(φ) such that the percentage difference between these two means dips below a given threshold p:

$$\begin{array}{c}\hfill \frac{⟨|V_n{(u,\phi )|⟩}_w-{⟨|V(u,\phi )|⟩}_w}{{⟨|V(u,\phi )|⟩}_w}<p⟹\begin{array}{c}n\text{th ring dominates}\\ \text{on}u\gtrsim b_n(\phi )\equiv u_w.\end{array}\end{array}$$(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 ≳ b₂, with the threshold b₂ ranging from ∼600 to 1400 Gλ. The corresponding visibility amplitude displays clear oscillations at levels varying between ∼10 and 50 μJy (see Table 2 and Fig. 7). Thus, these low-spin models match the idealized models considered in GLM 20 and Paugnat et al. (2022) reasonably well.

Fig. 7. Visibility amplitude profiles for all low-spin (a = 0.01) models. The full profile |V(u, φ)| is shown in black, the n = 1 profile |V1(u, φ)| in green, and the n = 2 profile |V2(u, φ)| in red. The four upper panels correspond to an orientation φ = 0° (perpendicular to the spin axis, ⊥), and the four lower panels to an orientation φ = 90° (parallel to the spin axis, ∥). The accretion flow model is specified in the top-right corner of each panel.

At high spin, by contrast, it is only on the baseline oriented parallel to the spin axis (φ = 90°) 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 ≈0.01 − 0.1 μ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 ≈100 − 1000.

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 φ = 0° (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.

On baselines perpendicular to the spin axis (φ = 0°), 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 t he 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₂(x) relative to the total flux in the full image I_ν(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 ν_obs = 230 GHz, we can also consider a higher frequency ν_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 ν_obs = 345 GHz. This higher frequency contains a clear n = 2 signal in all cases, and the associated visibility level is in the regime ∼100 μ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.

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).

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 of 2 across a decade of VLBI observations, which would influence 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.

Acknowledgments

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 helpful 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 horizon 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

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

$$\begin{array}{c}\hfill \frac{\mathrm{d}E}{\mathrm{d}t\mathrm{d}\nu }=\sqrt{3}\frac{e^{3}Bsin\theta }{m{c}^2}F\left(\frac{\nu }{{\nu }_{\mathrm{crit}}}\right),\end{array}$$(A.1)

where

$$\begin{array}{c}\hfill {\nu }_{\mathrm{crit}}=\frac{3}{2}{\gamma }^2{\nu }_{\mathrm{cyclo}}sin\theta ,{\nu }_{\mathrm{cyclo}}=\frac{\mathit{eB}}{2\pi mc}.\end{array}$$(A.2)

Here, γ denotes the Lorentz factor of the electron, θ the angle between the direction of emission and the magnetic field, and

$$\begin{array}{c}\hfill F(x)\equiv x{\displaystyle {\int }_x^{\infty }}{K}_{5/3}(u)\mathrm{d}u=\{\begin{array}{cc}{2}^{2/3}\mathrm{\Gamma }\left(\frac{2}{3}\right)x^{1/3}\hfill & \text{for}x\ll 1,\hfill \\ \sqrt{\frac{\pi x}{2}}{e}^{-x}\hfill & \text{for}x\gg 1,\hfill \end{array}\end{array}$$(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 θ, 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

$$\begin{array}{cc}\hfill j_{\nu }& \equiv \frac{\mathrm{d}E}{\mathrm{d}t\mathrm{d}\nu \mathrm{d}V\mathrm{d}\mathrm{\Omega }}=\frac{\sqrt{3}}{4\pi }\frac{e^{3}Bsin\theta }{m{c}^2}{\displaystyle {\int }_0^{\infty }}F\left(\frac{\nu }{{\nu }_{\mathrm{crit}}}\right)\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}\gamma }\mathrm{d}\gamma ,\hfill \end{array}$$(A.4)

where

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}\gamma }=\frac{n_{\mathrm{e}}}{{\mathrm{\Theta }}_{\mathrm{e}}}\frac{\gamma {({\gamma }^2-1)}^{1/2}}{K_2\left({\mathrm{\Theta }}_{\mathrm{e}}^{-1}\right)}exp(-\frac{\gamma }{{\mathrm{\Theta }}_{\mathrm{e}}})\end{array}$$(A.5)

is the relativistic thermal (Maxwell-Jüttner) distribution, K₂(u) is a modified Bessel function of the second kind, Θ_e = k_BT_e/mc² is the dimensionless electron temperature, and the factor of (4π)⁻¹ implicitly assumes an isotropic distribution of the electrons’ momenta. We note that the integral over γ in Eq. (A.4) is evaluated over the range [0, ∞), even though physically, γ ≥ 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 ≪ 1 and x ≫ 1, the emissivity (A.4) may be analytically expressed as (see, e.g., Leung et al. 2011)

$$\begin{array}{c}\hfill j_{\nu }\approx \{\begin{array}{cc}{\displaystyle \frac{2^{4/3}\pi }{3}\frac{n_{\mathrm{e}}{e}^2{\nu }_s}{c{\mathrm{\Theta }}_{\mathrm{e}}^2}{X}^{1/3}}\hfill & \text{for}X\ll 1,\hfill \\ {\displaystyle n_{\mathrm{e}}\frac{\sqrt{2}\pi e^2{\nu }_s}{6c{\mathrm{\Theta }}_{\mathrm{e}}^2}X{e}^{-X^{1/3}}}\hfill & \text{for}X\gg 1,\hfill \end{array}\end{array}$$(A.6)

where

$$\begin{array}{c}\hfill X=\frac{\nu }{{\nu }_s},{\nu }_s=\frac{2}{9}{\nu }_{\mathrm{cyclo}}{\mathrm{\Theta }}_{\mathrm{e}}^{2}sin\theta .\end{array}$$(A.7)

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

$$\begin{array}{c}\hfill j_{\nu }^{\mathrm{approx}}=n_{\mathrm{e}}\frac{\sqrt{2}\pi e^2{\nu }_s}{3c{K}_2\left({\mathrm{\Theta }}_{\mathrm{e}}^{-1}\right)}{(X^{1/2}+2^{11/12}{X}^{1/6})}^2{e}^{-X^{1/3}}.\end{array}$$(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 γ, 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 θ:

$$\begin{array}{c}\hfill \langle j_{\nu }^{\mathrm{approx}}\rangle =\frac{1}{4\pi }{\displaystyle \int }{j}_{\nu }^{\mathrm{approx}}\mathrm{d}\mathrm{\Omega }=\frac{1}{2}{\displaystyle {\int }_0^{\pi }}{j}_{\nu }^{\mathrm{approx}}sin\theta \mathrm{d}\theta ,\end{array}$$(A.9)

where dΩ = sin θdθdϕ and the emissivity is independent of the azimuthal angle ϕ.

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

$$\begin{array}{c}\hfill j_{\nu }^{\mathrm{exact}}=\frac{2\pi e^2\nu }{c|cos\theta |}{\displaystyle \sum _{n=n_{-}}^{\infty }{\int }_{{\gamma }_{-}(n)}^{{\gamma }_{+}(n)}\frac{1}{\beta }{\mathcal{J}}_n(\gamma )\frac{\mathrm{d}{n}_{\mathrm{e}}}{\mathrm{d}\gamma }\mathrm{d}\gamma ,}\end{array}$$(A.10)

where $n_{-}=\frac{\nu }{{\nu }_{\mathrm{cyclo}}}|sin\theta |$ and

$$\begin{array}{c}\hfill {\gamma }_{\pm }(n)=\frac{n\frac{{\nu }_{\mathrm{cyclo}}}{\nu }\pm |cos\theta |\sqrt{{\left(n\frac{{\nu }_{\mathrm{cyclo}}}{\nu }\right)}^2-{sin}^2\theta }}{{sin}^2\theta },\end{array}$$(A.11)

where n is the integer that is being summed over in Eq. (A.10). The dimensionless function 𝒥_n(γ) is given in terms of the Bessel function of the first kind J_n(z) and its derivative $J_n^{\prime }(z)$ by

$${\mathcal{J}}_n(\gamma )=\left[M{J}_n{(z)}^2\right]+\left[N{J}_n^{\prime }{(z)}^2\right],$$(A.12a)

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

$$N=\beta sin\alpha ,cos{\alpha }_0=\frac{1}{\beta cos\theta }1-\frac{n}{\gamma }\frac{{\nu }_{\text{cyclo}}}{\nu }.$$(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 θ = 30° and 80° is displayed as a function of Θ_e and ν/ν_cyclo. This error is typically on the order of a few percent for Θ_e ∈ [1, 10] and ν/ν_cyclo ∈ [10³, 10⁵], 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 ≈ 10¹¹ K, so Θ_e ≈ 15 (and the accuracy gets better with increasing Θ_e), while the magnetic field is B ≈ 10 G, which implies a cyclotron frequency ν_cyclo ≈ 10⁷ Hz, so that ν/ν_cyclo ≈ 10⁴. On the other hand, the error typically exceeds 10% when Θ_e ≲ 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 α_ν = j_ν/B_ν, where B_ν 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*

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\to 0}{\approx }2x^{-2}$. This substitution results in the further approximation

$$\begin{array}{c}\hfill j_{\nu }\stackrel{{\mathrm{\Theta }}_{\mathrm{e}}\gg 1}{\approx }{n}_{\mathrm{e}}\frac{\sqrt{2}\pi e^2{\nu }_{\mathrm{cyclo}}sin\theta }{27c}{(X^{1/2}+2^{11/12}{X}^{1/6})}^2{e}^{-X^{1/3}},\end{array}$$(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

$$\begin{array}{c}\hfill X\approx \frac{1600}{sin\theta }\left(\frac{\nu [\mathrm{Hz}]}{{10}^{12}}\right)\left(\frac{10}{B[\mathrm{G}]}\right){\left(\frac{10}{{\mathrm{\Theta }}_{\mathrm{e}}}\right)}^2,\end{array}$$(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 ≫ 1 for all typical values of M87* parameters, we can further simplify the approximation (B.1) for the emissivity by taking its X → ∞ limit. In the regime X ≫ 1, we may replace (X^1/2 + 2^11/12X^1/6)² by X to find

$$\begin{array}{c}\hfill j_{\nu }\stackrel{{\mathrm{\Theta }}_{\mathrm{e}},X\gg 1}{\approx }{n}_{\mathrm{e}}\frac{\sqrt{2}\pi e^2{\nu }_{\mathrm{cyclo}}sin\theta }{27c}X{e}^{-X^{1/3}},\end{array}$$(B.3)

which is a slight underestimation because (X^1/2 + 2^11/12X^1/6)² 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*.

Fig. B.1. Validity of the approximated emissivity given by Eq. (B.3). This figure shows the ratio of (X1/2 + 211/12X1/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.

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

$$\begin{array}{c}\hfill j_{\nu }\propto \frac{n_{\mathrm{e}}\nu }{{\mathrm{\Theta }}_{\mathrm{e}}^2}{e}^{-X^{1/3}}=\frac{n_{\mathrm{e}}\nu }{{\mathrm{\Theta }}_{\mathrm{e}}^2}exp[-\mathcal{C}{\left(\frac{\nu }{B{\mathrm{\Theta }}_{\mathrm{e}}^2}\right)}^{1/3}],\end{array}$$(B.4)

where 𝒞 = (9πmc/e)^1/3 is a constant.

Since we are assuming that the density and temperature fall off as r⁻² and r⁻¹, respectively, we may immediately conclude that the radial dependence of the emissivity is simply given by

$$\begin{array}{c}\hfill j_{\nu }(r)\propto exp[-X{(r)}^{1/3}],\end{array}$$(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 ∝ r⁻¹ [Eq. (4)], we may now write

$$\begin{array}{c}\hfill X\approx \frac{3.7\times {10}^5}{B_{\mathrm{inner}}{\mathrm{\Theta }}_{\mathrm{e};\mathrm{inner}}^{2}sin\theta }{\left(\frac{r}{r_{\mathrm{inner}}}\right)}^3,\end{array}$$(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,

$$\begin{array}{c}\hfill j_{\nu }(r)\propto exp[-\zeta \frac{r}{r_{\mathrm{inner}}}],\zeta ={\left(\frac{3.7\times {10}^5}{B_{\mathrm{inner}}{\mathrm{\Theta }}_{\mathrm{e};\mathrm{inner}}^{2}sin\theta }\right)}^{1/3}.\end{array}$$(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 ≫ 1 and X ≫ 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 ∝ r⁻², Θ_e ∝ r⁻¹, B ∝ r⁻¹). 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).

Fig. B.2. Values of the parameter ζ [Eq. (B.7)] as a function of the magnetic field strength B and dimensionless temperature Θe. Here, we fix sin θ = 1. Our models always stay in the part of the plane located above the dashed horizontal red line (Θ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 ζ.

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

$$\begin{array}{c}\hfill {\alpha }_{\nu }=\frac{j_{\nu }}{B_{\nu }(T_{\mathrm{e}})},\end{array}$$(B.8)

where B_ν(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ν ≪ k_BT_e, so that this relation simplifies to

$$\begin{array}{c}\hfill {\alpha }_{\nu }\approx \frac{c^2}{2k_B{T}_{\mathrm{e}}{\nu }^2}{j}_{\nu }.\end{array}$$(B.9)

Assuming as usual a temperature profile T ∝ r⁻¹, we obtain

$$\begin{array}{c}\hfill {\alpha }_{\nu }(r)\propto \frac{r}{r_{\mathrm{inner}}}exp[-\zeta \frac{r}{r_{\mathrm{inner}}}].\end{array}$$(B.10)

The highly lensed emission originates from small radii r ≳ r_inner, for which the exponential profiles (B.7) and (B.10) provide good approximations to the emissivity and absorptivity, respectively.

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

$$\begin{array}{cc}\hfill j_{\nu }^{\mathrm{PL}}& =n_{\mathrm{e}}\frac{e^2{\nu }_{\mathrm{cyclo}}}{c}\frac{3^{p/2}(p-1)sin\theta }{2(p+1)({\gamma }_{\min }^{1-p}-{\gamma }_{\max }^{1-p})}\hfill \\ \hfill & \times \mathrm{\Gamma }\left(\frac{3p-1}{12}\right)\mathrm{\Gamma }\left(\frac{3p+19}{12}\right){\left(\frac{\nu }{{\nu }_{\mathrm{cyclo}}}\right)}^{-(p-1)/2},\hfill \end{array}$$(B.11)

where γ_max/min are the maximum or minimum values of the Lorentz factor for electrons in the distribution and Γ(x) the gamma function. Following the same steps as in the previous section leads to

$$\begin{array}{c}\hfill j_{\nu }^{\mathrm{PL}}\propto {\left(\frac{r}{r_{\mathrm{inner}}}\right)}^{-(p+5)/2},\end{array}$$(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_{\nu }^{\text{PL}}\propto r^{-2}$.

The power-law absorptivity ${\alpha }_{\nu }^{\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 }_{\nu }^{\text{PL}}$ ∝ 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^μ and its associated 1-form u_μ take the general form

$$\begin{array}{cc}\hfill u_{\mathrm{circ}}^{\mu }{\mathit{\partial }}_{\mu }& =u^t({\mathit{\partial }}_t+\mathrm{\Omega }{\mathit{\partial }}_{\varphi }),\hfill \end{array}$$(C.1a)

$$\begin{array}{cc}\hfill u_{\mu }^{\mathrm{circ}}\mathrm{d}{x}^{\mu }& =-u_t(-\mathrm{d}t+\ell \mathrm{d}\varphi ),\hfill \end{array}$$(C.1b)

where Ω = u^ϕ/u^t denotes the angular velocity and ℓ = −u_ϕ/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

$$\begin{array}{c}\hfill -u_t=\frac{r-2}{\sqrt{r(r-3)}},\ell =\frac{r^{3/2}}{r-2},\mathrm{\Omega }=r^{-3/2}.\end{array}$$(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_ISCO = 6.

Below the ISCO, Cunningham (1975) prescribed that the geodesic constants of motion −u_t and u_ϕ (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

$$\begin{array}{c}\hfill u_{\mu }^{\mathrm{C}75}\mathrm{d}{x}^{\mu }=u_t^{\mathrm{ISCO}}\mathrm{d}t+u_r\mathrm{d}r+u_{\varphi }^{\mathrm{ISCO}}\mathrm{d}\varphi ,\end{array}$$(C.3)

where the (t, ϕ) components $u_{t,\varphi }^{\text{ISCO}}\equiv u_{t,\varphi }(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:

$$\begin{array}{c}\hfill u_r^2=g_{\mathit{rr}}[-1-g^{\mathit{tt}}{\left(u_t^{\mathrm{ISCO}}\right)}^2-g^{\varphi \varphi }{\left(u_{\varphi }^{\mathrm{ISCO}}\right)}^2].\end{array}$$(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 ≥ r_ISCO, and by Eq. (C.3) for r_H < r < r_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 (ρ > 6), we simply defined the specific angular momentum ℓ to take the same value at every height z above the equator as it does in the equatorial Keplerian disk model.

That is, for ρ > 6, we imposed the axially symmetric profile

$$\begin{array}{c}\hfill \ell (\rho ,z)=\frac{{\rho }^{3/2}}{\rho -2},\end{array}$$(C.5)

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

$$\begin{array}{c}\hfill u_{\mu }^{\mathrm{thick};>\mathrm{ISCO}}\mathrm{d}{x}^{\mu }=-u_t(\rho ,z)(-\mathrm{d}t+\frac{{\rho }^{3/2}}{\rho -2}\mathrm{d}\varphi ),\end{array}$$(C.6)

whose unit-normalization fixes the time component to be

$$\begin{array}{c}\hfill -u_t(\rho ,z)={(-g^{\mathit{tt}}-\frac{\rho }{{(\rho -2)}^2})}^{-1/2},\end{array}$$(C.7)

where we used the fact that g^ϕϕ = (r sin θ)⁻² = ρ⁻². 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 ρ > 6 (i.e., that the expression in parentheses remains positive).

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

$$\begin{array}{c}\hfill u_{\mu }^{\mathrm{thick};<\mathrm{ISCO}}\mathrm{d}{x}^{\mu }=u_t^{\mathrm{ISCO}}(z)\mathrm{d}t+u_r(\rho ,z)\mathrm{d}r+u_{\varphi }^{\mathrm{ISCO}}(z)\mathrm{d}\varphi ,\end{array}$$(C.8)

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

$$\begin{array}{c}\hfill {\left[u_r(\rho ,z)\right]}^2=g_{\mathit{rr}}(-1-g^{\mathit{tt}}{\left[u_t^{\mathrm{ISCO}}(z)\right]}^2-g^{\varphi \varphi }{\left[u_{\varphi }^{\mathrm{ISCO}}(z)\right]}^2).\end{array}$$(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),

$$\begin{array}{c}\hfill u_{\mu }^{\mathrm{circ}}(\rho ,z)\mathrm{d}{x}^{\mu }=-u_t(\rho ,z)(-\mathrm{d}t+\ell (\rho )\mathrm{d}\varphi ),\end{array}$$(C.10)

which assumes that the specific angular momentum depends only on cylindrical radius. The function ℓ(ρ) 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

$$\begin{array}{c}\hfill \mathrm{\Omega }\equiv \frac{u^{\varphi }}{u^t}=\frac{g^{\varphi \varphi }{u}_{\varphi }}{g^{\mathit{tt}}{u}_t}=\frac{\ell }{{\rho }^2}(1-\frac{2}{r}),\end{array}$$(C.11)

and unit-normalization fixes its time component to be

$$\begin{array}{c}\hfill -u_t={(-g^{\mathit{tt}}-g^{\varphi \varphi }{\ell }^2)}^{-1/2}=\sqrt{\frac{1-2/r}{1-\mathrm{\Omega }\ell }},\end{array}$$(C.12)

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

$$\begin{array}{c}\hfill \mathrm{\Omega }\ell <1,\end{array}$$(C.13)

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

$$\begin{array}{c}\hfill \frac{{\ell }^2}{{\rho }^2}(1-\frac{2}{r})<1.\end{array}$$(C.14)

Since ρ = r sin θ, we have ρ − 2 sin θ > 0 outside the horizon, where the condition (C.14) becomes

$$\begin{array}{c}\hfill |\ell (\rho )|<\frac{{\rho }^{3/2}}{\sqrt{\rho -2sin\theta }}.\end{array}$$(C.15)

We have thus far imposed no conditions on the function ℓ(ρ). At this point, we will assume that it converges to the Newtonian profile ℓ_Newton = ρ^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

$$\begin{array}{c}\hfill {\ell }_{\alpha }(\rho )=\frac{{\rho }^{3/2}}{\rho +\alpha },\end{array}$$(C.16)

where α is a real constant such that α = 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 α > 0, so that the denominator may never vanish. This also ensures that ℓ > 0, and reduces condition (C.15) to

$$\begin{array}{c}\hfill {(\rho +\alpha )}^2-\rho +2sin\theta >0.\end{array}$$(C.17)

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

$$\begin{array}{c}\hfill \alpha >\frac{1}{4},\end{array}$$(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 α = 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 ρ = 0. Taking the limit ρ → 0 keeping z > 2 fixed results in

$$\begin{array}{c}\hfill u^t\approx {(1-\frac{2}{r})}^{-1/2},u^{\varphi }\approx {\rho }^{-1/2}{(1-\frac{2}{r})}^{1/2},\end{array}$$(C.19)

so the ϕ component of the four-velocity, as well as the angular velocity Ω, diverge in the limit ρ → 0. Indeed, Eq. (C.11) shows that the angular velocity grows like Ω ≈ ℓρ⁻² ≈ ρ^−1/2 near the axis for the Keplerian ℓ, which behaves as ℓ ≈ ρ^3/2. Hence, any family ℓ(ρ) that includes the Keplerian profile as a member must have a divergent Ω near ρ = 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 }{\mathit{\partial }}_{\mu }=u^t{\mathit{\partial }}_t+u^r{\mathit{\partial }}_r+u^{\varphi }{\mathit{\partial }}_{\varphi },$$(C.20a)

$$u_t=-1,u_{\varphi }=0.$$(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, ϕ) contravariant components are

$$\begin{array}{c}\hfill u^t=g^{\mathit{tt}}{u}_t=-g^{\mathit{tt}},u^{\varphi }=g^{t\varphi }{u}_t=-g^{t\varphi }.\end{array}$$(C.21)

We note that u^ϕ 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

$$\begin{array}{c}\hfill u_r=-\sqrt{\frac{-1-g^{\mathit{tt}}}{g^{\mathit{rr}}}},\end{array}$$(C.22)

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

$$\begin{array}{c}\hfill u^r=g^{\mathit{rr}}{u}_r=-\sqrt{(-1-g^{\mathit{tt}})g^{\mathit{rr}}}.\end{array}$$(C.23)

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

$$\begin{array}{c}\hfill u^t={(1-\frac{2}{r})}^{-1},u^r=-\sqrt{\frac{2}{r}},u^{\varphi }=0.\end{array}$$(C.24)

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

$$\begin{array}{c}\hfill \mathbf{v}=v^r{\mathbf{e}}_r=\frac{\mathrm{d}r}{\mathrm{d}t}{\mathbf{e}}_r=\frac{u^r}{u^t}{\mathbf{e}}_r,\end{array}$$(C.25)

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

$$\begin{array}{c}\hfill v^r=-\sqrt{\frac{2}{r}}(1-\frac{2}{r})\stackrel{r\to \infty }{\approx }\sqrt{\frac{2}{r}},\end{array}$$(C.26)

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

$$\begin{array}{c}\hfill \frac{1}{2}{\left(v^r\right)}^2-\frac{1}{r}=0.\end{array}$$(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.

Fig. D.1. A null geodesic of order n = 2 (two θ turning points), in blue. The black disk is the event horizon and the yellow lines indicate the angular turning points θ− and θ+ = π − θ−, which are computed from the analytical expression for θ− 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.

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

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

where ℛ(r) and Θ(θ) are radial and angular geodesic potentials, and the integrals are evaluated along the entire trajectory with the signs ±_r, θ 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, θ_s) to an observer at (r_o, θ_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_θ|, as well as of the conserved energy E = −p_t and azimuthal angular momentum L = p_ϕ.

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 = −pt (blue) and L = pϕ (orange), as well as the quantity |Ir − Gθ| (green), as a function of the integration step. The peak of |Ir − Gθ| happens close to an angular turning point where the geodesic evolves in the Kerr spherical null geodesics region.

These quantities should remain constant along the geodesic, and the figure shows that their error is consistently below 10⁻¹³, except for a peak of 10⁻⁹ in |I_r − G_θ| 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 sn(x|k):

$$\begin{array}{c}\hfill r_{\mathrm{eq}}^{(m)}=\frac{r_4{r}_{31}-r_3{r}_{41}{\mathrm{sn}}^2\left(\frac{1}{2}\sqrt{r_{31}{r}_{42}}{\tau }_{\mathrm{eq}}^{(m)}-{\mathcal{F}}_o|k\right)}{r_{31}-r_{41}{\mathrm{sn}}^2\left(\frac{1}{2}\sqrt{r_{31}{r}_{42}}{\tau }_{\mathrm{eq}}^{(m)}-{\mathcal{F}}_o|k\right)}.\end{array}$$(D.2)

Here, r_ij = r_i − r_j with {r₁, r₂, r₃, r₄} the four roots of the quartic potential ℛ(r) (with exact expressions given in Eq. (A8) of Gralla & Lupsasca 2020a), and ℱ_o is the elliptic integral of the first kind

$$\begin{array}{c}\hfill {\mathcal{F}}_o=F\left(arcsin\sqrt{\frac{r_{31}}{r_{41}}}|k\right),k=\frac{r_{32}{r}_{41}}{r_{31}{r}_{42}}.\end{array}$$(D.3)

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

$$\begin{array}{c}\hfill {\tau }_{\mathrm{eq}}^{(m)}=\frac{2mK-\mathrm{sign}\left(p_o^{\theta }\right)F_o}{a\sqrt{-u_{-}}},\end{array}$$(D.4)

where $p_o^{\theta }$ is the polar photon momentum at the observer,

$$\begin{array}{c}\hfill u_{\pm }={\mathrm{\Delta }}_{\theta }\pm \sqrt{{\mathrm{\Delta }}_{\theta }^2+\frac{\eta }{a^2}},{\mathrm{\Delta }}_{\theta }=\frac{1}{2}(1-\frac{\eta +{\lambda }^2}{a^2}),\end{array}$$(D.5)

are the zeros u = cos²θ of Θ(u) in terms of the energy-rescaled angular momentum λ = L/E and Carter constant η = Q/E², and

$$\begin{array}{c}\hfill F_o=F\left(arcsin\left(\frac{cos{\theta }_o}{\sqrt{u_{+}}}\right)|\frac{u_{+}}{u_{-}}\right),K=F\left(\frac{\pi }{2}|\frac{u_{+}}{u_{-}}\right)\end{array}$$(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_{\text{eq}}^{(m)}$ agree to within 10⁻⁴ − 10⁻⁵M. We note that the equatorial crossing labeled τ^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_θ| at that crossing. Finally, we have compared the value of the radial turning point r₄ (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 ≲10⁻⁸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 δ = 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)

$$\begin{array}{c}\hfill \tau =I_r=G_{\theta },\end{array}$$(E.1)

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

$$\begin{array}{c}\hfill \tau (\theta )=\mathrm{sign}\left(p_s^{\theta }\right)({\mathcal{G}}_{\theta }^o-{\mathcal{G}}_{\theta }^s),\end{array}$$(E.2)

where $p_s^{\theta }$ is the polar photon momentum at θ = θ_s, while

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

and ${\mathcal{G}}_{\theta }^o$ and ${\mathcal{G}}_{\theta }^s$ are evaluated at the 𝒢_θ evaluated at θ = θ_o and θ = θ_s, respectively. We can determine the Mino time τ₁ elapsed at the first angular turning point by setting

$$\begin{array}{c}\hfill {\theta }_s=arccos(\pm \sqrt{u_{+}}),\pm =\mathrm{sign}\left(p_o^{\theta }\right),\end{array}$$(E.4)

while the Mino time elapsed between two successive angular is

$$\begin{array}{c}\hfill \mathrm{\Delta }\tau =\frac{2K}{a\sqrt{-u_{-}}}.\end{array}$$(E.5)

The order of any null geodesic is then defined as

• n = 0: between the observer and τ = τ₁,

• n = 1: between τ = τ₁ and τ = τ₁ + Δτ,

• n = 2: between τ = τ₁ + Δτ and τ = τ₁ + 2Δτ.

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 ≤ 2), but also individual images of each layer n ∈ {0, 1, 2} consisting of photons loaded onto each null geodesic at the corresponding order.

All Tables

All Figures

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.

In the text

Fig. 2. Density profiles of the “thick” (α = 1, left) and “thin” (α = 0.1, right) disk models we consider. The plots display the poloidal (ρ, z) plane, with all azimuthal angles ϕ 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 3sz above the equatorial plane, where sz 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 θ turning points (T1 and T2) along the geodesic. Very sharp changes of direction appear at the θ turning points. These are due to the 2D projection in the (ρ, 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 T1, the n = 1 between T1 and T2, and the n = 2 between T2 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.

In the text

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 × 1010 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 109 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.

In the text

	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 = 104 and r = 108). The right column has a different scale compared to Fig. 3.
In the text

Fig. 5. Origin of n = 2 emission. We show geodesics corresponding to the maximum intensity of the n = 2 ring of each model along a horizontal cut, on the left (blue) and right (red) sides of the image (see inset of the upper-left panel). The accretion flow and spin parameters are specified in the upper-left corner of each panel. The blue and red ellipses encircle the region emitting most of the n = 2 photons loaded onto each geodesic. The blue geodesic is abruptly cut in the lower-right panel because the medium becomes optically too thick (defined in the code as a transmission smaller than 10−6). For each ellipse, the local values of the emitted frequency (νem), self-absorbed emission (em, in cgs units), transmission (tr), and redshift factor (g3) are provided, as well as the resulting increment of specific intensity δIν (see Eq. (13) for the definition of these quantities).

In the text

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.

In the text

Fig. 7. Visibility amplitude profiles for all low-spin (a = 0.01) models. The full profile |V(u, φ)| is shown in black, the n = 1 profile |V1(u, φ)| in green, and the n = 2 profile |V2(u, φ)| in red. The four upper panels correspond to an orientation φ = 0° (perpendicular to the spin axis, ⊥), and the four lower panels to an orientation φ = 90° (parallel to the spin axis, ∥). The accretion flow model is specified in the top-right corner of each panel.

In the text

	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 φ = 0° (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.
In the text

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).

In the text

Fig. B.1. Validity of the approximated emissivity given by Eq. (B.3). This figure shows the ratio of (X1/2 + 211/12X1/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.

In the text

Fig. B.2. Values of the parameter ζ [Eq. (B.7)] as a function of the magnetic field strength B and dimensionless temperature Θe. Here, we fix sin θ = 1. Our models always stay in the part of the plane located above the dashed horizontal red line (Θ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 ζ.

In the text

Fig. D.1. A null geodesic of order n = 2 (two θ turning points), in blue. The black disk is the event horizon and the yellow lines indicate the angular turning points θ− and θ+ = π − θ−, which are computed from the analytical expression for θ− 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.

In the text

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 = −pt (blue) and L = pϕ (orange), as well as the quantity |Ir − Gθ| (green), as a function of the integration step. The peak of |Ir − Gθ| happens close to an angular turning point where the geodesic evolves in the Kerr spherical null geodesics region.

In the text