The Photon Ring in M87*

Radio emission from the immediate vicinity of the supermassive black hole M87* was used to reconstruct the first-ever image of a black hole, estimate its mass, and interpret its theoretical environment by the Event Horizon Telescope (EHT) Collaboration (Event Horizon Telescope Collaboration et al. 2019a, 2019b, 2019c, 2019d, 2019e, 2019f, hereafter Papers I-VI).

The analyses behind these findings included multiple image-reconstruction algorithms, model fitting in the visibility domain employing various geometric shapes, and an extensive investigation of physical emission models. The latter are based on a large library of simulations that model the source as a hot, magnetized accretion flow, that form the input for theoretical model images constructed via ray-tracing and solving the equations of radiative transfer for a thermal population of relativistic electrons emitting synchrotron radiation. Typically, these models predict images with both a sharp ring component—i.e., associated with the location of photon rings in the underlying spacetime—as well as a comparatively diffuse but still compact emission structure (see, e.g., Paper V).

Image reconstructions, model fitting to geometric shapes, and direct fitting to general relativistic magnetohydrodynamical (GRMHD) synchrotron models all yielded an inferred mass for M87* of (6.5 ± 0.7) × 10⁹ M_⊙, which agrees with stellar-dynamical mass inferences of $({6.14}_{-0.62}^{+1.07})\times {10}^{9}{M}_{\odot }$ (Gebhardt et al. 2011; Paper VI). The error in each method was estimated based on the variable emission structure from the GRMHD-based model images.

In the meantime, novel imaging schemes have been developed that can either approximate (Arras et al. 2019; Sun & Bouman 2020) or directly sample (Broderick et al. 2020b; Pesce 2021) the posterior distribution over possible image structures. The Bayesian nature of these schemes yield meaningful posterior distributions for the images. These posteriors permit a more rigorous characterization of the credibility of image features and a measure of image consistency. In addition, they permit a hybrid approach that combines image reconstruction with modeling specific expected features. This is demonstrated in Broderick et al. (2020b), where imaging is accomplished with a "nonparametric" model comprised of a rectilinear raster of control points, and additional geometric components (Gaussians, rings, etc.). Of particular interest here is the ability to reconstruct ring features in the data beyond the diffraction limit, provided the signal-to-noise ratios (S/Ns) of the measured complex visibilities are sufficiently high.

Such ring features are theoretically expected outcomes in black hole images due to the propagation of photons in close proximity to the black hole, and the associated strong gravitational lensing predicted by general relativity. It is useful to think of the resulting image that is measured at infinity to be composed of a direct emission component that is dominated by the typically nontrivial and uncertain astrophysical environment, plus a series of ring components that are confined to distinct narrow regions (far less influenced by the details of the overall flow structure), which arise from photons that traversed the black hole once, twice, and so on before reaching a distant observer (Bardeen 1973; Broderick et al. 2022). Thus far, the M87* images present the total image structure, i.e., the sum of all of those components; the diffraction-limited image reconstructions in Paper IV cannot distinguish the ring components from the rest of the emission directly.

In all of the theoretical models presented in Paper V that were consistent with the data, the direct emission provided the majority of the flux, with a substantial minority arising from the first higher-order image component; in what follows we will refer to these as the n = 0 and n = 1 "photon rings," respectively, despite the former often not forming a ring. Based on inspection of the GRMHD library presented in Paper V, the n = 1 ring typically contains 10%–30% of the total compact flux, depending on the model (e.g., Gralla et al. 2019; Johnson et al. 2020). Here, we demonstrate that statistically preferred reconstructions are achieved when a thin ring is included in addition to the standard nonparameteric image-reconstruction component. The thin ring is then identified with the n = 1 photon ring of the underlying spacetime. Interpreting the separate emission components as the n = 1 and n = 0 photon rings provides a powerful mapping of spacetime that is far less sensitive to the astrophysical processes in the emitting region than previous analyses (Broderick et al. 2022).

In this paper we report the results of applying the hybrid image model from Broderick et al. (2020b) to the M87* 2017 EHT data set. By fitting a geometric ring model component simultaneously with a flexible image component, we are able to isolate the n = 1 photon ring from the surrounding diffuse emission. We demonstrate that our reconstructed emission structure is in agreement with prior EHT imaging results, and the properties of the n = 1 ring yield novel constraints on the black hole mass with significantly reduced systematic uncertainties. Furthermore, because the bright ring emission is effectively extracted by the geometric model component, the image component is free to capture subtler details within the remaining diffuse emission than would otherwise be possible. We find that the removal of the bright foreground ring uncovers additional low-brightness image structures that are consistent with originating from the base of the forward jet; the observed brightness asymmetry in context with previous EHT results constraining the black hole spin orientation Paper V aligns with expectations for a jet driven by the Blandford–Znajek mechanism (Blandford & Znajek 1977).

The structure of the paper is as follows. In Section 2 we summarize the algorithm used to reconstruct the images. We present the reconstructed, structure-agnostic images in Section 3, before presenting the hybrid reconstructions in Section 4. In Section 5 we describe how our findings are related to the black hole parameters. The implications of the structure and evolution of the diffuse emission within the broader context of M87* are collected in Section 6. We summarize and conclude in Section 7.

We employ the forward-modeling Markov Chain Monte Carlo (MCMC) algorithm presented in Broderick et al. (2020b) and implemented in the Themis analysis framework (Broderick et al. 2020a). In this scheme, the image is forward-modeled by a rectilinear grid of control points, at which the intensity is a free parameter, and between which the intensity is modeled via an approximate cubic spline. Station gains are simultaneously modeled, and are assumed to be fixed in time over a single scan. For the current work, this algorithm has been supplemented in four ways.

First, we have included the field of view (FOV) as two additional parameters in the underlying image models (one for each dimension in the image plane). This change promotes the two FOVs from hyperparameters that must be manually surveyed to parameters that are continuously explored in a fashion more consistent with the general Bayesian approach employed. It also permits efficient exploration of asymmetric FOVs, i.e., different FOVs in the two directions of the rectilinear grid. While this flexibility makes only a small difference for the results presented here, it does enable analysis of highly asymmetric and extended systems (e.g., active galactic nuclei jets).

Second, we permit a rotation of the rectilinear grid. This change permits the control points to optimally arrange themselves within the confines of the grid. Typically, this freedom permits smaller grid dimensions than grids that are fixed to be oriented along the cardinal directions. Again, this additional flexibility makes only a modest difference in the applications here.

