Broadband Multi-wavelength Properties of M87 during the 2017 Event Horizon Telescope Campaign

In 2017, the Event Horizon Telescope (EHT) Collaboration succeeded in capturing the first direct image of the center of the M87 galaxy. The asymmetric ring morphology and size are consistent with theoretical expectations for a weakly accreting supermassive black hole of mass approximately 6.5 x 10^9 M_solar. The EHTC also partnered with several international facilities in space and on the ground, to arrange an extensive, quasi-simultaneous multi-wavelength campaign. This Letter presents the results and analysis of this campaign, as well as the multi-wavelength data as a legacy data repository. We captured M87 in a historically low state, and the core flux dominates over HST-1 at high energies, making it possible to combine core flux constraints with the more spatially precise very long baseline interferometry data. We present the most complete simultaneous multi-wavelength spectrum of the active nucleus to date, and discuss the complexity and caveats of combining data from different spatial scales into one broadband spectrum. We apply two heuristic, isotropic leptonic single-zone models to provide insight into the basic source properties, but conclude that a structured jet is necessary to explain M87's spectrum. We can exclude that the simultaneous gamma-ray emission is produced via inverse Compton emission in the same region producing the EHT mm-band emission, and further conclude that the gamma-rays can only be produced in the inner jets (inward of HST-1) if there are strongly particle-dominated regions. Direct synchrotron emission from accelerated protons and secondaries cannot yet be excluded.


Introduction
M87 is the most prominent elliptical galaxy within the Virgo Cluster, located just 16.8 ± 0.8 Mpc away (Blakeslee et al. 2009;Bird et al. 2010;Cantiello et al. 2018, and see also EHT Collaboration et al. 2019f).As one of our closest active galactic nuclei (AGNs), M87 also harbors the first example of an extragalactic jet to have been noticed by astronomers (Curtis 1918), well before these jets were understood to be a likely signature of black hole accretion.By now this famous one-sided jet has been well-studied in almost every wave band from radio (down to subparsec scales; e.g., Reid et al. 1989;Junor et al. 1999;Hada et al. 2011;Mertens et al. 2016;Kim et al. 2018b;Walker et al. 2018), optical (e.g., Biretta et al. 1999;Perlman et al. 2011), X-ray (e.g., Marshall et al. 2002;Snios et al. 2019), and γ-rays (e.g., Abdo et al. 2009a;Abramowski et al. 2012;MAGIC Collaboration et al. 2020).Extending over 60 kpc in length, the jet shows a system of multiple knot-like features, including an active feature HST-1 at a projected distance of ∼70 pc from the core, potentially marking the end of the black hole's sphere of gravitational influence (Asada & Nakamura 2012).In contrast, the weakly accreting supermassive black hole (SMBH) in the center of our own Milky Way galaxy, Sgr A * , does not show obvious signs of extended jets, although theory predicts that such outflows should be formed (e.g., Dibi et al. 2012;Mościbrodzka et al. 2014;Davelaar et al. 2018).The conditions under which such jets are launched is one of the enduring questions in astrophysics today (e.g., Blandford & Znajek 1977;Blandford & Payne 1982;Sikora & Begelman 2013).
In 2019 April, the Event Horizon Telescope (EHT) Collaboration presented the first direct image of an SMBH "shadow" in the center of M87 (EHT Collaboration et al. 2019a, 2019b, 2019c, 2019d, 2019e, 2019f).The key result in these papers was the detection of an asymmetric ring (crescent) of light around a darker circle, due to the presence of an event horizon, along with detailed explanations of all the ingredients necessary to obtain, analyze, and interpret this rich data set.The ring itself stems from a convolution of the light produced near the last unstable photon orbit, as it travels through the geometry of the production region with radiative transfer in the surrounding plasma, and further experiences bending and redshifting due to the effects of general relativity (GR).Photons that orbit, sometimes multiple times before escaping, trace out a sharp feature revealing the shape of the spacetime metric, the so-called "photon ring".From the size of the measured, blurry ring and three different modeling approaches (see EHT Collaboration et al. 2019e), the EHT Collaboration (EHTC) calibrated for these multiple effects, to derive a mass for M87's SMBH of (6.5 ± 0.7) × 10 9 M e .
One of the primary contributions to the ∼10% systematic error on this mass is due to uncertainties in the underlying accretion properties.As detailed in EHT Collaboration et al. (2019d) and Porth et al. (2019), the EHTC ran over 45 high-resolution simulations over a range of possible physical parameters in, e.g., spin, magnetic field configuration, and electron thermodynamics, using several different GR magnetohydrodynamic (GRMHD) codes.These outputs were then coupled to GR ray-tracing codes, to generate ∼60,000 images captured at different times during the simulation runs, which were verified in Gold et al. (2020).While the photon ring remains relatively robust to changes in spin, in part because of the small viewing angle (see, e.g., Johannsen & Psaltis 2010), the spreading of the light around this feature strongly depends on the plasma properties near the event horizon, introducing significant degeneracy.For instance, a smaller black hole produces a smaller photon ring, but in some emission models there is extended surrounding emission leading to a larger final blurry ring.Similarly, a larger black hole with significant emission produced along the line of sight will appear to have emission within the photon ring, and when convolved produces the appearance of a smaller blurry ring.The error in calibrating from a given image to a unique black hole mass is therefore a combination of image reconstruction limits, as well as our current level of uncertainty about the plasma properties and emission geometry very close to the black hole.
However, it is important to note that even in these first analyses, several of the models could already be ruled out using complementary information from observations with facilities at other wavelengths.For instance, the estimated minimum power in the jets, P jet 10 42 erg s −1 , from prior and recent multiwavelength studies (e.g., Reynolds et al. 1996;Stawarz et al. 2006;de Gasperin et al. 2012;Prieto et al. 2016), was already enough to rule out about half of the initial pool of models, including all models with zero spin.Furthermore, the X-ray fluxes from quasi-simultaneous observations with the Chandra X-ray Observatory and NuSTAR provided another benchmark that disfavored several models based on preliminary estimates of X-ray emission from the simulations.However, detailed fitting of these and other precision data sets were beyond the focus of the first round of papers and GRMHD model sophistication.
The current cutting edge in modeling accreting black holes, whether via GRMHD simulations or semi-analytical methods, focuses on introducing more physically self-consistent, reliable treatments of the radiating particles (electrons or electron-positron pairs).In particular, key questions remain about how the bulk plasma properties dictate the efficiency of heating, how many thermal particles are accelerated into a nonthermal population, and the dependence of nonthermal properties such as spectral index on plasma properties such as turbulence, magnetization, etc. (see, e.g., treatments in Howes 2010; Ressler et al. 2015;Mościbrodzka et al. 2016;Ball et al. 2018;Anantua et al. 2020).
To test the newest generation of models, it is important to have extensive, quasi-simultaneous or at least contemporaneous multi-wavelength monitoring of several AGNs, providing both spectral and imaging data (and ideally polarization where available) over a wide range of physical scales.These types of campaigns have been limited by the difficulty in obtaining time on multiple facilities and by scheduling challenges, so often data are combined from different epochs.However, the variability timescales for even SMBHs such as M87's are short enough (days to months) that combining data sets from different periods can skew modeling results for sensitive quantities such as radiating particle properties.
This Letter intends to be the first of a series, presenting the substantive multi-wavelength campaigns carried out by the EHT Multi-Wavelength Science Working Group (MWL WG), including EHTC members and partner facilities, for both our primary sources M87 and Sgr A * , as well as many other targets and calibrators such as 3C 279, 3C 273, OJ 287, Cen A, NGC 1052, NRAO 530, and J1924−2914.These legacy papers are meant to be companion papers to the EHT publications, and will be used for the detailed modeling efforts to come, both from the EHTC Theory and Simulations WG as well as from other groups.They serve as a resource for the entire community, to enable the best possible modeling outcomes and a benchmark for theory.Here we present the results from the 2017 EHT campaign on M87, combining very long baseline interferometry (VLBI) imaging and spectral index maps at longer wavelengths, with spectral data from submillimeter (submm) through TeV γ-rays (covering more than 17 decades in frequency).
In Section 2 we describe these observations, including images (a compilation of MWL images in one panel is shown in a later section), spectral energy distributions (SEDs) and, where relevant, lightcurves and comparisons to prior observations.In Section 3 we present a compiled SED together with a table of fluxes.We also fit this SED with a few phenomenological models and discuss the consequences for the emission geometry and high-energy properties.In Section 4 we give our conclusions.All data files and products are available for download, as described in Section 3.2.

Observations and Data Reduction
In Figure 1 we give a schematic overview of the 2017 MWL campaign coverage on M87.In the following subsections we provide detailed descriptions of the observations, data processing procedures, and band-specific analyses.To aid readability, all tabulated data are collected in Appendix A.

Radio Data
In this subsection we describe the observations and data reduction of radio/mm data obtained with various VLBI facilities and connected interferometers.Especially regarding VLBI data, here we introduce the term "radio core" to represent the innermost part of the radio jet.A radio core in a VLBI jet The top and bottom panels are for connected interferometers and VLBI, respectively.The corresponding beam sizes are indicated in Table A1.KVN data at 22 and 43 GHz are not shown here since KVN captures the data from the shortest baselines of EAVN.
image is conventionally defined as the most compact (often unresolved or partially resolved) feature seen at the apparent base of the radio jet in a given map (e.g., Lobanov 1998;Marscher 2008).For this reason, different angular resolutions by different VLBI instruments/frequencies, together with the frequency-dependent synchrotron optical depth (Marcaide & Shapiro 1984;Lobanov 1998), can make the identification of a radio core not exactly the same for each observation.See also Section 3.2 for related discussions.M87 was observed with the European VLBI Network (EVN) at 1.7 GHz on 2017 May 9.The observations were made as part of a long-term EVN monitoring program of activity in HST-1, located at a projected distance of ∼70 pc from the core (Cheung et al. 2007;Chang et al. 2010;Giroletti et al. 2012;Hada et al. 2015).A total of eight stations joined a 10 hr long session with baselines ranging from 600 km to 10 200 km, yielding a maximum angular resolution of ∼3 mas at 1.7 GHz.The data were recorded at 1 Gbps with dual-polarization (a total bandwidth of 256 MHz, 16 MHz × 8 subbands for each polarization), and the correlation was performed at the Joint Institute for VLBI ERIC (JIVE).Automated data flagging and initial amplitude and phase calibration were also carried out at JIVE using dedicated pipeline scripts.This step was followed by frequency averaging within each spectral band (IF) and in time to 8 s.The final image was produced using the Difmap software (Shepherd 1997) after several iterations of phase and amplitude self-calibration (see Giroletti et al. 2012, for more detail).Here we provide the peak flux of the radio core (see Figure 2 and Table A1) and a large-scale jet image (presented in Section 3), while a dedicated analysis on the HST-1 kinematics will be presented in a separate paper.M87 was observed with the High Sensitivity Array (HSA) at 8.4, 15, and 24 GHz on 2017 May 15, 16 and 20, respectively, which are roughly a month after the EHT-2017 observations.Each session was made with a 12 hr long continuous track and the phased Very Large Array and Effelsberg 100 m antenna participated in the observations along with the 10 stations of the NRAO Very Long Baseline Array (VLBA).The data were recorded at 2 Gbps with dual-polarization (a total bandwidth of 512 MHz, 32 MHz × 8 subbands for each polarization), and the correlation was done with the VLBA correlator in Socorro (Deller et al. 2011).The initial data calibration was performed using the NRAO Astronomical Image Processing System (AIPS; Greisen 2003) based on the standard VLBI data reduction procedures (Crossley et al. 2012;Walker 2014).Similar to other VLBI data, images were created using the Difmap software with iterative phase/amplitude self-calibration.
A detailed study of the parsec-scale structure of the M87 jet from this HSA program will be discussed in a separate paper.
Here we provide a core peak flux and VLBI-scale total flux at each frequency (Table A1 in Appendix A).We adopt 10% errors in flux estimate, which is typical for HSA.

