First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole

When surrounded by a transparent emission region, black holes are expected to reveal a dark shadow caused by gravitational light bending and photon capture at the event horizon. To image and study this phenomenon, we have assembled the Event Horizon Telescope, a global very long baseline interferometry array observing at a wavelength of 1.3 mm. This allows us to reconstruct event-horizon-scale images of the supermassive black hole candidate in the center of the giant elliptical galaxy M87. We have resolved the central compact radio source as an asymmetric bright emission ring with a diameter of 42+/-3 micro-as, which is circular and encompasses a central depression in brightness with a flux ratio ~10:1. The emission ring is recovered using different calibration and imaging schemes, with its diameter and width remaining stable over four different observations carried out in different days. Overall, the observed image is consistent with expectations for the shadow of a Kerr black hole as predicted by general relativity. The asymmetry in brightness in the ring can be explained in terms of relativistic beaming of the emission from a plasma rotating close to the speed of light around a black hole. We compare our images to an extensive library of ray-traced general-relativistic magnetohydrodynamic simulations of black holes and derive a central mass of M = (6.5+/-0.7) x 10^9 Msun. Our radio-wave observations thus provide powerful evidence for the presence of supermassive black holes in centers of galaxies and as the central engines of active galactic nuclei. They also present a new tool to explore gravity in its most extreme limit and on a mass scale that was so far not accessible.


Introduction
Black holes are a fundamental prediction of the theory of general relativity (GR; Einstein 1915).A defining feature of black holes is their event horizon, a one-way causal boundary in spacetime from which not even light can escape (Schwarzschild 1916).The production of black holes is generic in GR (Penrose 1965), and more than a century after Schwarzschild, they remain at the heart of fundamental questions in unifying GR with quantum physics (Hawking 1976;Giddings 2017).
Black holes are common in astrophysics and are found over a wide range of masses.Evidence for stellar-mass black holes comes from X-ray (Webster & Murdin 1972;Remillard & McClintock 2006) and gravitational-wave measurements (Abbott et al. 2016).Supermassive black holes, with masses from millions to tens of billions of solar masses, are thought to exist in the centers of nearly all galaxies(Lynden-Bell 1969; Kormendy & Richstone 1995;Miyoshi et al. 1995), including in the Galactic center (Eckart & Genzel 1997;Ghez et al. 1998;Gravity Collaboration et al. 2018a) and in the nucleus of the nearby elliptical galaxy M87 (Gebhardt et al. 2011;Walsh et al. 2013).
Active galactic nuclei (AGNs) are central bright regions that can outshine the entire stellar population of their host galaxy.Some of these objects, quasars, are the most luminous steady sources in the universe (Schmidt 1963;Sanders et al. 1989) and are thought to be powered by supermassive black holes accreting matter at very high rates through a geometrically thin, optically thick accretion disk (Shakura & Sunyaev 1973;Sun & Malkan 1989).In contrast, most AGNs in the local universe, including the Galactic center and M87, are associated with supermassive black holes fed by hot, tenuous accretion flows with much lower accretion rates (Ichimaru 1977;Narayan & Yi 1995;Blandford & Begelman 1999;Yuan & Narayan 2014).
In many AGNs, collimated relativistic plasma jets (Bridle & Perley 1984;Zensus 1997) launched by the central black hole contribute to the observed emission.These jets may be powered either by magnetic fields threading the event horizon, extracting the rotational energy from the black hole (Blandford & Znajek 1977), or from the accretion flow (Blandford & Payne 1982).The near-horizon emission from low-luminosity active galactic nuclei (LLAGNs; Ho 1999) is produced by synchrotron radiation that peaks from the radio through the farinfrared.This emission may be produced either in the accretion flow (Narayan et al. 1995), the jet (Falcke et al. 1993), or both (Yuan et al. 2002).
When viewed from infinity, a nonrotating Schwarzschild (1916) black hole has a photon capture radius R r 27 c g = , where r GM c g 2 º is the characteristic lengthscale of a black hole.The photon capture radius is larger than the Schwarzschild radius R S that marks the event horizon of a nonrotating black hole, R S ≡ 2 r g .Photons approaching the black hole with an impact parameter b<R c are captured and plunge into the black hole (Hilbert 1917); photons with b>R c escape to infinity; photons with b=R c are captured on an unstable circular orbit and produce what is commonly referred to as the lensed "photon ring."In the Kerr (1963) metric, which describes black holes with spin angular momentum, R c changes with the ray's orientation relative to the angular-momentum vector, and the black hole's cross section is not necessarily circular (Bardeen 1973).This change is small (4%), but potentially detectable (Takahashi 2004;Johannsen & Psaltis 2010).
The simulations of Luminet (1979) showed that for a black hole embedded in a geometrically thin, optically thick accretion disk, the photon capture radius would appear to a distant observer as a thin emission ring inside a lensed image of the accretion disk.For accreting black holes embedded in a geometrically thick, optically thin emission region, as in LLAGNs, the combination of an event horizon and light bending leads to the appearance of a dark "shadow" together with a bright emission ring that should be detectable through very long baseline interferometery (VLBI) experiments (Falcke et al. 2000a).Its shape can appear as a "crescent" because of fast rotation and relativistic beaming (Falcke et al. 2000b;Bromley et al. 2001;Noble et al. 2007;Broderick & Loeb 2009;Kamruddin & Dexter 2013;Lu et al. 2014).
The observed projected diameter of the emission ring, which contains radiation primarily from the gravitationally lensed photon ring, is proportional to R c and hence to the mass of the black hole, but also depends nontrivially on a number of factors: the observing resolution, the spin vector of the black hole and its inclination, as well as the size and structure of the emitting region.These factors are typically of order unity and can be calibrated using theoretical models.
Modern general-relativistic simulations of accretion flows and radiative transfer produce realistic images of black hole shadows and crescents for a wide range of near-horizon emission models (Broderick & Loeb 2006;Mościbrodzka et al. 2009;Dexter et al. 2012;Dibi et al. 2012;Chan et al. 2015;Mościbrodzka et al. 2016;Porth et al. 2017;Chael et al. 2018a;Ryan et al. 2018;Davelaar et al. 2019).These images can be used to test basic properties of black holes as predicted in GR (Johannsen & Psaltis 2010;Broderick et al. 2014;Psaltis et al. 2015), or in alternative theories of gravity (Grenzebach et al. 2014;Younsi et al. 2016;Mizuno et al. 2018).They can also be used to test alternatives to black holes (Bambi & Freese 2009;Vincent et al. 2016;Olivares et al. 2019).
VLBI at an observing wavelength of 1.3 mm (230 GHz) with Earth-diameter-scale baselines is required to resolve the shadows of the core of M87 (M87 * hereafter) and of the Galactic center of Sagittarius A * (Sgr A * , Balick & Brown 1974), the two supermassive black holes with the largest apparent angular sizes (Johannsen et al. 2012).At 1.3 mm and shorter wavelengths, Earthdiameter VLBI baselines achieve an angular resolution sufficient to resolve the shadow of both sources, while the spectra of both sources become optically thin, thus revealing the structure of the innermost emission region.Early pathfinder experiments (Padin et al. 1990;Krichbaum et al. 1998) demonstrated the feasibility of VLBI techniques at ∼1.3 mm wavelengths.Over the following decade, a program to improve sensitivity of 1.3 mm-VLBI through development of broadband instrumentation led to the detection of event-horizon-scale structures in both Sgr A * and M87 * (Doeleman et al. 2008(Doeleman et al. , 2012)).Building on these observations, the Event Horizon Telescope (EHT) collaboration was established to assemble a global VLBI array operating at a wavelength of 1.3 mm with the required angular resolution, sensitivity, and baseline coverage to image the shadows in M87 * and Sgr A * .
In this paper, we present and discuss the first event-horizonscale images of the supermassive black hole candidate M87

