Photon Ring Interferometric Signatures Beyond The Universal Regime

Introduction

The first, 230 GHz, Event Horizon Telescope (EHT) images of the supermassive black holes (SMBH) M87* [1] and Sgr A* [2] feature a bright, asymmetric ring, in complete accordance with ray traced images based on general-relativistic magnetohydrodynamic (GRMHD) simulations [3]. Theoretical models predict that, after enough time-averaging, the small-scale structure in such SMBH images comes to be dominated by a series of (not yet observed) narrow lensed photon rings [4, 5, 6, 7] exhibiting a universal interferometric signature from which one can infer the projected diameter of the ring. This observable diameter encodes information about the black hole mass and spin [4, 8], and can be measured using very-long-baseline interferometric (VLBI) observations on a single baseline extending from a ground-based telescope to one orbiting satellite, without requiring a dense sampling of the visibility plane to reconstruct the full image.

Prior studies [4, 9, 10] of the interferometric signature of a thin ring with diameter $d$ and width $w \ll d$ have shown that in the universal regime $1/d \ll u \ll1/w$, the visibility amplitude displays a ringing pattern of weakly damped oscillations with respect to the baseline length $u$. While the envelope of the oscillations depends on the detailed intensity profile of the ring, their frequency only depends on the ring diameter in the zero width limit, thus enabling a measurement of the ring diameter from the visibility amplitude alone. Unfortunately, since realistic photon rings have finite width, the condition $u \ll1/w$ may not hold in practice, especially for the $n=1$ photon ring, which will be the first target of future high-resolution space-VLBI observations.

In this work, we study the interferometric signatures of finite width photon rings via a perturbative approach. We show that, beyond the universal regime, the visibility amplitude of a finite-width, non-radially-symmetric ring displays a baseline-dependent oscillation frequency (the precise meaning of this will be clarified in what follows). In principle, this implies that observations on different baselines measure different ring diameters (due to finite ring-width effects), leading to some bias in the diameter inference when analyzing VLBI data (and hence, a bias in the inference of physical parameters from such data). We show that this bias can be modeled by including in the visibility phase correction terms proportional to $u^3$, $u^5$, etc. Our results also clarify what precise diameter is measured by VLBI observations of a finite-width ring: the center-of-light diameter defined in Sec. III below.

The remainder of this paper is organized as follows. In Sec. II, we review the visibility calculations of Ref. [9] for general convex thin rings within the universal regime. Next, in Sec. III, we introduce the notion of sliceable thick rings and describe an algorithm for decomposing a finite-width ring into a sum of thin rings. We then calculate the visibility of a sliceable ring in a perturbative expansion in $uw \lesssim1$, which allows us to push beyond the universal regime. We apply our method to ray traced GRMHD images in Sec. IV, and we find that keeping the higher-order corrections is necessary for accurately modeling the visibility for M87*. Finally, we conclude in Sec. V with a discussion of the implications of our results for spaceborne measurements of the photon ring, as well as potential directions for future research. We provide a python implementation of our method to slice a thick ring and calculate its visibility in the publicly available package hiring.¹

Thin ring visibility revisited

In this section, we briefly review the prior work [9] on the visibility of a thin, convex ring that forms the basis of our subsequent study of finite-width rings. In Sec. II A, we follow Gralla’s analysis [9] to obtain the visibility of a general thin ring in the large-$u$ limit, assuming a smooth angular profile for the ring. Alternatively, in Sec. II B, we derive the exact visibility of a circular, non-axisymmetric thin ring, before taking the same large-$u$ limit to show that the result agrees with that of Sec. II A.

General derivation for convex rings

An interferometer samples the radio visibility (the complex Fourier transform) of the sky intensity $I(x,y)$,⁴

where $\mathbf{r}=(x,y)$ is the sky position in radians while the baseline $\mathbf{u}=(u_x,u_y)=(u\cos\phi_{\mathrm{obs}},u\sin\phi_{\mathrm{obs}})$ is the distance (in observation wavelengths) between telescopes projected into the plane perpendicular to the line of sight. Without loss of generality, we can assume that $\mathbf{u}$ lies along the $x$-axis (i.e., that $\phi_{\mathrm{obs}}=0$), so that $V(u)$ only depends on the projected intensity $\mathcal{P}_x I(x)=\int I(\mathbf{r})\,\mathrm{d}y$,

For an infinitely thin ring parameterized by its arclength $s$ as $(x_s(s),y_s(s))$, the sky intensity is

and thus, the projected intensity is given by

where the sum runs over all $s$ such that $x=x_s(s)$. If for some $s_0$ in the sum, $x_s'(s_0)=0$ with $\mathcal{I}(s_0)\ne0$, then $\mathcal{P}_x I(x)$ diverges at $x_0=x_s(s_0)$, where it behaves as $\mathcal{I}(s)/|\mathrm{d}x_s/\mathrm{d}s|\approx\mathcal{I}(s_0)/\sqrt{|2(x-x_0)/\mathcal{R}_0|}$, with $\mathcal{R}_0$ the curvature radius of the ring at $x_0$ [9, 10].