VERA 22 GHz
The core of M87 was frequently monitored over the entire year of 2017 at 22 GHz with the VLBI Exploration of Radio Astrometry (VERA; Kobayashi et al. 2003), as part of a regular monitoring program of a sample of γ-ray bright AGNs (Nagai et al. 2013).A total of 17 epochs were obtained in 2017 (see Figure 2 and Table A1 in Appendix A).During each session, M87 was observed for 10-30 minutes with an allocated bandwidth of 16 MHz, sufficient to detect the bright core and create its light curves.All the data were analyzed in the standard VERA data reduction procedures (see Nagai et al. 2013;Hada et al. 2014, for more details).Note that VERA can recover only part of the extended jet emission due to the lack of short baselines, so the total fluxes of VERA listed in Table A1 in Appendix A significantly underestimate the actual total jet fluxes.Since 2013 a joint network of the Korean VLBI Network (KVN) and VERA (KaVA; Niinuma et al. 2014) has regularly been monitoring M87 to trace the structural evolution of the pcscale jet (Hada et al. 2017;Park et al. 2019).From 2017, the network was expanded to the East Asian VLBI Network (EAVN; Wajima et al. 2016;Asada et al. 2017;An et al. 2018) by adding more stations from East Asia, enhancing the array sensitivity and angular resolution.Between 2017 March and May, EAVN densely monitored M87 for a total of 14 epochs (five at 22 GHz, nine at 43 GHz; see Figure 2 and Table A1 in Appendix A).The default array configurations were KaVA +Tianma+Nanshan+Hitachi at 22 GHz and KaVA+Tianma at 43 GHz, respectively, while occasionally a few more stations (Sejong, Kashima, and Nobeyama) additionally joined if they were available.In addition, we also had four more sessions with KaVA alone (2+2 at 22/43 GHz) in 2017 January-February.
Each of the KaVA/EAVN sessions was made in a 5-7 hr continuous run at a data recording rate of 1 Gbps (2-bit sampling, a total bandwidth of 256 MHz was divided into 32 MHz × 8 IFs).All the data were correlated at the Daejeon hardware correlator installed at KASI.All the EAVN data were calibrated in the standard manner of VLBI data reduction procedures.We used the AIPS software package for the initial calibration of visibility amplitude, bandpass, and phase calibration.The imaging/CLEAN (Högbom 1974) and selfcalibration were performed with the Difmap software.In Section 3 we present one of the 22 GHz EAVN images (taken in 2017 March 18) where the KaVA, Tianma-65 m, Nanshan -26 m, and Hitachi-32 m radio telescopes participated.
2. 1.5. KVN 22,43,86,and 129 GHz The KVN regularly observes M87 at frequencies of 22, 43, 86, and 129 GHz simultaneously via the interferometric Monitoring of γ-ray Bright Active galactic nuclei (iMOGABA) program, starting in 2012 December and lasting until 2020 January.The total bandwidth of the observations at each frequency band is 64 MHz and the typical beam sizes of the observations are 6.1 × 3.1 mas at 22 GHz, 2.8 × 1.6 mas at 43 GHz, 1.5 × 0.8 mas at 86 GHz, and 1.1 × 0.5 mas at 129 GHz.Details of the scheduling, observations, data reduction including frequency phase transfer technique, analysis (i.e., imaging and model-fitting), and early results for M87 are shown in Lee et al. (2016) and Kim et al. (2018a).Despite the limited coverage of baselines and capability of the array to image the extended jet structure in M87, the flux density of the compact core can be rather reliably measured (Kim et al. 2018a).We extract the core flux densities at the four frequencies and obtain light curves spanning observing periods between 2017 March and December at seven epochs.While typical flux density uncertainties at 22-86 GHz are of order of ∼10%, residual phase rotations and larger thermal noise at 129 GHz often lead to uncertainties of ∼30%.Accordingly, we adopt these uncertainties for all KVN observing epochs in this Letter.
2.1.6.VLBA 24 and 43 GHz M87 was observed with the VLBA at central frequencies of 24 and 43 GHz on 2017 May 5.These observations were carried out in the framework of the long-term monitoring program toward M87, which was initiated in 2006 (Walker et al. 2018) and lasted until 2020.For the 2017 session, the total on-source time amounts to about 1.7 hr at 24 GHz and 6 hr at 43 GHz.The sources OJ 287 and 3C 279 were observed to use as fringe finders and bandpass calibrators.In each band, eight 32 MHz wide frequency channels were recorded in both right and left circular polarization at a rate of 2 Gbps, and correlated with the VLBA software correlator in Socorro.The initial data reduction was conducted within AIPS following the standard calibration procedures for VLBI data.Deconvolution and self-calibration algorithms, implemented in Difmap, were used for phase and amplitude calibration and for constructing the final images.Amplitude calibration accuracy of 10% is adopted for both frequencies.
The resulting total intensity 43 GHz image of M87 is presented in Section 3, with the details given in Table A1 in Appendix A. The synthesized beam size amounts to 0.76 × 0.40 mas at the position angle (PA) of the major axis of −8°at 24 GHz, and 0.41 × 0.23 mas at PA = 0°at 43 GHz.We note that these observations were used for the study of a linear polarization structure toward the M87 core, details of which can be found in a separate paper (Kravchenko et al. 2020).

GMVA 86 GHz
M87 was observed by the Global Millimeter-VLBI-Array (GMVA; Boccardi et al. 2017) on 2017 March 30 (project code MA 009).In total, 14 stations participated in the observations: VLBA (eight stations; without HN and SC), 100 m Green Bank Telescope, IRAM 30 m, Effelsberg 100 m, Yebes, Onsala, and Metsahovi.The observation was performed in full-track mode for a total of 14 hr.Nearby bright sources 3C 279 and 3C 273 were observed as calibrator targets.The raw data were correlated by using the DiFX correlator (Deller et al. 2011).247Further post-processing was performed using the AIPS software package, following typical VLBI data reduction procedures (see, e.g., Martí-Vidal et al. 2012).After the calibration, the data were frequency-averaged across the whole subchannels and IFs, and exported outside AIPS for imaging with the Difmap software.Within Difmap, the data were further time-averaged for 10 s, followed by flagging of outlying data points (e.g., scans with too low amplitudes due to pointing errors).Afterward, CLEAN and phase self-calibrations were iteratively performed near the peak of the intensity, but avoiding CLEANing of the faint counterjet side at the early stage.When no more significant flux remained for further CLEAN steps, a first amplitude self-calibration was performed using the entire time coverage as the solution interval, in order to find average station gain amplitude corrections.A similar procedure was repeated with progressively shorter selfcalibration solution intervals, and the final image was exported outside Difmap when the shortest possible solution interval was reached and no more significant emission was visible in the dirty map compared to the off-source rms levels.
Calibrated visibilities are shown in Figure 3, and the CLEAN image is given in Section 3. We note that the final image has an rms noise level of ∼0.5 mJy beam −1 .This noise level is a factor of a few higher than other 86 GHz images of M87 from previous observations, which were made with similar array configuration as the 2017 session (see Hada et al. 2016;Kim et al. 2018b).Therefore, we refer to this image as tentative, and it only reveals the compact core and faint base of the jet, mainly due to poor weather conditions during the observations.We consider an error budget of 30% for the flux estimate.The peak flux density on the resultant image amounts to ∼0.52 Jy beam −1 for the synthesized beam of 0.243 × 0.066 mas at PA = −9°.3(see Table A1 in Appendix A for other details).We note that the peak flux density as well as the flux, integrated over the core region, are comparable to their historical values (Hada et al. 2016;Kim et al. 2018b), except for the 2009 May epoch, when the GMVA observations revealed about two times brighter core region in both intensity and integrated flux (Kim et al. 2018b).The observations at Band 6 with phased-ALMA (Matthews et al. 2018) were conducted as part of the 2017 EHT campaign (Goddi et al. 2019).The VLBI observations were carried out while the array was in its most compact configurations (with longest baselines <0.5 km).The spectral setup at Band 6 includes four spectral windows (SPWs) of 1875 MHz, two in the lower and two in the upper sideband, correlated with 240 channels per SPW (corresponding to a spectral resolution of 7.8125 MHz).The central frequencies at this band are 213.1, 215.1, 227.1, and 229.1 GHz.Details about the ALMA observations and a full description of the data processing and calibration can be found in Goddi et al. (2019).
Imaging was performed with the Common Astronomy Software Applications (CASA; McMullin et al. 2007) package using the task tclean.Only phased antennas were used to produce the final images (with baselines <360 m), yielding synthesized beam sizes in the range 1 0-2 4 (depending on the observing band and date).We produced 256 × 256 pixel maps, with a cell size of 0 2 yielding a field of view of 51″ × 51″.
The main observational and imaging parameters are summarized in Table A1 in Appendix A. For each data set and corresponding image, the table reports the flux-density values of the central compact core and the overall flux, including the extended jet.In order to isolate the core emission from the jet, we compute the sum of the central nine pixels of the model (CLEAN component) map (an area of 3 × 3 pixels); 248 the contribution from the jet is accounted for by also summing the clean components along the jet.The extended emission accounts for less than 20% of the total emission at 1.3 mm.The large-scale jet image is presented in Section 3. Details about the imaging and flux extraction methods can be found in Goddi et al. (2021).
2.1.9.Submillimeter Array (SMA) 230 GHz The long term 1.3 mm band (230 GHz) flux density lightcurve for M87 shown in Figure 2 was obtained at the SMA near the summit of Maunakea (Hawaii).M87 is included in an ongoing monitoring program at the SMA to determine flux densities for compact extragalactic radio sources that can be used as calibrators at mm wavelengths (Gurwell et al. 2007).Available potential calibrators are occasionally observed for 3-5 minutes, and the measured source signal strength calibrated against known standards, typically solar system objects (Titan, Uranus, Neptune, or Callisto).Data from this program are updated regularly and are available at the SMA Observer Center website (SMAOC 249 ).Data were primarily obtained in a compact configuration (with baselines extending from 10 to 75 m) though a small number were obtained at longer baselines up to 210 m.The effective spatial resolution, therefore, was generally around 3″.The flux density was obtained by vector averaging of the calibrated visibilities from each observation.
Observations of M87 were additionally conducted as part of the 2017 EHT campaign, with the SMA running in phasedarray mode, operating at 230 GHz.All observations were conducted while the array was in compact configuration, with the interferometer operating in dual-polarization mode.The SMA correlator produces four separate but contiguous 2 GHz spectral windows per sideband, resulting in frequency coverage of 208-216 and 224-232 GHz.Data were both bandpass and amplitude calibrated using 3C 279, with flux calibration performed using either Callisto, Ganymede, or Titan.Due in part to poor phase stability at the time of observations, phase calibration is done through self-calibration of the M87 data itself, assuming a point-source model.Data are then imaged and deconvolved using the CLEAN algorithm.
A summary of the measurements made from these data are shown in Table A1 in Appendix A, along with measurements taken within a month of these observations from the SMAOC data mentioned above.The reported core flux for M87 is the flux measured in the center of the cleaned map.Combined images of all data show the same jet-like structure seen in the ALMA image in Figure 13, although recovery of the flux for individual days through imaging is limited by a lack of (u, v)coverage within individual tracks.Therefore, we estimate the total flux by measuring the mean flux density of all baselines within a (u, v)-angle of 110°± 5°, as these baselines are not expected to resolve the jet and central region.

Optical and Ultraviolet (UV) Data
We performed optical and UV observations of M87 with the Neil Gehrels Swift Observatory during the EHT campaign, and have also analyzed contemporaneous archival observations from the Hubble Space Telescope (HST).

