First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole

We present measurements of the properties of the central radio source in M87 using Event Horizon Telescope data obtained during the 2017 campaign. We develop and fit geometric crescent models (asymmetric rings with interior brightness depressions) using two independent sampling algorithms that consider distinct representations of the visibility data. We show that the crescent family of models is statistically preferred over other comparably complex geometric models that we explore. We calibrate the geometric model parameters using general relativistic magnetohydrodynamic (GRMHD) models of the emission region and estimate physical properties of the source. We further fit images generated from GRMHD models directly to the data. We compare the derived emission region and black hole parameters from these analyses with those recovered from reconstructed images. There is a remarkable consistency among all methods and data sets. We find that>50% of the total flux at arcsecond scales comes from near the horizon, and that the emission is dramatically suppressed interior to this region by a factor>10, providing direct evidence of the predicted shadow of a black hole. Across all methods, we measure a crescent diameter of 42+/-3 micro-as and constrain its fractional width to be<0.5. Associating the crescent feature with the emission surrounding the black hole shadow, we infer an angular gravitational radius of GM/Dc2 = 3.8+/- 0.4 micro-as. Folding in a distance measurement of 16.8(+0.8,-0.7) Mpc gives a black hole mass of M = 6.5 +/- 0.2(stat) +/-0.7(sys) 10^9 Msun. This measurement from lensed emission near the event horizon is consistent with the presence of a central Kerr black hole, as predicted by the general theory of relativity.


Introduction
Einstein's general theory of relativity not only predicts the existence of black holes, but also provides a means to directly observe them.Photons can escape from near the event horizon via an unstable circular orbit (von Laue 1921;Bardeen 1973), whose observational manifestation would be a bright ring of emission surrounding a dark interior black hole "shadow" (Luminet 1979;Falcke et al. 2000).The diameter of the shadow for a black hole of mass M as seen by a distant observer is predicted to be ;9.6-10.4GM/c 2 , which is larger than twice the coordinate radius of the event horizon due to light-bending effects (Takahashi 2004;Johannsen & Psaltis 2010).The range results from different values of black hole spin and observer inclination angle.The black hole shadow can only be seen if (i) there are a sufficient number of emitted photons to illuminate the black hole, (ii) the emission comes from close enough to the black hole to be gravitationally lensed around it, and (iii) the surrounding plasma is sufficiently transparent at the observed wavelength.For nearby low-luminosity black holes accreting via a radiatively inefficient flow, these conditions can be met at millimeter wavelengths(e.g., Özel et al. 2000;Ho 2008;Yuan & Narayan 2014).
The nearby massive elliptical galaxy M87 provides an ideal laboratory to search for such a black hole shadow.It is relatively nearby (   D 16.8 Mpc) and has long been known to host a bright, compact radio source at its center (Cohen et al. 1969).Starting with Young et al. (1978) and Sargent et al. (1978), several attempts have been made to "weigh" the supermassive black hole (SMBH) hypothesized to power the radio source.Recent stellardynamics observations by Gebhardt et al. (2011) found M= (6.6 ± 0.4)×10 9 M e , while the latest gas dynamics observations by Walsh et al. (2013)  Both values assumed a distance D=17.9Mpc.Strictly interpreted within the predictions of the general theory of relativity, these measurements make a strong case that M87 does harbor some sort of compact massive dark object at its center, but they have insufficient resolution to formally demonstrate that it is indeed an SMBH.For M87 the expected shadow angular diameter is ; 20 or ; 38μ as, which is now accessible using global very long baseline interferometry (VLBI) at millimeter wavelengths (EHT Collaboration et al. 2019a, hereafter Paper II).
Accreting black holes are powered by matter flowing in via an accretion disk that in many cases launches a powerful jet of magnetized, relativistic plasma (e.g., Blandford & Znajek 1977;Blandford & Payne 1982).M87 exhibits the characteristic flat-toinverted radio/millimeter (mm) synchrotron spectrum considered to be a hallmark of the compact innermost jet core in lowluminosity active galactic nuclei (AGNs; e.g., Blandford & Königl 1979;Ho 1999).In this picture, the jet photosphere moves inward with increasing frequency up to the spectral break, at which point the entire jet becomes optically thin.The average broadband spectrum of M87 (Reynolds et al. 1996;Di Matteo et al. 2003;Prieto et al. 2016) indicates that the mm-band should straddle this transition.While images at longer radio wavelengths reveal extended jet structure (e.g., Asada & Nakamura 2012;Hada et al. 2016;Mertens et al. 2016;Kim et al. 2018;Walker et al. 2018), both the observed core shift (Hada et al. 2011) and compact size of ;40 μas from past mm-VLBI (Doeleman et al. 2012;Akiyama et al. 2015) are consistent with an origin of the mm-band emission near the event horizon of the central black hole.
Although the photon ring and shadow predictions are clear, the image morphology will depend on the physical origin of the surrounding emission and spacetime of the black hole.If the observed synchrotron radiation at 1.3 mm originates far from the black hole, the forward jet will dominate the observed emission and the lensed emission and shadow feature should be weak (Broderick & Loeb 2009).If instead the emission comes from near the event horizon, either the counter-jet or the accretion flow can produce a compact ring-or crescent-like image surrounding the shadow (Dexter et al. 2012).This type of image is now known to be commonly produced in radiative models of M87 based on general relativistic magnetohydrodynamic (GRMHD) simulations (Dexter et al. 2012;Mościbrodzka et al. 2016;Ryan et al. 2018;Chael et al. 2019b; see EHT Collaboration et al. 2019d, hereafter Paper V).
The outline of the shadow is expected to be nearly circular, if the central object in M87 is a black hole described by the Kerr metric (Bardeen 1973;Takahashi 2004).Violations of the nohair theorem generically change the shadow shape and size (e.g., Bambi & Freese 2009;Johannsen & Psaltis 2010;Falcke & Markoff 2013;Broderick et al. 2014;Cunha & Herdeiro 2018).Therefore, detecting a shadow and extracting its characteristic properties, such as size and degree of asymmetry, offers a chance to constrain the spacetime metric.
The high resolution necessary to resolve horizon scales for M87 has been achieved for the first time by the Event Horizon Telescope (EHT) in April 2017, with an array that spanned eight stations in six sites across the globe (Paper II).The EHT observed M87 on four days (April 5, 6, 10, and 11) at 1.3 mm (EHT Collaboration et al. 2019b, hereafter Paper III), and imaging techniques applied to this data set reveal the presence of an asymmetric ring structure (EHT Collaboration et al. 2019c, hereafter Paper IV).A large library of model images generated from GRMHD simulations generically finds such features to arise from emission produced near the black hole (Paper V) which is strongly lensed around the shadow.
In this Letter we use three different methods to measure properties of the M87 230 GHz emission region using EHT 2017 observations.In Section 2, we describe the EHT data set used.We present in Section 3 a pedagogical description showing how within compact ring models the emission diameter and central flux depression (shadow) can be inferred directly from salient features of the visibility data.In Section 4 we describe the three analysis codes used to infer parameters from the data, and in Section 5 and Section 6 we fit both geometric and GRMHD-based models.In Section 7 we extract properties of the reconstructed images from Paper IV.
We show that asymmetric ring ("crescent") geometric source models with a substantial central brightness depression provide a better statistical description of the data than other comparably complex models (e.g., double Gaussians).We use Bayesian inference techniques to constrain the size and width of this crescent feature on the sky, as well as the brightness contrast of the depression at its center compared to the rim.We show that all measurements support a source structure dominated by lensed emission surrounding the black hole shadow.
To extract the physical scale of the black hole at the distance of M87, GM/Dc 2 , from the observed ring structure in geometric models and image reconstructions, we do not simply assume that the measured emission diameter is that of the photon ring itself.
We instead directly calibrate to the emission diameter found in model images from GRMHD simulations.The structure and extent of the emission preferentially from outside the photon ring leads to a 10% offset between the measured emission diameter in the model images and the size of the photon ring.The scatter over a large number of images, which constitutes a systematic uncertainty, is found to be of the same magnitude.
We use independent calibration factors obtained for the geometric models and reconstructed images (Paper IV), providing two estimates of GM/Dc 2 .We also fit the library of GRMHD images described in Paper V directly to the EHT data, which provides a third.All three methods are found to be in remarkable agreement.We consider prior dynamical measurements of M/D and D for M87 in Section 8.In Sections 9 and 10, we discuss the evidence for the detection of lensed emission surrounding the shadow of a black hole in EHT 2017 data.We use the prior distance information to convert the physical scale to a black hole mass and show that our result is consistent with prior stellar, but not gas, dynamical measurements.We further discuss the implications for the presence of an event horizon in the central object of M87.Further technical detail supporting the analyses presented here has been included as Appendices.

Observations and Data
Operating as an interferometer, the EHT measures complex visibilities on a variety of baselines b ij between stations i and j.A complex visibility is a Fourier component of the source brightness distribution I(x, y), u v e I x y dxdy , , , 1 i ux vy 2 where (x, y) are angular coordinates on the sky, and (u, v) are projected baseline coordinates measured in units of wavelengths (Thompson et al. 2017;hereafter TMS).
The 2017 EHT observations of M87 and their subsequent correlation, calibration, and validation are described in detail in Paper III.On each of the four days-April 5, 6, 10, and 11the EHT observed M87 in two 2 GHz frequency bands centered on 227.1 GHz (low-band; LO) and 229.1 GHz (high-band; HI); the baseline coverage for April 11 is shown in Figure 1.For the modeling results presented in this Letter we analyze all four observing days and both bands.We use Stokes I visibility data reduced via the EHT-Haystack Observatory Processing System (HOPS) pipeline (Blackburn et al. 2019), coherently averaged in time by scan.Scan averaging decreases the data volume with negligible coherence losses (see Paper III) and it further serves to increase the signal-to-noise ratio (S/N) of the data.Fewer timestamps correspond to fewer gain terms (see Section 2.1) and higher S/N improves the validity of our Gaussian likelihood functions (see Section 4.1).

Data Products
Visibility measurements are affected by a combination of thermal noise and systematic errors.The thermal noise ò ij is distributed as a zero-mean complex Gaussian random variable with variance determined by the radiometer equation (TMS), while the dominant systematic noise components are associated with station-based complex gains g i .The measured visibilities can thus be expressed as are the measured visibility amplitude and phase.Measured visibility amplitudes are biased upward by thermal noise, so we use A ij to denote debiased visibility amplitude measurements (see TMS).
All noise sources in Equation ( 2) are functions of time and frequency, but the gain phase variations are particularly important for EHT data.Characteristic atmospheric timescales at 230GHz are on the order of seconds, rendering visibility phase calibration unfeasible (Paper II; Paper III).We instead recover source phase information via the construction of closure phases ψ C , given by the argument of a product of visibilities around a triangle of baselines, Because each gain term in the triple product gets multiplied by its complex conjugate, closure phases are immune to gain phase corruptions (Rogers et al. 1974).
Visibility amplitudes suffer less severely than visibility phases from station gain noise but the use of gain-free amplitude quantities can still aid modeling efforts.Closure amplitudes A C are constructed from four visibilities on a quadrangle of baselines,

= ( )
The appearance of each station in both the numerator and denominator of this expression causes the station gain amplitudes to cancel out.Because the closure amplitude is constructed from products and ratios of visibility amplitudes, it is often convenient to work instead with the logarithm of the closure amplitude,

Data Selection and Preparation
From the scan-averaged visibility data, we increase the uncertainty associated with the debiased visibility amplitudes, A, by adding a 1% systematic uncertainty component in quadrature to the thermal noise (Paper III; Paper IV); we refer to this increased uncertainty as the "observational error."These debiased visibility amplitudes are then used to construct a set of logarithmic closure amplitudes, A ln C , per Equation (5).Closure amplitude measurements are generally not independent because a pair of quadrangles may have up to two baselines in common, and a choice must be made regarding which minimal (or "non-redundant") subset of closure amplitudes to use.We select the elements of our minimal set by starting with a maximal (i.e., redundant) set, from which we systematically remove the lowest-S/N quadrangles until the size of the reduced set is equal to the rank of the covariance matrix of the full set (see L. Blackburn et al. 2019, in preparation).This construction procedure serves to maximize the final S/N of the resulting closure amplitudes.
We construct closure phases, ψ C , from the visibilities using Equation (3), after first removing visibilities on the short intrasite baselines (James Clerk Maxwell Telescope-Submillimeter Array (JCMT-SMA), Atacama Large Millimeter/submillimeter Array-Atacama Pathfinder Experiment (ALMA-APEX)) which produce only "trivial" closure phases ;0°(Paper III; Paper IV).As with closure amplitudes, closure phase measurements are in general not independent of one another because a pair of triangles may share a baseline.However, a suitable choice of non-redundant closure phase subset can minimize the covariance between measurements.We select our subset such that the highest-S/N baselines are the most frequently shared across triangles.Such a subset can be obtained by selecting one station to be the reference and then choosing all triangles containing that station (TMS).Because ALMA is so much more sensitive than the other stations in the EHT 2017 array, this construction procedure using ALMA as a reference station ensures the near-diagonality of the closure phase covariance matrix (see Section 4.1).
We list the number of data products in each class, along with the number of station gains, in Table 1 for each observing day and band.Information about accessing SR1 data and the The (u, v)-coverage has two primary orientations, east-west in blue and north-south in red, with two diagonal fillers at large baselines in green and black.Note that the Large Millimeter Telescope (LMT) and the Submillimeter Telescope (SMT) participate in both orientations, and that the LMT amplitudes are subject to significant gain errors.There is evidence for substantial depressions in the visibility amplitudes at ∼3.4 Gλ and ∼8.3 Gλ.The various lines in the right panel show the expected behavior of (dotted line) a Gaussian, (dashed line) a filled disk, and (green area) a crescent shape along different orientations.The image of M87 does not appear to be consistent with a filled disk or a Gaussian.software used for analysis can be found on the Event Horizon Telescope website's data portal.107

