© The Authors 2025

1. Introduction

The first images by the Event Horizon Telescope (EHT) of the supermassive black hole (SMBH) in the heart of the M87 galaxy (M87^∗) revealed a ring structure with a diameter of approximately five times the projected Schwarzschild radius of a 6.5 × 10⁹ M_⊙ black hole (M87^* 2017 I; M87^* 2017 II; M87^* 2017 III; M87^* 2017 IV; M87^* 2017 V; M87^* 2017 VI). The EHT Collaboration recently published two new images of the M87^∗ ring based on EHT data that were collected in April 2018 (M87^* 2018 I, hereafter Paper I), almost exactly one year after the first observations. While the persistent ring structure revealed in the new images strongly supports the idea that the central depression of the M87^∗ image is indeed the shadow of an event horizon of a supermassive object (e.g., Hilbert 1917; Bardeen 2018; Luminet 1979; Jaroszynski & Kurpiewski 1997; Falcke et al. 2000, M87^* 2017 VI, Wielgus et al. 2020), the new data show a different brightness distribution in the ring. The new observations therefore place additional constraints on physical models of the emitting plasma close to the event horizon.

The conventional model for M87^∗ is a black hole surrounded by a magnetized, geometrically thick, optically thin, radiatively inefficient, advection dominated, rotating accretion disk (e.g., Ichimaru 1977; Rees et al. 1982; Narayan & Yi 1994, 1995; Reynolds et al. 1996) that launches a relativistic jet. There is no consensus model for the jet-launching mechanism, but the two main scenarios are that the jet is either dominated by Poynting flux and extracts rotational energy from the supermassive black hole (Blandford & Znajek 1977) and/or that the jet is a magnetohydrodynamically collimated wind from the accretion disk that is launched close to the event horizon (Blandford & Payne 1982; Lynden-Bell 2006). The jet from M87^∗ is clearly visible at longer radio wavelengths, in optical light, X-rays, and even γ-rays as an elongated feature that is resolved from submilliarcsecond to arcsecond scales (EHT MWL Science Working Group 2021; Lu et al. 2023).

The synchrotron emission from the disk-jet system close to the black hole is gravitationally lensed and Doppler boosted so that it appears to an external observer as an asymmetric ring structure. Because of astrophysical uncertainties and strong gravity, the physical interpretation of the ring emission is not straightforward and typically involves forward modeling using either semianalytic models (Broderick & Loeb 2009) or numerical simulations (e.g., Dexter et al. 2012; Mościbrodzka et al. 2016; Fromm et al. 2022).

The theoretical interpretations of the first 2017 EHT images of M87^∗ have been presented in three preceding works (M87^* 2017 V, M87^* 2017 VIII, M87^* 2017 IX). In M87^* 2017 V, we created a library of 60 000 mock black hole images based on general relativistic magnetohydrodynamic (GRMHD) simulations of black hole accretion. The library surveyed different black hole spins and electron temperature parameterizations, and two distinct accretion flow modes: the standard and normal evolution (SANE) models, and the magnetically arrested disk (MAD) models. The initial setup for all the simulations featured a magnetized torus of plasma orbiting in the equatorial plane of the black hole. The different types of magnetic field geometry in the initial torus evolve to either the SANE or the MAD state. In M87^* 2017 V, we constrained the models using M87 EHT 2017 data, observational jet power, and the M87 core X-ray luminosity. The models that passed these observational constraints included both SANEs and MADs, both with positive and negative spins. The sign of the spin is positive here when the accretion flow angular momentum and black hole spin are parallel, and it is negative when they are antiparallel (a counter-rotating, or retrograde, accretion flow). M87^* 2017 V reported two main findings. First, the southern brightness asymmetry in the image strongly supports the interpretation that the spin axis of the M87 black hole points away from Earth. Second, it was predicted that if the black hole spin axis is normal to the disk (parallel to the large-scale jet axis), then future observations would most often show the brightest part of the ring appearing counterclockwise from the position seen in the 2017 observations.

Polarimetric images are more sensitive than total intensity alone to the plasma properties around the supermassive black hole. In a subsequent publication (M87^* 2017 VIII), we compared the library of GRMHD simulations (72 000 images in total, extending the previous library with additional models for electron thermodynamics) to the 2017 linear polarimetric image of the M87^∗ ring (M87^* 2017 VII). We found that models that fell within the allowed ranges of the measured polarimetric characteristics are typically MAD. The observed, azimuthally dominated electric vector position angle (EVPA) pattern is usually inconsistent with SANE simulations. Although several SANE snapshots are consistent with EHT polarimetric observations, these models fail to produce sufficient jet power. The M87^* 2017 VIII constraints narrowed the range of allowed mass accretion rates onto the M87 black hole horizon, Ṁ, to (3 − 20)×10⁻⁴ M_⊙/yr. The near horizon jet power, P_jet, was narrowed down to 10⁴² − 10⁴³ ergs s⁻¹. The measurements of circular polarization on horizon scales of M87^* presented in M87^* 2017 IX did not change these conclusions.

The current paper is the second in a series of papers dedicated to the analysis of the EHT 2018 observations. In this paper, we focus on the theoretical interpretation of the observed total intensity of M87^∗. To do this, we prepared an updated image library that covers a wider range of possible synchrotron emission natures and parameters of the surrounding plasma, and our images are based on more advanced GRMHD simulations. Considering the constraints imposed by multi-epoch observations, we compare the model images separately with both the 2018 (Paper I) and the initial (2017) observations of M87^∗. The information derived from these comparisons is then combined to constrain the GRMHD models and draw new conclusions about the physical conditions around the supermassive black hole and to constrain the properties of the supermassive black hole itself, assuming that the 2017 and 2018 observations represent a typical state of the source.

The paper is structured as follows. In Sect. 2 we provide a list of the new observational constraints to motivate the more advanced library of GRMHD simulations. In Sect. 3 we describe our library. In Sect. 4 we describe our data-model comparison scheme and present the results of these comparisons. We discuss the new results in Sect. 5 and conclude in Sect. 6.

2. New observational constraints and expectations

Figure 1 (top and bottom left panels) shows representative reconstructed images of M87^∗ based on the 2017 and 2018 datasets. The 2018 image reconstructions of M87^∗ show the bright ring of emission with a diameter¹$d=43.3_{-3.1}^{+1.5}\mu$as; this angular size is consistent with the 2017 EHT measurements ($d={42}_{-3.0}^{+3.0}\mu$as, M87^* 2017 VI) and also with measurements from the less precise 2009-2013 proto-EHT data (Wielgus et al. 2020), but with a significantly different position angle of the brightest part of the ring (rotation by ∼30° from ∼180° as measured in 2017 to ∼210° in the 2018 data).

Fig. 1. Representative reconstructed images and model images of M87∗. Left panels: EHT images of M87∗ from the 2018 (M87* 2018 I; Paper I of this series) and 2017 (M87* 2017 IV) campaigns. Middle panels: Example GRMHD model images drawn from the same model at two different times. Right panels: Theoretical images convolved with a 20 μas FWHM Gaussian beam.

The shift in the observed brightness position angle is expected from both observational precedent and theoretical predictions. Observationally, an analysis of historical VLBI data at 230 GHz using a prototype EHT array indicates even more variability on the horizon scale between 2009 and 2017 (Wielgus et al. 2020). Theoretically, images from GRMHD simulations feature significant variability in their image features. The middle and right panel in Fig. 1 show two snapshots from an example movie based on the same GRMHD simulation and their corresponding image when convolved with a beam similar to EHT resolutions. These frames illustrate that the observed shift in the brightest location on the ring can naturally arise from the turbulent nature of the accretion environment around a black hole. The 2018 images of M87^∗ are consistent with the forecasts we made in M87^* 2017 V. Compared to the 2017 images, the 2018 observed ring brightness distribution is closer to the mean value expected from coaxial models in which the black hole spin axis is aligned with the orientation of the large-scale jet seen at lower radio frequencies (Walker et al. 2018).

Additional constraints come from the coordinated simultaneous broadband multiwavelength observations of M87 carried out during the EHT campaigns in April 2017 and 2018 by the EHT Multiwavelength (MWL) Science Working Group (EHT MWL Science Working Group 2024). These observations covered more than 17 decades in frequency from radio (cm wavelength) to very high-energy (VHE) γ-rays (> 100 GeV). In April 2017, the source was found to be in a historically low state at all frequencies (EHT MWL Science Working Group 2021), but in April 2018, M87 underwent a γ-ray flaring episode (EHT MWL Science Working Group 2024), the first detected episode since 2010, with the peak on April 21 (MJD 58229.07). The VHE γ-ray flux doubled within a period of 36 hours, which enabled us to constrain the size of the VHE γ-ray emitting region to $R_{\mathrm{HE}}\approx 8(GM/c^2)\delta \left(\frac{\mathrm{\Delta }t}{3\mathrm{days}}\right)$, where δ is the standard Doppler beaming factor (Aharonian et al. 2006). This timescale is consistent with the γ-ray variability arising from the EHT-emitting region, or from a larger region in a more relativistic part of the jet. While the radio and millimeter (mm) core fluxes are compatible with (or potentially lower than) the April 2017 emission, a likely longer-term core flux enhancement was also observed in the X-ray band.

For neither the 2017 nor the 2018 observation can heuristic single-zone models fit the entire SED, suggesting the need for a stratified model. More sophisticated modeling will be needed to localize the γ-ray flare, and nonthermal processes will need to be included. Intriguingly, the hint of a decline in the submm core flux observed in 2018 is a prediction of reconnection-powered VHE flare mechanisms near the event horizon, which may lead to mass ejections out of the EHT-resolved region (Ripperda et al. 2022; Gelles et al. 2022; Hakobyan et al. 2023; Jia et al. 2023). This dip was not a statistically significant detection in 2018, but the possibility argues for more coordinated MWL observations, particularly between EHT and VHE facilities, to determine the location of the flare events.

In this work, we extend the previously constructed GRMHD libraries that were used to interpret M87 EHT observations. The new models are motivated by (1) a more limited set of the old models, which were not in steady state and did not include radiation and tilt, and (2) lower-frequency (86 GHz) observations by Cui et al. (2023), which motivate the exploration of tilted models that may lead to jet precession. The theoretical models are compared to selected EHT observations. Selected best-bet models, according to EHT constraints and jet power constraints, are discussed in the context of the radio and high-energy constraints.

3. Image library of the extended GRMHD simulations

Our numerical model library consists of a number of 3D GRMHD simulations of gas from a magnetized torus accreting onto a black hole and producing a jet. Table 1 summarizes all GRMHD models used in this work, together with details of their parameters. The parameters that are common across different simulation codes include the dimensionless black hole spin (a_*), the adiabatic index in the adiabatic equation of state (Γ_ad), the final time of the simulation (t_f), the outer radius of the computational domain (r_out), and the numerical resolution of the simulation in three dimensions (3D). Compared to the previous image library (M87^* 2017 V), most of the images in our new library were computed based on GRMHD simulations evolved for a longer time (> 25 × 10³GM/c³ ≈ 25 years), with model images derived from an interval at the end of the simulation with a length of 5 × 10³GM/c³ ≈ 5 yr.

The black hole spin is described by the dimensionless spin parameter a_* ≡ Jc/GM² (e.g., Koide et al. 2000; De Villiers et al. 2003; Gammie et al. 2003; Porth et al. 2019), where J, M, G, and c are the black hole spin angular momentum, black hole mass, gravitational constant, and speed of light, respectively. The black hole spin is a free parameter within