UV/Optical Telescope (UVOT) Observations
The Neil Gehrels Swift Observatory (Gehrels et al. 2004) is equipped with UVOT (Roming et al. 2005) Poole et al. (2008).The reduction of the data followed the standard prescriptions of the instrument team at the University of Leicester. 251We checked the UVOT observations for small-scale sensitivity issues and found that none of our observations are affected by bad charge-coupled device (CCD) pixels.In 2020 November the Swift team released new UVOT calibration files, along with coefficients of the flux density correction as a function of time, for data reduced with the previous version of CALDB.For our observations the coefficients of the flux density correction (multiplicative factors) are very close to 1: 0.974 ± 0.001 (for the optical filters), 0.947 ± 0.001 (uvw1), 0.964 ± 0.001 (uvm2), and 0.958 ± 0.001 (uvw2).We have used these coefficients to correct our measurements for the UVOT sensitivity change. 252 We performed aperture photometry for each individual observation using the tool UVOTSOURCE, with a circular aperture of a radius of 5″ centered on the sky coordinates of M87 and detection significance σ 5.This aperture includes the M87 core, the knot HST-1, and some emission from the extended jet, which we cannot separate due to the size of the UVOT point-spread function (PSF; ∼2 5).Since there is contamination from the bright host galaxy surrounding the core of M87, we measure the background level in three circular regions, each of 30″ radius in a source-free area located outside the host galaxy radius (see below).We used the count-rate to magnitude and flux density conversion provided by Breeveld et al. (2011) and retained only those measurements with magnitude errors of σ mag < 0.2.This calibration of the UVOT broadband filters (Breeveld et al. 2011) includes additional calibration sources with a wider range of colors with respect to the one reported by Poole et al. (2008).
To estimate the contributions of the host galaxy to the derived flux densities, we modeled the UVOT images.First we combined individual images using the tool UVOTIMSUM to produce a stacked image in each filter.Ordinarily, this would have allowed us to obtain the best signal-to-noise ratio (S/N) for modeling.However, since the source was not centered in the different images at the same CCD position, the area of intersection between the images was too small.Hence, even in the outer regions the flux from the galaxy was still significant.Therefore, we made a set of models using individual images and then found a mean model for each band after a mask of background objects was used to exclude all objects (except for M87 itself).Since a unified mask was used, the images have the same pixels included and excluded from the analysis, so no bias is introduced.
After this procedure, the decomposition process was run for each image using a model consisting of three components: the core region, jet, and host galaxy.The core region includes the core and knot HST-1, which was modeled by a point source, while the jet was fitted by a highly elliptical Gaussian.The galaxy was modeled by a Sérsic function: , where n is the Sérsic parameter, b n = 2n − 1/3, R e is the halflight radius, and I e is the intensity at R e .All three structural features can be seen in the combined image in the uvw1 band presented in Figure 13.A point source model for the core + HST-1 region was convolved with the PSF determined for each filter using several isolated non-saturated stars in a number of images.The parameters of each PSF were averaged over stars and images.Before the decomposition, the images were normalized to a one-second exposure to make them uniform.The program imfit (Erwin 2015) was employed to perform decomposition for each image.Each model parameter for a given filter corresponds to the median value averaged over all images using the best-fit image parameters (according to χ 2 ), with >2σ outliers removed (number of outliers for each filter 4).