Heuristically, the visibility at large $u$ is dominated by the small-scale features in $\mathcal{P}_x I$. Assuming that $\mathcal{I}(s)$ is everywhere positive with characteristic correlation scale $\xi_{\mathcal{I}}$, the visibility at $u\gg\max(1/d,1/\xi_{\mathcal{I}})$ is dominated by the $|x-x_0|^{-1/2}$ divergences in $\mathcal{P}_x I$ introduced by the projection effects at the edge points where $\mathrm{d}x_s/\mathrm{d}s=0$.

Strictly speaking, $\xi_{\mathcal{I}}$ should also account for any “jaggedness” in the ring shape $(x_s(s),y_s(s))$, which can also introduce small-scale structure in the projection (4) beyond the $|x-x_0|^{-1/2}$ singularities near the edge points, even if $\mathcal{I}(s)$ is smooth. To our knowledge, the condition $u\gg1/\xi_{\mathcal{I}}$ has not been explicitly noted before.

For a convex shape, there are only two edge points, namely the leftmost ($x=x_L$) and rightmost ($x=x_R$) points in the projected intensity profile, which together define the projected diameter $d=x_R-x_L$ of the ring. The Heaviside step function $H(x)$ has Fourier transform

which lets us compute the visibility of an infinitely thin, convex ring in the limit of large $u\gg\max(1/d,1/\xi_{\mathcal{I}})$:

where $A_i=\mathcal{I}(x_i)\sqrt{\mathcal{R}_i}$. Having obtained this formula, we conclude that it continues to hold for a thin ring (of width $w\ll1/u$), which must produce the same interferometric response on baselines too short to resolve its profile.

In practice, because of the difficulties inherent in VLBI phase calibration, it is likely that in the near future we will only be able to measure the visibility amplitude

which gives the universal interferometric signature of a (thin) black hole photon ring. From (7), it is evident that the oscillating frequency of $|V(u)|^2$ with respect to $u$ is determined by the projected diameter $d$, while the angular asymmetry ($A_L\ne A_R$) of the ring only affects the amplitude of the oscillation.

To summarize, (7) gives the visibility of a (convex) thin ring in the universal regime

In this regime, the visibility amplitude takes the universal form (7), such that the locations of its peaks and valleys depend only on the projected diameter $d$ of the ring.

Finally, when $\mathcal{I}_L\mathcal{I}_R=0$, at least one of the left/right singularities in the projected intensity disappears, and hence (7) may not accurately represent the visibility (which may instead be dominated by other features in the intensity profile).⁵ In the next section, we show explicit examples where the condition $u\gg1/\xi_{\mathcal{I}}$ fails, so that the locations of the peaks and valleys shift as a result.

Fourier-based derivation for circular rings

In the last derivation, we first argued—following Gralla [9]—that in the large-$u$ limit, the visibility is dominated by the $|x-x_0|^{-1/2}$ divergences in the projection of the ring intensity, and then we approximated the visibility of the ring as the sum of the (shifted) visibilities of the leftmost and rightmost singularities. In this section, we first derive the exact visibility (valid for any $u \ll1/w$) of a circular yet non-axisymmetric thin ring, and then we demonstrate that in the $u \gg\max(1/d,1/\xi_I)$ limit, our result is consistent with (7).

First, a general circular thin ring may be decomposed into Fourier modes $\beta_m$ as

and its visibility is then given by [4]

where $J_m(\cdot)$ denotes the $m$th Bessel function. Again, we assume that $u$ is parallel to the $x$-axis; for the projection along a different direction, one can simply rotate the $\beta_m$ coefficients and our main conclusions here will still hold. As shown in App. A, in the long-baseline limit, we have

where the even and odd coefficients $C_e$ and $C_o$ are real, determined by the $\beta_m$, and satisfy

such that (10) is consistent with (7).

For a physical photon ring, $I_L>0$ and $I_R>0$, so that $C_e^2-C_o^2>0$. Hence, the peaks and valleys of $|V(u)|$ are always the same location for any ring of diameter $d$, regardless of its angular profile. The approximation (10) follows from the large-argument expansion of $J_m(\pi d u)$, which is valid only if $\pi d u \gg m$ [11]. This amounts to the condition $u \gg1/\xi_I$, since $\xi_I \approx d/m$ for the $m$th Fourier component: if $I(\phi)$ has substantial contribution from large-$m$ modes, then the universal interferometric signature of the ring must be sampled on longer baselines $u$ in order for the power due to $\xi_I$ to become negligible.