$$\begin{array}{c}\hfill -1<a_{\ast }<1.\end{array}$$(1)

We assumed unless stated otherwise that the angular momentum of the accretion flow, the spin of the black hole, and the large-scale jet are coaxial. The inclination is defined as the angle between the accretion flow angular momentum and the line of sight, so that i = 0° implies that the accretion flow angular momentum is pointed at the observer. Positive (negative) spin implies that the angular momentum of the accretion flow and the black hole spin are aligned (anti-aligned), and thus, that the accretion flow is prograde (retrograde). When the angle between the approaching jet and the line of sight is 17° (e.g., Walker et al. 2018) then the model inclination can be either 17° or 163°.

M87^* 2017 V showed that in coaxial models, the ring asymmetry follows the spin. In models in which the spin is pointed along the approaching jet (toward Earth), the mean position of the peak brightness on the ring is 90° counterclockwise from the jet, at a position angle of 20°. In models in which the spin is pointed away from Earth (along the counterjet), the mean position of the peak brightness on the ring is about 90° clockwise from the jet, at a position angle of 200°. Physically, this arises from a combination of effects (Doppler beaming, lensing, and aberration) that are difficult to separate in a relativistic context, but which all affect the fluid-frame frequency at which synchrotron emission is produced. Because GRMHD models are turbulent, the position of the peak brightness fluctuates around these expected values. M87^* 2017 V concluded that because the peak brightness on the ring in the 2017 data was closer to a PA = 200° than 20°, the spin direction was pointed away from Earth.

All GRMHD simulations in this work began with a hydrodynamic torus of gas (Fishbone & Moncrief 1976) with a constant adiabatic index Γ_ad seeded with a loop of poloidal magnetic fields, except for the radiative model. For the radiative GRMHD simulation, the initial torus followed Penna et al. (2013), and a variable adiabatic index (Sądowski et al. 2017) was used. These initial conditions for the simulations were reported in the papers related to the simulations (Prather 2022; Chatterjee & Narayan 2022; Cruz-Osorio et al. 2022; Fromm et al. 2022; Chatterjee et al. 2020; Chael et al. 2019a).

The accreting gas eventually becomes highly turbulent due to disk instabilities such as the magneto-rotational instability (MRI; Balbus & Hawley 1991), and occasionally, the magnetized Rayleigh–Taylor instabilities (e.g., Ripperda et al. 2022). The initial magnetic flux content of the disk determines the evolution of the accretion flow and therefore sets the magnetic flux at the black hole event horizon ϕ_BH. The horizon magnetic flux is an important quantity in our simulations as it determines the expected jet power according to the Blandford & Znajek (1977) effect.

As in M87^* 2017 V, we considered two different accretion states depending on the mean magnetic flux crossing one hemisphere of the event horizon, ϕ_BH: (i) a weakly magnetized accretion flow (“SANE”; ϕ_BH ≈ 1; Narayan et al. 2012; Porth et al. 2019), and (ii) a magnetically arrested disk (“MAD”; ϕ_BH ≈ 15; Bisnovatyi-Kogan & Ruzmaikin 1974; Igumenshchev et al. 2003; Narayan et al. 2003). ϕ_BH is defined² as the dimensionless form of the absolute magnetic flux (Φ_BH) threading a hemisphere at the horizon, that is, ${\varphi }_{\mathrm{BH}}\equiv {\mathrm{\Phi }}_{\mathrm{BH}}{(\dot{M}{r}_g^{2}c)}^{-1/2}$.

SANE accretion flows show MRI-driven disk turbulence and produce relativistic jets with powers P_jet ≲ 0.1 Ṁc². MADs are characterized by dynamically strong magnetic fields near the black hole. MADs occasionally exhibit eruptions of magnetic flux tubes that affect the disk evolution (Chatterjee & Narayan 2022). MADs also show highly efficient jets with powers that can exceed the available accretion power for high black hole spins, P_jet ≳ Ṁc² (Tchekhovskoy et al. 2011). Figure 2 displays the comparison for a typical MAD and SANE simulation. The top panel shows the dimensionless temperature of the fluid (≡ P_gas/(ρc²), where P_gas and ρ are the gas pressure and gas density, respectively), and the bottom panel shows the plasma β ≡ P_gas/P_mag (where P_mag ≡ B²/8π is the magnetic field pressure). Compared to a SANE flow, a MAD state has higher temperatures and more low-β regions.

The GRMHD library includes a set of fiducial models that span the parameters of the accretion state magnetization (SANE and MAD) and black hole spin a_*. Part of the library parameters overlap with the libraries created for the interpretation of the first M87 images (M87^* 2017 V; M87^* 2017 VI; M87^* 2017 VIII). Unlike the previous libraries, however, the fiducial GRMHD simulations in this work were run out to 3 × 10⁴ GM/c³ (cf. 10⁴ GM/c³ in M87^* 2017 V). GRMHD simulations at later times were preferred as the turbulent flow close to the black hole (the emission region) is fully relaxed, which minimizes any influence of the initial conditions on the resulting images (Narayan et al. 2022).

In addition to the fiducial simulations, we included the following exploratory sets: (i) tilted disks, where we relaxed the assumption that the angular momentum of the accreting gas and the angular momentum of the black hole are (anti-) aligned (Liska et al. 2018; Chatterjee et al. 2020), and (ii) radiative GRMHD simulations that self-consistently accounted for electron heating and cooling via radiative losses (Chael et al. 2019a). These exploratory models span a much smaller parameter space than the fiducial models. They mainly served as spot checks for physics that is not accounted for in the fiducial set, and they come from simulations that were on hand at the time of this analysis. This means that no additional exploratory models were run as part of this project.

A total of four GRMHD codes contributed to our model libray: BHAC (Porth et al. 2017), H-AMR (Liska et al. 2022), KHARMA (Prather 2022), and koral (Sadowski et al. 2013; Sądowski et al. 2017). These codes have been shown to provide results that are broadly consistent with each other, with some variation in the disk evolution or turbulence depending on the choices of the numerical implementation (Porth et al. 2019).

The GRMHD simulations were post-processed with three ray-tracing codes BHOSS (Younsi et al. 2012, 2020), grtrans (Dexter 2016), and ipole (Mościbrodzka & Gammie 2018 and see Wong et al. 2022 for pipeline description) that integrate the radiative transfer equations assuming synchrotron emission and self-absorption. The performance of these and other independent radiative transfer schemes was consistent given the same GRMHD simulation snapshot (Gold et al. 2020; Prather et al. 2023). For the sake of the computational efficiency, all images were computed using the so-called fast-light approximation, which assumes that the light has infinite speed. In this approximation, each GRMHD time-slice produces a single image. The alternative slow-light approach, which traces light rays through the evolving fluid simulation across a sequence of GRMHD time-slices, tends to have limited influence on the average images, but may alter the smoothness of the individual images and light curves (e.g. Dexter et al. 2010; Mościbrodzka et al. 2021).

Mock 228 GHz EHT images of M87^∗ were constructed by the ray-tracing techniques assuming a black hole mass M_BH = 6.5 × 10⁹ M_⊙, a distance D = 16.8 Mpc³, and inclination angles close to a face-on view of the accretion disk (i = 17° ,163°). The scaling of the model densities and magnetic field strength was determined by the averaged compact flux density of roughly 0.5 Jy at 228 GHz (Paper I; M87^* 2018 I). The field of view and resolution of the mock images usually were 200 μas × 200 μas and 400 × 400 pixels⁴.

The GRMHD simulations were used to model a collisionless plasma in which the ions and electrons were weakly coupled by Coulomb collisions, but were able to undergo partial relaxation due to wave-particle interactions (e.g., Kunz et al. 2014). It is computationally expensive to evolve ion and electron temperatures separately. Therefore, with the exception of the radiative models, most of our GRMHD simulations were single-fluid simulations, in which the electron temperature was not calculated directly. Instead, the electron temperatures, or the nonthermal electron distribution functions (eDF), were parameters of the post-processing radiative transfer model. In what follows, we describe the assumptions for the radiative transfer post-processing of the fiducial models and the exploratory models in more detail. The radiative transfer models and their parameters are summarized in Table 2.

All ray-traced images excise the jet region with high magnetization parameter, σ >  σ_cut (σ ≡ 2P_mag/ρc²) as GRMHD simulations do not accurately evolve the internal energy in these near-vacuum regions. The value of σ_cut varies between pipelines. For example, we present in Fig. 2 the representative GRMHD simulations from the KHARMA simulation. The contours for σ = σ_cut (σ_cut = 1 for KHARMA-ipolepipelines) are indicated by the solid black lines in the poloidal plane plots.

3.1. Fiducial thermal and nonthermal models

Our fiducial GRMHD simulations from the KHARMA-ipole, BHAC-BHOSS, and H-AMR-BHOSS pipelines assumed that radiating electrons have a thermal relativistic Maxwell-Jüttner distribution function. In these models, we assumed that the ion and electron temperatures are coupled as a function of the plasma-β using the R_high − R_low prescription of Mościbrodzka et al. (2016), which reads

$$\begin{array}{c}\hfill \frac{T_i}{T_e}=R_{\mathrm{low}}\frac{1}{1+{\beta }^2}+R_{\mathrm{high}}\frac{{\beta }^2}{1+{\beta }^2},\end{array}$$(2)

where R_low and R_high are model parameters (the surveyed values we used are listed in Table 2), and the ion temperatures are given by the GRMHD models (see also Fig. 2). Representative time-averaged images of the fiducial thermal models are shown in Fig. 3.

Fig. 2. Poloidal (azimuthally averaged) and equatorial plane cuts of selected KHARMA GRMHD simulation snapshots of the two considered accretion modes, both with a* = +0.5. The MAD (left panels) and SANE (right panels) snapshots feature different ion temperatures (displayed in dimensionless units in the upper panels) and the plasma β parameter (displayed in the bottom panels), resulting in different emission properties. The solid black lines represent the σ = 1 surface. The high temperature near the polar axes is a consequence of conservative codes. The temperatures in the polar regions are considered unphysical and are masked out in the radiative transfer calculations (regions with σ >  1 or σ >  3). Here, rg ≡ GM/c2, the gravitational radius.

Fig. 3. Time-averaged images of all the fiducial models with Rlow = 1 and Rhigh = 40 in the image library. The white arrows indicate the projected direction of the black hole spin for the simulation geometry. The lower left corner is blank because there is no H-AMR SANE a* = −0.94 case in our model library. All images are shown in linear scale.

Images that include radiation from nonthermal populations of electrons were constructed using the BHAC-BHOSS pipeline. These models have a nonthermal κ distribution, which has a thermal core with a width w that smoothly matches a power-law tail with dn_e/dγ ∝ γ_e^{−κ − 1} for γ_e ≫ w. The width w was partly determined from the R_low − R_high prescription above, but the model includes a parameter ϵ that represents the fraction of magnetic energy that contributes to the electron acceleration (first introduced by Davelaar et al. 2019; see also Cruz-Osorio et al. 2022; Fromm et al. 2022; Davelaar et al. 2023 for details) and contributes to the width w in a region near an injection radius that we took to be 10 GM/c². We considered models with ϵ = 0 and 0.5. The parameter κ was determined by β and σ following the prescription of Ball et al. (2018). Example nonthermal models are shown in the fourth and fifth panels in Fig. 3. The image morphology for models with ϵ = 0 (fourth panel) is similar to that of the thermal models (third panel) from the same pipeline. With increasing fraction of nonthermal emission (ϵ = 0.5; fifth panel), the resulting image is more extended, which agrees with Mao et al. (2017), for example.