Third, we make use of an updated set of samplers implemented within Themis. Details on the sampler, including demonstrations, may be found in Tiede (2021). In summary, the improved sampler makes use of a deterministic even–odd swap tempering scheme (Syed et al. 2019) using the Hamiltonian Monte Carlo sampling kernel from the Stan package (Carpenter et al. 2017). MCMC chain convergence is assessed using standard criteria, including integrated autocorrelation time, approximate split $\hat{R}$, and visual inspection of individual traces.

Fourth, we fit to complex visibilities rather than some combination of visibility or closure amplitudes and phases. This requires that we reconstruct the time-variable station gain phases in addition to their amplitudes. Both are stable across a 10 minute scan. Fitting the complex visibilities simplifies the treatment of the errors at low S/N: thermal errors are strictly Gaussian (Thompson et al. 2017). It also improves the structure of the likelihood surface to be explored, which is both smoother and has fewer modes in this case.

Fit quality is assessed using the χ² statistic, comparing log-likelihoods, and by inspecting the distribution of residuals.

Prior to applying the complete hybrid imaging model, we first perform nonhybrid image reconstructions to more easily enable comparison with previous imaging work and to provide context for subsequent interpretation of the hybrid images.

We reconstruct images for each of the four EHT observations of M87* taken in 2017 April. The observations were carried out using two frequency bands (Paper III), which we refer to as a high band and a low band, and we image both bands simultaneously. We assume that the image structure is shared across bands, such that a single image is produced for each day, but we permit the station gains to be completely independent among bands. Prior to fitting, we preprocess the visibility data as described in Broderick et al. (2020b): we coherently time-average the complex visibilities within each scan, and we add a systematic uncertainty of 1% in quadrature to the thermal uncertainties to account for residual calibration errors. This additional error budget is motivated by the analysis of nonclosing errors in Paper III and the same as that found in Papers V and VI sufficient to produce high-quality fits. The model we employ uses a 5 × 5 grid of control points to capture the image structure, and a large-scale asymmetric Gaussian (a major-axis FWHM above 0.2 mas; see Appendix A for details) to accommodate structure seen on only the shortest baselines (again utilizing the hybrid modeling+imaging approach).

Table 1 provides an accounting of the parameters and data quantities used for each of the fits, and it lists various fit statistics for each of our reconstructions. The quantities most comparable to the fit statistics reported in Table 5 of Paper IV and Table 2 of Paper VI would be the χ²/(N_data,HI + N_data,LO), which here range from 0.58 to 0.93 and are otherwise comparable to those reported previously. Direct comparisons with the complex visibility data and the corresponding residuals are shown for representative days in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Fit Quality Assessment

Day	Model a	Nparams	Ndata, HI b	Ndata, LO b	Ng,HI b	Ng,LO b	${\chi }_{\mathrm{HI}}^{2}$	${\chi }_{\mathrm{LO}}^{2}$	χ2 c	ΔBIC d	ΔAIC e
April 5	I5×5+A	34	336	336	162	162	295.2	246.4	628.0	⋯	⋯
	I5×5+A+X	41	336	336	162	162	219.6	196.2	454.7	−127.7	-107.4
April 6	I5×5+A	34	548	568	226	243	436.3	376.2	872.1	⋯	⋯
	I5×5+A+X	41	548	568	226	243	356.8	344.0	756.2	−66.8	-68.9
April 10	I5×5+A	34	182	192	79	73	87.9	96.7	216.5	⋯	⋯
	I5×5+A+X	41	182	192	79	73	80.5	93.2	194.0	18.9	35.5
April 11	I5×5+A	34	432	446	185	190	388.9	338.1	800.6	⋯	⋯
	I5×5+A+X	41	432	446	185	190	365.6	313.6	742.6	−10.5	-8.0

Notes.

^a Model components are as follows: an N × M dimensional image raster (I_N×M ), a large-scale asymmetric Gaussian (A), and a slashed ring (X). Detailed descriptions of the model components and priors can be found in the main text and Appendix A. ^b Each complex visibility is counted as two data points (N_data,HI/LO). Similarly, each complex gain is counted as two gain parameters (N_g,HI/LO). ^c Includes contributions from the gain priors. ^d Differences between the Bayesian information criterion (BIC) with and without the slashed ring. The BIC is defined by ${\chi }^{2}+k\mathrm{ln}(N),$ where k ≡ N_params + N_g,HI + N_g,LO is the total number of model parameters and N ≡ N_data,HI + N_data,LO. ^e Differences between the Akaike information criterion (AIC) with and without the slashed ring. The AIC is defined by χ² + 2k + 2k(k + 1)/(N − k − 1), where k and N are defined as they are for the BIC.

Download table as:  ASCIITypeset image

Figure 2 shows the resulting image reconstructions for each of the four days. The top row shows the maximum-likelihood samples from each chain, the middle row shows the posterior means, and the bottom row shows the posterior standard deviations. We find that the posterior mean images show a qualitatively similar ringlike structure to image reconstructions produced using regularized maximum-likelihood methods in Paper IV, though we note that we have not imposed any comparable regularization (e.g., maximum entropy, total variation) on our likelihood function. The general shape, size, and total flux of the emission structure are similar across all four days, as is the pronounced north–south asymmetry in the brightness distribution. We find that the image control point raster prefers a FOV of ∼50 μas in both axis directions and a modest rotation with respect to the equatorial coordinate system (see Appendix A).

Zoom In Zoom Out Reset image size

The bottom row of Figure 2 illustrates a measure of the uncertainty in the image reconstructions, as previously demonstrated in Broderick et al. (2020b). We find that the uncertainty is not uniform across the image, nor does it seem to be proportional to the image intensity. Rather, the uncertainty tends to be lowest within an approximately circular region running azimuthally around the ring, and it increases both radially inwards and outwards of this region. We note that this behavior contrasts with the appearance of the image uncertainty reported in Sun & Bouman (2020); we do not explore these differences in this paper, but we expect that they arise primarily from the large differences in likelihood and prior specification between these two algorithms.

The image model described in the previous section is incapable of reconstructing features on scales much smaller than the raster spacing, which is comparable to the nominal array resolution of ∼20 μas. However, the lowest-order lensed image around the black hole—i.e., the n = 1 photon ring—is expected to have a thickness of only ∼1 μas (Johnson et al. 2020). While we may not expect to be able to spatially resolve the thickness of this ring with the 2017 EHT array, its ∼40 μas diameter should still imprint itself on the visibility data. By enforcing the prior expectation that a putative n = 1 component of the observed emission structure should originate from a thin ring (i.e., one having a thickness that is much less than its diameter), Broderick et al. (2020b) showed that the diameter of the photon ring could be reliably recovered from EHT-like synthetic data sets generated from input images produced from GRMHD simulations.