In Figure 1, we compare some toy models of circular thin rings and their visibilities. In the universal regime $\max(1/d,1/\xi_I)\ll u\ll1/w$, the visibility amplitude is given by Eqs. (7), and the peaks and valleys align for all rings. For $I(\phi)=(1+\cos6\phi)I_0/2$, the condition $u\gg1/\xi_I$ is not met for $u\lesssim50$ G$\lambda$, as can be seen in the phase difference between the purple and blue/orange curves in the bottom panel of Figure 1. For $I(\phi)=(1+\cos\phi)I_0/2$, the oscillation in the visibility amplitude vanishes altogether because $I_L=C_e^2-C_o^2=0$.

Effects of finite ring width

We now study the visibility of rings of finite width $w$ in the region where the condition $u \ll1/w$ no longer holds. We find that the ring width may introduce additional $u$-dependent corrections to the oscillating phase of the visibility amplitude, which must be correctly modeled to ensure an unbiased estimate of the ring diameter.

General sliceable rings

The visibility of a general non-axisymmetric thick ring can only be computed after devising some criterion for deciding whether the image is indeed a “ring”; otherwise, the intensity $I(x,y)$ may represent anything and will not necessarily have a visibility whose perturbative expansion in small $uw \ll1$ reduces to the universal form (7) in the limit of zero width $w \to0$. For this reason, in this work, we focus on thick rings that are sliceable:

Definition 1 (Sliceable ring) A (possibly thick) ring is sliceable within the visibility baseline window $(u_-,u_+)$, if it can be decomposed (or “sliced”) into a series of thin convex rings, each of which satisfies the universal regime conditions$^{5}$ $\max(1/d,1/\xi_I) \ll u \ll1/w$ and $I_L I_R > 0$ for any $u \in(u_-,u_+)$.

If a ring is sliceable, then the intensity near the center of the ring should be (at least approximately) zero, or else $1/d \ll u$ cannot hold for the innermost thin ring slice. Definition 1 only gives the necessary condition for a ring to be sliceable. In this work, we do not address the question of what conditions on an image are sufficient for it to be sliceable. We do, however, present a numerical scheme in Sec. III D that shows that GRMHD-simulated models of the time-averaged $n \ge1$ subrings of M87* at 230 GHz are sliceable.$^{6}$

Technically, a slicing scheme is defined via the family of closed convex curves that decomposes the thick ring into thin rings, as illustrated in Figure 2. We take the origin to lie inside the ring and adopt the polar parameterization $r(\varphi,q)$, where $q \in[0,1]$ labels a family of nested curves with $r(\varphi,q_1) > r(\varphi,q_2)$ if $q_1 > q_2$. When projecting onto the axis $\hat{\boldsymbol{u}} \equiv(\cos\varphi_{\mathrm{obs}},\sin\varphi_{\mathrm{obs}})$, the left (right) edge need not be at $\varphi=\varphi_{\mathrm{obs}}+\pi$ ($\varphi=\varphi_{\mathrm{obs}}$), but rather corresponds to the angle $\varphi=\varphi_{L/R}$ that minimizes (maximizes) the projection $\boldsymbol{r}(\varphi,q)\cdot\hat{\boldsymbol{u}}$, which is also the angle where the normal vector to $\boldsymbol{r}(\varphi,q)$ is (anti-)parallel to $\hat{\boldsymbol{u}}$. If a thick ring has surface intensity $I(r,\varphi)$, then the thin ring slice between $q$ and $q+dq$ has an effective line intensity

At the edge points $\varphi=\varphi_{L/R}$, this becomes

where the edge point angles $\varphi_{L/R}$ need to be solved for in terms of the given $r(\varphi,q)$ and $\varphi_{\mathrm{obs}}$ as described above, and $\varphi-\varphi_{\mathrm{obs}}$ is the angle between $\hat{\boldsymbol{r}}$ and the vector normal to the ring. The curvature radius of a ring slice with polar parameterization $r_q(\varphi)=r(\varphi,q)$ is given by

Visibility of a sliceable ring

If a thick ring is sliceable, then each of its thin ring slices has a visibility of the form (6) in the universal regime (8). Summing over the contributions from all the ring slices then produces the visibility of the thick ring. Here, we assume that a viable slicing scheme exists and compute the visibility of the thick ring—we will explore possible strategies for slicing a ring in practice later.

In Sec. II A, we saw that the visibility of a convex thin ring is dominated by the $|x-x_0|^{-1/2}$ singularities in its projected intensity that arise from its left and right edges.

A similar conclusion holds for a sliceable (convex) thick ring, except that each peak in its projected intensity now consists of a sum over the $|x-x_0|^{-1/2}$ weighted by the radial intensity profile. Thus, the visibility of a sliceable thick ring can be calculated by adding the visibility of these two peaks, with each contributing a radial integral over the visibility (5) of a single $|x-x_0|^{-1/2}$ singularity.