The Radio Core in M87
In Curtis (1918), Heber Curtis detected a linear feature in M87, later called a "jet" by Baade & Minkowski (1954).The jet is seen as a bright radio source, VirgoA or 3C 274 (Bolton et al. 1949;Kassim et al. 1993;Owen et al. 2000), which extends out to 65 kpc with an age estimated at about 40 Myr and a kinetic power of about 10 42 to 10 45 erg s −1 (de Gasperin et al. 2012;Broderick et al. 2015).It is also well studied in the optical (Biretta et al. 1999;Perlman et al. 2011), X-ray (Marshall et al. 2002), and gamma-ray bands (Abramowski et al. 2012).The upstream end of the jet is marked by a compact radio source (Cohen et al. 1969).Such compact radio sources are ubiquitous in LLAGNs (Wrobel & Heeschen 1984;Nagar et al. 2005) and are believed to be signatures of supermassive black holes.
The radio structures of the large-scale jet (Owen et al. 1989;de Gasperin et al. 2012) and of the core of M87 (Reid et al. 1989;Junor et al. 1999;Hada et al. 2016;Mertens et al. 2016;Kim et al. 2018b;Walker et al. 2018) have been resolved in great detail and at multiple wavelengths.Furthermore, the leveling-off of the "core-shift" effect (Blandford & Königl 1979), where the apparent position of the radio core shifts in the upstream jet direction with decreasing wavelength from increased transparency to synchrotron self-absorption, indicates that at a wavelength of 1.3 mm M87 * is coincident with the supermassive black hole (Hada et al. 2011).The envelope of the jet limb maintains a quasi-parabolic shape over a wide range of distances from ∼10 5 r g down to ∼20 r g (Asada & Nakamura 2012;Hada et al. 2013;Nakamura & Asada 2013;Nakamura et al. 2018;Walker et al. 2018).
VLBI observations at 1.3 mm have revealed a diameter of the emission region of ∼40 μas, which is comparable to the expected horizon-scale structure (Doeleman et al. 2012;Akiyama et al. 2015).These observations, however, were not able to image the black hole shadow due to limited baseline coverage.
Based on three recent stellar population measurements, we here adopt a distance to M87 of 16.8±0.8Mpc (Blakeslee et al. 2009;Bird et al. 2010;Cantiello et al. 2018, see Paper VI).Using this distance and the modeling of surface brightness and stellar velocity dispersion at optical wavelengths (Gebhardt & Thomas 2009;Gebhardt et al. 2011), we infer the mass of M87 * to be M 6.2 10 0.6 1.1 9 = -+ M e (see Table 9 in Paper VI).On the other hand, mass measurements modeling the kinematic structure of the gas disk (Harms et al. 1994;Macchetto et al. 1997) yield M 3.5 10 0.3 0.9 9 = -+ M e (Walsh et al. 2013, Paper VI).These two mass estimates, from stellar and gas dynamics, predict a theoretical shadow diameter for a Schwarzschild black hole of 37.6 as , respectively.