Following Broderick et al. (2020b), we perform hybrid imaging of the four M87* data sets, in which we fit for a "slashed ring" model component alongside the image and large-scale Gaussian described in the previous section. For the M87* black hole, which has a spin-axis inclination of ≲20° with respect to the line of sight (Walker et al. 2018; Paper V), the photon ring is expected to have a nearly circular geometry with only small (≲2%; Johnson et al. 2020) deviations from circularity even for large spin values. We thus model the ring as a thin (fractional thickness less than 5%; see Appendix A) circular annulus with a linear brightness gradient ("slash"), and we permit the diameter, flux, and slash magnitude and orientation to be free parameters. The thickness parameter is restricted by a tight prior that forces it to be <5% of the diameter. Additionally, we permit the center coordinates of the ring model component to drift with respect to the center of the image model component. We describe the model and prior distribution specification in more detail in Appendix A. Direct comparisons with the complex visibility data and the corresponding fit residuals for representative days are presented in Figure 3, and indicate high-quality fits across the entirety of the range of baselines probed by the EHT.

Zoom In Zoom Out Reset image size

As tabulated in Table 1, we find lower χ² values when fitting a hybrid image to the data relative to fitting the image model described in Section 3. For comparison with Table 5 of Paper IV and Table 2 of Paper VI, the analogous χ² quantities to those reported there range from 0.52 to 0.86, again a roughly comparable fit quality. Such improved fit quality is expected given the increased complexity of the hybrid image model, but information criteria considerations indicate that the fit improvement outweighs the additional model complexity for all but the April 10 data set. The April 10 data set is the sparsest of the four, containing a factor of ∼2–3 fewer data points than any of the other data sets, and the information criteria indicate that the increased complexity of the hybrid images is not statistically necessary in this case. ⁴⁹ However, we note that the statistical preference for a ring component in the other three data sets—for which the coverage is more complete and the emission structure correspondingly better constrained—implies its presence in the April 10 data set, as well.

Our hybrid image reconstructions for each of the four days are shown in Figure 4. We find that after convolution with a 15 μas Gaussian beam, the gross structural properties of the emission qualitatively match those recovered from the imaging in Section 3 (see Figure 2) and image reconstructions in Paper IV, though with an evident preference for a smoother and more uniformly circular emission structure in the hybrid images. The ability of the thin ring model component to capture aspects of the source structure frees the remaining image model component to devote its flexibility to recovering fainter features. This increased focus on fainter structures can also be seen in the behavior of the image control point raster, which shows a broader east–west extent for the hybrid image fits than for the image-only fits (see Appendix A).

Zoom In Zoom Out Reset image size

Figure 5 shows the hybrid image reconstructions with the ring and image components separated. The slashed ring model component is shown at nearly its native resolution; it is smoothed by 0.5 μas for visualization purposes only. The diffuse emission map is displayed without any additional convolution (in contrast to Figures 4 and 2). In this figure, the slashed ring model component is visually distinct from the much-lower-brightness diffuse emission associated with the image model component.

Zoom In Zoom Out Reset image size

Though the model permits both the ring and image components to drift freely with respect to one another, we find that the data prefer reconstructions in which the two components are nearly concentric. We see that the diffuse emission is primarily concentrated along the southern portion of the ring, helping to define several "knots" of emission; a by-eye decomposition indicates that there are approximately two such knots on April 5 and 6, with a third knot appearing at the southernmost point of the ring on April 10 and 11. Additionally, we find that all four days show a feature that is significantly detected to the southwest of the ring. This southwestern component is present in some of the reconstructions from Paper IV and features more prominently in Arras et al. (2022) and Carilli & Thyagarajan (2022), but only the latter two works provided any comment; we discuss the potential origin of and implications for this feature in Section 6.2.

The reconstructed flux densities in both the image and ring model components are listed in Table 2, and we find that the ring component contains between ∼54%–64% of the total flux in the image. This range exceeds the fraction of flux contained in the narrow ringlike features in GRMHD simulations, from which we anticipate only ∼10%–30% of the total image flux to be contained in the narrow ring. The measured fluxes are, however, consistent with the results from Broderick et al. (2020b) for hybrid image reconstructions of simulated data. The excess ring flux appears to be a consequence of the absorption of a portion of the surrounding direct emission into the ring component. By virtue of their sparse (u,v) coverage and finite S/N, the EHT observations have an effective angular resolution limit of approximately ∼10 μas (Paper IV), which is smaller than the nominal beam size of ∼20 μas but still much larger than the anticipated n = 1 ring thickness of ∼1–2 μas. Image structures on scales smaller than this ∼10 μas threshold are effectively unresolved by the array. Our hybrid image model priors confine the ring thickness to be ≲5% of the ring radius (roughly 1 μas; see Appendix A), but source flux contained within an annulus of thickness ∼10 μas around the ring radius is structurally indistinguishable from flux residing within the ring itself. In GRMHD simulations, such an annulus contains roughly ∼50%–80% of the total source flux, consistent with the results in Table 2. Broderick et al. (2020b) have demonstrated that this excess flux capture does not appear to substantially bias the recovered n = 1 ring radius.

Table 2. Hybrid Image-ring Fit Parameters

Day	Idiff (Jy)	θdiff (μas)	Iring (Jy)	θring (μas)
April 5	${0.252}_{-0.017-0.044}^{+0.019+0.045}$	${17.0}_{-1.2-2.6}^{+1.1+3.4}$	${0.301}_{-0.007-0.017}^{+0.007+0.022}$	${21.88}_{-0.14-0.37}^{+0.13+0.36}$
April 6	${0.190}_{-0.012-0.041}^{+0.014+0.037}$	${20.2}_{-1.3-5.2}^{+1.0+2.6}$	${0.276}_{-0.008-0.019}^{+0.010+0.049}$	${21.44}_{-0.10-0.23}^{+0.09+0.21}$
April 10	${0.176}_{-0.019-0.052}^{+0.017+0.047}$	${20.4}_{-4.7-6.7}^{+4.2+7.0}$	${0.302}_{-0.021-0.054}^{+0.016+0.038}$	${21.89}_{-0.39-0.90}^{+0.27+0.51}$
April 11	${0.193}_{-0.018-0.070}^{+0.013+0.035}$	${20.8}_{-0.8-6.0}^{+0.7+2.8}$	${0.246}_{-0.013-0.026}^{+0.020+0.058}$	${22.51}_{-0.17-0.47}^{+0.16+0.37}$

Note. Values quoted are the median, 50-, and 90-percentile ranges.

Download table as:  ASCIITypeset image