3.2. Exploratory models

The tilted models were adopted from Chatterjee et al. (2020, see also Liska et al. 2018) and evolved using H-AMR code for a considerably longer time than the fiducial simulations. In tilted models, the accretion flow angular momentum and black hole spin are no longer coaxial, but are instead separated by a tilt angle. We considered tilt angles 0° (as a reference model), 30°, 60°, a single black hole spin (a_* = +0.94), and an initial magnetic field setup between the SANE and MAD setup with ϕ_BH ≈ 7 − 14 (see Chatterjee et al. 2020 for details). The tilted-disk runs adopted Eq. (2) as a prescription for the electron temperatures. The tilted-disk models imaged for one inclination and three azimuthal angles are listed in Table 2. Examples of time-averaged images of tilted disks are shown in Appendix B.

Our library also includes radiative two-temperature models from Chael et al. (2019a) that were simulated using the KORAL-grtrans pipeline. These simulations evolved a two-temperature magnetized fluid and a radiation field as a second fluid, coupled to the plasma using the M1 approximation. These radiative GRMHD (GRRMHD) simulations account for self-consistent radiation physics, incorporating both particle heating and radiative cooling via bremsstrahlung, synchrotron, Compton, and Coulomb losses. Both KORAL simulations in our library are MAD and assume a_* = +0.94. The two simulations feature identical initial conditions, but differ in their assumed prescription for subgrid electron heating: turbulent Landau damping (Howes 2010), or magnetic reconnection (Rowan et al. 2017). The radiative simulations were not scale-free. The gas density was scaled to physical units such that the compact emission at 230 GHz is roughly 0.98 Jy which is larger than the 0.5 Jy assumed in the fiducial simulations. This flux normalization was chosen based on pre-EHT constraints on the compact flux, but since the computational cost can be higher by an order of magnitude than a typical fiducial model, we chose not to produce new radiative simulations for this analysis and instead used what was already available. As we discuss in Sect. 4, we only scored on closure quantities, which do not contain information about the total flux. Examples of time-averaged images of radiative models are shown in Appendix B.

3.3. Model uncertainties and library summary

The simulations have several uncertainties that may impact the EHT data interpretation. If the ions and electrons are thermal and the electron temperature T_e is lower than the ion temperature T_i, then the ions dominate the pressure. Since the ions are nonrelativistic (kT_i/(m_ic²)< 1), this implies that Γ_ad = 5/3. If T_e = T_i, as in R_high = 1 models, and in R_low = 1 models when β ≲ 1, then both components contribute to the pressure, and we can show that Γ_ad = 13/9. Some of our models used Γ_ad = 4/3, but this was purely for numerical convenience because the Γ_ad = 4/3 models are more robust. A change of the adiabatic index in the GRMHD model changes the temperature profile and model images. For example, the KHARMA and H-AMR, MAD or SANE simulations all assumed a different Γ_ad index and therefore do not necessarily have to produce statistically consistent images. Moreover, even when the simulations assumed the same adiabatic index, the radiative transfer codes may interpret it in a different way when they calculate T_e (see Appendix H in Sgr A^* 2024 VIII). It is therefore possible that even simulations with the same Γ_ad, for instance, the KHARMA and BHAC MAD models, produce different images. Moreover, models with σ_cut = 3 have slightly lower optical depths (they are slightly brighter and therefore require a smaller density normalization) than the same models with σ_cut = 1, which introduces another discrepancy. Finally, the initial torus size or the grid resolution may affect the results also to some small degree. The effects of changing the equation of state, σ_cut, and other parameters in model-data comparisons have yet to be fully investigated, and we therefore treat the variations between the models that are produced by different simulation pipelines as a measure of the model uncertainty. We expect that the same accretion state models produced by different simulation pipelines may not always produce consistent results.

The total number of fiducial model images is ∼170 000 and the total number of images from tilted simulations is ∼10 000, including ∼2000 snapshots of the reference nontilted H-AMR simulation. Our extended image library consists of ∼184 000 snapshots, and it is three times larger than the library assembled in M87^* 2017 V.

4. Comparing GRMHD models with data

4.1. Overview of the data selection and comparison

The EHT collected multiple epochs of M87 data in two frequency bands in 2017 and four frequency bands in 2018. We focused our comparisons on the EHT observations centered at 229 GHz (high-band in 2017 and band 4 in 2018). To obtain uncorrelated constraints, we focused our comparisons on two selected days from each campaign: April 6, 2017, and April 21, 2018. April 6, 2017, was chosen following M87^* 2017 V based on the highest number of scans in the 2017 datasets and because of the four days of observations in 2017, April 5/6 differ most from the 2018 image. April 21, 2018, was chosen because it has the best (u, v)-coverage in the 2018 campaign. The dynamical timescale of M87^∗ near the event horizon (GM/c³ ∼ 9 h) is comparable to the single-epoch observation, and no intra-day variability of M87^∗ is therefore expected. The images are expected to be very strongly correlated below timescales of about 50 GM/c³ ∼ 20 days.

To be consistent with M87^* 2017 V and M87^* 2017 VI, we first compared the models to the observations using two previously developed algorithms: the snapshot scoring, and the average imaging scoring (hereafter, AIS). To effectively perform the scoring procedures on a significant number of model images, a few improvements were made to reduce the computational cost of the calculations. First, we applied an optimization approach instead of the MCMC method with gain fitting (Broderick et al. 2020) (also adopted in M87^* 2017 V; M87^* 2017 VI). Second, we fit the closure quantities of the VLBI data, such as the closure phases and the closure amplitudes. Third, the scoring was performed with an efficient Julia-based driver developed following Tiede (2022). The snapshot and AIS scoring procedures were performed separately in the 2017 and 2018 observations. In addition to this, we also introduced the multi-epoch scoring, which combines the scoring results for individual observations from different years. Finally, in addition to the EHT constraints of the three aforementioned scoring schemes, the M87 jet power in the models was required to exceed a lower limit of 10⁴² erg s⁻¹. The following subsections briefly describe each of the scoring schemes and their results. The AIS and multi-epoch scoring methods are described in depth in Appendix C.

4.2. Snapshot scoring constraints

The snapshot scoring approach determines how consistent each model snapshot is with the observation. In this procedure, we fit a given model with the EHT data by rotating and resizing every image in the model to find the best-fit values of the image orientation and the implied mass-to-distance ratio, M/D. The hundreds of images from each model are then used to estimate an ensemble-based posterior.

The mean values of the reduced χ², ${\chi }_{\nu }^2$, of the fiducial model fits are ∼10 (for the 2017 datasets) and ∼15–30 (for 2018 datasets). Values much higher than unity are expected because a single model image is unlikely to fit the observation perfectly, given the stochastic character of the turbulence and the finite simulation time. It is expected that the 2018 ${\chi }_{\nu }^2$ values are higher, since the newly added EHT baselines (formed by adding the Greenland Telescope to the EHT network) add more information to the (u, v) domain, while the single-snapshot model retains the same model complexity as in the 2017 analysis.

For exploratory models, the ${\chi }_{\nu }^2$ distributions for the radiative models and tilted-disk models with zero tilt angle are comparable to their distributions from fiducial models. We find that tilted-disk models with a tilt angle of 30° and 60° clearly have much higher values of ${\chi }_{\nu }^2$ than the other models. This is likely due to the more crescent-like, instead of ring-like, morphology of the tilted images (compare Fig. 3 with Fig. B.1).

Figure 4 shows the posterior distributions M/D for representative fiducial models fit to the 2017 (blue) and 2018 (red) datasets. The distributions shown in Fig. 4 include the results for all geometries available in the image library: for cases where the black hole spin axis pointed away from Earth, and for cases where it pointed toward Earth (see also Table 2). The distributions from different years overlap each other, but the 2018 distributions peak around slightly higher values of M/D for all simulation pipelines. Similar to the result in M87^* 2017 V, the distributions for both EHT 2017 and 2018 observations are more consistent with the mass estimation using stellar dynamics studies (e.g. Gebhardt et al. 2011; Liepold et al. 2023; Simon et al. 2024) than with the gas dynamical model (e.g. Walsh et al. 2013).

Fig. 4. Violin plots for M/D distributions from the snapshot scoring with 2017 and 2018 EHT observations for models with Rlow = 1 and Rhigh = 40 in the image library. The horizontal dashed black line and shaded gray region mark the ranges of M/D = 3.8 ± 0.4 μas reported in M87* 2017 I.

In Fig. 5 we rotate each snapshot image with its best-fit value of the PA and compare the distribution of the best-fit PAs with the observed forward-jet direction. This rotation is possible based on the prior knowledge of the black hole rotation axis and of the forward-jet axis of each image in our image library, as described in Sect. 3. Figure 5 shows the resulting distribution of the position angles of the models using only the 10% of images with the lowest χ_ν² from the KHARMA and H-AMR pipelines. We split the models into prograde/retrograde cases and considered both models with a spin axis pointing toward and pointing away from Earth (see discussion in Sect. 3).

Fig. 5. Normalized distributions of the fitted position angle of the forward jet in selected fiducial thermal models based on the snapshot scoring with 2017 and 2018 observations, using the best-fit 10% of the images. The PA = 288° ±10° of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023) is marked by the vertical dashed black line and shaded gray region. Left: Models with a black hole spin axis or a coordinate system axis (case a* = 0) pointing toward Earth. Prograde cases are shown in the top (a* >  0) and middle (a* = 0) panels, and the retrograde case (a* <  0) is shown in the bottom panel. Right: Models with a black hole spin axis or with a coordinate system axis (case a* = 0) pointing away from Earth. Prograde cases are shown in the top (a* >  0) and middle (for a* = 0 the matter angular momentum is always aligned with the coordinate axis) panels, and the retrograde case (a* <  0) is shown in the bottom panel. See text in this paper or Fig. 5 in M87* 2017 V for a description of the geometry of the black hole axis, flow rotational axis, and inclination angle.

Figure 5 shows that the jet directions in fiducial models with a black hole spin axis or angular momentum axis (in case of a_* = 0) that points away from Earth are significantly better aligned with the large-scale jet position angle (PA = 288° ±10° of the milliarcsecond-scale jet measured by Walker et al. (2018) and Cui et al. (2023), marked with a vertical dashed line) than in the opposite cases. This is true for the 2017 and 2018 scores, but compared to 2017, the 2018 preferred PA distributions are overall significantly better aligned with the large-scale jet. When we assume that the jet axis and black hole spin axis are aligned, our GRMHD simulations imply that the position angle of the brighter side of the ring structure in 2018 (southwest part of the ring) is roughly perpendicular to the jet direction. The distribution functions of the retrograde models are often broader than those of the prograde models, indicating that rings in retrograde models produce a greater angular variability.

The best-fit 2017 PA distributions presented here agree with those reported in M87^* 2017 V. The forward-jet PA distributions from the snapshot scoring of the previous smaller model image library from M87^* 2017 V is presented in Fig. 20 of Paper I (M87^* 2018 I).

4.3. Average imaging scoring constraints