Table A2 in Appendix
A lists the derived median values of the Sérsic model n-parameter and its uncertainty for different bands.The uncertainty corresponds to a scatter of models among the individual images.The derived values of the Sérsic n-parameter show a dependence on wavelength that is caused by a difference in the stellar population and, probably, the scattered light from the core, which is blue, leading to higher nvalues for blue filters.Note that there is significant scatter among values of the n-parameter reported in the literature, from n = 2.4 (Vika et al. 2012) to n = 3.0 (D'Onofrio 2001), n = 6.1 (Ferrarese et al. 2006), n = 6.9 (Graham & Driver 2007), and n = 11.8 (Kormendy et al. 2009).
Table A2 in Appendix A also gives the derived values of the effective radius of the galaxy in different bands, R e , ellipticity, ò, and PA, Φ, of the major axis of the galaxy calculated counter clockwise from north. Figure 4 shows an example of a comparison between the mean observed surface brightness profile along the semimajor axis of the host galaxy and the results of the modeling.The mean observed profile is obtained by azimuthally averaging along a set of concentric ellipses with the ellipticity and PA equal to those of the host galaxy.
According to our X-ray data analysis (see Section 2.3.3), the hydrogen column density corresponding to absorption in both our Galaxy and near M87 is equal to N H = -+ 0.050 10 0.002 0.003 22 cm −2 , while there is evidence from our X-ray data for additional X-ray absorption within the central 1″ around M87 with a column density of -+ 0.12 10 0.04 0.05 22 cm −2 .However, this latter value is most likely variable since much lower values of N H for M87 were detected previously, e.g., N H ≈ 0.01 × 10 22 cm −2 (Sabra et al. 2003) based on UV spectroscopy.Therefore, as discussed in Prieto et al. (2016) we assume that the additional X-ray absorber is dust-free and employ the extinction curve (R V = 3.1) and the extinction value (E B−V = 0.022) given by Schlegel et al. (1998) along with the formalism provided by Cardelli et al. (1989) to derive the  extinction values in different bands, which are listed in Table A2 in Appendix A.
Tables A3 and A4 in Appendix A give flux densities of the core region of M87 corrected for the host galaxy contamination and extinction in optical and UV bands, respectively.Figure 5 shows light curves of the core region during the campaign.According to Table A2 in Appendix A and as can be seen from Figure 5 the standard deviation of F core averaged over the EHT campaign is significantly less than the uncertainty of an individual measurement in all filters.The uncertainty is dominated by that of the host galaxy decomposition, which is added (in quadrature) to the photometric uncertainty of each measurement.Therefore, based on the UVOT observations we cannot detect variability in the core region of M87 during the EHT campaign, which is consistent with apparently a low activity state of both the core and HST-1 knot.Based on these results we have calculated the UVOT spectral index of the core region during the EHT campaign using the average values of the flux density in the UVOT bands as S ∝ ν − α , which results in α = 1.88 ± 0.55.Although α has a significant error, most likely connected with large uncertainties in the host galaxy decomposition in different bands, the spectral index is consistent with the optical/UV spectral index of the core of M87 reported by Perlman et al. (2011).This indicates that the core dominates the innermost region of M87 at UV/optical wavelengths during the campaign.

HST Observations
We have downloaded HST images of M87 from the HST archive obtained on 2017 April 7, 12, and 17 with the Wide Field Camera 3 (WFC3) camera in two wide bands, F275W and F606W (ID5o30010, ID5o30020, ID5o31010, ID5o31020, ID5o32010, and ID5o32020).We used fully calibrated and dither combined images.The image obtained on 2017 April 12 in the F606W filter is a part of a composite of M87 multiwavelength images presented in Section 3. The decomposition of the HST images was made using the same imfit package as for the images obtained with the UVOT (Section 2.2.1).The substantially higher spatial resolution of the HST images (1 pixel is 0 04), however, led to a somewhat different approach.First, we masked out the jet on the images: the HST images reveal its very complex structure, which cannot be approximated with a simple analytical model, as was the case for the low-resolution UVOT images.We also masked out knot HST-1 and the core, which are clearly resolved in the HST images (Figure 6).Second, we used a more complex model for the host galaxy: a sum of two Sérsic models (instead of one as in the case of the UVOT images), which gave us significantly better residual (observations-model) maps.Finally, we used the HST library of PSF images provided at the HST website253 to model the PSF in each filter.Figure 6 shows the resulting map of M87 on 2017 April 11 in F606W filter, with the host galaxy subtracted.
The decomposition analysis was limited to the central region of the images (16″ × 16″), since the main goal of the process was to subtract the host galaxy flux before the aperture photometry of the core and HST-1.Fitting the full image of the galaxy would have required a more complex model to have comparable quality of the fit for the central regions (see, e.g., Huang et al. 2013).After subtracting the best-fit model from the images, we performed the aperture photometry with a radius of 0 4 to find the flux densities from the core and HST-1 in each band.imfit allows one to perform a bootstrap method to derive estimates of uncertainties of the decomposition parameters.For each image, the program was run 100 times, each with same number of random pixels involved in the decomposition process.For each decomposition, a residual FITS file was constructed containing only the core and HST-1 (100 cases for each band and date).For each such residual file, we performed the aperture photometry with a radius of 0 4 and calculated the standard deviation of flux density measurements of the core and HST-1 over different decompositions.This standard deviation was added in quadrature to the uncertainty of the photometry from the decomposition of the original image in each band and date.Table A5 in Appendix A gives flux densities of the core and HST-1, measured in two different bands and at epochs contemporaneous with the EHT observations.The flux densities are corrected for the extinction in the same manner as described in Section 2.2.1.
According to Table A5 in Appendix A the core shows a slight increase in flux density over 10 days, while knot HST-1 has a constant flux density.We have estimated optical/UV spectral indices of the core and HST-1, which are 1.44 ± 0.09 and 0.60 ± 0.02, respectively (the spectral index is defined in the same way as in Section 2.2.1).These are in good agreement with those given by Perlman et al. (2011), ∼1.5 for the core and ∼0.5 for HST-1, which confirm a flatter spectral index of HST-1 with respect to that of the core at UV/optical wavelengths.
To determine the activity state of M87 during the 2017 campaign, we have constructed light curves of the core and HST-1 in two bands, F275W (from 1999 to 2017) and F606W (from 2002 to 2017).During the period 1999 to 2010 different instruments were used at HST: UV observations were performed with STIS/F250QTZ λ eff = 2365 Å, ACS/F220W λ eff = 2255.5Å, and ACS/F250W λ eff = 2716 Å (e.g., Madrid 2009).We have used the UV measurements presented in Madrid (2009) and translated them into WFC3/F275W using spectral indices reported by Perlman et al. (2011), who observed M87 in four UV/optical bands during the same period.In addition, Madrid (2009) performed photometry with an aperture of radius 0 25, so that to construct a uniform UV lightcurve, we have recalculated our measurements in F275W band using the same aperture.For the optical lightcurve in F606W band, we have used measurements provided by Figure 6.HST image of M87 in F606W filter with host galaxy subtracted; the core and HST-1 are designated, the distance between the features is 0 86 ± 0 04.Perlman et al. (2011) obtained with Advanced Camera for Surveys (ACS)/F606W and Wide Field Camera 2 (WFC2)/ F606W.Since λ eff of these instruments in F606W band is the same as for WFC3/F606W and photometry was performed with the same aperture of radius 0 4, no corrections to the measurements have been applied.Figure 7 shows the historical UV/optical light curves of M87.Independent of the filter used, HST-1 was in its lowest brightness state ever observed during the EHT campaign.Although the core was in its lowest brightness state at optical wavelengths in 2017, the bettersampled UV lightcurve suggests that the core was in an average quiescent state.

Chandra
We requested Director's Discretionary Time (DDT) observations of M87 with the Chandra X-ray Observatory to coordinate with the EHT campaign.The source was observed with Advanced CCD Imaging Spectrometer (ACIS)-S for 13.1 ks starting on 2017 April 11 23:46:58 UT (ObsID 20034) and again starting on 2017 April 14 02:00:28 UT (ObsID 20035).In order to perform spatially resolved spectral and variability analysis (see Section 2.3.3),we also analyzed several archival Chandra observations.Following Wong et al. (2017), for constraints on the intra-cluster medium (ICM), we included ObsIDs 352, 3717, and 2707, which were acquired on 2000 July 29, 2002 July 5, and 2004 July 6 and have good exposures  of 37.7 ks, 98.7 ks, and 20.6 ks, respectively.The Chandra observations were processed using standard data reduction procedures in CIAO v4.9. 254We focused on extracting spectra from the core, the knot HST-1, and the outer jet, along with instrumental response files for spectral analysis.We took positions for the core and HST-1 from Perlman & Wilson (2005).For the core, we used a circular source extraction region with a radius of 0 4 centered on M87 (the approximate FWHM of 0 8 is quoted in Table A8 in Appendix A).The core background region is a half annulus centered on M87 with inner and outer radii of 2″ and 3 5, respectively; we excluded the half of the annulus that is on the same side of the core as the extended X-ray jet.For HST-1, we used a similar circular source region, but for the background annulus we excluded ∼90°wedges containing the core on one side and the extended jet on the other.For the jet itself, we used a 19 5 × 3″ rectangular source region centered on the jet, with 19 5 × 1 5 rectangular regions on either side.To illustrate the relative brightness and variability of the core and HST-1, we show their lightcurves (1 ks bins) in Figure 8.Sun et al. (2018) presented a Chandra study spanning ∼8 yr from 2002 to 2010 with coverage of the core and HST-1 in low and high states (see their Figure 3).In their study, the core flux drops as low as ∼10 −12 erg s −1 cm −2 (the average is ∼4 × 10 −12 erg s −1 cm −2 ).In the 2017 Chandra observations, the unabsorbed core flux in the 0.3-7 keV band is 3 × 10 −12 erg s −1 cm −2 , and the absorbed flux is 2 × 10 −12 erg s −1 cm −2 -hence our observations show the core below the historical mean.
Because the ICM contributes significantly to the NuSTAR background, we used a single set of extraction regions to produce the Chandra ICM spectrum and the NuSTAR spectra (circular regions of radius 45″, see Section 2.3.2 for more details).In extracting the Chandra ICM spectrum, we excluded the source regions for the core, HST-1, and the extended jet, all of which fit well within the NuSTAR PSF (Section 2.3.2).

NuSTAR
We also requested two DDT observations of M87 with NuSTAR (Harrison et al. 2013) to coordinate with the EHT campaign in 2017 April.These observations are contemporaneous with the Chandra observations described in Section 2.3.1, and they are summarized in Table A6 in Appendix A. NuSTAR observations of M87 in 2017 February and April have been presented in Wong et al. (2017).
Raw data from NuSTAR observations were processed using standard procedures outlined in the NuSTAR data analysis software guide (Perri et al. 2017).We used data analysis software (NuSTARDAS, version 1.8.0), distributed by HEA-SARC/HEASOFT, version 6.23.Instrumental responses were calculated based on HEASARC CALDB version 20180312.We cleaned and filtered raw event data for South Atlantic Anomaly (SAA) passages using the nupipeline script; both minimal and maximal filtering levels were considered, with no substantial differences in the results.We used a source extraction region that is a circle with 45″ radius, while the background region was chosen as an identical region well separated from the peak of the hard X-ray emission (which is a point source above 12 keV, per Wong et al. 2017).Note that this does not include all of the low-energy extended emission resolved with Chandra (see Section 2.3.1),nor does it subtract it out.However, we modeled the low-energy (3 keV) X-ray and high-energy (8 keV) X-ray emission jointly, as described in Section 2.3.3.Source and background spectra were computed from the calibrated and cleaned event files using the nuproducts tool.

Joint Chandra and NuSTAR Spectral Analysis
Given the spatial complexity of the underlying X-ray emission, any detailed analysis of the X-ray spectrum of M87 must involve a joint treatment of Chandra and NuSTAR data.Here we discuss our strategy for this analysis.
To properly recover the intrinsic spectrum of the core in M87 up to 40-50 keV, it was necessary to either subtract or model the bright emission of the ICM.Given our interest in the core and HST-1, which are moderately bright point sources, we must correct our spectra for pileup (see Section 2.3.1).Since this correction depends on the total count rate per pixel per frame time (Davis 2001), it must be applied without background subtraction.Thus, we chose to model the ICM spectrum.
We also needed to account for variations between different detectors and epochs.It is common practice to include a crossnormalization constant between NuSTAR's FPMA and FPMB detectors.We opted to take the same approach when jointly modeling NuSTAR and Chandra data: we allowed additional cross-normalization constants between the Chandra spectrum and the NuSTAR spectrum.Because the archival Chandra observations we used to constrain the ICM spectrum are nearly 20 yr old (and because of Chandra's changing effective area at soft X-ray energies due to contaminant buildup255 ), we allowed those data to have a different cross-normalization constant than the 2017 data.
Schematically, then, we can represent the model for the NuSTAR spectra of M87 as follows: where A, B, C, and D are cross-normalization constants, F Chandra,ICM is the model for the Chandra ICM emission, and similar notation holds for the Chandra spectra of the core, the nearby knot HST-1, and the rest of the jet.Again, the purpose of this procedure is to use Chandra to account for all of the X-ray emission inside the NuSTAR extraction region, while allowing for time-dependent normalization constants between missions and detectors.
With the structure of our model defined, we proceeded to select models for each of the components.The ICM has a complex physical structure and three-dimensional temperature profile.Our purpose here was not to model the cluster gas physics but to reconstruct the ICM spectrum with sufficient accuracy to model the X-ray continuum from M87 itself.Following Wong et al. (2017), we adopted a two-temperature Astrophysical Plasma Emission Code (APEC) model (Smith et al. 2001) with variable abundances (vvapec).We allowed the Ne, Na, Mg, Al, Si, S, Ar, Ca, and Fe abundances to vary and included three Gaussian emission lines; note that this treatment differs slightly from Wong et al. (2017), who used solar metallicities and 5 Gaussians to account for their residuals relative to the APEC models.We also included a cutoff powerlaw to model the low-mass X-ray binary (LMXB) emission in the extraction region with photon index Γ = 0.5 and cutoff energy E c = 3 keV (Revnivtsev et al. 2014).This model was modified by interstellar absorption; we use the model tbnew (Wilms et al. 2000) with atomic cross-sections from Verner & Yakovlev (1995).
For the core, HST-1, and the remainder of the jet, we modeled their continuum emission using power laws (modified by the same interstellar absorption component as the ICM emission).However, inspection of the spectra in the top panel of Figure 9 revealed a steeper rollover in the soft X-ray spectrum of the core (orange) than in HST-1 or the outer jet.This is typical of ISM absorption, so we included a second tbnew component in the model for the spectrum of the core (see Wilson & Yang 2002;Perlman & Wilson 2005, but also Di Matteo et al. 2003).As noted above, these spectra required a pileup correction; the "grade migration" parameter α is tied between the two contemporaneous Chandra spectra of each spatial component.
While the Chandra data effectively disentangle the spatial components of M87, the NuSTAR data are superior when it comes to constraining the photon index, given their wider energy range.But because the count rates are low, we allowed the model X-ray flux to vary between observations while assuming the photon index was constant.
For fitting, the Chandra spectra were binned to a minimum S/N of 3 between 0.4 and 8 keV.To maximize the useful energy range of the NuSTAR data, we grouped the spectrum from each FPM to a minimum S/N of only 1.1; we also rebinned the spectra by factors of 3, 4, 5, 6, and 7 over the energy ranges 3-51 keV, 51-59 keV, 59-67 keV, 67-74.9keV, and 74.9-79 keV, respectively (this reduces oversampling of the energy resolution; J. Steiner 2021, private communication).Because very low count rates are involved for some energy bins, we used Cash statistics (Cash 1979).Our best fit gave a reduced Cash statistic of 1.26 (a Cash statistic of 1392 for 1139 data bins and 38 free parameters).Our power-law model does not rule out possible curvature in the high-energy X-ray spectrum.
Given the possibility of complex correlations between parameters in this model, we opted to determine confidence intervals for our parameters using Markov Chain Monte Carlo (MCMC) methods.We used an implementation of the affineinvariant sampler emcee, after Goodman & Weare (2010) and Foreman-Mackey et al. (2013), with 10 walkers for each of the 38 free parameters and let the sampler run for 20,000 steps.
At present, we focus on the spectral index, flux, and ISM absorption for each of our spatial components; additional parameters shall be described below as necessary.Unless otherwise noted, quoted uncertainties represent 90% credible intervals; these generally align well with our 90% confidence intervals from direct fitting.
Our results confirm our inspection of the Chandra spectrum of the core of M87: it is more highly absorbed than the rest of the jet.All components required a small absorbing column density = -  Wilson & Yang (2002) and Perlman & Wilson (2005).If it is real, it is most likely absorption in the galaxy M87 itself, though it is difficult to speculate on its properties without knowing more about its variability, which will be discussed in future work (J.Neilsen et al. 2021, in preparation).An alternative interpretation is that there is measurable curvature in the core spectrum at these energies.
The flux itself is not a parameter of our fits, but it is easily calculated from the normalization of the power law.We drew 1000 samples from our chain, excluding the first 4000 steps to account for the burn-in period.For each sample we used the power-law photon index and normalization for each spectral component to calculate the 2-10 keV luminosity, assuming a distance of 16.8 Mpc as in EHT Collaboration et al. (2019a).