Fixing a projection direction $\hat{u}$, for each thin ring slice labeled by $q$, we have an effective line intensity $I_q(q)\,dq$ [defined in Eqs. (12)–(13)], and two edge points $\varphi=\varphi_{L/R}$ where the curvature radius of the ring is $\mathcal{R}(q)=\mathcal{R}(\varphi,q)$ and its projected position is $x(q)$.⁸ We can then define the (one-side) center-of-light radius of the thick ring as

where $\alpha=+1$ for the right side and $\alpha=-1$ for the left side. The center-of-light diameter of the thick ring is then the sum of the left and right center-of-light radii,

With the local coordinate $r'$ defined by $x=\alpha(r_{\mathrm{COL}}+r')$ on each side, the visibility of the left or right peak in the projected intensity profile is (we use eq:36 of [9])

where in the last step, we introduced new functions

Using the moments $M_m$ of the ring profile,

the first of which vanishes ($M_1=0$) by definition of $r_{\mathrm{COL}}$, we can expand the functions $A_{L/R}(u)$ and $\phi_{L/R}(u)$ as polynomials in $u$ with coefficients $a_m$ given in App. B:

The total visibility of a sliceable thick ring is the sum of the left and right peaks,

In the zero-width limit, all the $M_m$ and $a_m$ with $m>0$ vanish and we have $A_{L/R}=I_{L/R}\sqrt{\mathcal{R}_{L/R}}$ and $\phi_{L/R}=0$, so that (22) reduces to (7). When the ring has finite width, we still have $M_1=a_1=0$ due to our choice of $r'$, such that the linear component of the oscillating phase only depends on the center-of-light diameter of the thick ring. However, there are additional corrections to the oscillating envelope involving even powers of $u$, with the leading-order term proportional to $u^2$. Similarly, the oscillating phase includes corrections in odd powers of $u$, with the leading-order term proportional to $u^3$.