Descriptive Features of the Visibility Data
In Paper IV, different image reconstruction methods all obtained similar looking images of M87 from the 2017 EHT observations, namely, a nearly circular ring with a dark center and azimuthally varying intensity.In this Letter, we consider a range of source models and calculate the corresponding visibilities as a function of the model parameters.We then employ statistical tools to select between models and to estimate model parameters by comparing the model visibilities to the observed ones.Because imaging and visibility domain analysis rely on the entire complex visibility data set and/or closure quantities, it is useful as a pedagogical guide to show first how some of the simple image characteristics are imprinted on the visibility data.We emphasize that a complete and accurate description of the source requires imaging analysis and visibility modeling, which we perform in Paper IV and in later sections of this Letter.
Here we show that the data are consistent with the presence of a ring structure with a characteristic emission diameter of ∼45 μas.These aspects also match the predictions for an image dominated by lensed emission near the photon ring surrounding the black hole shadow of M87 (Paper V).
The 2017 EHT observations of M87 have good (u, v)coverage, primarily along an east-west (blue) and a northsouth (red) orientation, with additional diagonal long baselines (green and black; see Figure 1  There is evidence for a minimum of the visibility amplitudes at baseline lengths of ∼3.4 Gλ, followed by a second peak around ∼6 Gλ.Such minima are often associated with edges or gaps in the image domain. Further, the visibility amplitudes are similar in the north-south and east-west directions, suggestive of a similar characteristic image size and shape in both directions (Figure 1).This is naturally accomplished if the image possesses a large degree of azimuthal symmetry, such as in a ring or disk.Differences in the visibilities as a function of baseline length for different orientations do exist, however, particularly in the depth of the first null, indicating that the source is not perfectly symmetric.
Next, we consider the presence of the central flux depression.For the case of a uniform disk model, the second visibility amplitude minimum occurs at 1.8 times the location of the first minimum, i.e., at ∼6.3-7Gλ, which is not seen in the visibility amplitudes.For a ring or annular model, however, the second minimum moves to longer baseline lengths, consistent with what is seen in the visibility amplitudes.
Indeed, the Fourier transform of an infinitesimally thin ring structure shows the first minimum in visibility amplitude at a baseline length b 1 for which the zeroth-order Bessel function is zero (see TMS).This allows us to estimate the source size for a ring model as In subsequent sections, we quantify our characterization of this model through fitting in the visibility and image domains.

Model Fitting to Interferometric Data
We utilize three independent algorithms for parameter space exploration to quantify the size, shape, and orientation of this asymmetric ring structure.We fit both geometric and GRMHD models to the 2017 EHT interferometric data.In this section we outline the modeling framework used to extract parameter values from the M87 data.We first detail the construction of the corresponding likelihood functions in Section 4.1, and we then describe in Section 4.2 the three different codes we have used to estimate model parameters.

Likelihood Construction
Our quantitative modeling approach seeks to estimate the posterior distribution Q ( | ) D P of some parameters Q within the context of a model and conditioned on some data D, is the prior probability of the model parameters, and is the Bayesian evidence.In this section we define our likelihood functions  Q ( ), which we note differ in detail from those adopted for the regularized maximum likelihood (RML) imaging procedures presented in Paper IV.
For each scan, the measured visibility amplitude, A, corresponds to the magnitude of a random variable distributed according to a symmetric bivariate normal distribution (TMS).This magnitude follows a Rice distribution, which in the Note.N A is the number of visibility amplitudes, N g is the number of gain terms, y N C is the number of closure phases, and N A ln C is the number of logarithmic closure amplitudes.We show counts for the visibility amplitudes both with and without the inclusion of short intra-site baselines (ALMA-APEX and JCMT-SMA); the visibility amplitudes including intra-site baselines are used in Section 5, while those without are used in Section 6.The closure phase count always excludes triangles containing intra-site baselines, while the logarithmic closure amplitude count always includes quadrangles containing intra-site baselines.Both closure phase and logarithmic closure amplitude counts are for minimal (nonredundant) sets.high-S/N limit reduces to a Gaussian near the mode, Here, s ij 2 is the variance of the visibility measurement, Âij is the model visibility amplitude of the source, and | | g i and | | g j are the gain amplitudes for stations i and j.Both the mean and standard deviation of Equation (9) are biased with respect to the true visibility amplitude distribution, but for A/σ2.0 these biases are below 10%; at least 94% of our visibility amplitude data for any day and band meet this criterion.
For scan-averaged EHT 2017 data, the gain amplitudes constitute on the order of 100 additional nuisance parameters per data set (see Table 1).These numerous additional parameters may be efficiently addressed by directly marginalizing the likelihood in Equation (9), a procedure detailed in Appendix A and A. E. Broderick et al. (2019, in preparation).Once the gain amplitudes have been reconstructed, the joint likelihood function for all visibility amplitude measurements within a data set is then given by the product over the individual likelihoods, where this product is taken over all baselines and scans.
The logarithm of the visibility amplitudes also follows a Gaussian distribution in the high-S/N limit, with an effective logarithmic uncertainty of The Gaussianity of the logarithmic visibility amplitudes implies that the logarithm of the closure amplitudes will similarly be Gaussian distributed in the same limit, with variances given by 12 This Gaussian approximation for logarithmic closure amplitudes holds well (i.e., the mean and standard deviation are biased by less than 10%) for  s 2.0; A ln C at least 87% of our logarithmic closure amplitude data for any day and band meet this criterion.
The likelihood function for a set of logarithmic closure amplitudes also depends on the covariances between individual measurements, which in general are not independent.We can construct a covariance matrix S A that captures the combined likelihood via where A is an ordered list of logarithmic closure amplitude residuals, and q is the number of non-redundant closure amplitudes; for a fully connected array with N el elements, q=N el (N el − 3)/2 (TMS).The covariance between logarithmic closure amplitude measurements A ln C,1234 and A ln C,1235 is (Lannes 1990a; L. Blackburn et al. 2019, in preparation) , ln 14 C,1235 ln ,12 2 ln ,13 2 in the Gaussian limit; here, we have used the fact that σ ij =σ ji to simplify notation.
The distribution of measured visibility phases, f, corresponds to the projection of a symmetric bivariate normal random variable onto the unit circle, which once again reduces to a Gaussian distribution in the high-S/N limit.Closure phases, ψ C , in this limit will also be Gaussian distributed with variances given by where s f ij , 2 is the variance in the visibility phase measurement, f ij .The Gaussian approximation for closure phases is unbiased in the mean with respect to the true closure phase distribution, and the standard deviation is biased by less than 10% for  s y 1.5; C this criterion is satisfied for at least 92% of our closure phase data on any day and band.
Closure phase measurements are also generally covariant, and for a covariance matrix S y the joint likelihood is given by where y is an ordered list of closure phase residuals, and t is the number of non-redundant closure phases; for a fully connected array containing el el (TMS).The covariance between two closure phase measurements, y C,123 and y C,124 , is given by (Kulkarni 1989;Lannes 1990b;L. Blackburn et al. 2019, in preparation) . 17 where y ˆijk C, is the modeled closure phase.The joint likelihood is then where the product is taken over all closure phases in the selected minimal subset.Because closure phases wrap around the unit circle, we always select the branch of y y -Ĉ C such that the difference lies between −180°and 180°.

Parameter Space Exploration Techniques
We utilize three independent algorithms for parameter space exploration.For the geometric crescent model fitting presented in Section 5, we use both Markov chain Monte Carlo (MCMC) and nested sampling (NS) algorithms.The MCMC modeling scheme explores model fits to the visibility amplitude and closure phase data, while the NS scheme fits to the closure phase and logarithmic closure amplitude data.For the GRMHD model fitting in Section 6, we use both MCMC and a genetic algorithm to fit the visibility amplitude and closure phase data.
THEMIS is written in C++ and parallelized via the Message Passing Interface (MPI) standard.THEMIS implements a differential evolution MCMC algorithm, and it utilizes parallel tempering based on the algorithm described in Nelson et al. (2014) and Braak (2006).In particular, THEMIS uses the adaptive temperature ladder prescription from Vousden et al. (2016).All sampling techniques, validation tests, and implementation details for THEMIS are described in detail in A. E. Broderick et al. (2019, in preparation).
In Sections 5 and 6, we use THEMIS to model the visibility amplitude and closure phase data, using the corresponding likelihood functions given in Equations (10) and (19), respectively.Gain amplitude terms are incorporated as model parameters (see Equation (9)) and are marginalized as described in Appendix A.

dynesty
In Section 5, we also use an NS technique, developed by Skilling (2006) primarily to evaluate Bayesian evidence integrals.We use the Python code dynesty (Speagle & Barbary 2018) as a sampler for the NS analyses presented in this Letter.The NS algorithm estimates the Bayesian evidence,  , by replacing the multidimensional integral over Q (see Equation ( 8)) with a 1D integral over the prior mass contained within nested isolikelihood contours.We permit dynesty to run until it estimates that less than 1% of the evidence is left unaccounted for.
Our NS analyses employ a likelihood function constructed exclusively from closure quantities; we account for data covariances in the likelihood function using Equation (16) for closure phases and Equation (13) for logarithmic closure amplitudes.Additionally, the use of logarithmic closure amplitudes in our NS fits removes information about the total flux density.

GENA
In addition to the above sampling algorithms, which seek to reconstruct a posterior distribution, we also employ an optimization procedure for comparing GRMHD simulations to data in Section 6.The optimization code, GENA (Fromm et al. 2019), is a genetic algorithm written in Python and parallelized using MPI.GENA minimizes a χ 2 statistic on visibility amplitudes and closure phases, using the gain calibration procedures in the eht-imaging (Chael et al. 2016(Chael et al. , 2018(Chael et al. , 2019a) ) Python package to solve for the gain amplitudes.GENA implements the Non-dominated Sorting Genetic Algorithm II (NSGA-II; Deb et al. 2002) for parameter exploration and the differential evolution algorithm from Storn & Price (1997) for constrained optimization.

Geometric Modeling
As detailed in Paper IV, images reconstructed from the M87 data show a prominent and asymmetric ring ("crescent") structure.In this section we use the techniques described in Section 4 to fit the M87 data sets with a specific class of geometric crescent models.
We first quantify the preference for crescent structure in Figure 2, which summarizes the results of fitting a series of increasingly complex geometric models to the M87 data.We Relative log-likelihood values for different geometric models fit to the M87 data as a function of nominal model complexity; the number of parameters is given in parenthesis for each model.April 5 is shown here, and all days and bands show the same trend.The models shown in this figure are strict subsets of the "generalized crescent model" (labeled here as model "n"; see Section 5.1), and they have been normalized such that the generalized crescent model has a value of  = 1; the reduced-χ 2 for the generalized crescent fit is 1.24 (see Table 2).We find that the data overwhelmingly prefer crescent models over, e.g., symmetric disk and ring models, and that additional Gaussian components lead to further substantial improvement.Note that a difference of ∼5 on the vertical axis in this plot is statistically significant.see that simple azimuthally symmetric models (e.g., uniform disks, rings) do a poor job of fitting the data; indeed, the strong detection of nonzero closure phases alone precludes such models.Models that allow for a central flux depression and a degree of asymmetry (e.g., double Gaussians or crescents) show significantly better performance.The most substantial gain in fit quality occurs for the crescent family of models.A top-hat crescent model (the difference of two uniform disks with the inner disk shifted, described by five parameters; see Kamruddin & Dexter 2013, Appendix B) performs vastly better than a top-hat ring model (three parameters).It also significantly outperforms the sum of two circular Gaussians (six parameters) and even the sum of two elliptical Gaussians (10 parameters).Adding parameters to the simplest crescent model continues to result in statistically significant, but comparatively modest, improvements.

Generalized Crescent Models
Among the large number of potential crescent-like models, we aim for one having the simplest geometry that is capable of both adequately fitting the M87 data and constraining several key observables.The geometric parameters of interest are the crescent diameter, its width and orientation, the sharpness of the inner edge, and the depth of any flux depression interior to the crescent.
With these key features in mind, we use an augmented version of the "slashed crescent" construction from Benkevitch et al. (2016) to provide the basis for a family of "generalized crescent" (GC) models.We refer to the two variants of the GC model that we use to fit the M87 data as xs-ring and xsringauss.Both GC models can be constructed in the image domain using the following procedure (see Figure 3).

Starting with a uniform circular disk of emission with
radius R out , we subtract a smaller uniform disk with radius R in that is offset from the first by an amount r 0 .The resulting geometry is that of a "top-hat crescent." 2. We apply a "slash" operation to the top-hat crescent, which imposes a linear brightness gradient along the symmetry axis.The brightness reaches a minimum of h 1 and a maximum of h 2 .3. We add a "floor" of brightness K to the central region of the crescent.For the xs-ring model this floor takes the form of a circular disk, while for the xs-ringauss model we use a circular Gaussian with flux density V F and width s F .The total flux density of the crescent plus floor component for the xs-ring model is denoted as V 0 .4. For the xs-ringauss model, we add an elliptical Gaussian component with flux density V 1 whose center is fixed to the inner edge of the crescent at the point where its width is largest (see the right panel of Figure 3), and whose orientation is set to align with that of the crescent.This fixed Gaussian component is inspired by the "xringaus" model from Benkevitch et al. (2016), which in turn sought to reproduce image structure seen in simulations from Broderick et al. (2014).5.The image is smoothed by a Gaussian kernel of FWHM σ * .6.The image is rotated such that the widest section of the crescent is oriented at an angle f in the counterclockwise direction (i.e., east of north).
The xs-ring model is described by eight parameters, while the xs-ringauss model is described by 11 parameters.Though it is useful to conceptualize the GC models via their image domain construction, in practice we fit the models using their analytic Fourier domain representations.The Fourier domain construction of both models is described in Appendix B, along with a table of priors used for the model parameters.Throughout this Letter we fit the xs-ring model using the dynesty sampling code, while the xs-ringauss model is implemented as part of THEMIS.Below, we define the various desired key quantities within the context of this GC model parameterization.We define the crescent diameter d to be twice the average of the inner and outer crescent radii, The fractional crescent width fw is defined in a similar manner, as the mean difference between the outer and inner radii (normalized to the diameter) plus a term to account for the FWHM, * s s º ( ) 2 2ln 2 , of the smoothing kernel: The sharpness ŝ is the ratio of the FWHM of the smoothing kernel to the crescent diameter, i.e., The fourth quantity of interest is the ratio, fc , between the brightness of the emission floor (interior to the crescent) and the mean brightness of the crescent, The different specifications for the xs-ring and xs-ringauss models (see Appendix B) means that these brightness ratios must be computed differently, , xs ringauss , xs ring. 24 Finally, we determine the orientation of the crescent directly from the f parameter, such that f f = ˆ.
In addition to the crescent component of the GC model, we also include a small number (two to three) of additional "nuisance" elliptical Gaussian components intended to capture extraneous emission around the primary ring and to mitigate other unmodeled systematics; the parameterization and behavior of these additional Gaussian components are described in Appendix B.2. GRMHD simulations of M87 often exhibit spiral emission structures in the region immediately interior and exterior to the photon ring (see Paper V), and M87 is known to have a prominent jet that extends down to scales of several Schwarzschild radii (e.g., Hada et al. 2016, Kim et al. 2018).Such extra emission is unlikely to be adequately captured by the crescent component alone, and the nuisance Gaussian components serve as a flexible way to model generic emission structures.We define the total compact flux density, ĈF, of the model to be equal to the summed contributions from the crescent and nuisance Gaussian components, where V 0 is the crescent flux density, V 1 is the flux density of the fixed Gaussian component, V F is the flux density of the central emission floor, and V g,i are the flux densities of the nuisance Gaussian components.

M87 Fit Results
We carry out independent fits to all days and bands using both the THEMIS-and dynesty-based codes.We show example GC model fits to the April 6 high-band data in Figure 4, with their corresponding image domain representations shown in Figure 5.There are no apparent systematic trends in the normalized residuals for either model, and we see similar behavior in the residuals from fits across all data sets.The corresponding THEMIS gain reconstructions are described in Appendix A. While our model-fitting procedures formally optimize a likelihood function modified by a prior (see Section 4.1), reduced-χ 2 values serve as general-purpose distance metrics between model and data.In our case, reduced-χ 2 values have the added benefit of enabling cross-comparisons with the images produced in Paper IV.The reduced-χ 2 expressions we use are detailed in Appendix C, and Table 2 lists their values for all M87 fits.
In general we find joint reduced-χ 2 values for the fits of ∼0.9-1.5.The number of degrees of freedom in the data ranges from ∼20 to 120, corresponding to an expected reduced-χ 2 deviation from unity of ∼0.06-0.16if our likelihoods follow a χ 2 -distribution.The implication is that while we see little evidence for overfitting, there are instances in which the residual values are distributed more broadly than the data error budget would nominally permit.The fact that such models sometimes underfit the data then indicates that we need to empirically determine our uncertainties, as posterior widths alone may not be reliable in the face of a statistically poor fit.We describe the empirical determination of the "observational uncertainties" in Appendix D.3, and we list the derived values in Table 3.
An illustrative pairwise parameter correlation diagram for both the xs-ring and xs-ringauss GC model fits to the April 5 high-band data set is shown in Figure 6, with single-parameter posteriors plotted along the diagonal.Figure 7 shows constraints on the GC model crescent component parameters, split up by data set.We also list the best-fit crescent parameters, as defined in Section 5.1, in Table 3.
Despite the differences in data products, sampling procedures, and model specifications, we find broad agreement between the derived posteriors from the two fitting codes.We note that the systematically wider posteriors from the xs-ring fits (by anywhere from a few to several tens of percent) are seen across all data sets and are an expected consequence of the use of closure amplitudes rather than visibility amplitudes.
Our primary parameter of interest is the diameter of the crescent, d.The weighted mean xs-ringauss value of 43.4 μas is in excellent agreement with the corresponding xs-ring value of 43.2 μas, with an rms scatter in the measurements of 0.64 μas and 0.69 μas, respectively, across all days and bands.This remarkable consistency provides evidence for the diameter measurement being robustly recoverable.The posterior widths of diameter measurements for individual data sets are typically at the ∼1% level, but our empirically determined uncertainties associated with the changing (u, v)-coverage and other  given by c y observational systematics (see Appendix D.3) place a more realistic uncertainty of ∼2% on this measurement.This more conservative uncertainty value is consistent with the magnitude of the scatter observed between measurements across different days and bands.
The fractional width fw of the crescent is considerably less well constrained than the diameter, with a scatter between data sets that is comparable to the magnitude of individual measurements.Furthermore, we see systematically larger fractional width measurements from the xs-ring model than from the xs-ringauss model (see Figure 7).Some of this offset is expected from differences in model specifications.A larger contributor to the discrepancy is likely to be the different data products being fit by the two models.For the data sets in which the xs-ring model prefers a significantly larger fw than the xsringauss model, we find that the smoothing kernel (described by ŝ) for the xs-ring model is also systematically wider.This smoothing by a Gaussian has no effect on the modeled closure phases.The additional gain amplitude degrees of freedom permitted in the xs-ringauss model fits can thus compensate for the smoothing in a manner that is not possible for the xs-ring model, which is fit only to closure amplitudes.Such smoothing affects the inferred diameter in a correlated manner that is well understood (see Section 7 and Paper IV, their Appendix G).
In both models, and across all data sets, we consistently measure a value for fw that is significantly smaller than unity.This rules out a filled-in disk structure at high confidence.We find instead that the emission must be concentrated in a relatively thin annulus, with a fractional width of 0.5, indicating the presence of a central flux depression.The brightness ratio fc in this hole is also well constrained: we consistently measure ).This value corresponds to a brightness contrast between the crescent and hole of at least a factor of 10.
We find a sharpness  ŝ 0.2, indicating that the smoothing kernel stays smaller than ∼20% of its diameter.The inner and outer edges of the crescent are therefore well defined, even if their locations are uncertain due to the large uncertainty in the width measurement.
We find that the crescent position angle f ˆconsistently confines the brightest portion of the crescent to be located in the southern half.We see some evidence for a net shift in orientation from ∼150°-160°(southeast) on April 5-6 to ∼180°-200°(south/southwest) on April 10-11, which amounts to a difference of ∼20-50 degrees between the two pairs of days.The direction of this shift is consistent with structural changes seen in the images from Paper IV, though the magnitude is a factor of ∼2 larger.
The xs-ringauss model fits find a typical compact flux density value of » ĈF 0.75 Jy.The inter-day measurement scatter is at the ∼50% level, which is consistent with the expected magnitude of the observational uncertainties as predicted by synthetic data in Appendix D.3.We find that the modeled ĈF value is less well-constrained than, but in good agreement with, the determined in Paper IV from consideration of both EHT and multi-wavelength constraints.

Calibrating the Crescent Diameter to a Physical Scale
Using GRMHD Simulations Though the GC models have been constructed to fit the M87 data well, the geometric parameters describing these models do not directly correspond to any physical quantities governing the underlying emission.Our primary physical parameter of interest is Contours enclose 68% and 95% of the posterior probability.Note that recovery of the total compact flux density ĈF is not possible for the dynesty-based fits, which only make use of closure quantities.the angular size corresponding to one gravitational radius, The gravitational radius sets the physical length scale of the emission region.Most of the observed 230 GHz emission is expected to originate near the photon ring (see, e.g., Dexter et al. 2012), whose scaling with θ g is known for a given black hole mass and spin a * (Bardeen 1973;Chandrasekhar 1983).The crescent component in the GC models does not necessarily correspond to the photon ring itself.If, however, the crescent component is formed by lensed emission near the horizon, then Posterior medians and 68% confidence intervals for selected parameters derived from GC model fitting for all observing days and bands.Blue circular points indicate xs-ring fits using the dynesty-based fitting scheme applied to individual data sets (i.e., a single band on a single day).Red square points indicate xs-ringauss fits using THEMIS applied to individual data sets, while orange square points show THEMIS-based xs-ringauss fits to data sets that have been band-combined.All plotted error bars include the systematic "observational uncertainties" estimated from simulated data in Appendix D; these uncertainties are listed in the bottom row of Table 3.Note that recovery of the total compact flux density ĈF is not possible for the dynesty-based fits, which use only closure quantities.The light purple band in the lower-right panel is the range inferred in Paper IV. its size should obey a similar scaling with θ g , g For emission at the photon ring, α;9.6-10.4depending on the black hole spin parameter.For a realistic source model, the emission is not restricted to lie exactly at the photon ring.The value of this scaling factor α and its uncertainty are therefore unknown a priori.We have measured α and its uncertainty for both GC models using a suite of synthetic data sets generated from snapshots of GRMHD simulations from the GRMHD image library (Paper V).The full calibration procedure, including properties of the selected GRMHD simulations, the synthetic data generation process, and the calibration uncertainty quantification, is detailed in Appendix D. By fitting each calibration image with a GC model, and then comparing the corresponding d measurement to the input value of θ g for the simulation that produced the image, we determine the value of α for that image.Combining the results of such fits from a large number of GRMHD simulations yields a calibration (and uncertainty) for α.For the xs-ring model we find a mean value of a = 11.55, while for the xs-ringauss model we find a nearly identical a = 11.50.Both of these values are somewhat larger than the a » 10 expected for the photon ring itself, indicating that the GC models are accounting for emission in the GRMHD model images that preferentially falls outside of the photon ring.
Figure 8 shows the θ g values obtained as a result of applying our calibrated scaling factor to the crescent diameter values measured for each day and band, and Table 4 lists the results from combining the measurements from all data sets.There is excellent agreement between the two GC models, resulting in an averaged value of q m = -+ .We note that the 12% uncertainty in the θ g measurement is dominated by the diversity of GRMHD models used in the primary calibration; the quantification of this "theoretical" uncertainty component from the GRMHD simulations is described in Appendix D.2.

Direct Comparison with GRMHD Models
EHT data have the power to directly constrain GRMHD simulation based models of M87 and to estimate the physical properties of the black hole and emitting plasma.Such a direct comparison is challenging due to stochastic structure in the models.
The EHT 2017 data span a very short time frame for the M87 source structure.Its characteristic variability timescale at 230 GHz is ;50 days (Bower et al. 2015), much longer than our observing run.Thus, although the entire ensemble of model snapshots taken from a given simulation captures both the persistent structure as well as the statistics of the stochastic components, we do not yet have enough time coverage for M87 itself to measure its structural variations.
Model images from GRMHD simulations show a dominant, compact, asymmetric ring structure resulting from strong gravitational lensing and relativistic gas motions (Paper V).Hence, they capture the qualitative features found by image reconstructions in Paper IV and by geometric crescent models in Section 5.This motivates a direct comparison of the GRMHD model images with the EHT data.
In this section we summarize the GRMHD image library (Section 6.1) and fit individual simulation snapshot images (Sections 6.2 and 6.3) in a similar fashion to past work on  Note.We use the angular diameter measurements ( d ) from Section 5.2 and combine them with the calibrated scaling factors (α) from Section 5.3 to arrive at our measurements of θ g .We list scaling factors calibrated using both magnetically arrested disks (MAD) and standard and normal evolution (SANE) GRMHD simulations (see Section 6.1), as well as ones calibrated using only MAD and only SANE simulations; for the final measurement presented in the text we have used the MAD +SANE calibration.The various uncertainty components are described in Appendix D.2.We quote median values for all measurements and the associated 68% confidence intervals for the different categories of uncertainty.
GRMHD fitting to mm-VLBI data of Sgr A * (e.g., Dexter et al. 2010;Kim et al. 2016).We further develop and apply a method for testing the consistency of the M87 data with the simulation models (Section 6.4) and find that the majority of the simulation library models is consistent with the data.We use the results to estimate physical parameters, including the black hole angular radius GM Dc 2 (Section 6.5).The implications of our results for the physical properties of the emission region are discussed in more detail in Paper V. From EHT data alone it is difficult to rule out many of the broad range of possible models for the black hole and plasma properties.However, in combination with other data (especially the observed jet power), ultimately more than half of the models can be excluded.

Summary of Simulations
As described in detail in Paper V, we have constructed a large image library of horizon-scale synchrotron emission images at 230 GHz computed from GRMHD simulations.We summarize the broad features of this library here and direct the reader to Paper V for more information.The GRMHD simulations cover a wide range of black hole spins as well as initial magnetic field geometries and fluxes.These result in images associated with a variety of accretion flow morphologies and degrees of variability.The magnetic flux controls the structure of the accretion flow near the black hole.Low magnetic fluxes produce the standard and normal evolution (SANE) disks characterized by low-efficiency jet production.In contrast, magnetically arrested disks (MAD) are characterized by large magnetic fluxes, set by the ram pressure of the confining accretion flow.
From these models, families of between 100 and 500 snapshot images were produced assuming synchrotron emission from an underlying thermal electron population (see Section 3.2 of Paper V).The snapshot image generation introduces additional astrophysical parameters associated with the intrinsic scales in the radiative transfer.These parameters include the black hole mass, the viewing inclination, i, and a model for the electron thermodynamics: T i /T e ≈R high in gas-pressure dominated regions and is unity otherwise (Mościbrodzka et al. 2016), where T i and T e are the ion and electron temperatures.
The number density of emitting electrons is scaled independently for each simulation such that the typical 230 GHz flux density is ∼0.5Jy.The temporal separation between snapshots is selected such that adjacent snapshots are weakly correlated.

Single Snapshot Model (SSM)
Each snapshot image generates a three-parameter SSM defined by the total compact flux (CF), angular scale (θ g ), and orientation (defined to be the position angle of the forward jet measured east of north, PA FJ ).Variations in these parameters approximately correspond to variations in the accretion rate, black hole mass, and orientation of the black hole spin, respectively.Variations in mass are associated with changes in the diameter of the photon ring, a generic feature found across all of the images in the GRMHD image library.
Each snapshot image is characterized by a nominally scaled, normalized intensity map of the image, ˆ( ) I x y , , with a corresponding nominal total intensity ĈF, gravitational radius q ˆg, and forward jet position angle PA FJ =0°; associated with the intensity map are complex visibilities ˆ( ) V u v , .The SSM is then generated by rescaling, stretching, and rotating ˆ( ) V u v , : where I º ĈF CF, m q q º ĝ g , and (u′, v′) are counter-rotated from (u, v) by the angle PA FJ .This procedure is illustrated in Figure 9.
We show in Paper V that these approximations generally hold for flux and mass for rescalings by factors of 2 from their fiducial values.

Fitting Single Snapshots to EHT Data
For both model selection and parameter estimation, the first step is fitting an SSM model to the EHT data set described in detail in Section 2.1.The only difference here is that intra-site baselines are excluded.These probe angular scales between 0 1 and 10″, at which unmodeled large-scale features, e.g., HST-1, contribute substantially (see Section4 of Paper IV).We verify after the fact that the reconstructed compact flux estimates are consistent with the upper limits necessarily implied by these baselines.The fitting process is complicated by large structural variations between snapshots resulting from turbulence in the simulations (see Section 6.4).We employ two independent methods to fit SSMs to the EHT data sets.Both employ the likelihoods constructed as described in Section 4.1 for visibility amplitudes and closure phases.The two methods, THEMIS and GENA, are utilized as described in Section 4, producing posterior estimates or best-fit estimates, respectively, for CF, θ g , and PA FJ .
We adopt a uniform prior on the total compact flux density between 0.1 and 10 Jy for THEMIS and between 0.1 Jy and 4 Jy for GENA; both of these priors cover ranges that substantially exceed the limits placed by Paper IV.While the position angle (PA) of the large-scale radio jet in M87 is well known, we permit the PA of the horizon-scale jet to be unconstrained, setting a uniform prior from [0°, 360°).Finally, we place a flat prior on θ g ranging from 0.1 to 100 μas for THEMIS and between 0 and 10 μas for GENA, again substantially exceeding the physically relevant ranges.

Model Selection and Average Image Scoring (AIS)
Quantitatively assessing the quality of SSM fits presents unique challenges because of the presence of stochastic image features due to turbulence in the underlying GRMHD simulations that produce large variations in image structure.The structural variability leads to changes in EHT observables that are much larger than the observational errors.It is therefore not feasible to generate a sufficient number of images from existing GRMHD simulations to have a significant probability of finding a formally adequate (i.e., χ 2 ≈1) fit to the data (see the discussion in Paper V). Nevertheless, the ensemble of snapshots from a given simulation provides a natural way to characterize the impact of these stochastic features on the inferred SSM parameters.
We can thus assess individual GRMHD simulations, effectively comprised of many (100-500) snapshot images, by comparing the quality of SSM fits to a numerically constructed distribution of reduced χ 2 values.Note that this procedure is conceptually identical to the normal fitting procedure, in which a χ 2 is interpreted relative to the standard χ 2 -distribution and thus a reduced χ 2 ≈1 is a "good" fit, with the exception that an additional source of noise has been effectively introduced as part of the underlying model, altering the anticipated distribution of χ 2 values.In practice, this comparison is executed via the AIS method described in Appendix F using the THEMIS SSM fitting pipeline.
We estimate the appropriate χ 2 -distribution by performing multiple fits of an SSM constructed from the arithmetic average snapshot image to simulated data generated from each snapshot image within the underlying GRMHD model.Finally, the average snapshot image is compared to the EHT data, and the resulting χ 2 is assessed.Thus, THEMIS-AIS is effectively determining if the EHT data are consistent with being drawn from the GRMHD simulation.The result is characterized by a simulation p-value, which we call p AIS .Because of the significant variations in data quality and baseline coverage across days, this procedure must be repeated independently for each day.
It is possible for a model to be ruled out by the AIS procedure both by the average image SSM having a χ 2 that is too high (the data is "further" from the average snapshot image than typical for the simulation) and by the average image SSM having a χ 2 that is too low (the data is "closer" to the average snapshot image than typical for the simulation, usually a consequence of too much variability in the GRMHD model).These are similar to finding a reduced χ 2 that is much larger and smaller than unity, respectively, in traditional fits.Both occur in practice.
The result of the THEMIS-AIS procedure is limited by the number of snapshots from each model, and thus is currently capable of excluding models only at the 99% level for a given day and band (i.e., < p 0.01

AIS
). Due to the long dynamical timescale of M87-GRMHD snapshots exhibit correlations extending up to 20GM/c 3 ≈1 week-this is not significantly improved by combining bands and/or days.Therefore, unless otherwise noted, we set a threshold of p AIS >0.01,below which we deem a GRMHD model unacceptable and exclude it from consideration.

Ensemble-based Parameter Estimation
The posteriors for the SSM parameters for a given GRMHD model are estimated using both the observational errors in the EHT data and the stochastic fluctuations in the snapshot images themselves; among these "noise" terms, the latter significantly dominates.This is evidenced by the variations among SSM fit parameters from snapshots within a GRMHD model, which typically exceed the formal fit errors from the SSM fit alone (see Appendix G).Some care must be taken in extracting and interpreting parameter estimates from these.
Only SSM fits with likelihoods among the highest 10% within each model are used for parameter estimation.The results for θ g are insensitive to the precise value of this fit-quality cut after it passes 50%.The measured values are consistent over the range of GRMHD models explored (see Figure 10 for five examples).The consistency is likely a result of the dominant strong lensing and relativistic motion that are common to all models.
The full posteriors are then obtained by marginalizing over all GRMHD models that are found to be acceptable via the THEMIS- AIS procedure.At this point we also include the ancillary priors on jet power, X-ray luminosity, and emission efficiency described in Paper V (see their Table2).This procedure corresponds to Figure 10.Distributions of recovered θ g from five representative GRMHD simulations as measured by the THEMIS (left; maroon) and GENA (right; green) pipelines.Only those snapshots for which the likelihood is above the median (THEMIS), or for which the combined χ 2 statistic is below the median (GENA) are included.For a wide range of simulations, the recovered θ g are consistent with each other.All of the simulations shown are deemed acceptable by AIS (p AIS >0.01;see Section 6.4 for details).
simply summing the posteriors obtained for each acceptable GRMHD model.
In principle, it is possible for the posteriors constructed in this fashion to suffer from biases induced by the treatment of the high values of χ 2 found when including only observational errors.We have attempted to estimate this bias using a large number of mock analyses.In those, simulated data were generated from a snapshot within a GRMHD model.Posteriors were then generated for that model given the simulated data and compared to their known SSM parameters; no significant biases were found.We have further conducted two posterior mock analyses for the full suite of GRMHD models, independently for the THEMIS and GENA pipelines.In these the impact of including "incorrect" GRMHD models was a key difference with the previous tests.See Appendix G for more details.
The results of both of these experiments only hold if the GRMHD models used as synthetic data provide a good description of M87.As multiple observation epochs become available, it will be possible to explicitly measure the statistics of the stochastic fluctuations, empirically addressing this assumption.It will also enable direct ensemble-to-ensemble comparison like that described in Kim et al. (2016).A promising, alternative approach would be to perform a principal component analysis (PCA) decomposition of the snapshots within each model  ( Medeiros et al. 2018), and fit images generated by varying the weights of the PCA components to the data to mimic the set of possible realizations of the turbulence.

M87 Fit Results
An example SSM fit is shown for both the THEMIS and GENA methods in Figure 11.As anticipated, the presence of stochastic features within the images results in a poor formal fit quality for all simulation snapshots, that is, reduced  c 2

2
. Nevertheless, broad features of the visibilities and closure phases are accurately reproduced.These include the deep visibility amplitude minimum near 3.4Gλ and amplitude of the following bump at 6.0Gλ, both of which are key constraints on the image structure (see Section 3).Significant structure in the residuals, as seen in Figure 11, is a natural consequence of the presence of stochastic model components, but may also signify that the underlying model is insufficient to fully explain the EHT data.Representative fits obtained via this process are shown in Figure 12 for the best SSM from a selection of models that are not ruled out by THEMIS-AIS.
Even aggressive model selection via THEMIS-AIS produces a very weak cut on the GRMHD models-upon choosing p AIS >0.01 and p AIS >0.1 only 9.7% and 28.5% of models are excluded.After applying the additional observational constraints, we discard 65.3% of the models (for details of the model selection, see Paper V).Here we focus on estimates of the compact flux density, θ g and PA FJ , after marginalizing over the acceptable models that remain.
The GRMHD-based SSM models produce compact flux estimates between 0.3 and 0.7Jy.This is at the lower end of the range reported in Paper IV ( ).The PA FJ obtained via comparison of the GRMHD snapshots and the f found from fitting the GC models are generally consistent after accounting for the average 90°offset between the location of the forward jet in the GRMHD simulation and the location of the brightest region in the crescent image (Paper V).This appears marginally consistent with the mas-scale forward-jet PA measured at 3 and 7 mm of 288°(e.g., Walker et al. 2018), though see Paper V for further discussion on this point.
The posteriors for θ g are broadly consistent among days and bands, as illustrated in Figure 13 and Table 5.The combined value for both the THEMIS and GENA analyses is q m = -+

as
. Finally, note that the process by which this estimate is arrived at differs qualitatively from the analysis presented in Section 5.Here GRMHD snapshot images are fit directly to data, whereas in Section 5 GRMHD snapshots were used to calibrate the geometric models.These subsets are independent-the set of images used to calibrate the geometric models are not used in the GRMHD analyses.Nevertheless, a systematic correlation between these estimates may remain as a result of the use of the same set of underlying GRMHD simulations.

Image Domain Feature Extraction
In the previous sections, we fit geometric and numerical models in the visibility domain to measure the properties of features in the models that give rise to the observed interferometric data.In Paper IV, we performed direct image reconstruction using two RML methods (eht-imaging and  We describe the method used to extract parameters (Section 7.1).We then show that the extracted ring diameter is consistent across imaging methods and compare the ring diameter and width with that measured from geometric modeling (Section 7.2).We convert the diameter measurements to the physical scale around the black hole, θ g , by calibrating with GRMHD simulations (Section 7.3).Finally, we measure the image circularity (Section 7.4).

Measuring the Parameters of Image Features
We describe the M87 image morphology using five quantities that are similar in spirit but not identical to those introduced in Equations ( 20)-( 23) in Section 5.1: the diameter d of the ring/ crescent, its fractional width f w , the relative brightness f c of its center with respect to its rim, and the PA f. Figure 14 shows a schematic diagram of these parameters.Here we focus on the first two parameters.Measurements of the remaining parameters in the various image reconstructions and more details on the method can be found in Paper IV.
The reconstructions obtained in Paper IV use a pixel size of 2 μas (;1/10 of the EHT 2017 beam), but the GC image models considered below often show narrower structures.We interpolate both the GC models and image reconstructions onto a grid with pixel size of 0.5 μas before applying the feature extraction methods.When necessary, we use a 2D linear interpolation between adjacent pixels for finer sampling.
The radial brightness profile is characterized by a peaked distribution that declines toward the center.We first locate the (arbitrary) center of the ring.We measure the position of the radial brightness profile maximum along different azimuthal angles.For a trial center position, we calculate the dispersion of the radii defined by these maxima.The center is then chosen iteratively as the location that minimizes this dispersion.Once a center position is chosen, the mean diameter d of the image is defined as twice the distance to the peak, averaged over azimuth.
We define the ring width to be the FWHM of the bright region along each radial profile.The fractional width f w of the image is the average of the FWHM over azimuthal directions divided by its mean diameter.In the following we first smooth the image with a 2μas Gaussian.In Appendix H we show that variations on this method produce small (sub-pixel) differences in the results.

M87 Image Ring Diameters and Widths
Figure 15 shows the mean diameters of images reconstructed using the low-band data during all four days of M87 observations with three image reconstruction methods (ehtimaging, SMILI, DIFMAP; see Paper IV). Image samples are reconstructed for each data set and for a wide range of weights of the regularizers (∼2000 images for eht-imaging and SMILI and ∼30 for DIFMAP, see Paper IV). Diameters measured from the full set of images that produce acceptable fits to the visibility data (the "Top Set") are shown here.The diameters of the ring features found across all methods and all days span the narrow range ∼38-44 μas.
Figure 16 shows the fractional width versus mean diameter for the Top Set images from eht-imaging and SMILI compared with those from geometric models for the low-band M87 data of April 6 (Section 5).For the two visibility domain methods, the points correspond to the diameters and widths obtained from a sample of 100 images from the xs-ring and xsringauss model, chosen randomly from their posteriors.Those model images have been analyzed in the same way as the reconstructed ones.
The fractional widths for the reconstructed images are 0.5, which are consistent with the results of geometric crescent models in Section 5.While both the image reconstructions and model fits show a large uncertainty in fractional width, the widths measured from the image reconstructions in the Top Sets are 10 μas.Images of rings with narrower widths are consistent with the data and can be produced by the imaging algorithms; however, the parameters in the imaging algorithms (e.g., regularizer weights) that determine the Top Set images were trained on synthetic data from sources smoothed with a 10 μas beam (Paper IV, their Section 6.1).This may help explain the differences in the fractional widths measured with the different techniques.
An anti-correlation between diameter and fractional width for these image domain results is clearly present in three of the four days, but is less clear in the data of April 11.We consider two manifestations for this anti-correlation, based on the location of the first visibility amplitude minimum in simple ring models.First, we use the geometric crescent model of Kamruddin & Dexter (2013), and derive the following approximate fitting formula to the exact expression relating the mean ring diameter d and ψ: where d 0 is the diameter of an infinitesimally thin ring that produces a null at baseline length b 1 (Equation ( 6)).The fractional width is . The blue shaded region in Figure 16 shows the expected diameter and fractional width anti-correlation for a visibility minimum occurring at a range of baseline lengths 3.5-4.0Gλ, similar to that seen in the M87 data (Section 3, Paper III).This corresponds to d 0 ;39.5-45μas.
We repeat this exercise for an infinitesimally thin ring convolved with a Gaussian.Here the mean diameter d can be approximated for FWHM  w d as (AppendixG of Paper IV): Here the fractional width The pink shaded area in Figure 16 shows the shaded region for this model, which follows a similar trend.
The measured properties of the images and source models inferred by all methods generally fall within the expected bands.At least part of the systematic differences in our diameter measurements may be attributed to the relatively large uncertainty in width, as a result of their weak anti-correlation.

From Image Diameter to Angular Gravitational Radius
We can convert the diameters measured from the reconstructed images to the black hole angular gravitational radius using a scaling factor, α (Equation ( 27)), following the same procedure used in the calibration of crescent model diameters to GRMHD images discussed in Section 5.3.We calculate a separate value of α for each image reconstruction method (eht-imaging, SMILI, DIFMAP) as described in Appendix E and listed in Table 6.The image domain methods do not report posteriors.We therefore estimate the statistical along with the observational component of the uncertainty in the calibration procedure using synthetic data.As found for the geometric crescent models, the theoretical uncertainty dominates the final error budget.The total observational uncertainty is similar to but slightly larger than the statistical spread of the median diameter measurement between days.
After applying the calibrated scaling factor to the image domain diameter measurements, we find consistent results across all observing days and reconstruction methods (Figure 17).This is further indication that the statistical component of the calibration error is sub-dominant.The combined value from all methods is q m = -+

as
. Despite small differences in the measured mean diameters, the physical scale θ g is remarkably consistent with that found earlier from both geometric and GRMHD model fitting.This is because the corresponding calibration factors α obtained from synthetic data show the same trends as the measured diameters.

The Circular Shapes of the M87 Images
Our reconstructed M87 images appear circular.We quantify their circularity by measuring the fractional spread in the inferred diameters measured along different orientations for each of the reconstructed images.Here we define the fractional spread as the standard deviation of the diameters measured along different orientations divided by the mean diameter.For each image, the diameters along different orientations measure  Note.We quote median values and 68% confidence intervals.The angular diameters d are averages across the four observing days, weighted by the range within each day.Their uncertainties s d are the standard deviation of the mean over the four days.We combine the angular diameter measurements with the calibrated scaling factors (α) to arrive at measurements of θ g .The various uncertainty components of θ g are obtained using synthetic data, as described in Appendix E. the distance of the location of peak brightness from the center of the image along each orientation.A small fractional spread suggests that the locus of peak brightness on each image has a circular shape.
Figure 18 shows the distribution of the fractional spread in diameters for the Top Set images of M87 reconstructed using two image domain methods (eht-imaging and SMILI) for the April5 low-band data set.The figure also shows the fractional spread in diameters of a subset of ∼300 images among those in the GRMHD image library discussed in Paper V that provide acceptable fits according to AIS.Note that these comparison images have a much higher intrinsic resolution than the reconstructed images of M87.The distributions of fractional spreads for the reconstructed images peak at a fractional spread of ∼0.05-0.06.This is within the range found in the GRMHD images.The GRMHD models mostly show circular image structure (peak 0.1 in blue distribution in Figure 18).This corresponds to the bulk of the images where the photon ring is the dominant feature (see Paper V).The tail to larger fractional diameter spread may in part be the result of occasional bright filaments that light up in a small azimuthal range and then quickly fade.During such snapshots, the locus of points with the brightest emission around the shadow becomes non-circular and the fractional spread increases.Arrows in Figure 18 indicate the expected fractional diameter spread for elliptical models.Our results suggest an image axial ratio of 4:3.

Prior Mass and Distance Estimates
In this section we briefly review prior estimates of the massto-distance ratio for M87ʼs central black hole for comparison with our new EHT results.For a more detailed discussion, see Appendix I.

Distance
M87 is the central galaxy of the Virgo cluster, the closest galaxy cluster to the Milky Way.This proximity allows a highprecision measurement of its distance using primary or secondary standard candles available only for nearby galaxies within ∼20 Mpc.Among numerous distance estimators obtained for M87, we base the distance that we use here on two estimators of the absolute distance modulus (μ).The first uses as a primary candle the tip of the red giant branch (TRGB method).The second uses a well-calibrated secondary standard candle, the surface brightness fluctuations (SBF method) from resolved stellar distributions observed at high angular resolution.
We combine one TRGB measurement (Bird et al. 2010) with two SBF measurements (Blakeslee et al. 2009;Cantiello et al. 2018a) to arrive at a combined distance measurement.We consider the three distance measurements to be independent and calculate the combined posterior as their product.Based on these, we adopt a distance of = -+ D 16.8 Mpc 0.7 0.8 to M87.A more detailed description, together with the summary of each measurement, can be found in Appendix I.

Mass-to-distance Ratio and Related Quantities
For the mass of the central black hole in M87 we focus on the two latest dynamical measurements of the kinematic properties of surrounding stars and gas (Gebhardt et al. 2011;Walsh et al. 2013) that have high-quality data and state-of-theart modeling approaches.
Stellar dynamical measurements.Gebhardt et al. (2011) used infrared spectroscopic observations of stellar absorption lines over the central 2″ (∼150 pc) of M87 and modeled the observed stellar kinematics to infer a mass of (6.6±0.4)× 10 9 M e , assuming a distance of 17.9 Mpc.They compared the full line-of-sight stellar velocity distributions at all points with the projected kinematics of axisymmetric models of the galaxy.The models allowed for adjustable but radially constant massto-light ratios of the central stellar population and the central effects of the M87 dark matter halo.The analysis demonstrates that kinematics cannot be fitted without the inclusion of a central compact non-luminous mass.
Gas dynamical measurements., assuming a distance of 17.9 Mpc, on the assumption that the velocity field is Keplerian and is not subject to non-gravitational forces such as shocks and winds.Several factors could account for the lower mass inferred via this technique (see, e.g., Kormendy & Ho 2013;Jeter et al. 2018).

Lensed Emission Around the Black Hole Shadow
The results of fitting analytic models (Section 5) and GRMHD image snapshots (Section 6) to visibility domain data, as well as image reconstructions (Paper IV, Section 7) strongly favor an asymmetric ring (crescent) source morphology with a deep central brightness depression.This morphology is robust among the different analysis methods and is consistent with characteristic features in the interferometric visibilities (Section 3).
Fitting geometric models to the data and extracting feature parameters in the image domain both allow us to quantify the following properties of the crescents: (i) a compact flux density Figure 18.Distributions of fractional spreads in the inferred diameters of the M87 images measured along different orientations, for images reconstructed using two methods for the April 5 data set.The blue histogram is the distribution of fractional spreads for a subset of images from the GRMHD library.The values for elliptical shapes with different axes ratios are indicated with arrows (a value of 0 is an axis ratio of 1:1).The distribution for reconstructed images peaks at values ∼0.05-0.06,which is consistent with the range found from GRMHD model images.
of ;0.7 Jy, which is 50% of the total at arcsecond scale (see also Paper III and Paper IV); (ii) a mean emission diameter of ∼42±3 μas; (iii) a fractional width of 0.5; (iv) a deep central brightness depression with a brightness contrast ratio of 10; and (v) a persistent asymmetry with the brightest region to the South (PA of 150°-200°).
All of these features support the interpretation that we are seeing emission from near the event horizon that is gravitationally lensed into a crescent shape near the photon ring.The central flux depression is the result of photons captured by the black hole: the black hole shadow (Falcke et al. 2000).The asymmetry is a result of Doppler beaming due to the rotation of plasma at relativistic speeds.The peak brightness location is expected to be oriented ;90°away from the jet PA (Dexter et al. 2012;Mościbrodzka et al. 2016;Ryan et al. 2018;Chael et al. 2019b, see Paper V), roughly north-south as we find here.
In many active galactic nuclei, the radio VLBI jet core on parsecs and larger scales is inferred to be significantly offset from the black hole itself(e.g., Marscher et al. 2008).For M87, this is disfavored by past VLBI observations.The jet and counter-jet angular separation inferred by VLBA observations (Ly et al. 2004(Ly et al. , 2007;;Kovalev et al. 2007) locates the radio core close to the black hole.VLBA astrometric observations further tightly constrain the jet apex to be ;41±12 μas upstream of the VLBI core at 43 GHz (Hada et al. 2011(Hada et al. , 2013;;Nakamura & Asada 2013), suggesting that the mm emission originates very close to the black hole.
The observed EHT image morphology also disfavors an origin offset from the black hole.Models of M87 where the emission predominantly arises in the forward jet have found a larger size and blob-like Gaussian structure (e.g., Broderick & Loeb 2009) due to Doppler beaming resulting from the acceleration of the jet along the line of sight.Additionally, the stable source size over several years of mm-VLBI (Doeleman et al. 2012;Akiyama et al. 2015), despite changes in compact flux density by factors of ;2-3 (Paper IV), favors the extreme gravitational lensing scenario.Future EHT observations will provide a sensitive additional test of this paradigm by measuring the source morphology on year timescales.

The Black Hole Mass and Comparison to Prior Dynamical Measurements
We have used several techniques to arrive at θ g estimates from the EHT 2017 M87 data, and in Table 7 we list the result for each technique after averaging over observing days and bands (see also Figure 19).We find a striking level of consistency across all measurement methods, with GC modeling, direct GRMHD fitting, and image domain feature extraction procedures converging on a characteristic value of θ g =3.8±0.4μas as the angular size subtended by one gravitational radius.The measurement techniques themselves are entirely independent of one another.The mutual agreement that we see among the multiple measurements indicates that the θ g value we have converged upon is robust to the many different choices that can be made regarding data products, model specifications, and parameter space exploration algorithms.
All of the individual θ g estimates use the GRMHD simulation library, either through directly fitting GRMHD snapshots to the data (Section 6) or through calibration of diameters resulting from geometric models or reconstructed images (Sections 5 and 7).A degree of caution is therefore warranted.The measurements rely on images generated from GRMHD simulations and should be understood within that context (see also Paper V).We have also used only a small subset of 100 such randomly chosen frames (with randomly assigned SSM parameters) out of a much larger library.A simple (but poorly motivated) alternative to the full calibration presented here would be to assign the mean emission diameter to the size of the photon ring itself (α;9.6-10.4 for all methods).That would give a nearly identical θ g result for the image domain estimates and only a ;10% increase for the geometric models, within its current systematic uncertainty.
The mass measurements that we have carried out with EHT data show a high level of consistency, converging to an average value of M=6.5×10 9 M e .All measurement techniques share a large source of systematic uncertainty arising from the GRMHD calibration, with an average value of σ sys =0.7×10 9 M e .The geometric crescent modeling and image domain measurements further enable an estimate of the systematic uncertainty associated   7).Gray horizontal bands around dashed lines correspond to the 68% confidence levels for the stellar (top) and gas (bottom) dynamical mass measurements (see Table 9).
with the angular diameter measurement, which manifests in the mass measurement with an average, rounded value of σ stat =0.2×10 9 M e .As can be readily seen in Figure 19, our mass measurements are also in excellent agreement with those made by Gebhardt et al. (2011) using stellar dynamics (Section 8.2).Our mass measurement is inconsistent with the gas dynamics measurement of Walsh et al. (2013) at the ∼99% confidence level.

Implications for the Black Hole Environment
The measurement of a large black hole mass in M87 has a number of implications for the relationship between the black hole and its immediate environment.The observed line-of-sight velocities of atomic line emission produced within the gas flow at 10-100 pc place strong constraints on the dynamical state of gas on these scales and its relationship with the horizon-scale and milliarcsecond-scale features.
In combination with the stellar dynamics mass estimate, the mass M reported here provides an upper limit on the extended dark mass enclosed within the smallest scale probed by Gebhardt et al. (2011), 0 25 or approximately 20 pc, set by the upper limit on the difference between the two masses: , where M G11 is the stellar dynamics mass in Table 7 and the 68% confidence intervals have simply been added in quadrature.Thus, at 95% confidence this comparison implies that inside 20 pc 2.3×10 9 M e , or 36% of the central black hole mass, can exist as a dark extended component.This limit is not yet sufficient to constrain the presence of a density "spike" of dark matter associated with the formation of the central black hole (see, e.g., Lacroix et al. 2015).

Evidence for a Horizon
The constraint on the size of the emission region, coupled with multi-wavelength limits on the emission from within the shadow region, provide evidence for the existence of event horizons.These limits arise from looking for the boundarylayer emission that should be present if a photosphere is visible, i.e., not hidden from view by an event horizon (Narayan & McClintock 2008;Broderick et al. 2015).The source of such emission is often called a "surface," though we note that it need not be so in practice.
The results reported in the previous section strengthen this argument.First, they provide an accurate reconstruction of the source structure, which now locates the central mass to within the photon orbit.Second, as discussed in Paper V, the observed image is broadly consistent with models of jet launching driven by electromagnetic processes.Jets are collimated by accreting material; thus a significant boundary-layer emission component would be expected if M87 is not a black hole (Broderick et al. 2015).
The spectrum of this additional boundary-layer emission component is generally thermal when the photon escape fraction is small, in which case this emission effectively couples the photosphere to itself on a timescale that depends only logarithmically on redshift (though see Lu et al. 2017 regarding the possibility of longer equilibrium times).The color temperature of the emission is set uniquely by the accretion power and apparent size of the emission region.The former can be related to within factors of a few to the total jet power, ∼10 42 erg s −1 (see the discussion in Paper V).Note that this lower limit is an order of magnitude smaller than the lower end of the range considered in Broderick et al. (2015).The resulting putative photospheric emission would have a color temperature of T C 10 4 K and peak in the near-infrared/ optical.
While the EHT is not sensitive to this component, its existence can be excluded by near-infrared, optical, and ultraviolet flux measurements (see, e.g., Broderick et al. 2015;Prieto et al. 2016).We show this excess in Figure 20 for the mass measurement reported here.

Implications for the Kerr Nature of the Black Hole
The detection of a substantial central brightness depression in the image of M87 that we identify with the shadow of the black hole allows us, in principle, to perform two tests of the nature of the compact object and its spacetime.One is related to the size of the shadow, and the other to its shape.
There are two properties of spinning black holes that, depending on spin and observer inclination, could alter the size and shape of the black hole shadow.The spacetime quadrupole q tends to lend an oblate shape to the shadow, while frame dragging due to the spin a acts to compress the shadow in the direction perpendicular to its spin.When the spacetime quadrupole satisfies the Kerr condition q=−a 2 , the two effects nearly cancel each other out, leaving a quasi-circular shadow with a mean diameter that lies in the very small range (9.6-10.4)GM/c 2 .This shape is expected to be nearly circular for the low inclination of M87 (e.g., Takahashi 2004;Chan et al. 2013).
For a black hole of known mass-to-distance ratio, the angular size of the shadow on the sky is, therefore, fixed and can be used as a null-hypothesis test of the Kerr nature of the compact object (Psaltis et al. 2015).Alternatively, we can infer the mass-to-distance ratio of the black hole using the size of its shadow image, as we do here, and compare it to the Figure 20.Spectral energy distribution of putative photosphere in comparison to EHT and multi-wavelength data.The lower limit for the thermal bump associated with boundary-layer emission in the absence of a horizon is shown for a compact object with the size reported here (dark green band and arrow) and a jet power 10 42 erg s −1 .The photospheric component has been added to a jet synchrotron emission model (blue lines, M. Lucchini et al. 2019, in preparation).EHT compact flux densities and limits (orange points) are insensitive to a putative boundary-layer component.Its presence is ruled out by near-infrared, optical, and ultraviolet measurements (dark red and black points, see Broderick et al. 2015;Prieto et al. 2016).mass-to-distance ratio measured dynamically at much larger distances.Demonstrating agreement between the two inferences results in an equivalent null-hypothesis test.
The null hypothesis in this case consists of three ingredients: that the dynamical measurements provide an accurate determination of the black hole mass, that the brightness depression we detect in the EHT image of M87 is indeed the black hole shadow, and that the spacetime of the black hole is described by the Kerr solution.Demonstrating a violation of this null hypothesis would imply that one or more of these assumptions is invalid.
We have used the 2017 EHT observations of M87 to infer the posterior probability P obs (θ g ) of the angular gravitational radius at the distance of M87.Earlier dynamical measurements also provide us with the prior probability P dyn (θ dyn ) for the same quantity.We can, therefore, measure the fractional difference δ between the predicted and measured mass-todistance ratio for the black hole as and its posterior as . 32 dyn dyn dyn obs dyn dyn Figure 21 shows the posterior on the fractional deviation between the EHT and the dynamical inferences of the mass-todistance ratio, when we assume the prior measurements based on stellar and gas dynamics.We find δ=−0.01±0.17(68% credible intervals) for the stellar and δ=0.78±0.3 for the gas dynamics priors.The fact that our measurement θ g is consistent with one of the prior measurements θ dyn allows us to conclude that our null hypothesis has not been violated.
The second property of the black hole shadow that is affected by potential deviations of its spacetime from the Kerr metric is its shape.Shadows of Kerr black holes are nearly circular with a maximum asymmetry of 10%, as defined in Johannsen & Psaltis (2010).For the case of the ∼17°i nclination of M87 (and assuming that the large-scale jet points along the spin axis), the maximum asymmetry is 2% for all values of black hole spin (e.g., Chan et al. 2013).Measuring a shadow shape that deviates from circular to a degree larger than this amount would require that the spacetime of the compact object is not described by the Kerr metric (Johannsen & Psaltis 2010).
Figure 18 shows that the circularity measured from reconstructed images of M87 is consistent with that found in model images from GRMHD simulations in the Kerr metric.The images used in our analysis in Section 7 measure the shape of the emitting region.In future work a similar analysis can be done using images from models that violate the no-hair theorem, where the photon ring/shadow shape can be highly prolate, oblate, or even amorphous (Broderick et al. 2014;Giddings & Psaltis 2018).

Conclusions
The horizon-scale emission of M87 at 1.3 mm exhibits a robust crescent-like structure.This is seen in image reconstructions (Paper IV) and supported by features in the visibility data from the 2017 EHT campaign (Paper III).Geometric crescent models are overwhelmingly preferred to similarly complex models that we have explored that do not have a central flux depression.
The crescent structure is well defined, with a rapid decline to a dark interior with a brightness 10× lower than the average of that in the annular region.The average diameter of the crescent is well constrained by multiple methods with a combined range of 42 ± 3 μas.This is consistent among all observation days and bands.
Based on geometric models, GRMHD simulations, and image domain feature extraction, the angular size of the gravitational radius is θ g =3.8±0.4μas.The uncertainty of this value is currently dominated by a systematic component associated with the location of the mean emission diameter in GRMHD simulations that is used to convert from an angular to physical scale.Accounting for this uncertainty explicitly, the resulting black hole mass is . This mass is consistent with that measured from stellar but not gas dynamics.
The crescent morphology, rapid drop to a deep interior flux depression, and broad consistency among days, methods, and the stellar dynamics measurement all point to the emission structure from M87 being due to strong gravitational lensing around a central black hole.The consistency between the 10 pc-scale stellar dynamics and EHT mass measurements provides a null test of the Kerr metric.The size constraint and qualitative support of the standard jet formation paradigm argues for the presence of a horizon.Together, our results strongly support the hypothesis that the central object in M87 is indeed a Kerr black hole, and provide new evidence for the long-believed connection between AGNs and SMBHs.
The authors of this Letter thank the following organizations and programs: the Academy of Finland (projects 274477, 284495, 312496); the Advanced European Network of E-infrastructures for Astronomy with the SKA (AENEAS) project, supported by the European Commission Framework Programme Horizon 2020 Research and Innovation action under grant agreement 731016; the Alexander von Humboldt Posterior in the fractional difference δ between the inferred angular gravitational radius of the black hole in M87 based on the EHT data and the values measured using stellar and gas dynamics.The shaded region corresponds to the maximum ±4% discrepancy between the measurements that can be due to the weak dependence of the black hole shadow on the spin magnitude and orientation of a Kerr black hole.The fact that at least one prior mass measurement is consistent with the EHT measurement reported here provides a null hypothesis test of the Kerr nature of the central black hole.
Stiftung; the Black Hole Initiative at Harvard University, through a grant (60477) from the John Templeton Foundation; for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada.The LMT project is a joint effort of 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.Partial SPT support is provided by the NSF Physics Frontier Center award (PHY-0114422) to the Kavli Institute of Cosmological Physics at the University of Chicago (USA), the Kavli Foundation, and the GBMF (GBMF-947).The SPT hydrogen maser was provided on loan from the GLT, courtesy of ASIAA.The SPT is supported by the National Science Foundation through grant PLR-1248097.Partial support is also provided by the NSF Physics Frontier Center grant PHY-1125897 to the Kavli Institute of Cosmological Physics at the University of Chicago, the Kavli Foundation and the Gordon and Betty Moore Foundation grant GBMF 947.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, both 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 acknowledge the significance that Maunakea, where the SMA and JCMT EHT stations are located, has for the indigenous Hawaiian people.

Appendix A THEMIS Station Gain Amplitude Reconstruction
Both THEMIS and eht-imaging provide tools for analyzing EHT data within the visibility domain.Both address the reconstruction of individual station gains, though in substantially different ways.Here, we briefly describe the procedure employed by THEMIS and compare the reconstructed gains with the results reported in Paper IV.
We reconstruct gains on a scan-by-scan basis.Thus, for EHT 2017 data, the gain amplitudes constitute between 40 and 143 additional nuisance parameters per data set; more if sub-scan variations are necessary.For the GC models this represents a proliferation of parameters by 200%-700%.
Even with efficient sampling techniques, treating these identically with the other model parameters would result in a substantial (super-linear) increase in the computation expense needed to adequately sample the likelihood.Thus, in THEMIS, these are efficiently addressed instead by directly marginalizing the likelihood in Equation (9).This is done by numerically maximizing  A ij , over the | | g i (subject to Gaussian priors on the | | g i ) for each scan and then marginalizing over an expansion around this maximum (see A. E. Broderick et al. 2019, in preparation, for more details).Within this procedure it is assumed that gains from neighboring scans are independent.For all analyses presented here, the standard deviation of these priors for the LMT and all other stations are 100% and 20%, respectively.While the latter substantially exceeds the anticipated uncertainty stated in Table 3 of Paper III, the reconstructed station gains are typically very well constrained and close to unity.
Typically, the station gains on a given scan are significantly overdetermined.Thus, calibrating the station gains in this way produces model-independent patterns, driven primarily by the need for consistency among a given observation.
In Figure 22 we show the reconstructed LMT gains for the best-fit GC model model on each day (low-band only) in comparison to those obtained by image-domain methods for the LMT.This station was selected because it exhibits large gain excursions; the other reconstructed station gains are generally much closer to unity.
The two different procedures based on different intrinsic source models find remarkable agreement over all days.Importantly, there are no large gain offsets that would produce a significant impact on the reconstruction of the total compact flux.

Appendix B Generalized Crescent Models
In this appendix we describe the GC models used in Section 5 to fit the M87 data.Descriptions of the model parameters and their associated priors are listed in Table 8.

B.1. Fourier-domain Construction of the GC models
Following Kamruddin & Dexter (2013), we start with a unit circular "top-hat" function circ(r), defined in the image domain to be A disk of radius R can then be expressed as D(r)=circ(r/R).
The Fourier transform of such a disk is given by is the Fourier-domain radial coordinate and J 1 (x) is the Bessel function of the first kind of order 1.We can construct a circular crescent C via the difference of two such disks after shifting the smaller, inner disk by r 0 in the negative x-direction prior to subtracting it from the outer disk.A spatial shift in the image domain corresponds to a frequency shift in the Fourier domain, such that r pr r pr where R out and R in are the outer and inner disk radii, respectively.Similar to Benkevitch et al. (2016), we define a "slash" operation s(x) that imposes a linear brightness gradient across the image.The slash operator can be expressed in the image domain as where h 1 and h 2 are the amplitudes of the function at x=R out and x=−R out , respectively.We apply s(x) to the crescent via multiplication in the image domain, which corresponds to taking a derivative in the Fourier domain Application of the slash is aided by the fact that the derivative of a Bessel function can be expressed in terms of other Bessel functions and we denote the slashed crescent model as sC.

B.1.1. xs-ring Model
We now seek to add an emission "floor," F, to the center of sC to account for nonzero central emission.In the xs-ring construction, F takes the form of a circular disk with amplitude K and radius R in that is shifted by r 0 along the x-axis, In the image domain, we can express the slashed crescent plus emission floor [sC+F] as The model thus corresponds to the difference of two spatially offset circular disks having a linear brightness gradient across them, to which a central emission floor is then added.We then further allow the entire model to be rotated by some PA f, which amounts to replacing cos .Finally, we convolve the model in the image domain with a circular Gaussian kernel of width σ and then scale it to have a total flux density of V 0 .The compact component of the xs-ring model can thus be written in the  ´-- and we have followed Kamruddin & Dexter (2013) The xs-ring model has eight free parameters describing the compact emission: V 0 , R out , f, σ, ψ, τ, γ, and β.The ψ, τ, and γ parameters are restricted to lie between 0 and 1.The ψ parameter describes the typical fractional thickness of the crescent; as y  1 the crescent reduces to a filled disk, and as y  0 the crescent reduces to an infinitesimally thin ring (i.e., delta function in radius).The τ parameter is a measure of the crescent's structural symmetry; as t  1 the crescent reduces to a circular ring, and t  0 when the inner and outer disks meet at one edge.The γ parameter describes the amplitude of the emission floor, with γ=0 indicating no central emission.The β parameter describes the strength of the emission gradient across the crescent, with β=1 indicating no gradient.
In addition to the compact component, we add a large-scale circular Gaussian, G, to the model This component is permitted to have very large widths, σ G ?1 arcsec, to account for short-spacing flux that might be missed by the bulk of the array.The free parameters in G are its flux density (V G ) and width (σ G ), bringing the total number of parameters for the xs-ring model to 10.An example of the xsring model is shown in the left panel of Figure 3.

B.1.2. xs-ringauss Model
The xs-ringauss model differs from the xs-ring model in two key respects.The first difference is that the xs-ringauss model follows Benkevitch et al. (2016) and includes a fixed elliptical Gaussian component on top of the crescent.This component is centered at (x, y)=(r 0 −R in , 0) prior to applying the PA rotation, and its orientation is fixed to match that of the crescent.The additional free parameters provided by this Gaussian component are thus its widths (σ x , σ y ) and flux density V 1 .The second difference is that the emission floor F in the xs-ringauss model takes the form of a circular Gaussian rather than a disk, In addition to the amplitude parameter K, this form of the emission floor has a width parameter σ F that also enters into the model.The remaining parameters in the xs-ringauss model match those in the xs-ring model, for a total of 14 free parameters: 10 of these parameters are shared with the xs-ring model, and the additional four are V 1 , σ x , σ y , and σ g .An example of the xsringauss model is shown in the right panel of Figure 3.

B.2. Nuisance Gaussian Components
In addition to the crescent components of the GC models, we also fit a small number (two to three) of additional "nuisance" Gaussian components (see Section 5.1).These additional components are intended primarily to capture extraneous amorphous emission surrounding the main ring (see, e.g., Figure 12) and secondarily to provide additional flexible degrees of freedom that allow the model to absorb systematic uncertainties.

B.2.1. Parameterization
We parameterize the Gaussian components in one of two ways.The first of these, used when fitting the xs-ring model, is given by x y 2 2 and Here, V g,1 is the total flux density of the Gaussian, (x 0 , y 0 ) are its central coordinates, (σ x , σ y ) are the Gaussian widths along the two principal axes, and θ is the PA.A circular Gaussian model is constructed by setting σ x =σ y and θ=0.
When fitting multiple Gaussian components in the context of the xs-ring model, we only parameterize the first component, as described above, with a central coordinate position of (x 0,0 , y 0,0 ).All additional components have positions that are referenced to that of the previous component, such that . Such a parameterization enables the imposition of a spatial ordering by treating the Δx or Δy parameters as slack variables, thereby breaking the multimodal labeling degeneracy otherwise caused by pairwise swaps of different Gaussian components.
The second Gaussian parameterization, used by the xsringauss model, follows the form described in Broderick et al. (2011).Specifically, where (ρ, ψ) are cylindrical polar coordinates in the (u, v)plane.Here, is a width parameter and is an anisotropy parameter.For small numbers (5) of Gaussian components, tempering is sufficient to ensure that the MCMC exploration can navigate all posterior modes; because the xs-ringauss model never utilizes more than two additional Gaussian components, no mode collapsing is necessary.

B.2.2. Behavior in Fits to the M87 Data
In principle, it is possible to add a large number of nuisance Gaussian components to a model and thereby to fit the data arbitrarily well.However, our primary goal with model fitting is not to exactly reproduce all details of the emission structure, but rather to make robust measurements of the crescent component model parameters relevant for constraining the lensed part of the image.Once a sufficient number of Gaussian components have been added to the model, further additions do not substantially modify the derived crescent parameters (see Figure 23).We find that two nuisance Gaussians for the xsringauss GC model (for a total of 26 free parameters), and three for the xs-ring GC model (for a total of 28 free parameters), are the threshold values that most generally satisfy this criterion.These values are thus used for all GC model fits presented in this Letter.
When fit alongside the GC models to the M87 data, we find that the nuisance Gaussian components account for ∼10%-70% of the compact flux density.Assessing consistency in the best-fit parameter values for the nuisance Gaussian components across models is complicated by the differing model We see that after a small number of components is reached (three for this data set, marked in green), adding additional components does not substantially modify the parameter posteriors.
specifications.In general we find less consistency in the properties of nuisance Gaussian components for a given day across models (within a single data set) than we do across bands (within a single model).
Within the context of a particular GC model, we find that the nuisance Gaussian components approximately fall into two categories: (1) components that reside near or on the crescent component and whose location and structure are consistent across bands within a single day, and (2) components that tend to reside far (40 μas) outside the crescent component and whose location and structure may vary across bands and days.We associate the first class of nuisance Gaussian components with attempts by the model to account for a real flux distribution that is too complex for the crescent component alone.The second class of nuisance Gaussians acts in a similar manner to the artifacts frequently seen in image reconstructions (see Paper IV).We thus associate this second class of components with the image-domain manifestation of systematic errors in the data products.We emphasize, however, that we have no mechanism for enforcing a categorization such as that described on the nuisance Gaussian components, and in general we expect that any single component will be influenced both by real source emission and by systematic uncertainties.

Appendix C χ 2 Statistics
In this appendix we define the various χ 2 statistics used as diagnostics of model fit quality in Section 5. We note that the expressions presented here do not have any bearing on the fitting process itself; the likelihood functions used during fitting are described in Section 4.1.
The visibility amplitude reduced-χ 2 is given by where N A is the number of visibility amplitudes (see Table 1), N p is the number of GC model parameters, N g is the number of independent gain terms, σ is the uncertainty in visibility amplitude, and Â and A are the modeled and measured visibility amplitudes, respectively.For the xs-ring GC model fits using dynesty we have N p =28 model parameters, while for the xs-ringauss GC model fits using THEMIS we have N p =26.Analogous to Equation (54), the closure phase reduced-χ 2 is where y N C is the number of closure phases, s y C is the closure phase uncertainty, and y ˆC and ψ C are the modeled and measured closure phases, respectively.As in Section 4.1, we avoid phase wrapping by selecting a branch of the closure phase space such that differences always fall between −180°a nd 180°.
For the logarithmic closure amplitude reduced-χ 2 , we have where N A ln C is the number of logarithmic closure amplitudes, s A ln C is the logarithmic closure amplitude uncertainty, and Â ln C and A ln C are the modeled and measured logarithmic closure amplitudes, respectively.
For the THEMIS fits, we report the joint visibility amplitude and closure phase reduced-χ 2 ´-+ - ´-+ - The values for each of these statistics are reported in Table 2 for both fitting codes and for all data sets.

Appendix D Using GRMHD Simulations to Calibrate the Generalized Crescent Model
To create a mapping between the GC model diameter and the angular gravitational radius θ g , we require data sets for which we know both quantities.To this end, we take the GRMHD simulation library described in Paper V to provide our "ground truth" images.These simulations are necessarily limited in the range and resolution of physical parameter space that they probe, but they critically provide a large sample of plausible emission structures for which we know the underlying physical parameters.By treating these GRMHD simulations as "ground truth" images, generating synthetic data from them, and then fitting that data with the GC models, we aim to (1) determine a calibration of the measured crescent diameter d in terms of θ g , and (2) estimate the magnitude of the uncertainty in this calibration.

D.1. GRMHD Image Selection and Synthetic Data Generation
For a given GRMHD simulation, individual "snapshots" are ray-traced to produce images of the 230 GHz emission as it would be seen from Earth.These images are then used as input for generating synthetic observational data with ehtimaging, which generates visibilities corresponding to the (u, v)-coverage from the 2017 EHT observations of M87 and corrupts this data with realistic levels of thermal noise and station-based systematics.A detailed description of the synthetic data generation pipeline is provided in Paper IV.
From the GRMHD simulation library we selected 103 separate snapshot images with which to generate synthetic data sets.These images were selected randomly from within the parameter space explored by the library, and thus included examples of emission structures that were incompatible with the M87 data.The images were split into two groups, from which were generated two sets of synthetic data targeting different aspects of the calibration.
The first group of synthetic data sets, Part I, was created using only three simulated GRMHD images (shown in the top three panels of Figure 24).From each of these images 10 synthetic data sets were generated, corresponding to 10 different realizations of the observational noise (i.e., thermal noise and known systematics such as station gain fluctuations and polarization leakage) spanning the full range of (u, v)coverage in the 2017 EHT M87 data and containing independent instances of thermal and systematic errors.The goal of the Part I data sets is to assess the impact of such observational corruptions on the recovered model parameters, which enters into the error budget for the calibration.
The second group of synthetic data sets, Part II, contains the remaining 100 GRMHD images (samples are shown in the bottom four panels of Figure 24).Only a single synthetic data set was generated for each of these images, using a randomly selected realization of the (u, v)-coverage from among the four days of 2017 M87 observations.The goal of the Part II data sets is to provide a sample of images for which we know both the underlying values of θ g and the best-fit values of the diameter d from the GC model fits.

D.2. θ g Calibration
We know the input value of the angular gravitational radius θ g for each simulation in Parts I and II of the GRMHD calibration data set.The basic calibration task is then to fit each of the corresponding synthetic data sets with a GC model and determine the scaling factor α (see Equation ( 27)) relating the best-fit crescent diameter to θ g for the underlying simulation.
The value of this scaling factor is equal to the diameter measured in units of θ g , which Figure 25 shows plotted for all simulations in the calibration data set.We obtain an estimate of the scaling factor from each individual fit in Parts I and II, and we use the mean of these fits as our final determination of α; the uncertainty in this mean value is roughly an order of magnitude smaller than any of the calibration uncertainties we consider below, and we thus neglect it in our error budget.Table 4 lists the calibrated α values, which we determine separately for THEMIS and dynesty.
We can see in Figure 25 that each set of 10 synthetic data sets corresponding to a single simulation from Part I shows excellent agreement across all measured d values.Comparing these to the measurements from Part II demonstrates that the inter-simulation scatter is substantially larger than the interobservation scatter, and that the former will be the limiting factor in measuring θ g from the M87 data.

D.3. Uncertainty Budgeting for All Parameters
In general, the error budget for any single measurement of θ g will have three contributions: (1) a "statistical uncertainty" associated with the width of the posterior, (2) an "observational uncertainty" associated with incomplete (u, v)-coverage and unmodeled systematics that change from one observation to the next, and (3) a "theoretical uncertainty" associated with the data being a single sample from a dynamic system whose specific physical properties we are largely ignorant of.For other measured parameters that lack a counterpart input to the GRMHD simulations, we only consider the first two uncertainty components.To quantify the magnitude of these different uncertainty components, we employ a generalized lambda distribution (GλD).The GλD is an extension of the Tukey λ distribution (Tukey 1962) that provides a compact parameterization for representing a diverse family of probability densities.The GλD is also flexible and convenient, permitting highly asymmetric and heavy-tailed distributions to be analytically represented in terms of a function that depends only on their quantiles.We use the FMKL parameterization (Freimer et al. 1988), defined in terms of the inverse cumulative distribution function (CDF), which can be expressed using four parameters (λ 1 , λ 2 , λ 3 , λ 4 ) as where q is the quantile of the distribution.The fits were performed using the GLDEX package in R (Su 2007a(Su , 2007b)).The GλD fits to all three components of the error budget are shown in Figure 26 for the θ g measurements from each of the two fitting codes.We find that the theoretical uncertainty is dominant, being larger by a factor of ∼4-5 than either the statistical or observational components.However, the magnitude of this theoretical uncertainty depends strongly on the category of simulation being used for calibration.We find that SANE simulations show nearly a factor of ∼2 larger theoretical uncertainty than MAD simulations; the increased SANE scatter is visually apparent in Figure 25, and the effective decrease in calibration uncertainty that would result from using only MAD simulations is listed in Table 4.However, because such a preference for MAD over SANE is not well motivated from a theoretical standpoint (Paper V), we use the calibration factor determined from MAD+SANE for the measurements presented in this Letter.
The output of a MCMC or NS fitting run is a "chain" of samples from the posterior distribution.We determined the statistical uncertainty component for each parameter by fitting a single GλD to the ensemble of chains from all 130 Part I and II data sets simultaneously.Each chain was first mean-subtracted so that the GλD characterizes only the shape of the typical posterior distribution for a single GC model fit, rather than capturing variation in mean value between different fits.
The observational uncertainty component was determined using Part I data sets only.For each of the three simulations in Part I, we combined the 10 chains corresponding to the individual realizations of observational uncertainty for that simulation.We then subtracted out the mean of each combined chain before further combining all three chains and fitting a GλD to the resulting stack.This fitted distribution then  4).
Figure 26.Sources and magnitudes random and systematic uncertainty on the measurement of θ g with GC models.These include the average posterior (blue), the impact of different realizations of unmodeled observational systematic errors (e.g., polarization leakage, (u, v)-coverage, etc.; green), and the impact of variations in the assumed underlying GRMHD simulation (red).(The latter two match the appropriately colored ranges shown in Figure 25.)These are similar for analyses that employ the xs-ringauss-based (THEMIS; top) and xs-ring-based (dynesty; bottom) GC models.The gray band shows the range of photon ring diameters for a Kerr black hole from spin 0 to 0.998 over all viewing inclination angles.
Figure 27.Illustration of "good" (top panels) and two "bad" (middle and bottom panels) models following the THEMIS-AIS procedure.Closure phases are shown for three triangles at 5 UTC for April 5 as proxies for the χ 2 .Blue dots indicate the closure phases for the best-fit SSM associated with each simulation snapshot image.The red diamond shows the same for the average snapshot image.The green triangle shows the observed values on April 5, high-band, near 5 UTC.The manner in which the two models shown in the middle and bottom rows are excluded differs: in the middle case the reduced χ 2 is too large, while in the bottom case the χ 2 is too small to be consistent with that anticipated by the individual GRMHD model snapshots.
which was not accounted for in either of the two measurements.A second one, considered in Blakeslee et al. (2009) but not in Cantiello et al. (2018a), arise from the intrinsic scatter in the absolute fluctuation magnitude attributed to stellar populations (cosmic scatter; Tonry et al. 2000) and is estimated to be 0.06 mag (Blakeslee et al. 2009).Adding these in quadrature to the formal errors quoted above, we obtain for the absolute distant moduli μ SBF1 =31.11±0.13(Blakeslee et al. 2009) and μ SBF2 =31.15±0.12(Cantiello et al. 2018a).Finally, we average the two measurements weighted by their errors to get μ=31.14±0.12.We report these measurements in Table 9, together with the combined posterior as their product in Figure 31.
Assuming Gaussian posteriors P(μ) for the distance modulus measurements with the means and standard deviations given in Table 9, we calculate the posteriors for the distance Finally, considering the three distance measurements to be independent, we calculate the combined posterior as their product.We show the posteriors in Figure 31 and quote their most likely values and uncertainties in Table 9.

I.2. Mass Determination
Both the stellar dynamical and the gas dynamical mass measurement techniques that we described in Section 8.2 measure, in reality, a mass-to-distance ratio rather than the mass of the central object directly.To convert these into posteriors over a mass, we first convert the measurements into posteriors over the mass-to-distance ratio  and then rescale them based on the posterior over distance that we described in the previous section.
We use the marginal χ 2 distributions for mass M BH shown in Figure 6 of Gebhardt et al. (2011) and Figure 6 of Walsh et al. (2013) and set M BH =×(17.9Mpc).We denote the posteriors for the stellar and gas dynamical measurements as  ( ) for each model (i=gas, stars).Table 9 summarizes the credible regions of these distributions.We can also convert the mass-to-distance ratio from dynamical measurements directly to an angular size of one gravitational radius at the distance of M87, θ dyn ≡G/c 2 .This posterior is given by We list the credible intervals for this quantity in Table 9.
We can now use the posterior over the mass-to-distance ratio and that over distance to calculate the posterior over the mass for the stellar and gas dynamics measurements as  Note.We report median values and 68% confidence intervals.
Figure 31.Normalized posteriors for the distance measurements to M87.Two colored lines show the posterior from the measurements using the TRGB method (Bird et al. 2010) and combined posterior from the SBF method (Blakeslee et al. 2009;Cantiello et al. 2018a).The black line shows the normalized product of three posteriors obtained by combining these measurements.

Figure 1 .
Figure1.(u, v)-coverage (left panel) and visibility amplitudes (right panel) of M87 for the high-band April 11 data.The (u, v)-coverage has two primary orientations, east-west in blue and north-south in red, with two diagonal fillers at large baselines in green and black.Note that the Large Millimeter Telescope (LMT) and the Submillimeter Telescope (SMT) participate in both orientations, and that the LMT amplitudes are subject to significant gain errors.There is evidence for substantial depressions in the visibility amplitudes at ∼3.4 Gλ and ∼8.3 Gλ.The various lines in the right panel show the expected behavior of (dotted line) a Gaussian, (dashed line) a filled disk, and (green area) a crescent shape along different orientations.The image of M87 does not appear to be consistent with a filled disk or a Gaussian.
and also Paper III).The right panel of this figure shows the visibility amplitudes observed on April 11 color coded by the orientation of the baselines.
Section 2.2, a non-redundant subset of closure phases can be selected to maximize independence and thus ensure the near-diagonality of S y .In this case the closure phase measurements can be treated as individually Gaussian,

Figure 2 .
Figure2.Relative log-likelihood values for different geometric models fit to the M87 data as a function of nominal model complexity; the number of parameters is given in parenthesis for each model.April 5 is shown here, and all days and bands show the same trend.The models shown in this figure are strict subsets of the "generalized crescent model" (labeled here as model "n"; see Section 5.1), and they have been normalized such that the generalized crescent model has a value of  = 1; the reduced-χ 2 for the generalized crescent fit is 1.24 (see Table2).We find that the data overwhelmingly prefer crescent models over, e.g., symmetric disk and ring models, and that additional Gaussian components lead to further substantial improvement.Note that a difference of ∼5 on the vertical axis in this plot is statistically significant.

Figure 3 .
Figure 3. Schematic diagrams illustrating the crescent components of the xs-ring (left panel) and xs-ringauss (right panel) models.Dashed lines outline the inner and outer circular disk components that are differenced to produce the crescent models, and for the xs-ringauss model the FWHM of the fixed Gaussian component is additionally traced as a dotted line.The red and green curves above and to the right of each panel show cross-sectional plots of the intensity through the corresponding horizontal and vertical slices overlaid on the images.The circular and square markers indicate the centers of the outer and inner disks, respectively.The labeled parameters correspond to those described in Section 5.1.Both crescents are shown at an orientation of f f = =  ˆ90 .

Figure 4 .
Figure 4. Modeled data (top panels) and residuals (bottom panels) for GC model fits to the April 6 high-band data set, with the data plotted in gray; we show results for the median posterior fit.The panels show visibility amplitude (left panels), closure phase (middle panels), and logarithmic closure amplitude (right panels) data.The xs-ringauss model, shown in red, is fit to the visibility amplitudes and closure phases using THEMIS; the dynesty-based xs-ring model, plotted in blue, comes from a fit to closure phases and logarithmic closure amplitude.Because both models fit to closure phases, the center panel shows two sets of models and residuals.All residuals are normalized by the associated observational noise values.

Figure 5 .
Figure 5. Image domain representations of a random posterior sample from the xs-ring (left panel) and xs-ringauss (right panel) model fits to the April 6 high-band data set; note that these are representative images drawn from the posteriors, and thus do not represent maximum likelihood or other "best-fit" equivalents.The xs-ring model fit uses only closure quantities, so we have scaled the total flux density to be equal to the 1.0 Jy flux density of the xs-ringauss model fit.

2 C 2 C 2 C
), with the joint visibility amplitude and closure phase reduced-χ 2 denoted as c y + A .The xs-ring values are from fits to closure phases and logarithmic closure amplitudes (reduced-χ 2 given by c A ln ), with the joint closure phases and logarithmic closure amplitude reduced-χ 2 denoted as c y + A