The Event Horizon Telescope
The EHT (Paper II) is a VLBI experiment that directly measures "visibilities," or Fourier components, of the radio brightness distribution on the sky.As the Earth rotates, each telescope pair in the network samples many spatial frequencies.
The array has a nominal angular resolution of λ/L, where λ is the observing wavelength and L is the maximum projected baseline length between telescopes in the array (Thompson et al. 2017).In this way, VLBI creates a virtual telescope that spans nearly the full diameter of the Earth.
To measure interferometric visibilities, the widely separated telescopes simultaneously sample and coherently record the radiation field from the source.Synchronization using the Global Positioning System typically achieves temporal alignment of these recordings within tens of nanoseconds.Each station is equipped with a hydrogen maser frequency standard.With the atmospheric conditions during our observations the coherent integration time was typically 10 s (see Figure 2 in Paper II). Use of hydrogen maser frequency standards at all EHT sites ensures coherence across the array over this timescale.After observations, recordings are staged at a central location, aligned in time, and signals from each telescope-pair are cross-correlated.
While VLBI is well established at centimeter and millimeter wavelengths (Boccardi et al. 2017;Thompson et al. 2017) and can be used to study the immediate environments of black holes (Krichbaum et al. 1993;Doeleman et al. 2001), the extension of VLBI to a wavelength of 1.3 mm has required long-term technical developments.Challenges at shorter wavelengths include increased noise in radio receiver electronics, higher atmospheric opacity, increased phase fluctuations caused by atmospheric turbulence, and decreased efficiency and size of radio telescopes in the millimeter and submillimeter observing bands.Started in 2009 (Doeleman et al. 2009a), the EHT began a program to address these challenges by increasing array sensitivity.Development and deployment of broadband VLBI systems (Whitney et al. 2013;Vertatschitsch et al. 2015) led to data recording rates that now exceed those of typical cm-VLBI arrays by more than an order of magnitude.Parallel efforts to support infrastructure upgrades at additional VLBI sites, including the Atacama Large Millimeter/submillimeter Array (ALMA; Matthews et al. 2018;Goddi et al. 2019) and the Atacama Pathfinder Experiment telescope (APEX) in Chile (Wagner et al. 2015), the Large Millimeter Telescope Alfonso Serrano (LMT) in Mexico (Ortiz-León et al. 2016), the IRAM 30 m telescope on Pico Veleta (PV) in Spain (Greve et al. 1995), the Submillimeter Telescope Observatory in Arizona (SMT; Baars et al. 1999), the James Clerk Maxwell Telescope (JCMT) and the Submillimeter Array (SMA) in Hawai'i (Doeleman et al. 2008;Primiani et al. 2016;Young et al. 2016), and the South Pole Telescope (SPT) in Antarctica (Kim et al. 2018a), extended the range of EHT baselines and coverage, and the overall collecting area of the array.These developments increased the sensitivity of the EHT by a factor of ∼30 over early experiments that confirmed horizon-scale structures in M87 * and Sgr A * (Doeleman et al. 2008(Doeleman et al. , 2012;;Akiyama et al. 2015;Johnson et al. 2015;Fish et al. 2016;Lu et al. 2018).
For the observations at a wavelength of 1.3 mm presented here, the EHT collaboration fielded a global VLBI array of eight stations over six geographical locations.Baseline lengths ranged from 160 m to 10,700 km toward M87 * , resulting in an array with a theoretical diffraction-limit resolution of ∼25 μas (see Figures 1 and 2, and Paper II).

Observations, Correlation, and Calibration
We observed M87 * on 2017 April 5, 6, 10, and 11 with the EHT.Weather was uniformly good to excellent with nightly median zenith atmospheric opacities at 230 GHz ranging from 0.03 to 0.28 over the different locations.The observations were scheduled as a series of scans of three to seven minutes in duration, with M87 * scans interleaved with those on the quasar 3C 279.The number of scans obtained on M87 * per night ranged from 7 (April 10) to 25 (April 6) as a result of different observing schedules.A description of the M87 * observations, their correlation, calibration, and validated final data products is presented in Paper III and briefly summarized here.
At each station, the astronomical signal in both polarizations and two adjacent 2 GHz wide frequency bands centered at 227.1 and 229.1 GHz were converted to baseband using standard heterodyne techniques, then digitized and recorded at atotal rate of 32 Gbps.Correlation of the data was carried out using a software correlator (Deller et al. 2007) at the MIT Haystack Observatory and at the Max-Planck-Institut für Radioastronomie, each handling one of the two frequency bands.Differences between the two independent correlators were shown to be negligible through the exchange of a few identical scans for cross comparison.At correlation, signals were aligned to a common time reference using an apriori Earth geometry and clock model.
A subsequent fringe-fitting step identified detections in correlated signal power while phase calibrating the data for residual delays and atmospheric effects.Using ALMA as a highly sensitive reference station enabled critical corrections for ionospheric and tropospheric distortions at the other sites.Fringe fitting was performed with three independent automated pipelines, each tailored to the specific characteristics of the EHT observations, such as the wide bandwidth, susceptibility to atmospheric turbulence, and array heterogeneity (Blackburn et al. 2019;Janssen et al. 2019, Paper III).The pipelines made use of standard software for the processing of radio-interferometric data  (Greisen 2003;Whitney et al. 2004;McMullin et al. 2007, I. M. van Bemmel et al. 2019, in preparation).
Data from the fringe-fitting pipelines were scaled from correlation coefficients to a uniform physical flux density scale (in Jansky) by using an independent apriori estimate of the sensitivity of each telescope.The accuracies of the derived station sensitivities were estimated to be 5%-10% in amplitude, although certain uncharacterized losses (e.g., from poor pointing or focus) can exceed the error budget.By assuming total flux density values derived from ALMA interferometric data (Goddi et al. 2019) and utilizing array redundancy via network calibration (Paper III), we refined the absolute amplitude calibration of telescopes that are colocated and have redundant baselines, i.e., ALMA/APEX and JCMT/SMA.
The median scan-averaged signal-to-noise ratio for M87 * was >10 on non-ALMA baselines and >100 on baselines to ALMA, leading to small statistical errors in visibility amplitude and phase.Comparisons between the three independent pipelines, the two polarizations, and the two frequency bands enabled estimation of systematic baseline errors of around 1°in visibility phase and 2% for visibility amplitudes.These small limiting errors remain after fitting station sensitivities and unknown station phases via self-calibration (Pearson & Readhead 1984) and affect interferometric closure quantities (Rogers et al. 1974;Readhead et al. 1980).Following data validation and pipeline comparisons, a single pipeline output was designated as the primary data set of the first EHT science data release and used for subsequent results, while the outputs of the other two pipelines offer supporting validation data sets.
The final calibrated complex visibilities V(u, v) correspond to the Fourier components of the brightness distribution on the sky at spatial frequency (u, v) determined by the projected baseline expressed in units of the observing wavelength (van Cittert 1934; Thompson et al. 2017).Figure 2 shows the (u, v) coverage and calibrated visibility amplitudes of M87 * for April11.The visibility amplitudes resemble those of a thin ring (i.e., a Bessel function J 0 ; see Figure 10.12 in Thompson et al. 2017).Such a ring model with diameter 46 μas has afirst null at 3.4 Gλ, matching the minimum in observed flux density and is consistent with a reduced flux density on the longest Hawai'i-Spain baseline (JCMT/SMA-PV) near 8 Gλ.This particular ring model, shown with a dashed line in the bottom panel of Figure 2, is only illustrative and does not fit all features in the data.First, visibility amplitudes on the shortest VLBI baselines suggest that about half of the compact flux density seen on the ∼2 km ALMA-APEX baseline is resolved out by the interferometer beam (Paper IV).Second, differences in the depth of the first minimum as a function of orientation, as well as highly nonzero measured closure phases, indicate some degree of asymmetry in the source (Papers III, VI).Finally, the visibility amplitudes represent only half of the information available to us.We will next explore images and more complex geometrical models that can fit the measured visibility amplitudes and phases.