The AIS method aims to determine how likely it is that the EHT data are consistent with a random draw from a specific model image distribution. As demonstrated in the previous subsection (and in M87^* 2017 V), individual model snapshots are unlikely to provide an acceptable ${\chi }_{\nu }^2$ fit to the observational data because of the finite number of snapshots and turbulence in the underlying GRMHD simulations.

In M87^* 2017 V, we used the AIS scheme to compare the GRMHD models with EHT data that overcome the fluctuations of individual snapshots. The scheme measures the distribution of χ_ν² distances between individual snapshots and the time-averaged image for each model. This distribution was compared to the distance between the data and the time-averaged image for the same model. The model was then assigned a probability based on where the data lay in the distribution. The model was disfavored when the data lay outside or on the tails of the model distribution.

AIS scoring was applied to all models using both 2017 and 2018 EHT data. We considered models whose 2017 AIS or 2018 AIS probabilities were higher than a certain threshold as equally good. We note that the identification of good models that passed the AIS scoring was based on only two observations, so caution should be taken when interpreting these results. The AIS results for the fiducial models are summarized in the Appendix D. As the threshold for AIS is arbitrary, we present the AIS results for which we chose a probability of 15% or 10% as the threshold⁵ for passing the models. The AIS scoring results, together with the snapshot scoring results, are necessary ingredients of the multi-epoch scoring procedure introduced next.

4.4. Multi-epoch scoring constraints

The multi-epoch scoring approach aims to provide an odds ratio that computes the probability of the model relative to the probability of the best-performing model. The typical correlated variability time-scale of all the GRMHD simulations (∼30 GM/c³, ∼10 days for M87^*) is much shorter than the observational cadence between the EHT 2017 and 2018 observations, and therefore, we can consider the 2017 and 2018 datasets as independent measurements.

In brief, the multi-epoch scoring post-processes the snapshot scoring and AIS scoring results following a Bayesian approach in which AIS results are treated as approximations for Bayesian evidence, and the snapshot scoring results are used as approximations of the ensemble-based posteriors⁶. The details of the procedure are presented in Appendix C. The odds ratio returned by the multi-epoch scoring provides a measure of the relative preference for two models. In our procedure, the returned odds ratio for a model would be relatively higher if the AIS scores are higher, and the distribution according to the snapshot scoring has a larger overlap (see Eq. (C.16) for a mathematical description). As a result, the multi-epoch scoring is a compromise of how a model could (i) satisfy the 2017 AIS scoring, (ii) satisfy the 2018 AIS scoring, and (iii) explain the 2017 and 2018 observations with same M/D ratio and the same direction of the black hole rotation axis.

An example is given in Fig. 6 as a demonstration of point (iii) in the multi-epoch scoring procedure. In the parameter space of M/D and the PA of forward-jet direction, the overlap of the ensemble posteriors for the 2017 observation and the 2018 observation (indicated by the black contours and profiles) represents how likely the independent observations are to be explained by the same parameters. Points (i) and (ii) mentioned above were also included as conditions when we computed the odds ratio for each model. We then compared the relative odds ratios of the models and selected the models that were preferred by the multi-epoch observations.

Fig. 6. Demonstration of the multi-epoch scoring based on the KHARMA SANE model with (a*, Rlow, Rhigh)=(−0.5, 10, 160). The 2017 and 2018 snapshot scoring results for the given model are shown in the M/D – PA plane with the blue and red distributions, respectively. The larger the overlap between the red and blue distributions, the better the fit of the multi-epoch observations with a single model. The black profile indicates the normalized combined distribution of the best-fit parameters computed according to the multi-epoch scoring procedure (see Sect. 4.4 for details). The dashed lines and shaded areas are the same as explained in Figs. 4 and 5.

Figure 7 compares the odds ratios for all the models with either black hole spin axes pointing away from Earth or a nonspinning black hole, but the accretion flow spin axis points away from Earth. The scores shown by the vertical axis indicate the relative preference for each model (i.e., the normalized odds ratio according to the specific model with the highest score; therefore, all but one model have relative scores lower than 100%). We selected 1% of the normalized odds ratio as the arbitrary threshold for the passing models. The low threshold value is motivated by all the uncertainties associated with the modeling of the images.

Fig. 7. Overview of the multi-epoch scoring results for all the fiducial models with black hole axes pointing away from Earth, and nonspinning black holes with a inclination angle i = 163°. The y-axis is the normalized odds ratio (see text for details), which indicates the relative preference by multi-epoch (2017 and 2018) observations. Models with different spins are shown in different panels. The green region marks passing models for which the relative odds ratio is ≥1%. Models with a relative odds ratio ≤0.001% are not shown.

Figure 7 shows models from all three pipelines with high relative scores, and MAD models reach higher scores overall than SANE models. Moreover, the multi-epoch scoring rules out more of the prograde models than the retrograde models. This can be explained by the more turbulent character of the retrograde models, which results in broader PA distributions, as shown in Fig. 5, and a better overlap between two observational epochs. Interestingly, tilted-disk models (which have positive black hole spins a_* = +0.94) have relatively low scores. The models that passed multi-epoch scoring were further combined with jet-power constraints.

4.5. Jet-power constraint

Previous studies of M87 estimated its jet power to be P_jet = 10⁴² − 10⁴⁵ erg s⁻¹ (Reynolds et al. 1996; Li et al. 2009; de Gasperin et al. 2012; Broderick & Tchekhovskoy 2015; Prieto et al. 2016). Following M87^* 2017 V, we computed the powers of outflows in our GRMHD simulations and compared them with observations. We adopted the same criterion as in previous EHT publications for the failing models with an insufficient jet power, that is, a model fails when its outward energy flux (outflow power; P_out) over the black hole polar regions is lower than 10⁴² erg s⁻¹.

In our models, the Δt-averaged P_out is defined as

$$\begin{array}{c}\hfill P_{\mathrm{out}}\equiv \frac{1}{\mathrm{\Delta }t}{\displaystyle \int d}t{\displaystyle \int d}\varphi {\displaystyle {\int }_{\mathrm{poles}}}d\theta \sqrt{-g}(-{T^r}_t-\rho u^r),\end{array}$$(3)

where T_t^r is a component of the stress-energy tensor associated with the outward radial energy flux, and ρu^r subtracts the rest-mass energy flux, g is the metric determinant, u^r is the gas four-velocity in radial direction, and the other variables have their usual meaning. In Eq. (3), ϕ and θ are the azimuthal and polar angles, respectively, and we only included those θ whose time- and azimuth-averaged energy flux was directed outward. There are many choices to define the cut on the polar region to integrate only sufficiently relativistic outflows, for example, a cut based on βγ, or the Bernoulli parameter values. Here, the polar angle integration includes 1 radian about the north or south poles and only those regions whose energy flux is directed outward. This includes both the relativistic narrow jet near the polar axis and a portion of the nonrelativistic outflows at higher inclinations. P_out can therefore be larger than the Blandford-Znajek jet power associated with the relativistic part of the jet, and models with zero spin may have P_out >  0. In some cases, for instance, SANEs, the P_out measurement can be dominated by the nonrelativistic outflows. P_out also depends on the selected mass unit (ℳ), which was used to set the ∼230 GHz total flux, overall density, and magnetic field strength in the simulations. Therefore, models with very cold electrons will have higher ℳ and proportionally higher P_out. However, the variations in P_out in the imaging parameters (inclination and R − β prescription) are smaller than the variations in the GRMHD simulation parameters (spin and accretion mode). For all of the fiducial models, P_out is a strong function of |a_*|.

The jet-power scoring results for the fiducial MAD and SANE models are shown in combination with the multi-epoch scoring in Table 3. For all codes, the a_* = 0 MAD models fail the power constraint most often (16 out of 20 fail). SANE models with R_high = 1 also often fail, regardless of spin. These results are consistent with those reported in M87^* 2017 V. P_out in the fiducial nonthermal MAD models show the same trend as in the fiducial thermal models, but all a_* = 0 MAD models fail. All spinning nonthermal MAD models pass the jet-power constraint.

For exploratory models, all tilted models pass the jet power constraint. This is an expected result as the tilted models assumed high spins (a_* = +0.94), for which most of the aligned fiducial models pass as well (with the exception of the few aforementioned SANEs). In addition, both high-spin radiative MAD models produce jet powers of P_jet ∼ 10⁴³ erg s⁻¹ (see Table 1 of Chael et al. 2019b), and thus, they pass the jet-power constraint⁷.

5. Discussion

5.1. Summary of the model constraints and best-bet models

Table 3 displays the pass and fail information for fiducial models in which their black hole spin axes point away from Earth, or, when the black hole is nonrotating, the flow spin axes point away from Earth. As discussed in Sect. 4.2, models oriented in this way are more probable given the direction of the large-scale jet in M87^∗, assuming a coaxial model. The pass and fail table combines two constraints: the multi-epoch scoring, and the jet power. Models that pass (or fail) these two constraints are indicated by dark green (or yellow), and models that pass only one constraint are indicated by other colors. We refer to the models that pass both constraints as the best-bet models. The exploratory models (all with a_* = 0.94) are not shown in Table 3: The two radiative models are also among the best-bet models, but all tilted disks have a relatively low multi-epoch score below the chosen cutoff. For reference, the M/D distributions for radiative models are presented in Appendix E. The M/D distributions for tilted-disk models are not shown since these models in general have higher ${\chi }_{\nu }^2$ (see Sect. 4.2), and they perform poorly in all other scoring constraints.

Most of the best-bet models pass both the 2017 AIS and 2018 AIS scoring, as expected from the multi-epoch scoring introduced in Sect. 4.4. Most models that fail both 2017 AIS and 2018 AIS also fail the multi-epoch scoring (see also Table D.1). However, some models pass both the 2017 AIS and 2018 AIS scoring, but fail the multi-epoch scoring because the multi-epoch scoring returns an odds ratio that depends not only on AIS scoring, but also on a constant M/D and PA to explain the observation. This demonstrates that our multi-epoch scoring procedure is more powerful than a logical AND applied to the 2017 AIS and 2018 AIS scoring results alone. For the jet-power constraints, models with a_* = 0 are generally less favored than spinning models.

The 68 best-bet models of the 194 fiducial models are shown in Table 3. The best-bet models include all types of fiducial models: thermal MAD and SANE models, and nonthermal MADs. However, the survival rate for MADs is higher than that for SANEs. For certain spins, the survival rate of the nonthermal MAD models does not strongly depend on the details of nonthermal physics (e.g., models with ϵ = 0 and ϵ = 0.5 score similarly for a_* = ( − 0.95, −0.5, 0)). Among all the best-bet models, retrograde models have a relatively higher survival rate than all other cases. We propose that turbulence might account for the variation in the brightest segment of the ring-like structure. This conjecture finds further support in the fact that retrograde models, characterized by higher levels of turbulence near the event horizon compared to prograde models (see also Fig. 5, where the PA distributions in retrograde models are often broader than in prograde models), match the observations better. In addition, the M/D distribution for 2017 and 2018 shown in Fig. 4 also overlaps better for negative spin cases.

Not all pipelines exhibit the same pass and fail results across the flow type, spin, and thermal electron parameters. For example, the KHARMA MAD, a_* = +0.94, R_low = 1, R_high = 1 model only passes the jet-power constraint, while the H-AMR model with the same parameters passes both the jet power and multi-epoch scoring constraints. The discrepancy may be due to differences in the numerical scheme in the code itself, to the simulation setup, and to post-processing details (see Tables 1 and 2). The models are also subject to sampling errors: the snapshots are drawn from a limited number of correlation times, meaning that acceptable models may be occasionally rejected.