Figure 6 shows the posterior distributions for the ring component radius parameter on all four days. The radius measurements are consistent with being constant across the week; the apparent evolution from April 6 to April 11 is a modest 2σ tension after inclusion of the appropriate trials factor. Combining the ring radius measurement across all days yields θ_ring = 21.74 ± 0.10 μas. Posteriors for other ring properties (flux, position angle, thickness) can be found in Appendix B, and are similarly consistent among days. The radius of the direct emission region, θ_diff, is measured from the ring center by computing the radial location of the brightness peak in the image component. Though we note that the direct emission does not necessarily form a closed ring on any day, this definition nevertheless presents the most consistent measure of direct emission ring size to that invoked in Broderick et al. (2022). We find that the direct emission radius is a factor of ∼10 more uncertain than that of the ring component; we generate posteriors for θ_diff on each day from the hybrid image reconstructions, and the resulting values are listed in Table 2.

Zoom In Zoom Out Reset image size

Most directly constrained, and most comparable to other measurements of the mass, is the angular scale θ_M ≡ GM/c² D, where D is the distance to M87*. This is related to the angular radius of the bright ring. However, it is subject to additional systematic uncertainties associated with the nature of this relationship, dominated by a systematic bias associated with the location of the emission region and the dependence on black hole spin.

Here we consider three methods for estimating θ_M and the corresponding mass, M = M₉ × 10⁹ M_⊙. All of these identify the bright ringlike structure in the image with the n = 1 photon ring, produced by photons that execute a half-orbit about the black hole prior to reaching the distant image plane. This is motivated both geometrically—higher-order photon rings are suppressed exponentially (Johnson et al. 2020)—and from astrophysical predictions ranging from semi-analytical modeling (Broderick et al. 2016) to GRMHD simulations (Paper V; Gammie et al. 2003; Abramowicz 2013; Porth et al. 2019). They differ in the manner in which they attempt to systematically address the dependence on spin and the relationship to the asymptotic (n = ∞ ) photon ring, which defines the edge of the black hole shadow.

In all cases, where we transform from an angular measurement of the mass, θ_M , to a physical measurement, we will assume a distance of D = 16.8 ± 0.8 Mpc (see Appendix I of Paper VI, and references therein). These measurements and relevant comparisons are collected in Table 3.

Table 3. M87* Mass Estimates

Method/Origin	θM (μas)	M (109 M⊙)
§5.1 Direct n = 1 photon ring	4.15 ± 0.74	7.06 ± 1.26
§5.2 Corrected n = ∞ photon ring	4.22 ± 0.17	7.18 ± 0.29
§5.3 Joint M/a reconstruction	${4.20}_{-0.06}^{+0.12}$	${7.13}_{-0.11}^{+0.20}$
Paper VI	3.8 ± 0.4	6.5 ± 0.7
Gebhardt et al. (2011)	${3.62}_{-0.34}^{+0.60}$	${6.14}_{-0.62}^{+1.07}$
Walsh et al. (2013)	${2.05}_{-0.16}^{+0.48}$	${3.45}_{-0.26}^{+0.85}$

Note. Errors indicate 1σ statistical and systematic errors, added in quadrature. See the relevant sections for more detailed error budgets.

Download table as:  ASCIITypeset image

5.1. Direct Mass Estimates

As described in Broderick et al. (2022), the sensitivity of the size of the n = 1 photon ring to the location of an equatorial emission region is bounded. This is in contrast to the n = 0 photon ring, which can grow arbitrarily large with more distant equatorial emission. As a result, with the detection of the n = 1 photon ring, it is now possible to place a limit on the mass that is weakly dependent on the physics of the emission region, subject to the assumption that this emission is confined to near-equatorial regions, e.g., as anticipated by Magnetically Arrested Disk (MAD) models. ⁵⁰

The size of the n = 1 photon ring from equatorial emission is given by

$$\begin{eqnarray}&&{{\theta }}_{n=1}={{\vartheta }}_{n=1}(a,{r}_{{\rm{e}}{\rm{m}}},i)\,{{\theta }}_{M},\end{eqnarray} \tag{ 1 }$$

where ϑ_n=1(a, r_em, i) is a dimensionless function that ranges from 4.30 and 6.17 for polar observers, depending on the black hole spin, a, and radius of the emission peak, r_em. Therefore, from the measurement of θ_n=1 it is possible to generate an astrophysics-independent limit on the mass of M87*:

$$\begin{eqnarray}\begin{array}{rcl}{{\theta }}_{M} & = & {{\vartheta }}_{n=1}^{-1}\,{\theta }\\ & = & 4.15\pm 0.02\pm {\left.0.74\right|}_{{{\vartheta }}_{n=1}}\,{\mu }{\rm{a}}{\rm{s}},\end{array}\end{eqnarray} \tag{ 2 }$$

where we have separately indicated statistical uncertainty and the systematic uncertainty in the relationship between the n = 1 ring and the mass. The corresponding mass estimate is

$$\begin{eqnarray}&&{M}_{9}\,=\,7.06\pm 0.03\pm {\left.1.26\right|}_{{{\vartheta }}_{n=1}}\pm {\left.0.34\right|}_{D},\end{eqnarray} \tag{ 3 }$$

where the systematic uncertainty associated with that on the distance is also separately stated.

5.2. Corrected Mass Estimates

The n = 1 photon ring is biased relative to the asymptotic (i.e., n = ∞ ) photon ring associated with the boundary of the black hole shadow. The degree to which it is biased depends on the spatial distribution of the emission region (Broderick et al. 2020b, 2022). In Appendix C we describe two attempts to estimate the degree of this bias. The first, described in detail in Appendix C.1, is based solely on geometric arguments, assumes a polar observer, and is subject to the size constraints reported in Paper VI. These are shown in Figure 7. These indicate that this bias is robustly limited to less than 1.3 μas, and typically significantly smaller.

Zoom In Zoom Out Reset image size

A potentially more relevant estimate based on GRMHD simulations is described in Figure 8. We make use of the GRMHD simulations reported in Paper V.

Zoom In Zoom Out Reset image size

These incorporate two important additional effects: the small but nonzero inclination (i ranges from 12° to 22°) and emission from above and below the equatorial plane. Both of these tend to shrink the size of the n = 1 photon ring, as seen in Figure 8.

Models that exhibit extended emission, e.g., the R_high = 1 simulations, can have biases that exceed 1 μas. However, these are excluded by the M87* size constraints Paper V. When only models that are consistent with the source size and ancillary limits described in Paper V are considered, the size of the bias is below 0.8 μas; the typical shift is Δθ_n=∞ = 0.56 ± 0.32 μas.

Like that of the n = 1 photon ring, the angular size of the asymptotic photon ring is given by

$$\begin{eqnarray}&&{\theta }_{n=\infty }={\vartheta }_{n=\infty }(a,i){\theta }_{M},\end{eqnarray} \tag{ 4 }$$