Figure 6 .
Figure 6.Joint posteriors for the key physical parameters derived from GC model fits to the April 5, high-band data set.Blue contours (upper-right triangle) show xsring posteriors obtained from the dynesty-based fitting scheme, while red contours (lower-left triangle) show xs-ringauss posteriors obtained using THEMIS.Contours enclose 68% and 95% of the posterior probability.Note that recovery of the total compact flux density ĈF is not possible for the dynesty-based fits, which only make use of closure quantities.

Figure 7 .
Figure 7. Posterior medians and 68% confidence intervals for selected parameters derived from GC model fitting for all observing days and bands.Blue circular points indicate xs-ring fits using the dynesty-based fitting scheme applied to individual data sets (i.e., a single band on a single day).Red square points indicate xs-ringauss fits using THEMIS applied to individual data sets, while orange square points show THEMIS-based xs-ringauss fits to data sets that have been band-combined.All plotted error bars include the systematic "observational uncertainties" estimated from simulated data in Appendix D; these uncertainties are listed in the bottom row of Table3.Note that recovery of the total compact flux density ĈF is not possible for the dynesty-based fits, which use only closure quantities.The light purple band in the lower-right panel is the range inferred in Paper IV.

Figure 8 .
Figure 8. Constraints on θ g arising from the GRMHD simulation calibrated GC model fits, by day and band.Solid error bars indicate 68% confidence intervals, while dashed error bars indicate the systematic uncertainty in the calibration procedure.Circular blue points indicate independent analyses for each band and on each day in the context of the xs-ring GC model fit using the dynestybased method.Square points indicate independent analyses for each band (red) and band-combined analyses (orange) for each day using the xs-ringauss GC model with THEMIS.Colored bands around dashed lines (right) indicate the combined constraint across both bands and all days, neglecting the systematic uncertainty in the calibration procedure.