Swift-X-Ray Telescope (XRT)
M87 was observed using Swift-XRT from 2017 March 27 to April 20.A total of 24 observations in the photon counting mode were conducted, out of which 12 were discarded due to the fact that M87 was centered on the bad columns, which resulted from a micrometeorite impact in 2005.For the purpose of this work, we analyzed data for the remaining 12 epochs.First, all the cleaned level 3 event files were generated using xrtpipeline version 0.13.5.For spectral extraction, a circular source region with a diameter of 15 pixels (35″) and an annular background region were chosen.If the source counts exceeded 0.5 cts s −1 , a pile-up correction was performed.In this case, an appropriate annular region was selected as the source region for the final spectrum extraction, ensuring that the piled-up region is removed from the final extraction process.The exposure maps created by the xrtpipeline were utilized to create the Ancillary Response File (arf) for each spectrum employing the xrtmkarf command.The Response Matrix File (rmf) used in this process was later used to group all these spectral files for each epoch, using the grppha command.The spectral resolution of Chandra provided well-constrained model parameters for M87 data obtained during the same time as Swift-XRT.The composite model contained a complex background model, which is a sum of the ISM, ICM, and any remaining background.All the details for this model, in addition to the details for the power-law models for the core, HST-1 and the jet, respectively, are described in the previous section.These Chandra derived parameters were used as a guide to fit the spectra obtained from Swift-XRT, freezing all parameters other than (1) the overall normalization of the power-law continuum components (spectral slopes and relative normalizations remained fixed) and (2) the normalization of the background component, vvapec.The overall variability seen in Figure 10 contains contributions from both the power-law continuum and background component variations.The variability of the background component may be an artifact of our assumptions about the cross-normalization of the Swift-XRT models relative to Chandra and NuSTAR.It could also indicate some residual systematic uncertainty in the spectral decomposition (background versus core, jet, and HST-1) that is not accounted for in our error bars.In any case, the non-ICM fluxes represent an upper limit on the core X-ray emission for our SED modeling, and so this precise decomposition does not have a significant impact on our conclusions about the SED.
The overall evolution of the lightcurve for M87 in 2017 is displayed in Figure 10.We derived an overall total unabsorbed flux using the spectral fits for each epoch.We also calculated the net flux resulting from the core, HST-1, and outer jet only from these fits by removing the background model from the final fit, which helped estimate the combined flux only from these three components of M87.Both these flux values are reported in Figure 10.A day to day variability of the order of about ∼5%-20% is seen during this time.However, an overall variation in the flux by about a factor of 3 was seen during 2017.Note that the spatial location of the variability can not be provided due to the lower angular resolution of Swift-XRT.The γ-ray emission at GeV energies from M87 was first reported in 2009, using the first 10 months of observations with Fermi-LAT (Abdo et al. 2009a).Fermi-LAT (Atwood et al. 2009) mainly operates in survey mode, observing the whole sky every three hours.Connected to the EHT campaign, during the period 2017 March 22 to April 20, pointing mode observations256 of M87 have been specially performed providing additional 500 ks of data on this source, allowing the source to be effectively observed >56% of the time (w.r.t. the standard ∼48%) during the EHT campaign.
We performed a dedicated analysis of the Fermi-LAT data of M87 using 11 yr of LAT observations taken between 2008 August 4 and 2019 August 4. Similarly, we repeated the analysis using three-month time bins centered on the EHT observation period (i.e., 2017 March 1-May 31).We selected P8R3 Source class events (Bruel et al. 2018), in the energy range between 100 MeV and 1 TeV, in a region of interest (ROI) of 20°radius centered on the M87 position.The lowenergy threshold is motivated by the large uncertainties in the arrival directions of the photons below 100 MeV, leading to a possible confusion between point-like sources and the Galactic diffuse component.See Principe et al. (2018Principe et al. ( , 2019) ) for a different analysis implementation to solve this and other issues at low energies with Fermi-LAT.
The analysis (which consists of model optimization, and localization, spectrum, and variability study) was performed with Fermipy257 (Wood et al. 2017), a Python package that facilitates analysis of LAT data with the Fermi Science Tools, of which the version 11-07-00 was used.The counts maps were created with a pixel size of 0°. 1.All γ-rays with zenith angle larger than 95°were excluded in order to limit the contamination from secondary γ-rays from the Earth's limb (Abdo et al. 2009b).We made a harder cut at low energies by reducing the maximum zenith angle and by selecting event types with the best PSFs. 258For energies below 300 MeV we excluded events with zenith angle larger than 85°, as well as photons from PSF0 event type, while above 300 MeV we use all event type.The P8R3_Source_V2 instrument response functions (IRFs) are used.The model used to describe the sky includes all point-like and extended LAT sources, located at a distance <25°from the source position, listed in the Fourth Fermi-LAT Source Catalog (4FGL; Abdollahi et al. 2020), as well as the Galactic diffuse and isotropic emission.For these two latter contributions, we made use of the same templates259 adopted to compile the 4FGL.For the analysis we first optimized the model for the ROI (fermipy.optimize),then we searched for the possible presence of new sources (fermipy.find_sources)and finally we re-localized the source (fermipy.localize).We investigated the possible presence of additional faint sources, not in 4FGL, by generating test statistic260 (TS) maps and we found two candidate new sources, detected at 5σ level, that we added into our model.The best-fit positions of these new sources are R.A., decl.= (184°.85, 5°.82) and (191°.99, 7°.20), with 95% confidence-level uncertainty R 95 = 0°.07.We left free to vary the diffuse background and the spectral parameters of the sources within 5°of our target.For the sources at a distance between 5°and 10°only the normalization was fitted, while we fixed the parameters of all the sources within the ROI at larger angular distances from our target.The spectral fit was performed over the energy range from 100 MeV to 1 TeV.To perform a study of the γ-ray emission variability of M87 we divided the Fermi-LAT data into time intervals of 4 months.For the lightcurve analysis we have fixed the photon index to the value obtained for 11 yr and left only the normalization free to vary.The 95% upper limit has been reported in each time interval with TS < 10.

H.E.S.S., MAGIC, VERITAS Observations
In 1998 the first strong hint of very-high energy (VHE; E > 100 GeV) γ-ray emission from M87 was measured by the High Energy Gamma Ray Astronomy (HEGRA) Collaboration (Aharonian et al. 2003).Since 2004, the source has been frequently monitored (Acciari et al. 2009;Abramowski et al. 2012) (Abramowski et al. 2012) that lasted one to a few days.A weak two-month VHE enhancement was also reported for 2012 (Beilicke & VERITAS Collaboration 2012), but since 2010, no major outburst has been detected in VHE γrays.During the observed high states the flux typically increased by factors of two to 10, while no significant spectral hardening was detected.In the following we describe the participating instruments, observations, and data analyses of the coordinated 2017 MWL campaign.
The H.E.S.S. observations presented here were performed using stereoscopic observations with the four 12 m CT1-4 telescopes (see Table A7).These observations of a total of 7.9 hr of live time were taken in so-called "wobble" mode with an offset of 0°. 5 from the position of M87, allowing simultaneous background estimation.The reconstruction of the Cherenkov shower properties was done using the ImPACT maximum likelihood-based technique (Parsons & Hinton 2014;Parsons et al. 2015) and hadron events were rejected with a boosted decision tree classification (Ohm et al. 2009).The background was estimated using the so-called "ring-background model" for the signal estimation and the "reflectedbackground model" for the calculation of the flux and the spectrum (Aharonian et al. 2006b).For the 7.9 hr of data a total statistical significance of 3.7σ was calculated.A source with 1% of the Crab Nebula's flux can be detected in ∼10 hr.Based on the estimated systematic uncertainties following Aharonian et al. (2006b), the systematic uncertainty of the flux has been adopted to be 20%.
For this study MAGIC observed M87 for a total of 27.2 h after quality cuts with the standard wobble offset of 0°.4.The data were analyzed with the standard MAGIC software framework MAGIC Analysis and Reconstruction Software (MARS; Zanin et al. 2013;Aleksić et al. 2016).We used the package SkyPrism, a spatial likelihood analysis, to determine the PSF in true energy (Vovk et al. 2018).6.7 hr of these data were collected in the presence of moonlight and analyzed according to Ahnen et al. (2017).All MAGIC observations together yielded a total significance of 4.6σ.MAGIC can detect a source with 1% of the Crab Nebula's flux in ∼26 hr above an energy threshold of 290 GeV.A detailed description of the various sources and estimates of systematic uncertainties in the MAGIC telescopes and analysis can be found in Aleksić et al. (2016).From there we estimate a systematic flux normalization uncertainty of 11%, a systematic uncertainty on the energy scale of 15% and a systematic uncertainty of ±0.15 on the reconstructed spectral slope for the MAGIC observations, which sums up to a total of ∼30% integral flux uncertainty.
VERITAS collected 15 hr of quality-selected observations of M87 during the 2017 EHT campaign window using a 0°.5 wobble.Data were analyzed using standard analysis tools (Cogan 2007;Daniel 2008;Maier & Holder 2017) with background-rejection cuts optimized for sources with spectral photon indices in the 2.5-3.5 range.The analysis yields an overall statistical significance of 3.8σ.In its current configuration, VERITAS can detect a source with a flux of 1% of the γray flux of the Crab Nebula within 25 hr of observation (Park 2015).We estimate a systematic uncertainty on the quoted photon flux of 30%.
Separate differential spectra and differential upper limits were produced assuming a power law of the form µ -G dN dE E for each instrument.For the different instruments the γ-ray flux is measured in five independent bins per energy decade, respectively, with flux points being quoted for bins with a significance larger than 2σ.A differential upper limit at the 95% confidence level following the method described in Rolke et al. (2005) is quoted otherwise.The bin edges vary among the different Imaging Atmospheric Cherenkov Telescope (IACT) measurements due to differences in instrumental and observational conditions.Night-wise integrated flux points and upper limits were calculated for the light curves above an energy threshold of 350 GeV assuming a differential spectrum that follows a simple power law with a spectral index of Γ = 2.3.Due to technological limits, and ultimately limited by fluctuations in the development of air showers, IACTs are not able to spatially resolve M87 and so the measurements are compatible with a point source (Hofmann 2006).
A summary of the individual observation nights can be found in Table A7 in Appendix A. Figure 11 shows the measured fluxes for each observation night with 1σ statistical uncertainty while the flux upper limits for nights with <2σ are given in Table A7 in Appendix A. All measurements are compatible with constant emission within uncertainties at an overall low state during the 2017 campaign (MAGIC Collaboration et al. 2020 and references therein).The resulting SEDs assuming that the intrinsic flux is constant for each instrument are shown in Figure 12 together with historical measured spectra.The index of a power-law fit to the MAGIC data is Γ = 2.17 ± 0.46.This is in good agreement with previous results obtained during low-flux states MAGIC Collaboration et al. (2020).The data from H.E.S.S. and VERITAS do not allow one to perform a measurement of the VHE γ-ray spectral shape of M87, but yield upper limits and energy flux measurements that are consistent with the spectrum determined with the MAGIC data.

Multi-wavelength Images and Core Shift
In Figure 13 we show a composite of M87 multi-wavelength images obtained by various instruments during the 2017 campaign, including the EHT image.The M87 jet is imaged at all scales from ∼1 kpc down to a few Schwarzschild radii, summarizing a contemporaneous MWL view of M87 during the 2017 EHT campaign.
The large-scale jet structure is well imaged with ALMA, HST, and Chandra.The observed kpc-scale jet morphology in 2017 is overall consistent with that known in the literature: a straight, highly collimated jet continues along PA ∼ 290°with several discrete knots (e.g., Owen et al. 1989;Biretta et al. 1999;Snios et al. 2019).The near-core knot HST-1 is well isolated from the nucleus by EVN, HST, and Chandra.The EVN at 1.7 GHz further resolves the inner jet at ∼10-100 pc scales, revealing a continuous, straight jet morphology.
At parsec scales, high-resolution VLBI observations reveal the detailed transverse structure of the jet, resolving a limbbrightened structure that follows a parabolic collimation profile (Asada & Nakamura 2012;Hada et al. 2013;Nakamura & Asada 2013;Kim et al. 2018b).Closer to the black hole, the jet morphology may change more rapidly than at large scales (e.g., Britzen et al. 2017;Walker et al. 2018).Walker et al. (2018) reported a significant oscillation of the inner-jet position angle in their VLBA 43 GHz images on timescales of several years.During the 2017 campaign, our VLBA and EAVN 43 GHz images (see also Figure 14) reveal a jet direction within 2 mas to be close to east-west (PA ∼ 270°), significantly offset from the large-scale jet direction (PA ∼ 290°).This inner-jet direction observed in 2017 seems to follow the long-term periodic oscillation found by Walker et al. (2018).A consistent jet direction can also be seen in the GMVA 86 GHz image obtained in 2017 March, although the image quality is poorer than that of VLBA at 43 GHz.If the jet is conically precessing, an apparent change of ∼20°in PA implies a variation of jet inclination by Δθ ∼ PA × q sin ∼ 6°(for θ ∼ 17°) on timescales of several years.Interestingly, where the jet follows the east-west direction (see Figure 14), the southern jet limb tends to be brighter than the northern limb, which is similar to the north-south asymmetric brightness pattern seen in the EHT-2017 image (EHT Collaboration et al. 2019a).We do not see any clear signature of prominent component ejection from the core during our campaign, although some motion in the underlying flow may exist near the core (e.g., Mertens et al. 2016;Walker et al. 2018;Park et al. 2019).
In Figure 15, we show two-frequency spectral index distribution maps within the mas-scale region of the M87 jet.
For the analysis, the (u, v)-coverage was matched for the corresponding pair of frequencies and images.The 22-43 GHz map was obtained by stacking EAVN images of the three pairs of quasi-simultaneous epochs around the dates of the EHT campaign (2017 March 18-April 18) at 22 and 43 GHz.The total intensity images at different frequencies were aligned by using the two-dimensional (2D) cross-correlation analysis (Walker et al. 2000) considering the optically thin regions of the jet.The 22-43 GHz spectral index map shows a flatspectrum radio core, while the extended jet regions become progressively optically thin, which is typical for relativistic jets in AGN (e.g., Hovatta et al. 2014).A comparison of the spectral index maps produced from stacked 22-43 GHz EAVN and 24-43 GHz VLBA observations, which are separated in time by roughly a month, shows no principal difference.
In the bottom panel of Figure 15 we also show a spectral index map zoomed into the central ∼1 mas of the radio core, which was produced using the 43 GHz VLBA and 86 GHz GMVA images, aligned at the position of the peak flux densities on corresponding images.The core shows a flat spectrum, which together with the core shift indicates that the inner jets are stratified and synchrotron self-absorbed at least up to 86 GHz, as predicted by (e.g., Blandford & Königl 1979).Spectra of the jet details downstream of the core may give unreliable results owing to the sparseness of the GMVA data and non-simultaneity of GMVA and VLBA observations (more than a month), thus we omit their discussion here.