Table 3 also demonstrates that the multi-epoch observations provide more constraining power, because more models fail only the multi-epoch-constraints than fail only the jet constraint. In M87^* 2017 V, AIS scoring based on the 2017 M87^∗ EHT observation was applied as one of the constraints for the model selection. This single-year EHT constraint has less constraining power than the jet-power constraint (see Table 2 of M87^* 2017 V). That is, the jet-power constraint can rule out more models than are ruled out by the 2017 AIS results. With the multi-year M87^∗ EHT observations, based on Stokes I properties of the model images, we were able to provide EHT constraints with a relative more constraining power than the jet-power constraints. The increased constraining power with multi-year observations is closely related to the fact that the ring in M87^∗ looks substantially different in 2017 and 2018.

Our exploratory models are sampled more sparsely than the fiducial models. There is a larger parameter space left to be explored for these models in the future (e.g., they currently all have spin a_* = +0.94). A more comprehensive radiative image library from two-temperature GRRMHD simulations with different black hole spins and flow types would be valuable considering the relatively good performance in this scoring. In the tilted models, even the reference zero tilt-angle model is not selected as one of the best-bet models. We recall that the tilted simulations have an intermediate magnetic flux between MAD and SANE, and therefore, the zero-tilt exploratory model is distinct from any of the fiducial models.

5.2. Multiwavelength spectra of the selected best-bet models

We determined whether the best-bet models are consistent with the multi-wavelength observations of M87^∗. We calculated broadband spectral energy distributions (SEDs) for a few selected best-bet thermal synchrotron models. In Fig. 8 we present the SEDs for three best-bet models along with radio and X-ray observations collected in 2018 (EHT MWL Science Working Group 2024). The selected models are taken from the KHARMA pipeline. They include SANE with (a_*, R_low, R_high)=(0.94,  10,  160), SANE with (a_*, R_low, R_high)=(−0.5,  1,  160), and MAD with (a_*, R_low, R_high)=(−0.5,  10,  40). The SEDs were generated using the Monte Carlo relativistic radiative transfer code grmonty (Dolence et al. 2009), which includes synchrotron (emission and self-absorption), inverse-Compton scattering, and bremsstrahlung emission. In the EHT bands, the emission is dominated by the synchrotron process. The X-ray part of the SED is dominated by either Compton upscattering of synchrotron photons or by bremsstrahulung emission. In general, both the Compton and bremsstrahulung components are relatively more prominent in SANEs than in MADs. The observed X-ray emission can be contributed by scales that far exceed the size of our model (EHT MWL Science Working Group 2024). When we treat the observed X-ray luminosity as an upper limit for the models, the best-bet models shown in Fig. 8 are all acceptable. Among the three SEDs, the time-averaged X-ray emission of one of the SANE models,(a_*, R_low, R_high)=(0.94, 10, 160), matches the observations, but in this case, the synchrotron peak is at much higher frequency than in the EHT bands. We also note that the flux in the optical and X-ray bands may increase when the nonthermal synchrotron component is considered as well.

Fig. 8. Time-averaged SEDs of selected best-bet models. The data points are taken from the multi-wavelength observations during the 2018 EHT campaign (EHT MWL Science Working Group 2024). The EHT observation is marked by the vertical red bar. The gray region of the SED indicates the variations of the SED for different snapshots. The colored histograms correspond to different radiative processes: synchrotron emission (blue), synchrotron photons scattered once (orange), synchrotron photons scattered twice (green), and bremsstrahlung (red). The total emission is displayed in black.

5.3. Mass-to-distance ratios and implications for gravitational physics

As discussed in Sect. 4.2, models in which the black-hole spin axis points away from Earth together with nonspinning black holes with an inclination angle i = 163° are preferred by both 2017 and 2018 observations. In Fig. 9 we show the mass-to-distance ratio (M/D) distributions for each GRMHD pipeline and for all pipelines combined. These distributions only include the top 10% best χ² images from the snapshot scoring procedure. The 2017 and 2018 constraints are represented with the blue and red distributions, respectively. For each model library, the 2018 results (red distribution) prefer slightly higher values for M/D than the results from the 2017 observations (blue distribution). This preference for a slightly higher value for M/D is consistent with the results from the image reconstructions in Paper I (M87^* 2018 I), which prefer slightly larger diameters and slightly smaller ring widths for the 2018 data. The 2018 data prefer $M/D=4.{17}_{-0.66}^{+0.41}\mu \mathrm{as}$, and the 2017 data prefer an $M/D=3.{70}_{-0.59}^{+0.51}\mu \mathrm{as}$. When we convert the above M/D values into units of solar masses at a distance of 16.8 Mpc, the mass for M87^∗ using the updated 2018 EHT observations is $M=7.{09}_{-1.11}^{+0.69}\times {10}^9{\mathrm{M}}_{\odot }$, and it is $M=6.{30}_{-1.01}^{+0.87}\times {10}^9{\mathrm{M}}_{\odot }$ for the 2017 observations.

Fig. 9. Distribution of M/D obtained by fitting the 2017 (blue) and 2018 (red) data with the fiducial models shown in Table 3. These distributions are generated from the top 10% best-fitting snapshots by χ2. The results for models prepared by the KHARMA, H-AMR, and BHACpipelines are shown separately, together with the result for all the models combined. The distribution when combining both the 2017 and 2018 result by the multi-epoch scoring algorithm is shown as the black profile. The vertical dashed black line and shaded gray region mark the ranges of M/D = 3.8 ± 0.4 μas reported in Event Horizon Telescope Collaboration (2019a).

Applying the multi-epoch scoring algorithm, we also constructed a combined M/D distribution (black distribution) for all the fiducial models (see also Fig. 6 for comparison). The combined distribution prefers an $M/D=3.{62}_{-0.34}^{+0.41}\mu \mathrm{as}$ ($M=6.{16}_{-0.57}^{+0.71}\times {10}^9{\mathrm{M}}_{\odot }$). We summarize the mass estimates from the individual years and the multi-epoch results in Table 4. The slightly lower mass estimate when incorporating the multi-epoch constraints can be attributed to the general preference for retrograde models in the scoring results. The retrograde models generally produce less compact images, so that the snapshot scoring procedure tends to shrink the simulated images to match the size of the emission region in the real data, leading to lower mass estimates. This bias may disappear in the future when we also incorporate additional conditions in the image structure in the snapshot scoring procedure beyond rescaling and rotation. Dedicated image reconstruction techniques may also be able to incorporate a more complicated image structure and differentiate between models with less compact emission. These new M/D values are statistically consistent with the result published in 2019, that is, $M=6.{50}_{-0.7}^{+0.7}\times {10}^9{\mathrm{M}}_{\odot }$ (M87^* 2017 I; M87^* 2017 VI), and it differs by about 1.5 − 2σ from the mass estimates made in Liepold et al. (2023) and by slightly more than 1σ from the mass estimates found in Simon et al. (2024).

Notwithstanding the slight differences between the new mass estimate and that published in 2019, it is also important not to overinterpret them. Numerous changes have been introduced to the GRMHD pipelines since the 2019 analysis, as described in Sect. 3, and it is expected that the new library produces some differences from the previous library in the mass estimate. Even within the new model library, significant differences remain between the model pipelines beyond the black hole spin, electron temperature prescription, and magnetic field configuration that change the M/D distributions. In addition, the coarse sampling of these parameters within a given pipeline can produce significant spikes in the multi-epoch M/D distributions that come from a small number of high-performing models (see, e.g., the black KHARMA curve). Nevertheless, when we consider either the new individual or the combined mass estimates, the differences between the M87^* 2017 VI result and the mass results presented in this work are not statistically significant. In view of this, we regard the new M/D and mass estimates to be consistent with those published in the 2019 papers (M87^* 2017 I; M87^* 2017 VI). As anticipated above, this conclusion is also expected when we assume that M87^∗ is a supermassive black hole described by general relativity, whose corresponding spacetime is either static or stationary, and hence cannot exhibit a significant difference in its properties over a timescale of only one year.

6. Conclusions

In M87^* 2018 I, we presented the results of the 2018 EHT observations of M87^∗ at 1.3 mm. Leveraging the first 2017 and the newest 2018 EHT observations of M87^∗, we here compared the source Stokes I observations to the upgraded library of GRMHD simulations of black hole accretion. The new model library consists of fiducial thermal synchrotron emission models and compares for the first time tilted-disk models and radiative GRMHD models.

First, treating 2017 and 2018 observations as independent experiments, we presented the results of the snapshot scoring and average image scoring (AIS) procedures. Assuming that the black hole angular momentum axis is aligned with the jet axis, we find that the observations prefer a spin axis of the black hole that points away from Earth.

Next, we presented a Bayesian-based multi-epoch scoring procedure, within which the AIS results were used to approximate the Bayesian evidence, and the snapshot scoring results were used to approximate the ensemble-based posterior construction. In the multi-epoch scoring algorithm, the success of a model is also evaluated based on whether the same M/D ratio and black hole rotation axis in the model is consistent with both years of data. As a result, the multi-epoch scoring was performed as a post-processing step using the results of both the snapshot and AIS procedures. The multi-epoch scoring produces better constraints than the single-year scoring, for instance, the AIS scoring alone. The multi-epoch scoring finds that retrograde models have a higher survival rate than prograde systems, which might be related to their comparatively more turbulent nature. The main cause of the stronger variability of the retrograde flows may require further examinations. For example, the variability in the PA in the retrograde flows does not seem to be correlated with the accretion rate onto the black hole (Narayan et al. 2022).

The current analysis does not allow us to explain the M87 flaring in VHE reported in EHT MWL Science Working Group (2024) with the near horizon emission because passing purely thermal models do not produce emission at these energies. Future simultaneous observations at higher frequencies (∼THz) could help us to distinguish an accretion state or black hole spin in M87^∗ probably only for the nonthermal models.

The recent reanalysis of historical radio (86 GHz) observations of the M87 core indicated a quasi-periodic change in the 288° position angle of the large-scale jet by about 10° degrees with a period of 11 years (Cui et al. 2023). The tilted models were disfavored in this work, but current results do not yet have significant implications for the jet precession. A correlation of the jet motion with the variation in the horizon-scale structure would require many more EHT epochs, together with low-frequency VLBI monitoring, preferably sampling timescales shorter than one year. Nevertheless, for future tilted model simulations, it would be useful to explore the effect of a tilted accretion flow around black holes with various spins and with a wider range of tilt angles.

The sensitivities of the 2017 and 2018 EHT arrays are insufficient to reveal jet-related extended emission near the black hole with confidence (e.g. Broderick et al. 2022), although future EHT observations for M87^∗ with adequate dynamical range may show further details of an extended component in addition to the near-horizon ring structure. In this case, several emission features related to the jet physics may also be tested and constrained (e.g. Broderick & Tchekhovskoy 2015; Pu et al. 2017; Jeter et al. 2020). We also note that including the morphology constraint of M87^∗ at lower frequencies can aid in identifying the relative contributions of the emission from the accretion and jet components of M87^∗ (Lu et al. 2023), which would otherwise be challenging if based on the SED alone (Dexter et al. 2012).