Figure 9 .
Figure 9. Illustration of the parameters of the SSM described in Section 6.2.Both the original GRMHD simulation (left) and the corresponding SSM for an arbitrary set of parameter values, (flux rescaling, stretching of the image, and rotation; right) are shown.In both panels, the gray arrow indicates the orientation of the forward jet.

Figure 11 .
Figure11.Visibility amplitude and closure phase residuals for an SSM fit to the April 5, high-band data for a "good" snapshot image frame from a MAD simulation with a * =0, i=167°, and R high =160.The reduced-χ 2 values for the fits are 5.9 (THEMIS) and 7.3 (GENA).All residuals are normalized by their corresponding estimated observational errors.

Figure 12 .
Figure 12.Three sample SSM fits to April 6 high-band data.All are from models with THEMIS-AIS p AIS >0.1.In panels (a) and (b), the prominent photon ring places a strong constraint on θ g .In panel (c), the extended disk emission results in a smaller θ g estimate.

Figure 13 .
Figure 13.Constraints on θ g arising from the GRMHD model fitting procedure, by day and band.The maroon squares and green triangles are the constraints arising form the THEMIS and GENA pipelines.Solid error bars indicate the 68% confidence levels about the median.The maroon colored band indicates the combined constraint across both bands and all days.

Figure 15 .
Figure 15.Image domain measurements of the ring diameter of M87 over all observing days comparing three image reconstruction methods.The error bars show the full range of results for the Top Set images using low-band data.The gray dashed line and band show the weighted average and uncertainty across methods and observing days.