Multi-wavelength SED of the M87 Core and High-energy Jet Emission
In Figures 2,5,7,8,10,and 11 we show the lightcurves from the radio to the γ-ray bands, taken contemporaneously with the 2017 April 5-11 EHT observations.Decades of longterm flux-monitoring programs and other previous work on M87 allow us to place these multi-wavelength fluxes in a historical context.From this analysis we conclude that during the 2017 EHT campaign, the emission from the innermost regions of M87, to the extent to which they can be resolved with our multi-wavelength observations, is consistent with being in a historically low/quiescent state.Specifically, HST-1 remained in a typical low state at radio (15 mJy at 1.7 GHz).In the optical and X-rays, HST-1 seems to be in the lowest state observed so far, following the continued flux decaying trend over the last decade (e.g., Abramowski et al. 2012;Sun et al. 2018).While the GeV/TeV facilities cannot resolve the source of the emission, they also find M87 to be in a typical low state in 2017.The fact that HST-1 is so sub-dominant to the core in the higher energy bands lends confidence to the association of the total fluxes from NuSTAR and Swift with the core emission, while the localization of the VHE emission is less constrained.
As described in Section 2, many different facilities carried out observations of M87 during or around the EHT observing campaign.Although these observations were not all strictly simultaneous, within the context of the typical variability timescale of the core of several days to weeks, they effectively provide a snapshot of the broadband SED of the M87 core in a quiescent, low-flux state.In addition to near simultaneity of the multi-wavelength observations, our data processing included careful analysis of non-core emission components in order to isolate the core emission as much as possible.Therefore, this legacy data set comprises the definitive target constraints for any model seeking to explain the source behavior at the time of the EHT 2017 image acquisition.
The list of fluxes over a broad range between frequency of ;1 GHz and photon energy ;1 TeV is given in Table A8 in Appendix A. We provide both f ν and νf ν values, although the differing data analysis procedures typically yield only one of these values.The other one is calculated by multiplication or division with the representative frequency.In the high-energy regime (>0.1 keV), we adopt the geometrical mean of the band's limits as appropriate even for relatively wide bands sampling power-law spectra.As described in Section 2, all flux points are based on statistically significant detections (>3σ for all instruments except >2σ for IACTs) and provided with uncertainties equivalent to the 1σ confidence level.Upper limits on flux are given at the 2σ confidence level.
Table A8 in Appendix A also contains band labels typical for the given energy or frequency range, as well as representative angular scales and the time ranges over which the data were acquired.While most observations were performed within a few days of the EHT observations, some of the data included in our broadband SED were acquired in single observations up to one month later, e.g., EVN, HSA, VLBA.For the majority of instruments we have multiple observations; in those cases we average the fluxes over observations taken within a week of the EHT campaign.In the case of Fermi-LAT, the fluxes provided for the SED are integrated over a period of three months, chosen as a compromise between simultaneity and photonlimited data quality.
It is very important to note that the emission measured by the instruments contributing to this broadband SED spans a factor of 10 8 in angular scale: from the effective spatial resolution of the EHT, ;20 μas, to the resolving power of Fermi-LAT in low-energy γ-rays, ;2°.In the radio bands, we list unresolved peak fluxes to represent the core.One exception is the EHT, for which we list the estimated total flux within a 60 × 60 μas 2 region, as described in EHT Collaboration et al. (2019f, their Appendix B).For elliptical beams, we list the average of the axes as a representative angular scale.At higher energies, the scale typically corresponds to the FWHM of the point-spread function within a given band or the diameter of the region used for data extraction.For Chandra and NuSTAR data analyzed jointly, the scale is tied to Chandra's greater resolving power, even though NuSTAR can only separate out the core emission spectroscopically (see Section 2.3.3).Because Swift-XRT cannot resolve the core from other components detected by Chandra, we treat its unresolved flux as an upper limit.
The data listed in Table A8 in Appendix A are plotted in Figure 16 with labels highlighting different observatories and instruments.For clarity, we omitted every other point in the X-ray spectrum and one high upper limit from the figure.In the figure we also highlight the upper limits on emission region size for the VLBI measurements.The machine-readable form of the table is provided to the community via the EHT Collaboration Data Webpage261 along with associated documentation and all data files needed for spectral modeling of the Chandra, NuSTAR, and Swift-XRT data.The list of all supplementary material published along with this paper is provided in Appendix B.