where ϑ(a, i) is another dimensionless function, ranging from 4.90 to 5.20 depending on spin for the inclinations relevant for M87*. This range is considerably smaller than that for ϑ_n=1, corresponding to a significantly reduced dependence on the details of the emission region, which is otherwise encoded in Δθ. Thus, the uncertain spin introduces only an additional 3% systematic uncertainty in the relationship between the asymptotic photon ring radius and the mass.

Therefore, the mass of M87* in angular units is given by

$$\begin{eqnarray}\begin{array}{rcl}{\theta }_{M} & = & {\vartheta }_{n=\infty }{(a,i)}^{-1}\left({\theta }_{n=1}-{\rm{\Delta }}{\theta }_{n=\infty }\right)\\ & = & 4.22\pm 0.02\pm {\left.0.06\right|}_{{\rm{\Delta }}\theta }\pm {\left.0.16\right|}_{\vartheta }\,\mu \mathrm{as},\end{array}\end{eqnarray} \tag{ 5 }$$

where again we have separately listed the random measurement error and the systematic errors associated with the emission region (Δθ) and spin (ϑ). This corresponds to a mass estimate of

$$\begin{eqnarray}&&{M}_{9}\,=\,7.18\pm 0.04\pm {\left.0.11\right|}_{{\rm{\Delta }}\theta }\pm {\left.0.27\right|}_{\vartheta }\pm {\left.0.34\right|}_{D},\end{eqnarray} \tag{ 6 }$$

5.3. Joint Mass/Spin Estimate

Finally, following Broderick et al. (2022), we attempt to jointly reconstruct the size of the n = 0 and n = 1 photon rings, as listed in Table 2. This leverages the additional information presented by the diffuse emission to constrain the location of the emission region. However, because there is limited evolution in the maps during the 2017 EHT observing campaign and the simulated analyses in Broderick et al. (2022), the constraint on spin will be weak, rendering this primarily a demonstration in principle.

Details for the model, likelihood, and method of sampling are collected in Appendix D. The resulting joint posterior on θ_M and a is shown in Figure 9. Marginalizing over spin, we find for the mass

$$\begin{eqnarray}&&{\theta }_{M}={4.20}_{-0.06-0.10}^{+0.12+0.23}\,\mu \mathrm{as},\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{M}_{9}\,=\,{7.13}_{-0.11-0.17}^{+0.20+0.39}\pm {\left.0.34\right|}_{D},\end{eqnarray} \tag{ 8 }$$

where the indicated errors correspond to 1σ and 2σ. Note that because the spin is simultaneously reconstructed, these include what has previously been identified as the systematic uncertainties associated with the n = 1 photon ring size bias and spin.

Zoom In Zoom Out Reset image size

After marginalizing over mass, the spin estimate is

$$\begin{eqnarray}&&a\,=\,{1.0}_{-0.5-1}^{+0+0},\end{eqnarray} \tag{ 9 }$$

where again the indicated errors correspond to 1σ and 2σ. As anticipated, the spin is effectively unconstrained, although there appears to be a very weak preference for high spin.

To demonstrate this we repeated the analysis with synthetic size measurements constructed from the equatorial emission model with two sets of (θ_M , a). The first was set to the average values from the joint posteriors, (θ_M , a) = (4.32 μas, 0.63). The second was to set to (θ_M , a) = (4.2 μas, 0.0), exploring the posterior associated with a truth value of zero spin. These are shown in Figure 9 by the dashed red and green lines, respectively, in the one-dimensional, marginalized posteriors for θ_M and a. Both exhibit the same posterior excess near a = 1, implying that it is not significant.

5.4. Synthesis and Discussion

The detection of a bright ringlike feature, and its separation from the diffuse component, has a number of immediate implications for the origin of the emission and the properties of the central black hole. As shown explicitly in Broderick et al. (2020b), the morphology of the diffuse emission is accurately recovered despite the absorption of excess flux into the thin ring component. Thus, here we discuss the structure the diffuse component in the broader context of the environment and properties of M87*.

6.1. Evolution of M87 from 2017 April 5–11

6.2. Extended Southwestern Emission

In all of the diffuse emission maps shown in Figure 5, an extension to the southwest is visible. The ability to produce a statistically meaningful image posterior enables the assignment of a significance to these features, which ranges from 4σ to 12σ on the individual days. On three of the four days there is a matching northwestern extension, though this is less significant (1σ–2σ).

While such features appear similar to the dirty beams seen in some M87* images produced by other algorithms (see Figures 7 and 8 of Paper IV), none have been seen at statistically significant levels (as characterized by the image posteriors) in the various simulated data tests performed with the Bayesian scheme employed here (see, e.g., Figure 4 of Broderick et al. 2020b). Similar claims of a statistically significant detection have been made by other groups (Arras et al. 2022; Carilli & Thyagarajan 2022).

It is suggestive that the orientation of these diffuse extensions align with the limb-brightened jet seen at 3 mm, shown in Figure 10. We align the center of light of the mean 1.3 mm maps, averaged over the four observation days in 2017, and the centroid of the core component of the 3 mm maps. The ridgeline fits from Kim et al. (2018) are also shown in Figure 10, assuming a jet position angle of 69° east of north and a width, W ∝ z^−0.498, where z is the projected core distance. Additional core shifts along the jet of 25 μas and transverse to the jet of −10 μas are applied. Extrapolating the core shift power law determined from longer wavelengths by Hada et al. (2011) places the anticipated 1.3 mm core on top of the brightest diffuse component. ⁵¹ With these shifts, the ridgelines connect from the southwestern and northwestern extensions. We caution that the apparent structure of the jet base on horizon scales may depart substantially from the power-law behavior at large scales due to the combined effects of inclination, gravitation lensing at small radii, inhomogeneous evolution in the optical depth across the image, and relativistic motion (see, e.g., Broderick & Loeb 2009; Mościbrodzka et al. 2016; Chael et al. 2019; Davelaar et al. 2019).

Zoom In Zoom Out Reset image size

The dominance of the southwestern extension, in comparison to the more marginal northwestern extension, is consistent with a rapidly rotating jet structure near the black hole, aligned with the black hole spin. Brightness temperature profiles are shown in Figure 11 along the corresponding paths shown in Figure 10. The brightness ratio between these two components reaches values as low as ${ \mathcal R }\approx 0.2$ at similar projected distances from the black hole. The uncertainty of these profiles is a combination of the intrinsic uncertainty in the reconstructions and of the ridgelines themselves; we show the uncertainty associated with a Gaussian error in the position of the parabola apex (the rightmost open green point in Figure 10) of 2 μas for the EHT profiles and 10 μas for the 3 mm profiles.

Zoom In Zoom Out Reset image size