Figure 16 .
Figure16.Diameters and fractional widths inferred from image (black and green) and visibility (red and blue) domain measurements on April 6.The visibility domain measurements are from GC model fitting (see Section 5), while the reconstructed images are from Paper IV.The filled regions show diameter and width anti-correlations expected in simple ring models (Section 7.2).The anti-correlation between mean diameter and width helps to explain the small (;5%) offset in mean diameter found between methods.

Figure 17 .
Figure 17.Constraints on θ g from GRMHD simulation calibrated image domain feature extraction, by day and band.The solid lines show the full range of diameter measurements from the Top Set images.The dotted lines show the systematic calibration uncertainty.The shaded regions show the weighted average and uncertainty over observing days for each method.
Walsh et al. (2013) used Hubble Space Telescope (HST) Space Telescope Imaging Spectrograph observations of Hα and [N II] emission lines from the central gas-disk within the M87 nucleus to map its kinematic field within 40 pc of the black hole.The velocity field implies a compact mass of -

Figure 19 .
Figure19.Estimates of the mass of the central black hole in M87 for the three different measurement techniques employed in this Letter (see alsoTable 7).Gray horizontal bands around dashed lines correspond to the 68% confidence levels for the stellar (top) and gas (bottom) dynamical mass measurements (see Table9).