Single-zone Heuristic Modeling and Interpretation
As is the case for many nearby low-luminosity AGNs where broadband SEDs are available, M87's radio core and SED indicate a multi-scale, stratified, and self-absorbed jet up to at least 86 GHz (see Section 3.1 and Figure 16).As this requires detailed dynamical modeling at a level that goes well beyond the goal of this Letter, we have opted to use the simpler singlezone approach common to many AGN blazar papers, in order to highlight some baseline properties of the peak emission regions.Furthermore, this approach will enable relative comparisons to other M87 epochs, as well as to the other AGNs observed by the EHT, which we plan to present in future papers.
A structured, multi-zone jet can to first order be considered as a summation of many single-zones.We consider two different single-zones, one representing the launchpoint of the jets and the EHT-detected emission region, and one roughly 100 times larger and thus probing a region further down in the inner jets.For each of these size scales we also explore a different approach.The first seeks to maximize the contribution to the entire broadband SED from the most compact regions imaged by EHT, while the second focuses on statistical fitting of the X-rays, allowing an exploration of the degeneracies inherent to such modeling.Both classes of model share a spherical geometry, assume isotropy, and calculate synchrotron and synchrotron self-Compton (SSC) emission, but otherwise have somewhat different parameters and approaches.However, in all cases the predicted radio emission must be at or below the flux measured for the associated size scale, to avoid violating the radio constraints.We have labeled these constraints for clarity in Figure 16 as well as provide a quick-reference in Table 1.It is important to emphasize that any single-zone model that dominates over a large range of radio/mm frequency will be in conflict with the core size constraints.
For all models, we take the mass to be M BH = 6.5 × 10 9 M e , the distance to be 16.8 Mpc, and a source inclination (viewing angle of the emitting region with respect to the line of sight) of 17°, as adopted by EHT Collaboration et al. (2019a).A8 in Appendix A) with fluxes measured by various instruments highlighted with different colors and markers.Note that only every other point in the X-ray spectrum is plotted here and that one Fermi-LAT upper limit is missing, off to the upper right of the figure.For the mm-radio VLBI, the upper limits on emission size for several representative frequencies are labeled to clarify the constraints used in Section 3.3.An illustration of the resolved flux differences depending on spatial resolution is shown by the comparison of the differing EHT and ALMA-only 230 GHz fluxes and size limits.First, we consider a single-zone model that aims to provide a straightforward description of the flux and compact emission region size measured by EHT and other VLBI facilities.Even in the framework of a single-zone model, the predicted spectra depend strongly on how nonthermal electron (and positron in principle) distribution functions (eDF) are prescribed.In order to take into account such uncertainties, we therefore explore two different model scenarios with different eDF treatments.The first model, hereafter referred to as model 1a, includes the effects of radiative cooling on the initial single power-law eDF (see, e.g., Kino et al. 2002, and references therein).We also apply a broken power-law model to allow additional degrees of freedom in the eDF shape, but without directly calculating radiative cooling (RAIKOU code, see Kawashima et al. 2019, and also T. Kawashima et al. 2021, in preparation), hereafter referred to as model 1b.In this way we explore a large range of eDF parameter space for the single-zone approach.
Models 1a/1b share the following parameters determining the macroscopic characteristics of the spherical emission region: the radius (R), the Doppler factor of the bulk motion (δ), and the global magnetic field strength ( ¢ B ).The power law of nonthermal electrons injected in the emission region is characterized by the following quantities: the total (energyintegrated) number density of non-thermal electrons (n e ), the spectral indices of the injected power law N(γ) ∼ γ − p , and the minimum/maximum Lorentz factors (g min and g max ).Model 1a includes the feedback of radiative cooling on the eDF, thus it self-consistently calculates the break Lorentz factor (γ br ), which is not included in Table 2 as a free parameter.On the other hand, model 1b treats the γ br as a free parameter and thus it needs the second power-law index p 2 above γ br .Thus, model 1b includes additional degrees of freedom in the eDF.
Because models 1a/1b focus on describing an emission region with size and flux as determined from the EHT observations in 2017 April (EHT Collaboration et al. 2019a), we fix the emission region radius to around what is expected theoretically for the source of the observed blurred ring image in M87 (see Table 2).Since the bulk motion of the emission region has likely not yet reached a relativistic speed, we fix the Doppler factor to δ = 1 and this is consistent with the brightness temperature measured at 230 GHz to be between 10 9 and 10 10 K (EHT Collaboration et al. 2019a).We also set the minimum Lorentz factor of the injected nonthermal electrons to g = 1 min -2, which has little impact on the resultant spectra but will affect the derived total power for the singlezone.To avoid over-producing the UV flux density, we fix the injected power-law spectral index of nonthermal electrons to be steeper than p 1 = 2.
Once R and δ are fixed to the above explained values (see Table 2), the only properties that can change the normalization of the synchrotron flux at 230 GHz are the magnetic field strength ¢ B and the eDF (particularly the particle distribution normalization) and the synchrotron-self-absorption (SSA) turnover frequency (ν SSA ; the SSA frequency for the most compact region that we are modeling).As discussed above, the radio core must be optically thick up to at least 86 GHz, and we are "tuning" this component to account for the EHT flux at 230 GHz.In order to not underpredict this flux we find that ν SSA must in fact be rather close to 230 GHz.
Then, together with the standard theory of the synchrotron self-absorption process, one can obtain an order of magnitude estimation of the magnetic field strength as follows: g min γ br (10 4 ) g max (10 Notes.All values above in italics are not fitted parameters, they are either fixed before modeling, or derived for comparison between models (see the text).For model 2, the errors for ¢ n e and ¢ B are calculated from the 1σ range for all the model parameters, using Equations ( 2) and (3). a The total power for model 1a/1b is calculated in the same way as for model 2 assuming one proton per electron, but as protons are unconstrained this is solely for comparison purposes.b For M87 r g = 9.8 × 10 14 cm = 65.51Astronomical Units = 9.08 lt-hr.c The number density of total nonthermal electrons or electron/positron pairs.d Parameter pegged to the upper limit of the allowed interval.does not affect the estimation of the magnetic field strength shown here.It is thus clear that too small a value for ¢ B (∼mG) cannot consistently reproduce the observed EHT synchrotron flux and optical depth constraints.
With this approximate scale for ¢ B , we adjust the eDF to optimize the model, in order to have the SSC component also contribute to the X-ray flux.For model 1a, the resultant eDF steepens compared to its injected value of p 1 = 2.2, due to the rapid cooling of the injected electrons caused by the large ¢ B .We then tune the model parameters in order to maximize the γray flux, primarily by increasing n e to the maximum value that is still consistent with the VLBI observations.For model 1b, taking advantage of the flexibility of the double power law, we have tried to choose as large a p 1 value as possible within the range consistent with the UVOT and HST flux densities.
Our resulting best fits are shown in Figure 17.Both models 1a/1b show that the radio core emission peaking in the submm, as measured by VLBI, can be well explained by the optically thick part of the nonthermal synchrotron radiation.The model fitting gives the SSA turnover frequency ∼230 GHz which is consistent with the detection of the core shift constraints above, and in Hada et al. (2011).
The key difference between models 1a/1b can be seen by the spectra in the infrared (IR) band.The spectral break seen in model 1b is the result of tuning γ br in order to increase the γray flux.The optical/UV data impose the strongest limits on the synchrotron model spectra.Both models avoid overproducing the Swift-UVOT flux in the optically thin part of the synchrotron spectrum, although model 1b is higher than the HST points, which may argue for even less γ-ray production than in this maximal case.
For both models, we find that synchrotron emission dominates over the SSC contribution in the X-rays.A key result is that neither EHT-oriented model can produce sufficient GeV and TeV γ-ray emission.We attempted to maximize the γ-ray flux, so these results indicate the upper limit possible without violating the low-energy constraints.Therefore, our results support the idea that the γ-rays detected in 2017 likely originate from a region that is more extended than the EHTobserved region (R 10 r g ).

Model 2: High-energy Oriented Model
As it is clear that the γ-ray data cannot be explained simultaneously with the high-frequency/mm-radio emission via a single-zone model, we focus here on whether the non-VLBI (O/UV, X-ray and γ-ray) data could be consistent with originating in a single emission region.Due to the difficulties of disentangling the various jet and ICM components in the X-ray observations, because of their spatial resolution particularly in the NuSTAR band, it is important to model the data in detector space.We do this by importing a single-zone model into ISIS.We model all components except the core X-ray emission with the same models and parameters as described in Section 2.3.3,except swap in the single-zone model for the core power law.
The model is similar to models 1a/1b overall, but with somewhat different choices in fixed and fitted parameters.We assume that the emitting region is a sphere of radius R moving with a bulk Lorentz factor Γ j and corresponding speed βc, carrying a power in the comoving frame, divided between relativistic electrons, non-relativistic protons and a global magnetic field.We assume one cold proton per electron (no positrons), and that the electron eDF is described by a power law with slope p 2 between a minimum and maximum Lorentz factors g min and g max .We assign the index to p 2 to allow better comparison to model 1b, as our steeper distribution is likely in the radiatively cooled regime, although we do not calculate γ br explicitly.
The ratio between radiating particle and magnetic energy density is another fitted parameter, ¢ ¢ U U e B , where 2 , where 〈γ〉 is the average Lorenz factor of the electrons, and . We calculate the particle number density ¢ n e and magnetic field strength ¢ B from ¢ ¢ U U e B and L j : We take Γ j = 3, which for our assumed viewing angle of 17°r esults in a Doppler factor δ ≈ 3.4, consistent with being inwards of the region constrained to Γ  6 (Snios et al. 2019).
We fit the model to the data via the spectral modeling tool ISIS (Houck & Denicola 2000), using its implementation of the emcee algorithm (Foreman-Mackey et al. 2013) in order to thoroughly explore the parameter space, as well as evaluate model degeneracies and parameter uncertainties, using Cash statistics. 262We run the algorithm with 20 walkers, and for 20,000 loops, conservatively taking the initial 14,000 iterations as the chain "burn-in" period in order to ensure convergence.We define the best-fitting parameters as the median of the posterior distribution in the final 6000 loops, and the 1σ error bars as the intervals that contain 68% of the walkers.
The best-fitting parameters are shown in Table 2, and fits to the broadband SED and X-ray spectra of the core components are shown in Figure 18.Without a strong prior constraint on the size as in models 1a/1b, we found that most model parameters are highly degenerate and cannot be constrained satisfactorily by our fit.This behavior is shown in Figure 19, which shows the 2D posterior distributions where this degeneracy is most noticeable.This degeneracy is slightly reduced for models 1a/1b because of the additional constraints from optical depth effects and because g min is free in model 2 and frozen at 1 in models 1a/1b.

Conclusions from Single-zone Fitting
As expected from the VLBI constraints, a stratified outflow model is necessary in order to capture the primary features of the entire broadband M87 SED.A compact, isotropic, singlezone model that is mildly magnetically dominated can be "maximally tuned" to capture both the EHT-mm flux as well as the X-ray power, but falls significantly below the γ-rays.This model does not, however, account for the fact that the synchrotron cooling time for the X-ray producing electrons is ∼30 s in a ∼5 G magnetic field, which is in tension with the modest X-ray variability.A larger single-zone model that is significantly particle-dominated can describe the X-rays well statistically, but cannot provide a compelling fit to the radiomm VLBI core nor the VHE emission particularly in the GeV range, and even then requires rather extreme parameter values (discussed further below).These two approaches effectively approximate regions within either the base and inner jet of an expanding outflow (longitudinally structured), or an inner spine and outer sheath in a radially structured outflow.Therefore, we conclude that a structured jet is necessary to explain M87ʼs observational properties, and that an additional emission component other than the EHT-core seems necessary to explain the γ-ray emission observed in 2017.
The physical conditions in the jet base of M87 are critically important for understanding the jet power and jet-launching mechanism (see the discussion below regarding the jet Lorentz factor).Recently there have been divergent views based on more sophisticated modeling approaches applied to nonsimultaneous data sets, with Kino et al. (2015)   4 .Interestingly, our results fall between these two extremes, consistent with mildly magnetically dominated conditions in the jet launching region probed by EHT, evolving to moderate particle domination either further out along the jet axis, or radially (sheath).
This picture is somewhat consistent with the MAGIC Collaboration et al. (2020) results as the jet base in their model is already over 400 r g , and thus similar to our model 2 scenario.This larger size scale, however, is in conflict with the 262 For the bands outside of the X-ray, we use artificial diagonal response matrices with a combination of area/exposure times so that we fit a "counts" in each frequency bin with a Cash statistic in its χ 2 limit, with the resulting error corresponding to the Gaussian error of that frequency bin.Thus the 1σ Gaussian error corresponds to the square root of these "counts," which are in the χ 2 limit of the Cash statistic.Only the X-ray band is treated in the "low counts" regime where the Cash statistics deviate significantly from χ 2 statistics.VLBI size constraints above ∼5 GHz, and the core shift/ optical depth requirements described above.However, we note that their predicted mm flux is a factor of a few lower than that measured by EHT, and thus emission from this component could in principle dominate the VHE emission while contributing only slightly to the radio/mm core.
The bigger issue pinpointed by MAGIC Collaboration et al. (2020) and common to many VHE-focused models is the need for extreme particle domination (  ¢ ¢ U U 10 e B 3 ) on scales that correspond to M87's inner jets (z 10 4 r g ).The extreme values obtained by this and other models may originate in the nonsimultaneity of earlier data sets, and/or the non-inclusion of the VLBI constraints on geometry.VLBI images at cm wavelengths (Asada & Nakamura 2012;Hada et al. 2013;Nakamura & Asada 2013) show a roughly parabolic radial jet profile on these scales, usually interpreted as a magnetically accelerated bulk flow being collimated by the external medium (e.g., Nakamura et al. 2018).The fact that GRMHD simulations of large-scale jets can self-consistently obtain such a profile further supports this argument (e.g., Chatterjee et al. 2019, Figure 13).Upcoming EHTC results on the polarization properties of M87 will provide even more stringent constraints (EHT Collaboration 2021a, 2021b).
Furthermore, recent results constraining the innermost knot HST-1 to be moving superluminally outwards at an apparent speed requiring a Lorentz factor of ∼6 (Snios et al. 2019, and references therein) indicate a moderately relativistic bulk flow velocity.Even disregarding the above arguments, such a speed would be challenging to obtain for extreme particle domination in the regions probed by EHT.In the ideal MHD scenario, the bulk flow velocity at large scales reflects the magnetization (the ratio of magnetic to gas enthalpy) at the jet base (see, e.g., Komissarov et al. 2007;Tchekhovskoy et al. 2009), which is now also reproduced in large scale numerical GRMHD simulations (Chatterjee et al. 2019).If the base of M87's jets are extremely magnetically dominated, then instabilities leading to mass entrainment may be required very close in, to reduce the final bulk flow velocity.
Interestingly, despite being in a historically low state for the high-energy emission, our conclusions regarding M87's mm/ radio core properties and power are very similar to those presented in Reynolds et al. (1996).Based on earlier 5 GHz VLBI measurements, the authors concluded that the jets are more likely to be pair-dominated because of power requirements.Similarly, when we estimate the power assuming (somewhat arbitrarily) an electron-proton plasma with charge neutrality, we find models 1a/b to be in tension with estimates of the total jet power of L j ∼ 10 −4 L Edd , based on the pressure arguments for the knots.As Reynolds et al. also discuss, this power issue can be mitigated by raising g min , as we also found for model 2. This overall consistency between power and inferred properties for the inner jets over decades may provide further support for the scenario where the changes seen at highenergy originate further out in the jets.
Constraining the jet composition near the launch point will have key consequences for not only the power requirements, but potentially the high-energy emission.If hadrons are indeed present and accelerated, a population of >100 TeV protons in the jet base could potentially produce γ-rays via secondary neutral π 0 decay or via inverse Comptonization from secondary leptons injected at similar energies.As a consequence, the models as presented would need to be renormalized in order to avoid overproducing the lower frequency emission, in turn requiring even more particle domination over the magnetic fields.As described above, this energy partition would further exacerbate the problems explaining the jet acceleration and the core shift, as well as increasing the jet power, particularly for accelerated proton power laws with index p < 2. Another possibility is direct proton synchrotron from >EeV protons, as well as secondary muon synchrotron, given the potential for strong magnetic fields in the EHT-emission region (see, e.g., Reimer et al. 2004).We cannot rule this scenario out without more detailed modeling, however it may be difficult to reconcile with the accompanying low-energy emission if leptons are assumed accelerated at the same rate.Regardless of mechanism, all models seeking to produce the observed ∼10 42 erg s −1 in γ-rays within the innermost 10r g of M87 would have to find away around the extreme attenuation expected due to pair production off the core optical/UV photons.
In conclusion, while we cannot constrain the composition of the jets, power arguments favor lighter jets at least in the jet base.Our results also favor moderate magnetic domination at the launch region of M87's jets, at least where the EHT and possibly some X-ray emission are produced, while a larger and more particle-dominated region than the "EHT-zone" is likely responsible for the VHE emission.Either this region is not in the accelerating part of the jet flow, or it could be interacting with the surrounding ICM and may require an additional source of seed photons to the jet emission alone.The VHE emission could also be coming from the extended jet and knots outward of HST-1, not included in our radio through X-ray images but still unresolved in the γ-rays.Correlated multi-wavelength variability in future campaigns will be the ultimate test of these scenarios.We will return to these questions with a more time-dependent analysis in the EHT 2018 multi-wavelength campaign paper.

Summary
We present the most extensive, quasi-simultaneous, broadband spectrum of M87 taken yet, together with the highest ever resolution mm-VLBI images using the Event Horizon Telescope from its 2017 April campaign.The primary result of this Letter is the presentation of this legacy data set, for which all data and some analysis scripts are available to the community via a Cyverse repository; see Appendix B. From the observational side, our main conclusions are that the M87 core was in a relatively low state compared to historical observations, but clearly still dominating over the nearest knot HST-1, which was seemingly at its lowest historical brightness state.This provides the ideal observing conditions for a multiwavelength campaign combining data over a large range of spatial resolution, because M87's core was dominating the total flux in radio through X-ray bands.
While we defer detailed modeling to further work by ourselves and the broader community, we can already make some baseline conclusions based on simple single-zone models.The first conclusion is that M87's complex, broadband spectral energy distribution cannot be modeled by a single zone.The stratified nature of the radio through millimeter bands is clearly indicated by the high dynamic-range images that place significant size constraints on the emitting geometry per frequency band.However by selecting two representative regions for the EHT-observed region and the inner jets, respectively, we can already come to a few solid conclusions.First, it is not yet clear where the VHE γ-rays originate, but we can robustly rule out that they coincide with the EHT region for leptonic processes.The energetics and small size of this region, together with the VLBI constraints on optical depth and thus strong magnetic fields, cannot be reconciled with significant γray production via SSC emission.A ∼100 times larger region can provide VHE γ-ray fluxes but requires an uncomfortably high particle domination compared to magnetic fields, which may be inconsistent with the observed jet velocity and parabolic geometrical profile.Direct proton and muon synchrotron emission from the EHT-emission region contributing to the GeV/TeV range cannot be ruled out at this time, but can be tested in the future using the upcoming improved constraints on the magnetic fields expected from EHT polarization results.We emphasize in any case that a structured jet model including timedependence will be important for any detailed interpretation.Additional information will be presented in the upcoming paper on the 2018 campaign, where there will be new constraints on the multi-wavelength variability properties.
This "quiescent active" core spectrum for M87 is also important for comparison with the upcoming analysis on the Galactic center supermassive black hole Sgr A * .Sgr A * has also been in an extended, much weaker quiescent period since (at least) its first identification as a high-energy source (Baganoff et al. 2001(Baganoff et al. , 2003)), punctuated only by moderate flaring (e.g., Neilsen et al. 2015;Witzel et al. 2018;Do et al. 2019;Haggard et al. 2019).M87's low state thus provides the ideal comparison spectrum to understand how sources evolve from quiescence into low-luminosity AGN and launch more extended jet outflows.
With results from the 2017 EHT campaign for 3C 279 already published (Kim et al. 2020) 2021, in preparation), and many more in the pipeline, we expect the context and impact of the multi-wavelength results to significantly broaden in the coming years.The precision EHT mm-VLBI images together with broadband SED information will be essential ingredients for developing a fuller understanding of the AGN phenomenon.Beyond simple accretion rate changes, this new body of work will allow us to directly study the interplay between strong gravity and magnetic fields and to trace those effects to the larger jet structures and their feedback onto the external environment.With our legacy sample, enhanced by multi-year campaigns with the steadily enhanced EHT array, as well as eventually with the future next-generation EHT, we aim to provide a rich, public data set to the community for decades of projects to come.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 EHTspecific 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.
The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes.Scientific results from data presented in this publica-Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.