Images and Features
We reconstructed images from the calibrated EHT visibilities, which provide results that are independent of models (Paper IV).However, there are two major challenges in reconstructing images from EHT data.First, EHT baselines sample a limited range of spatial frequencies, corresponding to angular scales between 25 and 160 μas.Because the (u, v) plane is only sparsely sampled (Figure 2), the inverse problem is under-constrained.Second, the measured visibilities lack absolute phase calibration and can have large amplitude calibration uncertainties.
To address these challenges, imaging algorithms incorporate additional assumptions and constraints that are designed to produce images that are physically plausible (e.g., positive and compact) or conservative (e.g., smooth), while remaining consistent with the data.We explored two classes of algorithms for reconstructing images from EHT data.The first class of algorithms is the traditional CLEAN approach used in radio interferometry (e.g., Högbom 1974;Clark 1980).CLEAN is an inverse-modeling approach that deconvolves the interferometer point-spread function from the Fourier-transformed visibilities.When applying CLEAN, it is necessary to iteratively self-calibrate the data between rounds of imaging to solve for time-variable phase and amplitude errors in the data.The second class of algorithms is the so-called regularized maximum likelihood (RML; e.g., Narayan & Nityananda 1986;Wiaux et al. 2009;Thiébaut 2013).RML is a forward-modeling approach that searches for an image that is not only consistent with the observed data but also favors specified image properties (e.g., smoothness or compactness).As with CLEAN, RML methods typically iterate between imaging and self-calibration, although they can also be used to image directly on robust closure quantities immune to station-based calibration errors.RML methods have been extensively developed for the EHT (e.g., Honma et al. 2014;Bouman et al. 2016;Akiyama et al. 2017;Chael et al. 2018b; see also Paper IV).
Every imaging algorithm has a variety of free parameters that can significantly affect the final image.We adopted a twostage imaging approach to control and evaluate biases in the reconstructions from our choices of these parameters.In the first stage, four teams worked independently to reconstruct the first EHT images of M87 * using an early engineering data release.The teams worked without interaction to minimize shared bias, yet each produced an image with a similar prominent feature: a ring of diameter ∼38-44 μas with enhanced brightness to the south (see Figure 4 in Paper IV).
In the second imaging stage, we developed three imaging pipelines, each using a different software package and associated methodology.Each pipeline surveyed a range of imaging parameters, producing between ∼10 3 and 10 4 images from different parameter combinations.We determined a "Top-Set" of parameter combinations that both produced images of M87 * that were consistent with the observed data and that reconstructed accurate images from synthetic data sets corresponding to four known geometric models (ring, crescent, filled disk, and asymmetric double source).For all pipelines, the Top-Set images showed an asymmetric ring with a diameter of ∼40 μas, with differences arising primarily in the effective angular resolutions achieved by different methods.
For each pipeline, we determined the single combination of fiducial imaging parameters out of the Top-Set that performed best across all the synthetic data sets and for each associated imaging methodology (see Figure 11 in Paper IV).Because the angular resolutions of the reconstructed images vary among the pipelines, we blurred each image with a circular Gaussian to a common, conservative angular resolution of 20 μas.The top part of Figure 3 shows an image of M87 * on April11 obtained by averaging the three pipelines' blurred fiducial images.The image is dominated by a ring with an asymmetric azimuthal profile that is oriented at a position angle ∼170°east of north.Although the measured position angle increases by ∼20°between the first two days and the last two days, the image features are broadly consistent across the different imaging methods and across all four observing days.This is shown in the bottom part of Figure 3, which reports the images on different days (see also Figure 15 in Paper IV).These results are also consistent with those obtained from visibility-domain fitting of geometric and general-relativistic magnetohydrodynamics (GRMHD) models (Paper VI).