When analyzing measured visibility data, these higher-order corrections can be separated from the linear term based on their different $u$-dependence. If $I_q\sqrt{\mathcal{R}}\frac{\partial q}{\partial r'}$ is an even function of $r'$, then $\phi_{L/R}=0$ regardless of the ring width.⁹ In other words, higher-order phase corrections only appear if the ring profile is radially nonsymmetric.

In Eq. (21), we expanded $A_{L/R}$ and $\phi_{L/R}$ as power series in $u$. In practice, we also truncate these infinite series after a few terms, which amounts to approximating the $\sin2\pi r'u$ and $\cos2\pi r'u$ terms in Eqs. (18)–(19) with polynomials. Roughly speaking, for a thick ring of width $w$, this is a good approximation on baselines $u\lesssim1/w$, as opposed to the stricter condition $u\ll1/w$ required by the universal regime. For $u\gg1/w$, the power series in Eq. (21) diverge, although one can still evaluate $A_{L/R}$ and $\phi_{L/R}$ using the original integrals in Eqs. (18)–(19).

Circular slicing of a finite-width ring

If the thick ring can be sliced into circular thin rings, all within the universal regime, then the results from the previous subsection can be recovered from the circular thin ring calculations in Sec. II B. Indeed, let $C_{\mathrm{e/o}}(\rho)$ denote the “radial density” of $C_{\mathrm{e/o}}$, so $C_{\mathrm{e/o}}=C_{\mathrm{e/o}}(\rho)\,d\rho$ for the annulus between $\rho$ and $\rho+d\rho$. By Eqs. (10)–(11), we can then calculate the visibility of the thick ring by summing up the visibility of circular thin rings,

which agrees with Eq. (17).

While the circular slicing works for axisymmetric rings with $I(\rho,\varphi)=I_{\rho}(\rho)$, as well as for separable rings with $I(\rho,\varphi)=I_{\rho}(\rho)I_{\varphi}(\varphi)$ if $I_{\varphi}(\varphi)$ is sufficiently smooth, it is not always a viable strategy for generic rings since the circular thin ring slices may not satisfy the universal regime conditions in Eq. (8). As an example, consider a thick ring with nonzero intensity within the elliptical band $\gamma_{1}^{2}\le x^{2}/a^{2}+y^{2}/b^{2}\le\gamma_{2}^{2}$ ($a>b>0$, $\gamma_{2}>\gamma_{1}>0$) and zero otherwise. Then for $\gamma_{1}b<r<\gamma_{1}a$, the circle of radius $r$ will not lie entirely within the thick ring. It will thus have $I_{\mathrm{L}}=I_{\mathrm{R}}=0$ as well as discontinuities in $I(s)$, and hence its visibility cannot be evaluated using the universal regime solution. However, such a thick ring is actually sliceable, as one can slice it into elliptical thin rings $x^{2}/a^{2}+y^{2}/b^{2}=\gamma^{2}$ with $\gamma_{1}\le\gamma\le\gamma_{2}$, for which the universal regime conditions hold.

Non-circular slicing using radial intensity quantiles

The elliptical band example suggests that, in order to ensure the universal regime criteria are met by each thin ring slice, the slicing should follow the intensity profile of the thick ring. This observation motivates the following slicing scheme based on radial intensity quantiles. Along each ray of fixed $\varphi$, we compute the $q^{\mathrm{th}}$ quantiles $r_q(\varphi)$ of the radial intensity $I(r,\varphi)$ by solving the equation

The $q^{\mathrm{th}}$ thin ring slice then consists of the $q^{\mathrm{th}}$ quantiles at every $\varphi$, or in other words:

Such a construction guarantees that each thin ring slice has a smooth angular profile, as long as the intensity profile $I(r,\varphi)$ of the thick ring is sufficiently smooth.

With this quantile-based scheme, the exact slicing of the thick ring does depend on the choice of origin for the coordinate system, which we typically take to be the origin of the Bardeen coordinates [12]. We stress that the slicing scheme presented here is specifically designed for simulated images, for which the Bardeen coordinates of each pixel are always known.

By contrast, in real observations, the “true” black hole position is not known, but this should not be a concern since one works directly in visibility space. Nevertheless, our experiments indicate that with a slight shift of the origin—e.g., to Bardeen coordinates $(x,y)=(0.5M,0)$—the visibility amplitude predicted by our perturbative calculation remains unchanged.

Numerical experiments

We test our scheme against two numerical simulations of supermassive black hole accretion to gauge how well it performs in the face of complicated accretion physics.

We use the athenak code [13] for the fluid simulations. This code solves the hyperbolic, conservative equations of ideal GRMHD on a logically Cartesian, statically refined mesh in Cartesian Kerr-Schild coordinates. The domains in the simulations extend to $\pm1024\,GM/c^{2}$, have eight total refinement layers, and feature a resolution of 16 zones per $GM/c^{2}$ within a $\pm8\,GM/c^{2}$ cube centered on the black hole.

We perform two simulations in order to cover both the magnetically arrested disk (MAD) as well as the standard and normal evolution (SANE) accretion states. Both simulations assume a black hole spin of $a/M=93.75\%$ and are initialized from a Fishbone-Moncrief [14] torus configuration, with the inner edge at $r=20\,GM/c^{2}$ and the pressure maximum at $r=41\,GM/c^{2}$. We set the fluid adiabatic index to $13/9$ globally across the flow.

We raytrace the simulations with a modified version of the blacklight code [15] that allocates more pixels within the lensing bands where small scale demagnified structures are expected. The details of our nested ray tracing scheme can be found in [16]. For post-processing, we make use of the $(R_{\mathrm{low}},R_{\mathrm{high}})$ model [17] to compute the electron temperature $T_{\mathrm{e}}$ from GRMHD fluid pressure $p$, density $\rho$ and plasma $\beta$,

where

The camera lies at $r_0 = 100 r_g$, $\theta_0 = 163^\circ$ and $\phi_0 = 0^\circ$. We set the black hole mass to $M = 6.5 \times10^9 M_\odot$ and use $1M = 3.8\,\mu\mathrm{as}$ to convert from Bardeen’s coordinates to angular distances, corresponding to the mass-to-distance ratio for M87*. We get the image of the $n^{\mathrm{th}}$ photon ring by taking the difference between the ray tracing starting from the $(n+1)^{\mathrm{th}}$-to-last turning point¹⁰ on the geodesic, and the ray tracing starting from the $n^{\mathrm{th}}$-to-last turning point. Then, we average over 500 snapshots at a cadence of $10M$ and calculate the “true” visibility using the Fast Fourier Transform (FFT).

The averaged images and projected intensity profiles corresponding to $(R_{\mathrm{low}}, R_{\mathrm{high}})$ equal to $(1, 10)$, $(1, 100)$, $(10, 10)$ and $(10, 100)$ are shown in Figure 3 and Figure 9 for the MAD and SANE simulations, respectively. Though the exact electron temperature prescription for M87* is still uncertain, we note that the two intermediate models with $(1, 100)$ and $(10, 10)$ are in generally better agreement with the current EHT observations (see, e.g., [3]).

Rather than treating the $a_m$ as free parameters and fitting the visibility amplitude with (22), we directly apply our slicing scheme in Sec. III D to the coherently averaged images, and evaluate the $a_m$ from the moments $M_m$ of the ring profile. As further discussed in App. F, instead of calculating the moments (20) with the integral running over the whole range $[0, 1]$ for $q$, we truncate it to $[0.1\%, 99.9\%]$ so as to mitigate the instabilities caused by the tail of the integrand. We then compare the $|V(u)|$ calculated via (22) up to the $m^{\mathrm{th}}$ order and the actual visibility amplitude of the image in Figure 4 for the MAD simulation, with $m = 1$ representing the universal regime formula. Since in practice, we use the perturbative phase $\phi$ for the analysis of visibility data, as discussed in Sec. V, we include a case where both the envelope and phase are calculated with perturbation theory, and also a case with only the phase computed perturbatively. For the latter, we use the method in Sec. D to numerically decompose the true visibility into $A(u)$, $B(u)$ and $\phi(u)$ as follows:

The envelope of $|V(u)|^2$ can then be obtained from $A(u)$ and $B(u)$. As shown in Figure 4, the peaks fade when the perturbative $m = 1$ approximation to $|V(u)|$ misalign with the true visibility due to the error in the phase at $u \sim100\,\mathrm{G}\lambda$, which is improved by the $m = 5$ model. Moreover, the phase is more amenable than the envelope to perturbative treatment, since for $m = 5$, the envelope deviates from the truth at around $40$-$50\,\mathrm{G}\lambda$, even though the phase is still accurate up to about $80$-$100\,\mathrm{G}\lambda$.

We also compare the true $\phi(u)$ to the perturbatively predicted phase at $m^{\mathrm{th}}$ order with $m = 1$ and $m = 5$, and plot their difference in Fig. 5 for the MAD simulation, which clearly shows that $m = 5$ improves the accuracy of the $\phi$ modeling. An error in the center-of-light diameter leads to a phase error $\Delta\phi\equiv\phi_m - \phi_{\mathrm{true}}$ that is linearly dependent on $u$ [see Eq. (22)], and the boundary of the green bands in Fig. 5 represents the $\Delta\phi$ corresponding to a $0.2\,\mu\mathrm{as}$ difference in the diameter. Therefore, should there be no other forms of noise, the systematic error of the center-of-light diameter measurement will be less than $0.2\,\mu\mathrm{as}$ if (and only if) $\Delta\phi$ falls within the green bands. We find that the phase from the $m = 5$ model closely matches the truth at shorter baselines $u \lesssim50\,\mathrm{G}\lambda$, and the expected diameter measurement error remains less than $0.2\,\mu\mathrm{as}$ at $u \lesssim80\,\mathrm{G}\lambda$ for all four of our electron temperature models. By contrast, with the $m = 1$ model, we expect a systematic error of $0.2\,\mu\mathrm{as}$ in the center-of-light diameter, even when the measurement is conducted on shorter baselines $u \sim40\,\mathrm{G}\lambda$.

In Fig. 6, we show the analogous results for the SANE simulation, in which we find that the ring is generally thicker (see Fig. 7), therefore becoming non-perturbative on shorter baselines. More detailed discussion regarding the SANE results can be found in App. E.

Discussion

In this paper, we investigate the visibility of black hole photon rings beyond the universal regime by accounting perturbatively for the effects of finite ring width. Here, we focus on sliceable thick rings (Definition 1), which can be decomposed into multiple thin ring slices such that the visibility of the thick ring is the sum of the visibilities of the thin rings. In this paper, we do not determine the sufficient conditions on a ring image for it to be sliceable. However, we do show that photon ring images ray traced from GRMHD simulations can be sliced using the scheme in Sec. III D based on radial intensity quantiles.

The main contributions of this work are that:

1. We show that the exact diameter measured for sliceable thick rings from the oscillation of their visibility amplitude is precisely the center-of-light diameter defined in (16).

2. We compute higher-order corrections (polynomial in $u$) to the visibility beyond the universal regime. We obtain baseline-dependent corrections $\propto u^3$, $u^5$, etc., to the oscillation in the visibility amplitude. Incorporating these higher-order terms is crucial for accurately inferring the part of the oscillation phase that is linear in the baseline $u$. It is this linear term that encodes the diameter of the ring, which in turn depends on the black hole mass and spin.

The panels in Figs. 5 and 6 give an intuitive estimate of the systematic error introduced by the thickness of the $n=1$ photon ring in the inference of its center-of-light diameter. We find that this error stays below $0.2\,\mu\mathrm{as}$ for $u\lesssim80\,\mathrm{G}\lambda$ in MAD simulations and for $u\lesssim40\,\mathrm{G}\lambda$ in SANE simulations, provided one includes the additional perturbative corrections to the oscillating phase of the visibility amplitude calculated here. As a comparison, the typical width $w$ of the first lensing band for M87* is about $6\,\mu\mathrm{as}$ given our assumed mass-to-distance ratio, corresponding to a baseline length of $u\sim1/w=34.4\,\mathrm{G}\lambda$.

This roughly sets a lower bound (among different plasma prescriptions) on the baseline $u$ past which we expect our perturbative approach to breakdown. We find that the MAD models remain perturbative up to significantly longer baselines because the $n = 1$ image is significantly narrower than the total width of the first lensing band, which is not as true for the SANE models.

Recent work [18] exploring diameter-fitting for thick rings from their interferometric visibility, and applied the fitting procedure previously derived for a zero-width ring [8] without theoretically justifying the unbiasedness of the method when applied to thick rings. This approach results in different diameters when fitting within different baseline windows, likely because it does not account for the higher-order corrections to the phase of the visibility amplitude discussed herein.

The visibility of a ring fails to display the universal regime behavior (7) if any of the conditions in Eq. (8) are violated. Our analysis only addresses the change in the visibility of the ring due to finite ring width when the condition $u \ll1/w$ is not met. Our method may still fail (1) if the thick ring is not sliceable, that is, if it cannot be decomposed into multiple thin rings for which Eq. (8) holds, or (2) if the ring is sliceable but the slicing is not perturbative, meaning that the polynomial expansion in Eq. (21) fails at large $u$ of interest. It is not clear how to extend our perturbative calculation of the visibility when $u \gg1/d$ or $u \gg1/\xi_I$ does not hold, since the visibility then depends on the global intensity profile of the ring (rather than being dominated by its left and right edges).

We find that the inclusion of polynomial corrections up to fifth order to the oscillating visibility amplitude is sufficient for measurements on baseline lengths $u \lesssim40$G$\lambda$ likely achievable even with space-VLBI in the near term. However, our perturbative approach inevitably breaks down on yet-longer baselines $u \gtrsim1/w$, as seen in Figs. 5, 6 and 10.

It may be possible to further extend the convergence interval by, e.g., employing a different slicing scheme or rearranging higher-order terms via Padé resummation (see, e.g., Ref. [19]). Empirically, we find that if the ring is either nonsliceable or nonperturbative, then it is typically not possible to fit the oscillating phase of its visibility amplitude with a simple polynomial in $u$. This observation suggests that, in practice, if one can measure the visibility within a baseline window that is sufficiently large to separate the different $u$-dependent terms, and one finds that the model achieves a good fit to the data, then it is likely that the ring diameter inferred from such a fit is correct and unbiased. This conjecture will be the subject of more detailed tests on simulated data.

There are additional effects that should be accounted for before the models proposed herein can be applied to future observations; several of these will be addressed in a companion paper [16].

First, while we calculated the visibility amplitude of a single sliceable thick ring in Sec. III, realistic black hole images include multiple nested photon rings due to the lensing of Kerr spacetime [5], which in principle should be modeled as the sum of multiple sliceable rings. But, since these rings dominate the visibility on different baselines, the total visibility on a given range of baselines can typically be well-approximated with that of a single ring [4], provided one excludes the transition regions between regimes dominated by different rings [8, 20].

Second, due to the difficulties inherent in VLBI phase calibration, we can only measure the averaged visibility amplitude of the image (the incoherent average), and not the visibility amplitude of the time-averaged image (the coherent average) studied in this work. Even though an incoherent average changes the envelope of the visibility amplitude (relative to a coherent average), nevertheless the shift of the phase becomes negligible after averaging over sufficiently long timescales, such that one can still infer the ring diameter from the phase.

In addition, to obtain the visibility across different baselines, one needs to vary the orbit of a space-orbiting satellite or sample the visibility at different frequencies. In the latter case, one would need to model the frequency dependence of the image. Encouragingly, however, our analysis suggests that within the frequency window likely to be explored by forthcoming space-VLBI missions (e.g., between 230 GHz and 345 GHz), the black hole image only exhibits minor variation, while the center-of-light diameter $d_{\mathrm{COL}}$ and the coefficients $a_m$ can be treated as either constant or linearly dependent on the observing frequency, thereby enabling a straightforward modeling of the frequency dependence of the image.

Acknowledgements

We thank Alejandro Cárdenas-Avendaño, Andrew Chael, and Sihao Cheng for helpful discussions. The work presented in this article was performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. This work was supported in part by a Simons Investigator Grant from the Simons Foundation (EQ) and by NSF grants AST-2307888 and PHY-2340457 (AL). It greatly benefited from AL and EQ’s stay at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-15 2210452, and by EQ’s stay at the Kavli Institute for Theoretical Physics, supported by NSF PHY-2309135.

References

Visibility of a circular thin ring

In this section, we demonstrate how to derive (10) from (9). When $\pi d u \gg m$, the $J_m(\pi d u)$ terms in (9) can be approximated with trigonometric functions using the large-argument approximation of Bessel functions, namely

Note that $(-i)^{-m}=(-1)^m \times(-i)^m$, $J_{-m}(\pi d u)=(-1)^m \times J_m(\pi d u)$, and for real ring images, $\operatorname{Re}\beta_{-m}=\operatorname{Re}\beta_m$, $\operatorname{Im}\beta_{-m}=(-1)\times\operatorname{Im}\beta_m$. Therefore, we have¹¹

The $\operatorname{Im}\beta_m$ terms vanish in the summation, since the projection of a sine ring along the $x$-axis is zero. Taking the amplitude of (28) then yields (10).

Coefficients $a(m)$

Here, we list the expressions for the expansion coefficient $a_m$ introduced (21), up to $m=13$:

Comparison between different perturbative orders

In Figure 8, we compare the perturbative solution to the visibility amplitude up to different orders. Although $m\geq8$ improves the envelope prediction up to about $80\,G\lambda$, $m=5$ offers the most accurate results for the phase. Since in practice, the envelope will be modeled with smooth spline functions instead of perturbative polynomials due to the noise from incoherent averaging [16], we recommend choosing $m=5$ for a better modeling of the phase.

The envelope and phase of $|V(u)|^2$

We use the following approach to numerically decompose $|V(u)|^2$ into $A(u)$, $B(u)$ and $\phi(u)$ as in (27), which assumes that $A(u)$, $B(u)$ and $\phi(u)$ are all smooth functions of $u$ and $\phi(0)=0$. We first find the local maxima and minima of $|V(u)|^2$; at these points, $\phi$ should be close to, but slightly smaller than, $\left(\frac{1}{2}+2i\right)\pi$ and $\left(\frac{3}{2}+2i\right)\pi$, if $A(u)$ and $B(u)$ monotonically decrease with $u$. We obtain an initial decomposition by setting $\phi$ to $\left(\frac{1}{2}+2i\right)\pi$ or $\left(\frac{3}{2}+2i\right)\pi$ at the local extrema of $|V(u)|^2$, solving for envelope of $|V(u)|^2$ at these extrema, and then interpolating $A(u)$, $B(u)$ and $\phi(u)$ between the extrema using cubic splines. We then iteratively refine the decomposition by finding the local extrema of $\left[|V(u)|^2-A(u)\right]/B(u)$ to repeat this process, which corrects the aforementioned bias if one naively sets the phase to $\left(\frac{1}{2}+2i\right)\pi$ or $\left(\frac{3}{2}+2i\right)\pi$ at the local extrema of $|V(u)|^2$, and typically converges within about 5 iterations.

Results for the SANE simulation

We show the results for the SANE simulation in Figure 9 and Figure 6: while the $m = 5$ model still consistently outperforms the $m = 1$ model for the visibility modeling, they both diverge on shorter baselines than in the MAD simulation with the same electron temperature model. This is because the $n = 1$ subimages generally have “heavier tails” (i.e., more intensity well away from the peak of the ring subimage) in SANE models compared with MAD models, as demonstrated in Figure 7. In other words, the effective width of the $n = 1$ ring is larger in the SANE models, making the perturbative approach break down on smaller baselines.

However, we note that the subimage being non-perturbative does not necessarily mean that the total visibility does not follow our polynomial model, as the approach to decompose the total image into subimages is not unique. At baselines where the $n = 1$ subimage dominates the total visibility, while the $n = 1$ subimage and $n \le1$ subimage (i.e., the sum of $n = 0$ and $n = 1$) have the same visibility, the former will be more perturbative than the latter, as the latter also includes the long tail due to the $n = 0$ component. While this $n = 0$ tail has negligible contribution to the total visibility on longer baselines, it can still render the total ring non-perturbative since one can no longer expand the trigonometric functions as polynomials in the calculation. Therefore, throughout this paper, we apply our perturbative method to the $n = 1$ subimage, defined via the $\mathrm{d}z/\mathrm{d}\lambda= 0$ turning points. As the decomposition is not unique, it may be possible to define an alternative decomposition scheme, with which the $n' = 1$ subimage has the same visibility (within some desired baseline window) but even shorter tails than the standard $n = 1$ subimage, so that the ring remains perturbative on longer baselines. Nevertheless, since the standard subimage decomposition seems sufficient for the visibility modeling below $40\,G\lambda$ for all the plasma models considered in this work, we leave such an alternative decomposition scheme for future research.

Convergence with respect to the range of radial integration

In Figure 10, we compare the visibility phase calculations with the integral in (20) running over $0 \leq q \leq1$ (the whole ring) versus $0.1\% \leq q \leq99.9\%$; the latter effectively truncates the ring at the 0.1% radial intensity quantiles on each side. We find that while the $m = 1$ results remains indistinguishable, the $m = 5$ results on longer baselines do depend on the range of the radial integral: due to the $(r^{\prime})^{m}$ factor in (20), the higher order moments will be more sensitive to the tails of the ring intensity profile. Throughout this paper, we adopt the $0.1\% \leq q \leq99.9\%$ truncation as our default choice, which slightly improves the accuracy of the fifth order visibility modeling.