Figure 21 .
Figure 21.Posterior in the fractional difference δ between the inferred angular gravitational radius of the black hole in M87 based on the EHT data and the values measured using stellar and gas dynamics.The shaded region corresponds to the maximum ±4% discrepancy between the measurements that can be due to the weak dependence of the black hole shadow on the spin magnitude and orientation of a Kerr black hole.The fact that at least one prior mass measurement is consistent with the EHT measurement reported here provides a null hypothesis test of the Kerr nature of the central black hole.

Figure 22 .
Figure 22.Reconstructed LMT station gains by THEMIS from GC model fits to the April 5, 6, 10, and 11 low-band data sets.All other station gains are very close to unity, and are consistent with those from image reconstruction.The inferred LMT gains are consistent with those found in Paper IV across all days.

Figure 23 .
Figure23.Measured crescent parameter values from fitting the xs-ring model with an increasing number of nuisance Gaussian components, shown for the April 5 high-band data set.We see that after a small number of components is reached (three for this data set, marked in green), adding additional components does not substantially modify the parameter posteriors.

Figure 24 .
Figure 24.Examples of GRMHD snapshots used for calibration (see Appendix D).The top three panels show the images used to calibrate the "observational uncertainty" associated with changing (u, v)-sampling and systematics across days (Part I).The bottom four panels show example images taken from the pool of 100 that were used to calibrate the "theoretical uncertainty" caused by different physical parameters and turbulence realizations in the GRMHD simulations (Part II).40 μas scale bars are shown in white.