Theoretical Modeling
The appearance of M87 * has been modeled successfully using GRMHD simulations, which describe a turbulent, hot, magnetized disk orbiting a Kerr black hole.They naturally produce a powerful jet and can explain the broadband spectral energy distribution observed in LLAGNs.At a wavelength of 1.3 mm, and as observed here, the simulations also predict a shadow and an asymmetric emission ring.The latter does not necessarily coincide with the innermost stable circular orbit, or ISCO, and is instead related to the lensed photon ring.To explore this scenario in great detail, we have built a library of synthetic images (Image Library) describing magnetized accretion flows onto black holes in GR145 (Paper V).The images themselves are produced from a library of simulations (Simulation Library) collecting the results of four codes solving the equations of GRMHD (Gammie et al. 2003;Saḑowski et al. 2014;Porth et al. 2017;Liska et al. 2018).The elements of the Simulation Library have been coupled to three different general-relativistic ray-tracing and radiative-transfer codes (GRRT, Bronzwaer et al. 2018;Mościbrodzka & Gammie 2018;Z. Younsi et al. 2019, in preparation).We limit ourselves to providing here a brief description of the initial setups and the physical scenarios explored in the simulations; see Paper V for details on both the GRMHD and GRRT codes, which have been cross-validated for accuracy and consistency (Gold et al. 2019;Porth et al. 2019).
A typical GRMHD simulation in the library is characterized by two parameters: the dimensionless spin a Jc GM 2 * º , where J and M are, respectively, the spin angular momentum and mass of the black hole, and the net dimensionless magnetic flux over the event horizon MR g 2 1 2 f º F ( ˙) , where Φ and M ˙are the magnetic flux and mass flux (or accretion rate) across the horizon, respectively.Since the GRMHD simulations scale with the black hole mass, M is set only at the time of producing the synthetic images with the GRRT codes.The magnetic flux is generally nonzero because magnetic field is trapped in the black hole by the accretion flow and sustained by currents in the surrounding plasma.
These two parameters allow us to describe accretion disks that are either prograde (a * 0) or retrograde (a * <0) with respect to the black hole spin axis, and whose accretion flows are either "SANE" (from "Standard and Normal Evolution," Narayan et al. 2012) with f∼1, or "MAD" (from "Magnetically Arrested Disk," Narayan et al. 2003) with f∼15. 146In essence, SANE accretion flows are characterized by moderate dimensionless magnetic flux and result from initial magnetic fields that are smaller than those in MAD flows.Furthermore, the opening angles of the magnetic funnel in SANE flows are generically smaller than those in MAD flows.Varying a * and f, we have performed 43 high-resolution, three-dimensional and long-term simulations covering well the physical properties of magnetized accretion flows onto Kerr black holes.
All GRMHD simulations were initialized with a weakly magnetized torus orbiting around the black hole and driven into a turbulent state by instabilities, including the magnetorotational instability (Balbus & Hawley 1991), rapidly reaching a quasistationary state.Once a simulation was completed, the relevant flow properties at different times are collected to be employed for the further post-processing of the GRRT codes.The generation of synthetic images requires, besides the properties of the fluid (magnetic field, velocity field, and rest-mass density), also the emission and absorption coefficients, the inclination i (the angle between the accretion-flow angular-momentum vector and the line of sight), the position angle PA (the angle east of north, i.e., counter-clockwise on our images, of the projection on the sky of the accretion-flow angular momentum), the black hole mass M and distance D to the observer.
Because the photons at 1.3 mm wavelength observed by the EHT are believed to be produced by synchrotron emission, whose absorption and emission coefficients depend on the electron distribution function, we consider the plasma to be composed of electrons and ions that have the same temperature in the magnetically dominated regions of the flow (funnel), but have a substantially different temperature in the gas dominated regions (disk midplane).In particular, we consider the plasma to be composed of nonrelativistic ions with temperature T i and relativistic electrons with temperature T e .A simple prescription for the ratio of the temperatures of the two species can then be imposed in terms of a single parameter (Mościbrodzka et al. 2016), such that the bulk of the emission comes either from weakly magnetized (small R high , Te ; Ti/R high ) or strongly magnetized (large R high , Te ; Ti) regions.In SANE models, the disk (jet) is weakly (strongly) magnetized, so that low (high) R high models produce most of the emission in the disk (jet).In MAD models, there are strongly magnetized regions everywhere and the emission is mostly from the disk midplane.While this prescription is not the only one possible, it has the advantage of being simple, sufficiently generic, and robust (see Paper V for a discussion of nonthermal particles and radiative cooling).
Since each GRMHD simulation can be used to describe several different physical scenarios by changing the prescribed electron distribution function, we have used the Simulation Library to generate more than 420 different physical scenarios.Each scenario is then used to generate hundreds of snapshots at different times in the simulation leading to more than 62,000 objects in the Image Library.From the images we have created model visibilities that correspond to the EHT observing schedule and compared them to the measured VLBI visibilities as detailed in Paper VI.
As an example, the top row of Figure 4 shows three GRMHD model snapshots from the Image Library with different spins and flow type, which fitted closure phases and amplitudes of the April 11 data best.For these models we produced simulated visibilities for the April 11 schedule and weather parameters and calibrated them with a synthetic data generation and calibration pipeline (Blecher et al. 2017;Janssen et al. 2019;Roelofs et al. 2019a).The simulated data were then imaged with the same pipeline used for the observed images.The similarities between the simulated images (bottom row of Figure 4) and the observed images (Figure 3) are remarkable.
Overall, when combining all the information contained in both the Simulation Library and Image Library, the physical origin of the emission features of the image observed in M87 * can be summarized as follows.
First, the observed image is consistent with the hypothesis that it is produced by a magnetized accretion flow orbiting within a few r g of the event horizon of a Kerr black hole.The asymmetric ring is produced by a combination of strong gravitational lensing and relativistic beaming, while the central flux depression is the observational signature of the black hole shadow.Interestingly, all of the accretion models are consistent with the EHT image, except for the a * =−0.94MAD models, which fail to produce images that are sufficiently stable (i.e., the variance among snapshots is too large to be statistically consistent with the observed image).
Second, the north-south asymmetry in the emission ring is controlled by the black hole spin and can be used to deduce its orientation.In corotating disk models (where the angular momentum of the matter and of the black hole are aligned), the funnel wall, or jet sheath, rotates with the disk and the black hole; in counterrotating disk models, instead, the luminous funnel wall rotates with the black hole but against the disk.The north-south asymmetry is consistent with models in which the black hole spin is pointing away from Earth and inconsistent with models in which the spin points toward Earth.
Third, adopting an inclination of 17°between the approaching jet and the line of sight (Walker et al. 2018), the west orientation of the jet, and a corotating disk model, matter in the bottom part of the image is moving toward the observer (clockwise rotation as seen from Earth).This is consistent with the rotation of the ionized gas on scales of 20 pc, i.e., 7000 r g (Ford et al. 1994;Walsh et al. 2013) and with the inferred sense of rotation from VLBI observations at 7 mm (Walker et al. 2018).
Finally, models with a * =0 are disfavored by the very conservative observational requirement that the jet power be P jet >10 42 erg s −1 .Furthermore, in those models that produce a sufficiently powerful jet, it is powered by extraction of black hole spin energy through mechanisms akin to the Blandford-Znajek process.