Appendix A Tabulated Data
This section contains all tabulated data mentioned in Sections 2 and 3.2, as follows.Note.F uvw1 , F uvm2 , F uvw2 -flux densities within a radius of 5″ in bands uvw1, uvm2, and uvw2, respectively, corrected for the extinction and host galaxy contamination.

Note.
a Only statistical 1σ errors on the fluxes are given.95% confidence level upper limits are given for flux measurements with less than 2σ significance.

Figure 3 .
Figure 3. Visibility amplitude vs. uv-distance plot of GMVA data on 2017 March 30.In this display the visibilities are binned in 30 s intervals for clarity.

Figure 4 .
Figure 4. Mean surface brightness in uvw1 band as a function of the semimajor axes of concentric ellipses with ellipticity and PA equal to those of the host galaxy; the red, blue, and green curves show the core + HST-1, jet, and host galaxy profiles, respectively; the solid black curve represents the observed flux, and the black dotted curve gives the total flux of the model; all of the curves show the average profiles over all images in this band.

Figure 5 .
Figure 5. Optical and UV light curves of the core region of M87 in all six UVOT filters; the red lines show the average flux density in each filter during the EHT campaign and their standard deviations (red dotted lines) along with the time and duration of the campaign.

Figure 8 .
Figure 8. Chandra X-ray lightcurve of the core of M87 (black) and HST-1 (blue) in 2017 April, showing a small amount of variability between observations 20034 (left panel) and 20035 (right panel).

Figure 7 .
Figure 7. Historical optical and UV light curves of M87's core and the HST-1 knot (see also Madrid 2009; Perlman et al. 2011).Horizontal dashed lines indicate the lowest flux density level of the core in the F606W (red) and F275W (black) filters, the vertical green dotted line marks the time of the EHT campaign.

Figure 9 .
Figure 9. Chandra and NuSTAR spectra of M87.Top panel: count rate spectra from both observatories from 2017 April.Spectra have been rebinned for plotting purposes.In both panels, the Chandra spectrum of the underlying ICM is shown in green, while the NuSTAR FPMA and FPMB are shown in black and blue, respectively; the red curve is the total model spectrum for each data set.In the top panel, the spectrum of the core, HST-1, and the outer jet are shown in orange, magenta, and cyan, respectively.Bottom panel: unfolded spectra from Chandra and NuSTAR.Because the Chandra core/HST-1/jet spectra have the pileup kernel applied (Section 2.3.1),they are excluded from the flux plot.

Figure 10 .
Figure 10.Variability in the Swift-XRT unabsorbed flux (2-10 keV) from 2017 March 27 to April 20.The Total Flux refers to the overall unabsorbed flux derived from the spectral fitting to the composite model.The net flux from only the three components, i.e., core, HST-1 and outer jet is referred to as the Flux(core+HST-1+jet), which was calculated from their respective spectral fitted parameters.

Figure 11 .
Figure 11.Flux measurements of M87 above 350 GeV with 1σ uncertainties obtained with H.E.S.S., MAGIC, and VERITAS during the coordinated MWL campaign in 2017.Upper limits for flux points with a significance below 2σ are provided in Table A7 in Appendix A.

Figure 13 .
Figure 13.Compilation of the quasi-simultaneous M87 jet images at various scales during the 2017 campaign.The instrument, observing wavelength, and scale are shown on the top-left side of each image.Note that the color scale has been chosen to highlight the observed features for each scale, and should not be used for rms or flux density calculation purposes.Location of the Knot A (far beyond the core and HST-1) is shown in the top figures for visual aid.

Figure 14 .
Figure 14.Innermost region of the M87 jet observed with VLBA at 43 GHz.To better describe the transverse jet structure, the image is convolved with a circular beam of 0.31 mas (shown in the bottom-left corner of the figure), which is slightly superresolved along the north-south direction.Black dots indicate the ridges of the northern and southern limbs beyond 0.2 mas from the core.For reference, a polar axis is overlaid to show corresponding PA.Contour levels are scaled as (1, 1.4, 2, 2.8...) × 3 mJy beam −1 .

Figure 15 .
Figure 15.(a) Spectral-index map of the stacked EAVN image including three pairs of quasi-simultaneous epochs at 22 and 43 GHz.A common beam size of 1.2 mas × 1.2 mas, corresponding to the 22 GHz beam, is used for both frequencies.The images at different frequencies were aligned by doing the 2D cross-correlation analysis using optically thin regions of the jet.The dashed box indicates the region of the bottom panel.(b) Spectral-index map from VLBA 43 GHz on 2017 May 5 and GMVA 86 GHz on 2017 March 30.A common beam size of 0.2 mas × 0.4 mas at PA = 1°.2, corresponding to the 43 GHz beam, is used for both frequencies.We note that the small patch of high opacity is likely due to poor quality and large uncertainty in the flux of the 86 GHz data (see also Section 2.1.7 for details).Both of these two images are rotated clockwise by 18°.Contours in (a) indicate the EAVN 22 GHz total intensity map, while those in (b) are the VLBA 43 GHz total intensity map.Contour levels for both (a) and (b) are scaled as (1, 1.4, 2, 2.8...) × 1.9 mJy beam −1 .A restoring beam is indicated in the bottom-right corner of each panel.The clipping levels of both images are 10 mJy beam −1 , which corresponds to over 10 times the image rms.

Figure 16 .
Figure16.Observed broadband SED of M87 quasi-simultaneous with the EHT campaign in 2017 April (see TableA8in Appendix A) with fluxes measured by various instruments highlighted with different colors and markers.Note that only every other point in the X-ray spectrum is plotted here and that one Fermi-LAT upper limit is missing, off to the upper right of the figure.For the mm-radio VLBI, the upper limits on emission size for several representative frequencies are labeled to clarify the constraints used in Section 3.3.An illustration of the resolved flux differences depending on spatial resolution is shown by the comparison of the differing EHT and ALMA-only 230 GHz fluxes and size limits.
factor (conventionally denoted as b(p)) related to synchrotron absorption has been described in the literature(Marscher 1983;Hirotani 2005;Kino et al. 2014) in detail and the slight variation in b(p) in those literature values

Figure 17 .
Figure17.SED fit focusing on the EHT data.Blue and green lines display the resulting SEDs for models 1a and 1b, respectively.At γ-ray energy bands, one can see the SSC emission components although they underestimate the observed γ-ray flux density.

Figure 18 .
Figure18.Top panel: same data as Figure17, showing the model 2 fit focusing on the higher energy data.Bottom panel: statistical model 2 fit to the X-ray data.The orange points are the Chandra data, the blue and black points the NuSTAR data, and the red line shows the total model (which for NuSTAR includes the contribution of the ICM, HST-1 and kpc-scale jet).Note that when performing this fit, we modeled the X-ray spectra in detector space, simultaneously with the multiwavelength data shown in the top panel.
Collaboration et al. (2020), focusing on the MAGIC/VHE bands, prefer the opposite extreme of particle dominance with

Figure 19 .
Figure 19.2D posterior distributions for the parameters that show degeneracies in the high-energy single-zone SSC model.Left to right, top row: injected power L j (always given in units of L Edd ) vs. emitting region radius r, injected power L j vs. ¢ ¢ U U e B , injected power L j vs. minimum electron Lorentz factor g ; min left to right, bottom row: maximum electron Lorentz factor g max vs. ¢ ¢ U U e B , electron power-law index s vs. ¢ ¢ U U e B .The red lines represent the location of 68, 95 and 98 percent of the walkers.
represents the unresolved flux of the radio core (the most compact and brightest feature at the jet base) within the beam size.b Total flux density is measured by integrating the radio fluxes over the entire jet structure within the field of view of each observation.

,
with b n = 2n − 1/3, R e -half-light radius, and I e -intensity at R e .Columns are: (1)-UVOT band; (2) -total extinction; (3)-value of Sérsic parameter n; (4)-half-light radius of host galaxy; (5)-its ellipticity; (6)-PA of the major axis of host galaxy; (7)-host galaxy flux density in a circle aperture of a radius of 5″; (8)-flux density of the core in a circle aperture of a radius of 5″, averaged over the period of the EHT campaign and its standard deviation.Flux densities are corrected for the extinction.
flux represents the unabsorbed flux derived from the Swift-XRT spectral fits.See Table A8 for joint Chandra and NuSTAR constraints on the flux.b Net flux corresponds to the unabsorbed flux from core, HST-1, outer jet only in the 2-10 keV band.The flux units are 10 −12 erg cm −2 s −1 Matteo et al. (2003)over and above the column toward the rest of the jet.The origin of this excess absorption is not known.It was not detected by DiMatteo et al. (2003), but was mentioned by 22 cm −2 , though this is larger than the Galactic value from H I studies (0.025 × 10 22 cm −2 ; Dickey & Lockman 1990, 0.025 × 10 22 cm −2 ; HI4PI Collaboration et al. 2016).The core, however, required an additional = - (Aharonian et al. 2006aari et (Albert et al. 2008 al. 2012;.E.S.S.;Aharonian et al. 2006a), the Major Atmospheric Gamma Imaging Cherenkov (MAGIC;Albert et al. 2008;Aleksić et al. 2012; MAGIC Collaboration et al. 2020)telescopes, and the Very Energetic Radiation Imaging Telescope Array System (VERITAS;Acciari et al. 2008Acciari et al. , 2010;;Aliu et al. 2012; Beilicke & VERITAS Collaboration 2012).Over the course of its monitoring, M87 exhibited high VHE emission states in 2005(Aharonian et al. 2006a), 2008(Albert et al. 2008), and 2010

Table 1
Spatial Extent (Diameter of a Circle) within which Compact Radio Fluxes were Measured with VLBI Observations at Each Frequency

Table A1 ,
2. UV/optical host galaxy and core region parameters in Table A2, 3. Optical flux densities from Swift-UVOT observations in Table A3,

Table A2
UV/Optical Host Galaxy and Core Region Parameters Note.F v , F b , F u -flux densities within a radius of 5″ in bands v, b, and u, respectively, corrected for the extinction and host galaxy contamination.