Regardless of all these limitations, our current models can be still further constrained by their polarimetric properties (see, e.g., M87^* 2017 VII; M87^* 2017 VIII or M87^* 2017 IX). We postpone a full discussion of the multi-epoch polarimetric scoring to a forthcoming publication.

When we assume that the black hole spin axis is closely aligned with the large-scale jet axis, images of the black hole shadow from GRMHD simulations produce a distribution of brightness position angles centered roughly 90° offset from the jet position angle. If this is the case, then the brightness position angle of the 2017 images should be relatively rare when making many independent observations. For this reason, we predicted in the previous paper series (M87^* 2017 V) that future observations would most likely show a brightness position angle counterclockwise from the 2017 value. The 2018 observations appear to follow this prediction, but only represent one additional independent realization. The current multi-epoch analysis was performed on just two years of data, so that any differences between the years may be subject to concerns about small-number statistics. Future observations can help us to build confidence in our results. The analysis of the EHT data from later years (2021 and 2022) are already underway and will help us understand whether the brightness position angle of the 2017 observations is rare, following our predictions. In addition, a proper distribution of the brightness position angles in the observations can be a powerful tool for constraining the simulations by helping us reject simulations that return distributions of the position angle that are too wide or too narrow.

While the year-to-year analysis of EHT images is critical for understanding the uncorrelated variability of M87^∗, a dedicated monitoring campaign consisting of observations every few days over the course of one to two months is necessary to investigate the correlated variability around M87^∗. Observations with this high cadence would provide correlated constraints, unlike the independent constraints from the year-scale observation cadence considered in this paper. The EHT observations already revealed differences between the correlated variability seen in current GRMHD simulations and the variability seen in Sgr A* (the simulated variability is typically stronger than the observed variability in the total-intensity light curves; see Wielgus et al. (2022) and Sgr A^* 2017 V for further discussion), and a dedicated effort to investigate the correlated variability in M87^∗ will provide an independent constraint and valuable cross-check on models for the magnetic field and turbulent plasma environment around the black hole.

Acknowledgments

The full acknowledgements are available in Appendix A.

References

Appendix A: Acknowledgements

We thank the anonymous referee for helpful comments that improved the paper.