Model Comparison and Parameter Estimation
In Paper VI, the black hole mass is derived from fitting to the visibility data of geometric and GRMHD models, as well as from measurements of the ring diameter in the image domain.Our measurements remain consistent across methodologies, algorithms, data representations, and observed data sets.
Motivated by the asymmetric emission ring structures seen in the reconstructed images (Section 5) and the similar emission structures seen in the images from GRMHD simulations (Section 6), we developed a family of geometric crescent models(see, e.g., Kamruddin & Dexter 2013) to compare directly to the visibility data.We used two distinct Bayesian-inference algorithms and demonstrate that such crescent models are statistically preferred over other comparably complex geometric models that we have explored.We find that the crescent models provide fits to the data that are statistically comparable to those of the reconstructed images presented in Section 5, allowing us to determine the basic parameters of the crescents.The best-fit models for the asymmetric emission ring have diameters of 43±0.9μas and fractional widths relative to the diameter of <0.5.The emission drops sharply interior to the asymmetric emission ring with the central depression having a brightness <10% of the average brightness in the ring.
The diameters of the geometric crescent models measure the characteristic sizes of the emitting regions that surround the shadows and not the sizes of the shadows themselves (see, e.g., Psaltis et al. 2015;Johannsen et al. 2016;Kuramochi et al. 2018, for potential biases).
We model the crescent angular diameter d in terms of the gravitational radius and distance, GM c D g 2 q º , as d=αθ g , where α is a function of spin, inclination, and R high (α;9.6-10.4corresponds to emission from the lensed photon ring only).We calibrate α by fitting the geometric crescent models to a large number of visibility data generated from the Image Library.We can also fit the model visibilities generated from the Image Library to the M87 * data, which allows us to measure θ g directly.However, such a procedure is complicated by the stochastic nature of the emission in the accretion flow(see, e.g., Kim et al. 2016).To account for this turbulent structure, we developed a formalism and multiple algorithms that estimate the statistics of the stochastic components using ensembles of images from individual GRMHD simulations.We find that the visibility data are not inconsistent with being a realization of many of the GRMHD simulations.We conclude that the recovered model parameters are consistent across algorithms.
Finally, we extract ring diameter, width, and shape directly from reconstructed images (see Section 5).The results are consistent with the parameter estimates from geometric crescent models.Following the same GRMHD calibration procedure, we infer values of θ g and α for regularized maximum likelihood and CLEAN reconstructed images.
Combining results from all methods, we measure emission region diameters of 42±3 μas, angular sizes of the gravitational radius θ g =3.8±0.4μas, and scaling factors in the range α=10.7-11.5, with associated errors of ∼10%.For the distance of 16.8±0.8Mpc adopted here, the black hole mass is M=(6.5±0.7)×10 9M e ; the systematic error refers to the 68% confidence level and is much larger than the statistical error of 0.2×10 9 M e .Moreover, by tracing the peak of the emission in the ring we can determine the shape of the image and obtain a ratio between major and minor axis of the ring that is 4:3; this corresponds to a 10% deviation from circularity in terms of root-mean-square distance from an average radius.
Table 1 summarizes the measured parameters of the image features and the inferred black hole properties based on data from all bands and all days combined.The inferred black hole mass strongly favors the measurement based on stellar dynamics (Gebhardt et al. 2011).The size, asymmetry, brightness contrast, and circularity of the reconstructed images and geometric models, as well as the success of the GRMHD simulations in describing the interferometric data, are consistent with the EHT images of M87 * being associated with strongly lensed emission from the vicinity of a Kerr black hole.