In Figure 11 the difference in the brightness temperature profiles is naturally explained by black hole spin-driven jets like those first described in Blandford & Znajek (1977). Relativistic rotational motions in the emitting plasma near the black hole, with velocities that could reach as high as c/2, are a consequences of the twisted magnetic field lines that penetrate the horizon and extract the rotational energy of the black hole. These motions are seen explicitly in GRMHD simulations (see, e.g., Wong et al. 2021). At larger distances, this rotation is suppressed due to angular momentum conservation, at which point the jet plasma motions become predominantly poloidal. As a result, even only very modest increases in the jet height result in drastic reductions in the beaming-induced brightness asymmetry, naturally explaining its absence at 3 mm and longer wavelengths.

In both Figure 10 and Paper V the implied spin of the black hole is oriented into the sky, i.e., the black hole rotates clockwise, dragging the emitting plasma with it. As such, this provides a striking, direct confirmation of black hole spin as the driver of the jet in M87*. In practice, the angular scale over which the limbs become symmetric depends on the detailed jet structure, and thus black hole spin (Takahashi et al. 2018); however, we leave further astrophysical interpretation to future work.

We have applied the hybrid imaging algorithm outlined in Broderick et al. (2020b) to the 2017 EHT observations of M87*. This method is a Bayesian imaging and modeling scheme that reconstructs the brightness map from visibility data, accounting for station gains and atmospheric phase delays, and produces statistically meaningful posteriors for both components. We considered both imaging and imaging with a narrow ring, finding that information criteria indicate a preference for the latter on three of four observation days.

We demonstrate that the EHT observations of the horizon scale emission of M87* support the presence of a narrow ringlike feature. Its radius is consistent across the seven days of the 2017 EHT observing campaign (2017 April 5–11). The size and structure of this ring is consistent with the prominent lensed structure anticipated in horizon-resolving images of M87*. We associate the ring emission with the strong gravitational lensing that produces the n = 1 photon ring; this decomposition thus represents the first direct detection of the "back of the emission region." It also provides an important confirmation of the key role played by strong lensing in the formation of the images of M87* presented in Paper I.

Evolving, extended diffuse emission is clearly present in addition to the bright, narrow ring. The extended image is consistent among neighboring days, as anticipated by the long dynamical timescales in M87*. However, it differs from the beginning of the week to the end of the week. This may be due either to the appearance of a new southern component or the shearing of a western component southeastward. The latter is consistent with the direction of motion expected from the black hole spin orientation presented in Paper V.

The diffuse emission is dominated by compact components that surround the bright ring, similar to the morphology in many GRMHD simulation images. This is more compact than many such images, suggesting that additional constraints on the library in Paper V can be made based on the compactness of this portion of the diffuse emission.

Extended components within the diffuse emission are also detected at statistically significant levels. These include a southwestern extension, detected at between 4σ and 8σ across all days, and a northwestern extension that is only marginally detected on two days (April 5 and 11). Both of these components are reminiscent of the larger-scale limb-brightened features seen at 3 mm; their orientation matches those extrapolated from longer wavelengths. These may be tentatively identified with the rapidly rotating jet footprint. The difference in the luminosity of the southern and northern components would then be naturally explained by relativistic jet rotation at the jet base. Both of these would support the conclusion that the jet in M87* is driven by the black hole rotation, as described in Blandford & Znajek (1977).

The size of the bright ring and its relation to the diffuse emission presents a number of ways to estimate the black hole mass-to-distance ratio with varying degrees of astrophysical inputs. All of these estimates are consistent, with a θ_M = 4.15 ± 0.74 μas, independent of the extent of the emission region. Our most precise estimate arises from simultaneously reconstructing the ring and diffuse emission scales, and obtains ${\theta }_{M}={4.20}_{-0.06}^{+0.12}\,\mu \mathrm{as}$, which improves on the fractional uncertainty presented in Paper VI by a factor of more than four. The resulting mass estimate after folding in a distance of 16.8 Mpc is ${7.13}_{-0.11}^{+0.20}\times {10}^{9}\,{M}_{\odot }$. The uncertainty in the mass estimate is now dominated by the systematic uncertainty in the distance of 0.35 × 10⁹ M_⊙.

These mass-to-distance ratio estimates are consistent with those from the variety of methods presented in Paper VI and from stellar dynamics (Gebhardt et al. 2011). Because the latter is estimated from the dynamics of what are effectively test particles at distances four orders of magnitude larger than the photon orbit, the comparison of these mass measurements provide a direct test of general relativity, as described in Paper VI. In practice, this test remains limited by the uncertainty in the stellar dynamics measurement.

Nevertheless, the mass estimates presented here lie at the high end of the ranges presented in Paper VI. This suggests that the set of GRMHD simulations used to calibrate the mass estimates in Paper VI were more extended than the observed emission, biasing the calibration factor α toward high values and thus the mass toward low values. This is only partly ameliorated by selecting only MAD models in the calibration process. In comparison, no such calibration is required for two of the three mass estimates made here; in the one instance where calibration from simulations is performed, the measured systematic modification is small due to the much more robust size of the bright rings in simulated images.

Additional epochs of horizon-resolving observations will prove particularly useful in confirming the existence and nature of the bright ring. While small variations in its location are anticipated, associated with (potentially large) variations in the location of the emission region, if the ring structure detected here is, in fact, identified with the n = 1 photon ring, it should persist.

The evolution of the diffuse emission map will also be diagnostic of the location and origin of the emission in M87*. Similar quality observations that extend over observation epochs as short as two weeks will permit the conclusive differentiation between orbiting and outflowing features (Jeter et al. 2020). The ability to resolve dim, variable structures thus motivates longer-duration observing campaigns.

The constraints on black hole spin obtained here are inconclusive. However, the ability to constrain (θ_M , a) to a band in the mass-spin parameter space indicates that even a single, fortuitous future EHT observation of M87* may provide a measurement of black hole spin from gravitational lensing alone (Broderick et al. 2022). The strength of potential spin constraints depends on the degree to which the emission region location differs during future observations. Multiple additional measurements provides a lensing-only test of general relativity (Broderick et al. 2022). These provide a strong motivation for including M87* in future EHT campaigns and those of subsequent instruments.

It also suggests that future space-based millimeter-very long baseline interferometry experiments may be able to detect the next-order lensed image, i.e., that associated with the n = 2 ring, via a similar method to that presented here (Broderick et al. 2022). The astrophysics-independent mass measurements become substantially better constrained in this instance, and immediate tests of general relativity become possible.

Finally, this provides a direct demonstration of the ability to leverage high S/N data to estimate image features with precisions that significantly exceed those implied by the ostensible observing beam. This effective super-resolution may suggest that a similar effort applied to more distant sources may yield practical mass estimates even in the absence of resolved ring structures in images.