Figure 25 .
Figure 25.Reconstructed diameters using the GC models from synthetic data generated from GRMHD simulations.Green bands show the one-sided 68-percentile ranges within the Part I reconstructions (left).The red bands show the one-sided 68-percentile range of the combined Part I and II reconstructions (right).These are separated by GRMHD simulation type: MAD models are collected in white panels, while SANE models are collected in gray panels.The dotted lines indicate the MAD+SANE calibration factor, α for the xs-ringauss model (see Table4).

Table 1
The Number of Data Product and Gain Terms

Table 2
Reduced-chi2 Statistics for the GC Model Fits Note.Reduced-χ 2 values corresponding to the maximum likelihood posterior sample for individual M87 data sets from both fitting codes, split by data type and calculated as described in Appendix C. The xs-ringauss values are from fits to visibility amplitudes (reduced-χ 2 given by χ A2 ) and closure phases (reduced-χ 2

Table 3
Best-fit GC Model Parameters for All Data Sets and Both Models, as Defined in Section 5.1; Median Posterior Values are Quoted with 68% Confidence Intervals LNote.Observational uncertainties are listed at the bottom of the table and have been determined as described in Appendix D.3.Note that recovery of the total compact flux density ĈF is not possible for the dynesty-based fits, which only make use of closure quantities.

Table 4
Calibrated Scaling Factors, α, and Corresponding θ g Measurements

Table 5
Best-fit Parameters for All Data Sets from Direct GRMHD Simulation Fitting, as Defined in Section 6.2 SMILI;Akiyama et al.Akiyama et al. , 2017b) )17b) )and one CLEAN method (DIFMAP; Shepherd 1997).Here we present a first analysis of ring properties seen in reconstructed M87 images.

Table 6
Measured Diameters, d, Calibrated Scaling Factors, α, andCorresponding Gravitational Radii, θ g , for the Image Domain Analysis Presented in Section 7

Table 7
Summary of θ g and M Measurements Note.Measurements made in this work (top) and from the literature (middle) are quoted as median values with 68% confidence intervals.The final combined measurements from this work are listed in bottom row of the table, rounded to two significant figures.The distance used to compute M from θ g is  16.8 0.8 Mpc (see Section 8).

Table 8
Parameters and Priors for the GC Models

Table 9
Relevant Literature and Derived Quantities