The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19_058 (TITANs), Fondecyt 1221421 and BASAL FB210003; the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the ALMA North America Development Fund; the Astrophysics and High Energy Physics programme by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1); the Black Hole Initiative, which is funded by grants from the John Templeton Foundation (60477, 61497, 62286) and the Gordon and Betty Moore Foundation (Grant GBMF-8273) - although the opinions expressed in this work are those of the author and do not necessarily reflect the views of these Foundations; the Brinson Foundation; “la Caixa” Foundation (ID 100010434) through fellowship codes LCF/BQ/DI22/11940027 and LCF/BQ/DI22/11940030; Chandra DD7-18089X and TM6-17006X; the China Scholarship Council; the China Postdoctoral Science Foundation fellowships (2020M671266, 2022M712084); Conicyt through Fondecyt Postdoctorado (project 3220195); Consejo Nacional de Humanidades, Ciencia y Tecnología (CONAHCYT, Mexico, projects U0004-246083, U0004-259839, F0003-272050, M0037-279006, F0003-281692, 104497, 275201, 263356, CBF2023-2024-1102, 257435); the Colfuturo Scholarship; the Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de Andalucía (grant P18-FR-1769), the Consejo Superior de Investigaciones Científicas (grant 2019AEP112); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Autónoma de México (DGAPA-UNAM, projects IN112820 and IN108324); the Dutch Research Council (NWO) for the VICI award (grant 639.043.513), the grant OCENW.KLEIN.113, and the Dutch Black Hole Consortium (with project No. NWA 1292.19.202) of the research programme the National Science Agenda; the Dutch National Supercomputers, Cartesius and Snellius (NWO grant 2021.013); the EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, Center for Astronomical Mega-Science, Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute; the European Research Council (ERC) Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (grant 610058) and Synergy Grant “BlackHolistic: Colour Movies of Black Holes: Understanding Black Hole Astrophysics from the Event Horizon to Galactic Scales” (grant 10107164); the European Union Horizon 2020 research and innovation programme under grant agreements RadioNet (No. 730562), M2FINDERS (No. 101018682) and FunFiCO (No. 777740); the European Research Council for advanced grant “JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales” (grant No. 884631); the European Horizon Europe staff exchange (SE) programme HORIZON-MSCA-2021-SE-01 grant NewFunFiCO (No. 10108625); the Horizon ERC Grants 2021 programme under grant agreement No. 101040021; the FAPESP (Fundação de Amparo á Pesquisa do Estado de São Paulo) under grant 2021/01183-8; the Fondo CAS-ANID folio CAS220010; the Generalitat Valenciana (grants APOSTD/2018/177 and ASFAE/2022/018) and GenT Program (project CIDEGENT/2018/021); the Gordon and Betty Moore Foundation (GBMF-3561, GBMF-5278, GBMF-10423); the Institute for Advanced Study; the ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU; the Istituto Nazionale di Fisica Nucleare (INFN) sezione di Napoli, iniziative specifiche TEONGRAV; the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne; the Italian Ministry of University and Research (MUR) – Project CUP F53D23001260001, funded by the European Union – NextGenerationEU; DFG research grant “Jet physics on horizon scales and beyond” (grant No. 443220636); Joint Columbia/Flatiron Postdoctoral Fellowship (research at the Flatiron Institute is supported by the Simons Foundation); the Japan Ministry of Education, Culture, Sports, Science and Technology (MEXT; grant JPMXP1020200109); the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Joint Institute for Computational Fundamental Science, Japan; the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJSSW-SYS008, ZDBS-LY-SLH011); the Leverhulme Trust Early Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max Planck Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP21H01137, JP18H03721, JP18K13594, 18K03709, JP19K14761, 18H01245, 25120007, 19H01943, 21H01137, 21H04488, 22H00157, 23K03453); the MICINN Research Projects PID2019-108995GB-C22, PID2022-140888NB-C22; the MIT International Science and Technology Initiatives (MISTI) Funds; the Ministry of Science and Technology (MOST) of Taiwan (103-2119-M-001-010-MY2, 105-2112-M-001-025-MY3, 105-2119-M-001-042, 106-2112-M-001-011, 106-2119-M-001-013, 106-2119-M-001-027, 106-2923-M-001-005, 107-2119-M-001-017, 107-2119-M-001-020, 107-2119-M-001-041, 107-2119-M-110-005, 107-2923-M-001-009, 108-2112-M-001-048, 108-2112-M-001-051, 108-2923-M-001-002, 109-2112-M-001-025, 109-2124-M-001-005, 109-2923-M-001-001, 110-2112-M-001-033, 110-2124-M-001-007 and 110-2923-M-001-001); the National Science and Technology Council (NSTC) of Taiwan (111-2124-M-001-005, 112-2124-M-001-014 and 112-2112-M-003-010-MY3); the Ministry of Education (MoE) of Taiwan Yushan Young Scholar Program; the Physics Division, National Center for Theoretical Sciences of Taiwan; the National Aeronautics and Space Administration (NASA, Fermi Guest Investigator grant 80NSSC23K1508, NASA Astrophysics Theory Program grant 80NSSC20K0527, NASA NuSTAR award 80NSSC20K0645); NASA Hubble Fellowship Program Einstein Fellowship; NASA Hubble Fellowship grants HST-HF2-51431.001-A, HST-HF2-51482.001-A, HST-HF2-51539.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555; the National Institute of Natural Sciences (NINS) of Japan; the National Key Research and Development Program of China (grant 2016YFA0400704, 2017YFA0402703, 2016YFA0400702); the National Science and Technology Council (NSTC, grants NSTC 111-2112-M-001 -041, NSTC 111-2124-M-001-005, NSTC 112-2124-M-001-014); the US National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, OIA-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1440254, AST-1555365, AST-1614868, AST-1615796, AST-1715061, AST-1716327, AST-1726637, OISE-1743747, AST-1743747, AST-1816420, AST-1935980, AST-1952099, AST-2034306, AST-2205908, AST-2307887); NSF Astronomy and Astrophysics Postdoctoral Fellowship (AST-1903847); the Natural Science Foundation of China (grants 11650110427, 10625314, 11721303, 11725312, 11873028, 11933007, 11991052, 11991053, 12192220, 12192223, 12273022, 12325302, 12303021); the Natural Sciences and Engineering Research Council of Canada (NSERC, including a Discovery Grant and the NSERC Alexander Graham Bell Canada Graduate Scholarships-Doctoral Program); the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752; the Korea Research Fellowship Program: NRF-2015H1D3A1066561; Brain Pool Program: RS-2024-00407499; Basic Research Support Grant 2019R1F1A1059721, 2021R1A6A3A01086420, 2022R1C1C1005255, 2022R1F1A1075115); Netherlands Research School for Astronomy (NOVA) Virtual Institute of Accretion (VIA) postdoctoral fellowships; NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation; Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017-00648); 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 and by the Province of Ontario through the Ministry of Research, Innovation and Science); the Portuguese Foundation for Science and Technology (FCT) grants (Individual CEEC program - 5th edition, https://doi.org/10.54499/UIDB/04106/2020, https://doi.org/10.54499/UIDP/04106/2020, PTDC/FIS-AST/3041/2020, CERN/FIS-PAR/0024/2021, 2022.04560.PTDC); the Princeton Gravity Initiative; the Spanish Ministerio de Ciencia e Innovación (grants PGC2018-098915-B-C21, AYA2016-80889-P, PID2019-108995GB-C21, PID2020-117404GB-C21); the University of Pretoria for financial aid in the provision of the new Cluster Server nodes and SuperMicro (USA) for a SEEDING GRANT approved toward these nodes in 2020; the Shanghai Municipality orientation program of basic research for international scientists (grant no. 22JC1410600); the Shanghai Pilot Program for Basic Research, Chinese Academy of Science, Shanghai Branch (JCYJ-SHFY-2021-013); the Simons Foundation (grant 00001470); the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017- 0709); the Spanish Ministry for Science and Innovation grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033; the Spinoza Prize SPI 78-409; the South African Research Chairs Initiative, through the South African Radio Astronomy Observatory (SARAO, grant ID 77948), which is a facility of the National Research Foundation (NRF), an agency of the Department of Science and Innovation (DSI) of South Africa; the Swedish Research Council (VR); the Taplin Fellowship; the Toray Science Foundation; the UK Science and Technology Facilities Council (grant no. ST/X508329/1); the US Department of Energy (USDOE) through the Los Alamos National Laboratory (operated by Triad National Security, LLC, for the National Nuclear Security Administration of the USDOE, contract 89233218CNA000001); and the YCAA Prize Postdoctoral Fellowship. This work was also supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (RS-2024-00449206).

We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00841.V and ADS/JAO.ALMA#2016.1.01154.V. ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council (Canada), Ministry of Science and Technology (MOST; Taiwan), Academia Sinica Institute of Astronomy and Astrophysics (ASIAA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc. (AUI)/NRAO, and the National Astronomical Observatory of Japan (NAOJ). The NRAO is a facility of the NSF operated under cooperative agreement by AUI. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under contract No. DE-AC05-00OR22725; the ASTROVIVES FEDER infrastructure, with project code IDIFEDER-2021-086; the computing cluster of Shanghai VLBI correlator supported by the Special Fund for Astronomy from the Ministry of Finance in China; We also thank the Center for Computational Astrophysics, National Astronomical Observatory of Japan. This work was supported by FAPESP (Fundacao de Amparo a Pesquisa do Estado de Sao Paulo) under grant 2021/01183-8.

APEX is a collaboration between the Max-Planck-Institut für Radioastronomie (Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key Research and Development Program (No. 2017YFA0402700) of China and Natural Science Foundation of China grant 11873028. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica, Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA). The IRAM 30-m telescope on Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (Max-Planck-Gesellschaft, Germany), and IGN (Instituto Geográfico Nacional, Spain). The SMT is operated by the Arizona Radio Observatory, a part of the Steward Observatory of the University of Arizona, with financial support of operations from the State of Arizona and financial support for instrumentation development from the NSF. Support for SPT participation in the EHT is provided by the National Science Foundation through award OPP-1852617 to the University of Chicago. Partial support is also provided by the Kavli Institute of Cosmological Physics at the University of Chicago. The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA.

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. This research is part of the Frontera computing project at the Texas Advanced Computing Center through the Frontera Large-Scale Community Partnerships allocation AST20023. Frontera is made possible by National Science Foundation award OAC-1818253. This research was done using services provided by the OSG Consortium (Pordes et al. 2007; Sfiligoi et al. 2009), which is supported by the National Science Foundation award Nos. 2030508 and 1836650. Additional work used ABACUS2.0, which is part of the eScience center at Southern Denmark University, and the Kultrun Astronomy Hybrid Cluster (projects Conicyt Programa de Astronomia Fondo Quimal QUIMAL170001, Conicyt PIA ACT172033, Fondecyt Iniciacion 11170268, Quimal 220002). Simulations were also performed on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, on the HazelHen cluster at the HLRS in Stuttgart, and on the Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University. The computer resources of the Finnish IT Center for Science (CSC) and the Finnish Computing Competence Infrastructure (FCCI) project are acknowledged. This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and Compute Canada (http://www.computecanada.ca).

The EHTC has received generous donations of FPGA chips from Xilinx Inc., under the Xilinx University Program. The EHTC has benefited from technology shared under open-source license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER). The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers. This research has made use of NASA’s Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for EHT-specific support with the use of DiFX. We thank Martin Shepherd for the addition of extra features in the Difmap software that were used for the CLEAN imaging results presented in this paper. We acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

Appendix B: Average images of exploratory models

For reference, we present the average images for the exploratory models (see the main text for detailed descriptions).

Figure B.1 displays the selected images of the tilted disk models. Note that the morphology of tilted disk images (the middle and right panels) are more crescent-like, compared to the reference, un-tiled images (left panel).

Figure B.2 displays the averaged images of the radiative models for two electron heating models (Howes and Rowan electron heating; see (Chael et al. 2019c) for more details). The two different types of electron heating make only a small difference in the averaged images.

Fig. B.1. Time-averaged images of untilted and tilted disk models computed using H-AMR pipeline, for azimuthal viewing angle 0° and Rhigh = (1, 160). The white arrows indicate the projected direction of black hole spin, which is pointing away from Earth. The images are shown in linear scale.

Fig. B.2. Time-averaged images of radiative models computed using koralpipelines. The images are shown in linear scale.

Appendix C: Bayesian inference for combining multi-year results

Here we present the details of the multi-epoch scoring procedure.

C.1. Overview: a Bayesian view

We begin by making the following approximations and identifications:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \text{Likelihood:}{L}_{Y_i}(F_{Y_i},\mathrm{PA},M)\equiv P(D_i|F_{Y_i},\mathrm{PA},M)\\ \hfill \text{Bayesian evidence:}{Z}_{Y_i}\equiv {\displaystyle \int P}(D_{Y_i}|F_{Y_i},\mathrm{PA},M)\\ \hfill \text{Posterior:}\wp (F_{Y_i},\mathrm{PA},M)\equiv P(F_{Y_i},\mathrm{PA},M|D_{Y_i}),\end{array}\end{array}$$(C.1)

where D represents observed data, and (F, PA, M) respectively represents the total flux, the position angle of the black hole spin, and the mass of the source. In the above description, the subscript Y_i denotes the value observed at different years.

Bayes theorem,

$$\begin{array}{c}\hfill P(A|B)=P(B|A)P(A)/P(B)=P(B|A)P(A)/{\displaystyle \int P}(B|A)P(A),\end{array}$$(C.2)

here implies that

$$\begin{array}{c}\hfill \wp (F_{Y_i},\mathrm{PA},M)=L_{Y_i}(F_{Y_i},\mathrm{PA},M)/Z_{Y_i}.\end{array}$$(C.3)

We assume that the total flux is independent from year to year, and that the PA and mass are fixed. For two observations at different years, keeping these definitions for three sets of data: that in Y₁, that in Y₂, and that for the combined data set Y₁₂, we note that if the data sets are independent as assumed throughout our analysis, then

$$\begin{array}{c}\hfill L_{Y_{12}}(F_{Y_1},F_{Y_2},\mathrm{PA},M)=L_{Y_1}(F_{Y_1},\mathrm{PA},M)L_{Y_2}(F_{Y_2},\mathrm{PA},M),\end{array}$$(C.4)

and therefore,

$$\begin{array}{c}\hfill \wp (F_{Y_1},F_{Y_2},\mathrm{PA},M)Z_{Y_{12}}=\wp (F_{Y_1},\mathrm{PA},M)\wp (F_{Y_2},\mathrm{PA},M)Z_{Y_1}{Z}_{Y_2}.\end{array}$$(C.5)

Finally, we note that by definition,

$$\begin{array}{c}\hfill {\displaystyle \int \wp }(F_{Y_1},F_{Y_2},\mathrm{PA},M)d\mathrm{PA}dM=1,\end{array}$$(C.6)

and therefore, we trivially have,

$$\begin{array}{c}\hfill Z_{Y_{12}}=Z_{Y_1}{Z}_{Y_2}{\displaystyle \int \wp }(F_{Y_1},\mathrm{PA},M)\wp (F_{Y_2},\mathrm{PA},M)d\mathrm{PA}dM.\end{array}$$(C.7)

This approach can be generalized for additional years of data by, e.g., adding more liklihood L_Yi and Bayesian evidence Z_Yi into the right hand sides of Eqs. (C.4) and (C.5).

C.2. Ensemble-based posterior construction

Within a given model, the posterior probability distribution on rotation (PA), and mass, is a natural output of any fitting procedure, such as the snapshot scoring applied in this paper. Typically, the posteriors of fitting individual snapshots to the EHT data do not overlap, with the intrinsic posterior widths being small in comparison to the distances between those from different snapshots. For the 2017 EHT observations we solved this with what we called “Ensemble-based Posterior Construction” (Appendix G of M87^* 2017 VI, here EBPC). This is an example of Approximate Bayesian Computation (ABC) methods.

ABC is a method of posterior approximation in which forward modeled realizations of a random process are compared against data to generate a posterior on a subset of the model parameters. It is especially useful when the construction of realizations is expensive and the nature of the random process is not of intrinsic interest. In our case, the random process is the turbulence, and the model parameters of interest are the flux, PA, and mass. The procedure is to define a cut on a data/model statistic, here chosen to be χ_ν², keep fits that pass the cut, and use the resulting distribution of the parameters of interest. Demonstration that the resulting distributions faithfully capture truth values is presented in Appendices G.2 and G.3 of M87^* 2017 VI, and used to generate parameter constraints in M87^* 2017 V.

To produce ABC Posteriors (or EBPC posteriors) for a single simulation, the procedure is:

• 1.

  Fit each snapshot to the data (N fits).

• 2.

  Continuous posteriors are generated from kernal density estimations (KDEs) of the sampled values after applying desired χ_ν² cut.

Validation requires additional iterations using synthetic data sets drawn from simulation images. This procedure was already validated on the 2017 data.

Posteriors from the full set of models can be constructed by summing over posteriors from each model after selecting a model quality factor, e.g., AIS.

C.3. AIS – average image scoring

Average image scoring (AIS) is a procedure for assessing the likelihood that the data is drawn from the distribution of random processes associated with a given GRMHD simulation. This method is a modest extension of the standard χ² test. What differs explicitly in this case is the probability distribution of the relevant fit statistic, i.e. the analog of the χ²-distribution, which is otherwise used to assess the test. This distinction is a consequence of the presence of stochastic, or “turbulent”, structures in the image that must be marginalized over. The primary complication is that these structures are neither sufficiently sampled in the GRMHD simulation image sets for a formal fitting process nor have a priori known statistics.

The solution is to assess a measure of the “distance” between the data and a fiducial image, here chosen to be the average image from a given simulation⁸ though such a choice is not unique, in light of the typical range of distances we expect based on the simulation itself. The distance measure we adopt is the reduced χ², χ_ν², obtained by fitting the average image to the observational data. The value for the data is then evaluated against the value obtained by fitting the average image to synthetic data sets generated from individual snapshot images. A modest optimization is found by fitting the snapshot images to data generated from the average image, which is exactly equivalent in the absence of telescope gains due to the quadratic nature of the χ_ν² statistic, and is equivalent in practice in the presence of systematic uncertainties.

In principle, a Kolmogorov-Smirnov (KS) test may be used to formally assess probability that the data χ_ν² is drawn from the distribution of values from the simulation itself. In the past (M87^* 2017 V) we have made the very conservative binary assessment of “inside” or “outside” the distribution, which for one value is equivalent. The outcome of AIS is a probability that the measured χ_ν² is drawn from the distribution.

C.4. Multi-epoch scoring

C.4.1. Statistical considerations

In principle, for each year Y and each simulation S, we will have a Bayesian evidence Z_SY resulting from some overall fit quality assessment and a normalized posterior ℘_SY(F, PA, M). The likelihood of a given parameter set is then given by

$$\begin{array}{c}\hfill L_{\mathit{SY}}(F,\mathrm{PA},M)={\wp }_{\mathit{SY}}(F,\mathrm{PA},M)Z_{\mathit{SY}},\end{array}$$(C.8)

from which we identify

$$\begin{array}{c}\hfill {\wp }_{\mathit{SY}}(F,\mathrm{PA},M)=\frac{L_{\mathit{SY}}(F,\mathrm{PA},M)}{Z_{\mathit{SY}}}\end{array}$$(C.9)

and

$$\begin{array}{c}\hfill Z_{\mathit{SY}}={\displaystyle \int d}Fd\mathrm{PA}dM{L}_{\mathit{SY}}(F,\mathrm{PA},M).\end{array}$$(C.10)

The odds ratios between two simulations, S and S′, from a given year is simply

$$\begin{array}{c}\hfill O_{S{S}^{\prime },Y}=\frac{{\displaystyle \int L_{\mathit{SY}}}(F,\mathrm{PA},M)}{{\displaystyle \int L_{S^{\prime }Y}}(F,\mathrm{PA},M)}=\frac{Z_{\mathit{SY}}}{Z_{S^{\prime }Y}}.\end{array}$$(C.11)

The full parameter posterior in a given year is by Bayes theorem directly,

$$\begin{array}{c}\hfill {\wp }_Y(F,\mathrm{PA},M)={\displaystyle \sum _S}{\wp }_{\mathit{SY}}(F,\mathrm{PA},M)Z_{\mathit{SY}}/{\displaystyle \sum _S}{Z}_{\mathit{SY}},\end{array}$$(C.12)

where now the simulation-specific prior is the Z_SY.

Because the measurements and the turbulent realizations are expected to be uncorrelated across years, combining multiple years is straightforward in principle. For a given simulation, the combined likelihood is

$$\begin{array}{c}\hfill L_S(F_{Y_1},F_{Y_2},\mathrm{PA},M)=L_{S{Y}_1}(F_{Y_1},\mathrm{PA},M)L_{S{Y}_2}(F_{Y_2},\mathrm{PA},M).\end{array}$$(C.13)

This must be supplemented with (possibly joint) priors on F_Y1, F_Y2 (hereafter the joint prior is denoted as π). Under the assumption that total flux is independent from year to year, and that the PA is fixed (as is mass), π(F_Y1, F_Y2)=const and π(PA_Y1, PA_Y2)=δ(PA_Y1 − PA_Y2). We may then relate the posteriors and Bayesian evidence from the independent analyses on each year to those from the combined data set:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Z_S& ={\displaystyle \int d}{F}_{Y_1}d{F}_{Y_2}d\mathrm{PA}dM{L}_S(F_{Y_1},F_{Y_2},\mathrm{PA},M)\hfill \\ \hfill & =Z_{S{Y}_1}{Z}_{S{Y}_2}\times \hfill \\ \hfill & {\displaystyle \int d}{F}_{Y_1}d{F}_{Y_2}d\mathrm{PA}dM{\wp }_{S{Y}_1}(F_{Y_1},\mathrm{PA},M){\wp }_{S{Y}_{Y_2}}(F_2,\mathrm{PA},M).\hfill \end{array}\end{array}$$(C.14)

and, again, from Bayes theorem directly,

$$\begin{array}{c}\hfill \wp (F_{Y_1},F_{Y_2},\mathrm{PA},M)=\frac{Z_{S{Y}_1}{Z}_{S{Y}_2}}{Z_S}{\wp }_{S{Y}_1}(F_{Y_1},\mathrm{PA},M){\wp }_{S{Y}_2}(F_{Y_2},\mathrm{PA},M).\end{array}$$(C.15)

With the Z_S for each model, we can construct the odds ratios immediately for model selection, and with the ℘_S(F_Y1, F_Y2, PA, M) we can assess parameter constraints.

To implement the foregoing, we require estimates of the following functions:

• Z_SYi – The Bayesian evidence for a simulation in a given year Y_i.

• ℘_SYi}(F_Yi, PA,M) – The posterior probability distribution of the snapshot fit parameters in a given year Y_i.