This work was made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET: www.sharcnet.ca) and Compute/Calcul Canada (www.computecanada.ca). Computations were made on the supercomputer Mammouth Parallèle 2 from the University of Sherbrooke, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the ministère de l'Économie, de la science et de l'innovation du Québec (MESI), and the Fonds de recherche du Québec—Nature et technologies (FRQ-NT). This work was supported in part by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade. A.E.B. thanks the Delaney Family for their generous financial support via the Delaney Family John A. Wheeler Chair at the Perimeter Institute. A.E.B. and P.T. receive additional financial support from the Natural Sciences and Engineering Research Council of Canada through a Discovery grant. R.G receives additional support from the ERC synergy grant "BlackHoleCam: Imaging the Event Horizon of Black Holes" (grant No. 610058). D.W.P. is supported by the NSF through grant Nos. AST-1952099, AST-1935980, AST-1828513, and AST-1440254; by the Gordon and Betty Moore Foundation through grant No. GBMF-5278; and in part by the black hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation to Harvard University. Furthermore, we thank Ivar Coulson for useful contributions during this project. H.-Y.P. acknowledges the support of the Ministry of Education (MoE) Yushan Young Scholar Program, the Ministry of Science and Technology (MOST) under grant No. 110-2112-M-003-007-MY2, and National Taiwan Normal University. IMV acknowledges support from Research Project PID2019-108995GB-C22 of Ministerio de Ciencia e Innovación (Spain) and from the GenT Project CIDEGENT/2018/021 of Generalitat Valenciana (Spain).

We employ three model components in the comparisons made here. These are the "nonparameteric" raster splined raster model, the large-scale Gaussian, and a slashed crescent. Here we collect details of the implementation and assumed priors on each component.

The splined raster model is described in Broderick et al. (2020b), to which we direct the reader for more detailed information. It is characterized by a rectilinear set of control points, I_M,N , which may vary independently. An example with the control points explicitly marked is shown in Figure 12. The intensity map is interpolated between these control points using an approximate cubic spline. We modify this by permitting a variable field of view, (FOV_x ,FOV_y ), orientation, ϕ, and shift in the center of the raster, (x,y). Parameters and priors are listed in Table 4. Note that the ring radius in angular units is given by θ = R(1 − ψ/2).

Zoom In Zoom Out Reset image size

Table 4. Hybrid Image-ring Fit Priors

Comp.	Param.	Units	Prior a
Raster	IM,N	Jy μas−2	${ \mathcal L }({10}^{-14},3\times {10}^{3})$
IM×N	FOVx	μas	${ \mathcal U }(0,200)$
	FOVy	μas	${ \mathcal U }(0,200)$
	ϕ	rad	${ \mathcal U }(-0.25\pi ,0.25\pi )$
	x	μas	${ \mathcal U }(-40,40)$
	y	μas	${ \mathcal U }(-40,40)$
Gaussian	I0	Jy	${ \mathcal U }(0,10)$
A	σ	mas	${ \mathcal U }(0.1,{10}^{4})$
	A	⋯	${ \mathcal U }(0,1)$
	ϕ	rad	${ \mathcal U }(0,\pi )$
	x	mas	${ \mathcal U }(-2,2)$
	y	mas	${ \mathcal U }(-2,2)$
Ring	I0	Jy	${ \mathcal U }(0,2)$
X	R	μas	${ \mathcal U }(0,{10}^{2})$
	ψ	⋯	${ \mathcal U }(0,0.05)$
	f	⋯	${ \mathcal U }(0,1)$
	ϕ	rad	${ \mathcal U }(-\pi ,\pi )$
	x	μas	${ \mathcal U }(-40,40)$
	y	μas	${ \mathcal U }(-40,40)$

Note.

^a Linear priors from a to b are represented by ${ \mathcal U }(a,b)$, logarithmic priors from a to b are represented by ${ \mathcal L }(a,b)$.

Download table as:  ASCIITypeset image

The large-scale Gaussian is the asymmetric Gaussian model described in Broderick et al. (2020a). This is characterized by a total flux, I₀, a symmetrized standard deviation, σ, asymmetry parameter, A, position angle, ϕ, and location (x,y). Parameters and priors are also listed in Table 4.

Finally, we apply the slashed ring using the xs-ringauss model described in Broderick et al. (2020a). We force the asymmetry parameter and additional Gaussian to vanish via the choice of appropriate priors. This leaves the total flux, I₀, outer radius, R, fractional width, ψ, linear brightness gradient, f, brightness gradient position angle, PA, and location, (x,y). Of these, we restrict in particular the fractional width to ensure that the ring thickness remains well below the scale of the spacings within the splined raster.

These are supplemented with a model for the scan-specific complex station gains. The logarithmic amplitudes of the gains are assumed to have a Gaussian prior with width 20% for all stations except the Large Millimeter Telescope (LMT), for which we assume a prior with width 100%. The phases priors are flat and uninformative. The gain phases are trivially degenerate with an overall shift in the image location, and thus we remove two parameters from the total parameter count in Table 1.

The contribution of the individual model components are illustrated in Figure 13. The observed decomposition matches that anticipated by the GRMHD models explored in Broderick et al. (2020b). The large-scale Gaussian only impacts the intrasite baselines, the Atacama Large Millimeter/submillimeter Array (ALMA)-Atacama Pathfinder Experiment (APEX) telescope (APEX) and the James Clerk Maxwell Telescope (JCMT)-Submillimeter Array (SMA). The raster image and narrow ring both contribute at longer baselines, with the latter dominating at the longest baselines. The null in the visibility amplitudes near 3.5 Gλ is reconstructed by a combination of these components, both of which exhibit nulls shifted from the observed baseline length. In practice, Figure 13 shows the range of the visibility amplitudes at each baseline length, and thus the distinguishing power of the data exceeds that illustrated in Figure 13.

Zoom In Zoom Out Reset image size

The result of the image-reconstruction process is a posterior distribution of brightness maps, corresponding to the posterior distribution of the underlying model parameters. As shown in Figure 14 for representative days, the posterior range of visibility amplitudes associated with the image posterior is small at (u,v) positions near the EHT measurements and largest in the (u,v) holes, prominently in the north/south near 4 Gλ. Thus, as expected, the posterior encompasses the uncertainty associated with the unknown visibilities away from the locations of the EHT data, subject to the constraints of the underlying image model.

Zoom In Zoom Out Reset image size

In addition to the diameter of the bright ringlike feature, a number of additional model parameters are constrained. In Figure 15 we collect some of the parameter constraints that may be of interest.

Zoom In Zoom Out Reset image size