Discussion
A number of elements reinforce the robustness of our image and the conclusion that it is consistent with the shadow of a black hole as predicted by GR.First, our analysis has used multiple independent calibration and imaging techniques, as well as four independent data sets taken on four different days in two separate frequency bands.Second, the image structure matches previous predictions well (Dexter et al. 2012;Mościbrodzka et al. 2016) and is well reproduced by our extensive modeling effort presented in Section 6.Third, because our measurement of the black hole mass in M87 * is not inconsistent with all of the prior mass measurements, this allows us to conclude that the null hypothesis of the Kerr metric (Psaltis et al. 2015;Johannsen et al. 2016), namely, the assumption that the black hole is described by the Kerr metric, has not been violated.Fourth, the observed emission ring reconstructed in our images is close to circular with an axial ratio 4:3; similarly, the time average images from our GRMHD simulations also show a circular shape.After associating to the shape of the shadow a deviation from the circularity-measured in terms of root-mean-square distance from an average radius in the image-that is 10%, we can set an initial limit of order four on relative deviations of the quadrupole moment from the Kerr value (Johannsen & Psaltis 2010).Stated differently, if Q is the quadrupole moment of a Kerr black hole and ΔQ the deviation as deduced from circularity, our measurement-and the fact that the inclination angle is assumed to be small-implies that ΔQ/Q4 (ΔQ/Q=ε in Johannsen & Psaltis 2010).
Finally, when comparing the visibility amplitudes of M87 * with 2009 and 2012 data (Doeleman et al. 2012;Akiyama et al. 2015), the overall radio core size at a wavelength of 1.3 mm has not changed appreciably, despite variability in total flux density.This stability is consistent with the expectation that the size of the shadow is a feature tied to the mass of the black hole and not to properties of a variable plasma flow.
It is also straightforward to reject some alternative astrophysical interpretations.For instance, the image is unlikely to be produced by a jet-feature as multi-epoch VLBI observations of the plasma jet in M87 (Walker et al. 2018) on scales outside the horizon do not show circular rings.The same is typically true for AGN jets in large VLBI surveys (Lister et al. 2018).Similarly, were the apparent ring a random alignment of emission blobs, they should also have moved away at relativistic speeds, i.e., at ∼5 μas day −1 (Kim et al. 2018b), leading to measurable structural changes and sizes.GRMHD models of hollow jet cones could show under extreme conditions stable ring features (Pu et al. 2017), but this effect is included to a certain extent in our Simulation Library for models with R high >10.Finally, an Einstein ring formed by gravitational lensing of a bright region in the counter-jet would require a fine-tuned alignment and a size larger than that measured in 2012 and 2009.
At the same time, it is more difficult to rule out alternatives to black holes in GR, because a shadow can be produced by any compact object with a spacetime characterized by unstable circular photon orbits (Mizuno et al. 2018).Indeed, while the Kerr metric remains a solution in some alternative theories of gravity (Barausse & Sotiriou 2008;Psaltis et al. 2008), non-Kerr black hole solutions do exist in a variety of such modified theories (Berti et al. 2015).Furthermore, exotic alternatives to black holes, such as naked singularities (Shaikh et al. 2019), boson stars (Kaup 1968;Liebling & Palenzuela 2012), and gravastars (Mazur & Mottola 2004;Chirenti & Rezzolla 2007), are admissible solutions within GR and provide concrete, albeit contrived, models.Some of such exotic compact objects can already be shown to be incompatible with our observations given our maximum mass prior.For example, the shadows of naked singularities associated with Kerr spacetimes with are substantially smaller and very asymmetric compared to those of Kerr black holes (Bambi & Freese 2009).Also, some commonly used types of wormholes (Bambi 2013) predict much smaller shadows than we have measured.However, other compact-object candidates need to be analyzed with more care.Boson stars are an example of compact objects having circular photon orbits but without a surface or an event horizon.In such a spacetime, null geodesics are redirected outwards toward distant observers (Cunha et al. 2016), so that the shadow can in principle be filled with emission from lensed images of distant radio sources generating a complex mirror image of the sky.More importantly, accretion flows onto boson stars behave differently as they do not produce jets but stalled accretion tori that make them distinguishable from black holes (Olivares et al. 2019).Gravastars provide examples of compact objects having unstable photon orbits and a hard surface, but not an event horizon.In this case, while a single image of the accretion flow could in principle be very similar to that coming from a black hole, differences of the flow dynamics at the stellar surface (H.Olivares et al. 2019, in preparation), strong magnetic fields (Lobanov 2017), or excess radiation in the infrared (Broderick & Narayan 2006) would allow one to distinguish a gravastar from a black hole.
Altogether, the results derived here provide a new way to study compact-object spacetimes and are complementary to the detection of gravitational waves from coalescing stellar-mass black holes with LIGO/Virgo (Abbott et al. 2016).Our constraints on deviations from the Kerr geometry rely only on the validity of the equivalence principle and are agnostic about the underlying theory of gravity, but can be used to measure, with ever improved precision, the parameters of the background metric.On the other hand, current gravitational-wave observations of mergers probe the dynamics of the underlying theory, but cannot rely on the possibility of multiple and repeated measurements of the same source.
To underline the complementarity of gravitational-wave and electromagnetic observations of black holes, we note that a basic feature of black holes in GR is that their size scales linearly with mass.Combining our constraints on the supermassive black hole in M87 with those on the stellar-mass black holes detected via gravitational waves we can infer that this property holds over eight orders of magnitude.While this invariance is a basic prediction of GR, which is a scale-free theory, it need not be satisfied in other theories that introduce a scale to the gravitational field.
Finally, the radio core in M87 is quite typical for powerful radio jets in general.It falls on the fundamental plane of black hole activity for radio cores (Falcke et al. 2004), connecting via simple scaling laws the radio and X-ray properties of lowluminosity black hole candidates across vastly different mass and accretion rate scales.This suggests that they are powered by a scale-invariant common object.Therefore, establishing the black hole nature for M87 * also supports the general paradigm that black holes are the power source for active galaxies.