• π(F_Y1,F_Y2) – The joint prior on the fluxes from year Y₁ to year Y₂, which we assume is constant, but really must merely be independent.

Now we turn to evaluation of these quantities.

C.4.2. Approximate integration

Given the natural expectations for the joint priors on flux, PA, and mass already described, we can combine the constraints from multiple years. Adopting AIS as the Z_SY and the EBPC posterior as the ℘_SY, the integrals in Eqs. C.14 and C.15 may be numerically performed explicitly.

Explicitly, upon choosing a flux prior that is completely uninformative, i.e., π(F_Y1, F_Y2)=1, and a fixed PA, i.e., π(PA_Y1, PA_Y2)=δ(PA_Y1 − PA_Y2), the updated Bayesian evidence for a single simulation obtained from combining the two years is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Z_S& =Z_{S{Y}_1}{Z}_{S{Y}_2}\times \hfill \\ \hfill & {\displaystyle \int d}{F}_{Y_1}d{F}_{Y_2}d\mathrm{PA}dM{\wp }_{S{Y}_1}(F_{Y_1},\mathrm{PA},M){\wp }_{S{Y}_2}(F_{Y_2},\mathrm{PA},M),\hfill \end{array}\end{array}$$(C.16)

and the posterior on PA and M is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\wp }_S(\mathrm{PA},M)& =\frac{Z_{S{Y}_1}{Z}_{S{Y}_2}}{Z_S}\times \hfill \\ \hfill & {\displaystyle \int d}{F}_{Y_1}d{F}_{Y_2}{\wp }_{S{Y}_1}(F_{Y_1},\mathrm{PA},M){\wp }_{S{Y}_2}(F_{Y_2},\mathrm{PA},M).\hfill \end{array}\end{array}$$(C.17)

Combining multiple simulations gives

$$\begin{array}{c}\hfill \wp (\mathrm{PA},M)={\displaystyle \sum _S}{\wp }_S(\mathrm{PA},M)Z_S/{\displaystyle \sum _S}{Z}_S.\end{array}$$(C.18)

This is attractive for a number of reasons. First, it is done as soon as scoring is finished. Second, the Bayesian odds ratios can place significantly improved cuts on the GRMHD library. For example, in our analysis, we compared normalized odd ratio for different models, and applied an arbitrary cut (see Fig. 7). Third, it has the virtue of being formally sensible.

The AIS score is not, however, exactly a Bayesian evidence, despite being used as such. In the absence of a full description of the turbulence and the capacity to effectively parameterize it and/or sample it, we do not have the capacity to generate true Bayesian evidences. However, generally, the AIS score is a more permissive, largely because we conservatively treat poorer-than-expected fits as acceptable. Nevertheless, there remains a possibility that a much better fit exists associated with a turbulent realization that was simply never sampled by the simulation. Similarly, the EBPC is not, exactly, a true posterior; it is typically broader. Again, this is unlikely to reject the true model since it serves to contaminate the acceptable parameter sets with formally bad fits.

Appendix D: AIS result summary

In Table D.1, we summarize the 2017 and 2018 AIS results for the fiducial models.

Appendix E: M/D distribution for exploratory, radiative models

Figure E.1 displays the posterior distributions of M/D for the radiative models fitted to 2017 (blue) and 2018 (red) datasets. The range of M/D is similar to some fiducial models shown in Fig. 4. As seen in most of the fiducuial models, snapshots fit to the 2018 data produce a distribution with a slightly higher peak than the 2017 data.

Fig. E.1. Violin plots of the M/D distributions from the snapshot scoring method against the 2017 and 2018 EHT observations, for the exploratory radiative models.

All Tables

All Figures

	Fig. 1. Representative reconstructed images and model images of M87∗. Left panels: EHT images of M87∗ from the 2018 (M87* 2018 I; Paper I of this series) and 2017 (M87* 2017 IV) campaigns. Middle panels: Example GRMHD model images drawn from the same model at two different times. Right panels: Theoretical images convolved with a 20 μas FWHM Gaussian beam.
In the text

Fig. 2. Poloidal (azimuthally averaged) and equatorial plane cuts of selected KHARMA GRMHD simulation snapshots of the two considered accretion modes, both with a* = +0.5. The MAD (left panels) and SANE (right panels) snapshots feature different ion temperatures (displayed in dimensionless units in the upper panels) and the plasma β parameter (displayed in the bottom panels), resulting in different emission properties. The solid black lines represent the σ = 1 surface. The high temperature near the polar axes is a consequence of conservative codes. The temperatures in the polar regions are considered unphysical and are masked out in the radiative transfer calculations (regions with σ >  1 or σ >  3). Here, rg ≡ GM/c2, the gravitational radius.

In the text

	Fig. 3. Time-averaged images of all the fiducial models with Rlow = 1 and Rhigh = 40 in the image library. The white arrows indicate the projected direction of the black hole spin for the simulation geometry. The lower left corner is blank because there is no H-AMR SANE a* = −0.94 case in our model library. All images are shown in linear scale.
In the text

	Fig. 4. Violin plots for M/D distributions from the snapshot scoring with 2017 and 2018 EHT observations for models with Rlow = 1 and Rhigh = 40 in the image library. The horizontal dashed black line and shaded gray region mark the ranges of M/D = 3.8 ± 0.4 μas reported in M87* 2017 I.
In the text

Fig. 5. Normalized distributions of the fitted position angle of the forward jet in selected fiducial thermal models based on the snapshot scoring with 2017 and 2018 observations, using the best-fit 10% of the images. The PA = 288° ±10° of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023) is marked by the vertical dashed black line and shaded gray region. Left: Models with a black hole spin axis or a coordinate system axis (case a* = 0) pointing toward Earth. Prograde cases are shown in the top (a* >  0) and middle (a* = 0) panels, and the retrograde case (a* <  0) is shown in the bottom panel. Right: Models with a black hole spin axis or with a coordinate system axis (case a* = 0) pointing away from Earth. Prograde cases are shown in the top (a* >  0) and middle (for a* = 0 the matter angular momentum is always aligned with the coordinate axis) panels, and the retrograde case (a* <  0) is shown in the bottom panel. See text in this paper or Fig. 5 in M87* 2017 V for a description of the geometry of the black hole axis, flow rotational axis, and inclination angle.

In the text

Fig. 6. Demonstration of the multi-epoch scoring based on the KHARMA SANE model with (a*, Rlow, Rhigh)=(−0.5, 10, 160). The 2017 and 2018 snapshot scoring results for the given model are shown in the M/D – PA plane with the blue and red distributions, respectively. The larger the overlap between the red and blue distributions, the better the fit of the multi-epoch observations with a single model. The black profile indicates the normalized combined distribution of the best-fit parameters computed according to the multi-epoch scoring procedure (see Sect. 4.4 for details). The dashed lines and shaded areas are the same as explained in Figs. 4 and 5.

In the text

Fig. 7. Overview of the multi-epoch scoring results for all the fiducial models with black hole axes pointing away from Earth, and nonspinning black holes with a inclination angle i = 163°. The y-axis is the normalized odds ratio (see text for details), which indicates the relative preference by multi-epoch (2017 and 2018) observations. Models with different spins are shown in different panels. The green region marks passing models for which the relative odds ratio is ≥1%. Models with a relative odds ratio ≤0.001% are not shown.

In the text

Fig. 8. Time-averaged SEDs of selected best-bet models. The data points are taken from the multi-wavelength observations during the 2018 EHT campaign (EHT MWL Science Working Group 2024). The EHT observation is marked by the vertical red bar. The gray region of the SED indicates the variations of the SED for different snapshots. The colored histograms correspond to different radiative processes: synchrotron emission (blue), synchrotron photons scattered once (orange), synchrotron photons scattered twice (green), and bremsstrahlung (red). The total emission is displayed in black.

In the text

Fig. 9. Distribution of M/D obtained by fitting the 2017 (blue) and 2018 (red) data with the fiducial models shown in Table 3. These distributions are generated from the top 10% best-fitting snapshots by χ2. The results for models prepared by the KHARMA, H-AMR, and BHACpipelines are shown separately, together with the result for all the models combined. The distribution when combining both the 2017 and 2018 result by the multi-epoch scoring algorithm is shown as the black profile. The vertical dashed black line and shaded gray region mark the ranges of M/D = 3.8 ± 0.4 μas reported in Event Horizon Telescope Collaboration (2019a).

In the text

	Fig. B.1. Time-averaged images of untilted and tilted disk models computed using H-AMR pipeline, for azimuthal viewing angle 0° and Rhigh = (1, 160). The white arrows indicate the projected direction of black hole spin, which is pointing away from Earth. The images are shown in linear scale.
In the text

	Fig. B.2. Time-averaged images of radiative models computed using koralpipelines. The images are shown in linear scale.
In the text

	Fig. E.1. Violin plots of the M/D distributions from the snapshot scoring method against the 2017 and 2018 EHT observations, for the exploratory radiative models.
In the text