The ring flux posteriors present the same information reported in Table 2 and associated discussion, to which we direct the reader for more information.

The ring position angle estimates are consistent among days, in contrast to the reconstructions in Paper IV and Paper VI, which exhibited an evolution from April 5 and 6 to 10 and 11. The latter change is indicated in the second column of Figure 15 by the orange shaded region. In our analysis, this variation is entirely confined to evolution within the reconstructed diffuse emission.

The ring fractional width, ψ, fills the posterior, indicating that we have no traction on the width itself.

The orientation of the large-scale Gaussian is weakly constrained by the two sets of intrasite baselines, ALMA-APEX and JCMT-SMA. The posterior of the major-axis position angle is shown in the rightmost column of Figure 15. These are only very weakly nonuniform, though do peak in the octant consistent with the orientation of the large-scale radio jet. We do not attempt to interpret the large-scale Gaussian component further.

Generally, the n = 1 photon ring is biased relative to the asymptotic n = ∞ photon ring, associated with the boundary of the "shadow." To estimate the magnitude of this bias and its sensitivity to astrophysical uncertainties we explore both general geometric arguments and GRMHD modeling of M87* as described in Paper V.

C.1. Geometric Modeling

As described in Broderick et al. (2020b), the size of the emission region places a constraint on the magnitude of the difference between the n = 1 and n = ∞ photon rings. We compute this here by explicitly tracing rays back toward the black hole launched from a polar observer and categorize trajectories by the number of equatorial plane crossings. Relevant examples are shown in Figure 16.

Zoom In Zoom Out Reset image size

The diameter of the observed bright ring of 42 ± 3 reported in Paper VI places a limit on the location of an equatorial emission region. Assuming that this emission is dominated by the direct (i.e., n = 0) photon trajectories, this mapping may be made explicit. The population of orbits that terminate within the observed ring diameter is shown by the n = 0 range in each panel of Figure 16, and the region where they cross the equatorial plane is indicated by the thick bright line segments.

The n = 1 photon ring is then formed by those photon trajectories that intersect this portion of the equatorial plane after executing an additional half-orbit about the black hole. These may be traced forward to obtain the radial position on the image plane to which these trajectories contribute, θ_n=1. These are also shown by the n = 1 range in each panel of Figure 16. A similar range of trajectories for n = 2 is shown for completeness.

The n = ∞ envelope is also computed, providing an estimate of the size of the photon ring, θ_n=∞. These may then be compared, yielding an estimate of the bias, θ_n=1 − θ_n=∞. These limits are summarized for a handful of black hole spins in Figure 7.

In practice, due to the finite resolution of the instrument, the diameter is weakly degenerate with the subbeam substructure of the emission region. This manifests, for example, as a dependence on the width of the ring (see Appendix G of Paper IV). Therefore, these estimates of the range of the n = 1 photon ring bias are overestimates of the actual bias, and thus systematic uncertainty.

C.2. GRMHD Simulation

Because the bias between the n = 1 and n = ∞ photon ring radii depends on the emission distribution, we explored the relationship between the bright ring and the asymptotic photon ring in GRMHD simulation images. We employed 198 GRMHD simulation image libraries, corresponding to those reported in Paper V supplemented by a set of a = 0.75 simulations. These span MAD and Standard and Normal Evolution (SANE) morphologies, a ∈ [ − 0.9375, − 0.5, 0, 0.5, 0.75, 0.9375], i ∈ [12°, 17°, 22°], and R_high ∈ [1, 10, 20, 40, 80, 160], where negative spins correspond to counterrotating accretion flows, R_high describes the relationship between the electron and ion temperatures, and the a = 0.75 simulations were performed only for the MAD morphology.

For each simulation, we time-averaged the simulated images to improve smoothness and reduce the impact of turbulent fluctuations. Generally, the image center will be shifted from the image origin, i.e., (x, y) = (0, 0), as a result of the small nonzero inclinations and spin. Therefore, to determine the center we minimize the variance in the radial location of the image brightness maximum between 4GM/c² and 7GM/c². Finally, the brightness was computed along 36 radial chords originating from the image center, normalized by its maximum value between 2GM/c² and 8GM/c², and then averaged, generating I_avg. The radius of the bright ring was then taken to be the position of the peak of the I_avg. This procedure is shown for three examples in Figure 17.

Zoom In Zoom Out Reset image size

Each time-averaged simulation image-ring fit was visually inspected; there are no significant failures. As anticipated by the geometric calculation in Appendix C.1, the outward bias in the location of the bright ring is dependent on the extent of the emission region, growing with the second moment of the image intensity map. The distribution of the shifts are shown in Figure 8 when the results for all simulations are included and when only those that were deemed acceptable in Paper V.

Following Broderick et al. (2022), we model the position of the peak of the n = 0 and n = 1 photon rings. The angular radii of the former are determined relative to the center of the fitted ring. These are listed in Table 2.

From these measurements, we construct a Gaussian log-likelihood from the measurements listed in Table 2:

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }({\theta }_{M},a,{r}_{\mathrm{em},5},{r}_{\mathrm{em},6},{r}_{\mathrm{em},10},{r}_{\mathrm{em},11})\\ =\sum _{d=\mathrm{day}}\sum _{m\,=\,0,1}\displaystyle \frac{{\left[{\theta }_{n=m,d}-{\vartheta }_{n=m}(a,{r}_{\mathrm{em},{\rm{d}}},i){\theta }_{M}\right]}^{2}}{2{\sigma }_{n\,=\,m,d}^{2}}.\end{array}\end{eqnarray} \tag{ D1 }$$

The six parameters are the mass, spin, and the four radii at which the emission peaks on each day. In all cases we assume a polar observer (i = 0); due to the weak dependence of the n = 0 and n = 1 photon rings on inclination near the pole and the large uncertainties on the measured ring radii, this makes a small difference. The uncertainties σ_n=0,d and σ_n=1,d are computed directly from the standard deviations of the θ_n=0 and θ_n=1 estimates.

Priors on all parameters are adopted either on ranges large in comparison to their reconstructed values or (for spin) that encompass all physical values. Uniform priors are adopted on the mass and spin within the intervals [0, 10 μas] and [0, 1), respectively. Logarithmic priors are adopted for each r_em,d on the interval [r_h , 25r_h ], where ${r}_{h}=1+\sqrt{1-{a}^{2}}$ is the horizon radius in units of GM/c².

We sample the corresponding posterior using the ensemble MCMC method provided by the emcee python package (Foreman-Mackey et al. 2013). We employ 64 independent walkers and run for 10⁵ steps, discarding the first half of each chain. Explorations with simulated data sets and inspection of the resulting MCMC chains indicate that by this time they are well converged.