Conclusion and Outlook
We have assembled the EHT, a global VLBI array operating at a wavelength of 1.3 mm and imaged horizon-scale structures around the supermassive black hole candidate in M87.Using multiple independent calibration, imaging, and analysis methods, we find the image to be dominated by a ring structure of 42±3 μas diameter that is brighter in the south.This structure has a central brightness depression with a contrast of >10:1, which we identify with the black hole shadow.Comparing the data with an extensive library of synthetic images obtained from GRMHD simulations covering different physical scenarios and plasma conditions reveals that the basic features of our image are relatively independent of the detailed astrophysical model.This allows us to derive an estimate for the black hole mass of M=(6.5±0.7)× 10 9 M e .Based on our modeling and information on the inclination angle, we derive the sense of rotation of the black hole to be in the clockwise direction, i.e., the spin of the black hole points away from us.The brightness excess in the south part of the emission ring is explained as relativistic beaming of material rotating in the clockwise direction as seen by the observer, i.e., the bottom part of the emission region is moving toward the observer.
Future observations and further analysis will test the stability, shape, and depth of the shadow more accurately.One of its key features is that it should remain largely constant with time as the mass of M87 * is not expected to change measurably on human timescales.Polarimetric analysis of the images, which we will report in the future, will provide information on the accretion rate via Faraday rotation (Bower et al. 2003;Marrone et al. 2007;Kuo et al. 2014;Mościbrodzka et al. 2017) and on the magnetic flux.Higher-resolution images can be achieved by going to a shorter wavelength, i.e., 0.8 mm (345 GHz), by adding more telescopes and, in a more distant future, with space-based interferometry (Kardashev et al. 2014;Fish et al. 2019;Palumbo et al. 2019;F. Roelofs et al. 2019b, in preparation).
Another primary EHT source, Sgr A * , has a precisely measured mass three orders of magnitude smaller than that of M87 * , with dynamical timescales of minutes instead of days.Observing the shadow of Sgr A * will require accounting for this variability and mitigation of scattering effects caused by the interstellar medium (Johnson 2016;Lu et al. 2016;Bouman et al. 2018).Time dependent nonimaging analysis can be used to potentially track the motion of emitting matter near the black hole, as reported recently through interferometric observations in the near-infrared (Gravity Collaboration et al. 2018b).Such observations provide separate tests and probes of GR on yet another mass scale (Broderick & Loeb 2005;Doeleman et al. 2009b;Roelofs et al. 2017;Medeiros et al. 2017).
In conclusion, we have shown that direct studies of the event horizon shadow of supermassive black hole candidates are now possible via electromagnetic waves, thus transforming this elusive boundary from a mathematical concept to a physical entity that can be studied and tested via repeated astronomical observations.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 authors of this
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 opensource 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.

Figure 1 .
Figure 1.Eight stations of the EHT 2017 campaign over six geographic locations as viewed from the equatorial plane.Solid baselines represent mutual visibility on M87 * (+12°declination).The dashed baselines were used for the calibration source 3C279 (see Papers III and IV).

Figure 2 .
Figure 2. Top: (u, v) coverage for M87 * , aggregated over all four days of the observations.(u, v) coordinates for each antenna pair are the source-projected baseline length in units of the observing wavelength λ and are given for conjugate pairs.Baselines to ALMA/APEX and to JCMT/SMA are redundant.Dotted circular lines indicate baseline lengths corresponding to fringe spacings of 50 and 25 μas.Bottom:final calibrated visibility amplitudes of M87 * as a function of projected baseline length on April 11.Redundant baselines to APEX and JCMT are plotted as diamonds.Error bars correspond to thermal (statistical) uncertainties.The Fourier transform of an azimuthally symmetric thin ring model with diameter 46 μas is also shown with a dashed line for comparison.

Figure 3 .
Figure 3. Top: EHT image of M87 * from observations on 2017 April 11 as a representative example of the images collected in the 2017 campaign.The image is the average of three different imaging methods after convolving each with a circular Gaussian kernel to give matched resolutions.The largest of the three kernels (20 μas FWHM) is shown in the lower right.The image is shown in units of brightness temperature, T S k 2 b 2 B l = W, where S is the flux density, λ is the observing wavelength, k B is the Boltzmann constant, and Ω is the solid angle of the resolution element.Bottom: similar images taken over different days showing the stability of the basic image structure and the equivalence among different days.North is up and east is to the left.

Figure 4 .
Figure 4. Top: three example models of some of the best-fitting snapshots from the image library of GRMHD simulations for April 11 corresponding to different spin parameters and accretion flows.Bottom: the same theoretical models, processed through a VLBI simulation pipeline with the same schedule, telescope characteristics, and weather parameters as in the April 11 run and imaged in the same way as Figure 3.Note that although the fit to the observations is equally good in the three cases, they refer to radically different physical scenarios; this highlights that a single good fit does not imply that a model is preferred over others (see Paper V).
the image domain.b Derived from crescent model fitting.c The mass and systematic errors are averages of the three methods (geometric models, GRMHD models, and image domain ring extraction).d The exact value depends on the method used to extract d, which is reflected in the range given.e Rederived from likelihood distributions (Paper VI).
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 Stiftung; the Black Hole Initiative at Harvard University, through a grant (60477) from the John Templeton Foundation; the China Scholarship Council; Comisión Nacional de Investigación Científica y Tecnológica (CONICYT, Chile, via PIA ACT172033, Fondecyt 1171506, BASAL AFB-170002, ALMA-conicyt 31140007); Consejo Nacional de Ciencia y Tecnología (CONACYT, Mexico, projects 104497, 275201, 279006, 281692); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Planck-Gesellschaft, Germany) and IGN (Instituto Geográfico Nacional, Spain).
from an EHT VLBI campaign conducted in 2017 April at a wavelength of 1.3 mm.The accompanying papers give a more extensive description of the instrument (EHT Collaboration et al. 2019a, Paper II), data reduction (EHT Collaoration et al. 2019b, hereafter Paper III), imaging of the M87 shadow (EHT Collaboration et al. 2019c, hereafter Paper IV), theoretical models (EHT Collaboration et al. 2019d, hereafter Paper V), and the black hole mass estimate (EHT Collaboration et al. 2019e, hereafter Paper VI). *