© The Authors 2024

1. Introduction

The M87 elliptical galaxy in the Virgo Cluster is one of the nearest active galactic nuclei (AGN), located at a distance of 16.8 ± 0.8 Mpc (Blakeslee et al. 2009; Bird et al. 2010; Cantiello et al. 2018, and Event Horizon Telescope Collaboration 2019d) It also harbours one of the largest known supermassive black holes (SMBHs; ∼6.5 × 10⁹ M_⊙ determined by Gebhardt et al. (2011), but see also, e.g., Walsh et al. (2013), Liepold et al. (2023), Osorno et al. (2023), Simon et al. (2024) finding a slightly smaller mass). It is one of only two SMBHs where a ring of light emitted from near the event horizon has been directly imaged at 230 GHz, via the global Very Long Baseline Interferometry (VLBI) project, the Event Horizon Telescope (EHT). The ring’s size and shape are directly predicted by general relativity (GR) and depend on the black hole (BH) mass and spin and the orientation of the spin axis with respect to our line of sight (see, e.g., Bardeen 1973; Johannsen & Psaltis 2010).

In April 2019, the EHT Collaboration (EHTC) presented the first M87* images (we note that we use M87* to refer to the SMBH rather than its host galaxy) from the 2017 campaign along with modelling and interpretation of the broadband emission (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f). A primary result of this sequence of works by the EHTC was an independent determination of M87*’s mass: (6.5 ± 0.7)×10⁹ M_⊙, consistent with the previous measurements of Gebhardt et al. (2011)

For its modelling, EHTC uses state-of-the-art techniques (see, e.g., Event Horizon Telescope Collaboration 2019e; Porth et al. 2019; Gold et al. 2020, M87 2018 II) that combine ideal GR magnetohydrodynamics (GRMHD) with GR Ray Tracing (GRRT). GRRT is typically conducted as a post-processing step and is highly dependent on the assumed electron distributions that are not determined from first-principles methods. For example, the electrons emitting synchrotron radiation are accelerated by and change the local electromagnetic fields, and they may also interact with the radiation field. The addition of constraints that can guide microphysical kinetic models of particle acceleration is a clear way to reduce the current degeneracy in interpretation.

At the same time, the presence of strong, ordered magnetic fields near the SMBH is required to explain the detected linearly polarised radio synchrotron radiation as well as the angular momentum transport in the inflow and the launching of powerful jets of plasma (see, e.g., Blandford et al. 2019; Event Horizon Telescope Collaboration 2021a). M87* currently expels vast bipolar jets. However, in most observations, only the one that is relativistically beamed towards Earth at approximately 20° from the line of sight is visible (Mertens et al. 2016), as the jets are accelerated to apparent velocities of up to 6c (where c is the speed of light) on sub-arcsecond scales (Snios et al. 2019). Previous studies have found that the width and the outflow velocity of the approaching jet vary as a function of the distance from the SMBH (Asada et al. 2014), and other physical parameters possibly vary in a similar manner (Asada & Nakamura 2012). These jets radiate across the entire electromagnetic spectrum, with contributions to different parts of the spectral energy distribution (SED) arising at different places along their elongated structure. A complete model ultimately needs to include both the inflow and outflow components and be able to not only match the EHT image properties but also reproduce the inner jet dynamical properties as well as the overall multi-wavelength (MWL) spectrum

For the 2017 EHT papers on M87*, complementary MWL information, such as the estimated minimum jet power (P_jet ≥ 10⁴² erg s⁻¹; e.g. Reynolds et al. 1996; Stawarz et al. 2006; de Gasperin et al. 2012; Prieto et al. 2016), allowed roughly half of the ∼45 initial GRMHD models to be ruled out, including all models with zero spin. However, MWL constraints on the SED such as the X-ray luminosity were only applied as an upper limit and thus did not exclude many additional models. In stark contrast, the inclusion of MWL data to compare against SEDs for the second 2017 EHTC image, that was of Sagittarius A* (Sgr A*; the SMBH in the centre of the Milky Way), provided a significantly augmented constraining power (Event Horizon Telescope Collaboration 2022a,b,c,d,e,f; Wielgus et al. 2022). Specifically, the infrared and X-ray observations on Sgr A*’s emission from the Keck Observatory, the Very Large Telescope (VLT), the Chandra X-ray Observatory, and Nuclear Spectroscopic Telescope Array (NuSTAR) provided particularly tight constraints and allowed us to downselect from a large number of original synthetic images that were generated from models covering a much wider parameter space range than the original set for M87*.

We have acquired a full set of complementary, contemporaneous MWL observations for M87* to support modelling and interpretation as well as to provide a legacy dataset for the community. Together with the EHT images, these datasets provide rich input to help distinguish between the currently large range of theoretical scenarios for accretion and jet launching as well as for particle acceleration, that powers the emitted radiation. The first paper in this series was published as EHT MWL Science Working Group (2021) (hereafter M87 MWL2017), and it presents the extensive MWL observations carried out alongside the EHT 2017 campaign on M87 by the EHT-MWL Science Working Group. The working group included EHTC members as well members of more than a dozen other partner facilities on the ground and in space. M87 MWL2017 includes VLBI images and spectral index maps at several frequencies lower than the 230 GHz band probed by EHT together with an SED spanning centimetre-band radio through TeV γ-rays, a range of more than 17 decades in frequency. In 2017, this source was found to be in a historically low state at all frequencies, with the core flux dominating over the HST-1 knot at high energies, making it possible to combine the high energy core flux constraints with the more spatially precise VLBI data.

In M87 MWL2017 and the present work, we merge the VLBI and spectral constraints by ‘tagging’ the radio flux points with VLBI constraints on the emitting region size, a constraint that is often missed by single-zone SED modelling approaches. In both papers we also explore M87*’s jet properties in an order-of-magnitude approach via two heuristic single-zone models meant to help reveal general trends and allow rough comparisons between basic properties year to year. Single-zone models are a common first approach that often successfully account for most of the optical and higher energy emission in many AGN, particularly blazars. The success of these techniques implies that a significant fraction of energy in the jets is converted into particle energy in a relatively localised region. How this process occurs is one of the most urgent and outstanding questions in the field today, relating to identification of the source of very high energy (VHE) particles detected on Earth, such as cosmic rays (CRs) and neutrinos, as well as understanding the electromagnetic counterparts of explosive transient events, such as compact object mergers involving a neutron star. Localisation is thus a key step in constraining the particle acceleration mechanism as well as in investigating the hadronic content of the jet.

M87 MWL2017, however, demonstrated that M87* requires a stratified jet model, meaning that a single-zone leptonic model cannot adequately describe the sub-millimetre and high-energy SED simultaneously. In particular, we showed that the γ-rays cannot be fit by the same single-zone model that can match the EHT flux and size constraints. While this would seem to exclude scenarios where the γ-rays detected in 2017 in the quiescent (i.e. non-flaring) state are emitted from the same region as the sub-millimetre radiation in the EHT image, we discuss new advances in the modelling that may indicate otherwise.

Meanwhile, a new constraint that has yet to be incorporated into the EHT modelling is the source emission variability. In the EHTC Sgr A* papers, it was shown that the root-mean-square variability is one of the most difficult constraints to satisfy, and in fact, almost all of our GRMHD models are somewhat too variable compared to the data (see Event Horizon Telescope Collaboration 2022g). Based only on the 2017 observations, it was not possible to draw conclusions regarding the variability of M87*. For prograde orbits, the dynamical timescale for light to circle M87*’s innermost stable circular orbit ranges from just under a week to approximately a month, depending on the spin, that is currently unconstrained. Therefore, daily images taken over a single week can be considered to correspond effectively to a single snapshot. However, one year represents a significant jump in time for M87*, roughly ten to 70 dynamical timescales, providing a chance to temporally resolve its uncorrelated states.

The EHTC has recently published the first results from its April 2018 campaign on M87* (M87 2018 I). The 2018 EHT images of M87* reveal an asymmetric ring structure that is brighter in the southeast and with a diameter of ∼43 μas. This is remarkably consistent with the EHT images obtained from the 2017 observations, strongly confirming the presence of a supermassive black hole with a mass M_BH ∼ 6.5 × 10⁹ M_⊙. The 2018 images show a significant shift in the position angle (by ∼30°) of the ring brightness asymmetry with respect to that of the 2017 images, suggesting the presence of variations on yearly timescales in the event-horizon-scale structures. However, the compact flux density in the 2018 EHT data is less well constrained than in 2017 due to the more limited coverage of short-to-intermediate baselines. Bayesian image reconstruction methods yield a relatively large range of ∼0.5–1.0 Jy within a scale of ≲100 μ (see M87 2018 I for a full description of the 2018 EHT observations).

During the accompanying 2018 MWL campaign, all currently operating ground-based imaging atmospheric Cherenkov telescopes (IACTs; H.E.S.S., MAGIC, and VERITAS, see Sect. 2.4.2) as well as Fermi-LAT measured the first γ-ray flare seen in M87* in eight years. In this paper, we present the results of the 2018 MWL campaign as well as a comparison with the 2017 results. In Sect. 2, we describe the new MWL observations with band-specific results and, where relevant, light curves and comparisons to prior observations, while we report in Appendix B further details on the data processing. In Sect. 3, we present the compiled nearly simultaneous SED of M87* for the 2018 MWL campaign. In Sect. 4, we present the best fits using phenomenological single-zone models similar to those used in M87 MWL2017 and discuss indications for structural changes between 2017 and 2018. In Sect. 5, we give our conclusions. All data files and products are available for download, as described in Sect. 3 and in appendix. Data availability.

2. Observations and data reduction

The 2018 EHT observations of M87 were made in April 2018 using a total of eight stations (in six geographical locations) across the globe. While observations were performed over four nights (21, 22, 25, and 28 April), the second and last sessions suffered from various issues (bad weather conditions across the array, poor baseline coverage, etc.). Therefore, reliable EHT images of M87* could only be obtained from the data taken on 21 and 25 April

In the following subsections we present brief descriptions of the 2018 MWL campaign on M87* – more detailed descriptions including data processing procedures and band-specific analyses appear in Appendix B. To aid readability, tabulated data are collected in Appendix C. Figure 1 shows a schematic overview of the observations carried between March and May 2018 for the second EHT-MWL campaign on M87. In addition, to observations carried during March–May 2018 we considered for our study also global VLBI (carried on 1 February 2018) and HST (carried on 30 July 2018) observations

Fig. 1. Instrument coverage summary of the 2018 M87 MWL campaign covering MJD range 58198–58270. The grey band indicates the time interval of the VHE γ-ray flare episode.

2.1. Radio observations

Here we outline the radio-millimeter observations obtained during the 2018 campaign with various VLBI facilities and connected interferometers. Following M87 MWL2017, we use the term “radio core” to represent the innermost part of the radio jet imaged by VLBI. A radio core in a VLBI jet image is conventionally defined as the most compact (often unresolved or partially resolved) feature seen at the apparent base of the radio jet (e.g. Lobanov 1998; Marscher 2008). For this reason, different angular resolutions by different VLBI instruments and frequencies, together with the frequency-dependent synchrotron optical depth, can lead to differences in the identification of a radio core in each observation (see also M87 MWL2017, for a related discussion).

2.1.1. VERA 22 GHz

The nucleus of M87 was densely monitored over 2017 and 2018 at 22 GHz with the VLBI Exploration of Radio Astrometry project (VERA; Kobayashi et al. 2003), as part of a regular monitoring program of a sample of γ-ray bright AGN (Nagai et al. 2013). A total of 17 and 13 epochs were obtained in 2017 and 2018, respectively (see Table C.1). 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 (additional details appear in Appendix B.1.1). Peak flux density light curves obtained with VERA at 22/43 GHz during 2017–2018 are included in Fig. 2

Fig. 2. Radio light curves of the M87 core in 2017 and 2018 at multiple bands. The upper and lower panels are for connected interferometers and VLBI, respectively. The corresponding beam sizes are indicated in Table C.1. KVN data at 22 and 43 GHz are not shown here since KVN captures the data from the shortest baselines of EAVN.

2.1.2. EAVN/KaVA 22 and 43 GHz

Between January 2017 and June 2018, M87 was regularly monitored at 22 and 43 GHz with the East Asian VLBI Network (EAVN; Cui et al. 2021; Akiyama et al. 2022), a joint VLBI array of radio telescopes in East Asia. A total of 22 (9 at 22 GHz and 13 at 43 GHz) and 20 (9 at 22 GHz and 11 at 43 GHz) epochs were obtained in 2017 and 2018, respectively (details in Appendix B.1.2 and Table C.1). Peak flux density light curves obtained with EAVN at 22 and 43 GHz during 2017–2018 are included in Fig. 2.

2.1.3. VLBA 24 and 43 GHz

The program targeting M87 using the National Radio Astronomy Observatory (NRAO) Very Long Baseline Array (VLBA) at central frequencies of 24 and 43 GHz was initiated in 2006 (Walker et al. 2018) and lasted until 2020. Presented here are 10 epochs of observations that were obtained in 2018 and include two imaging (i.e. full-track) and eight short (snapshot) sessions on M87 (see summary in Appendix B.1.3 and Table C.1). Peak flux density light curves obtained with VLBA at 22 and 43 GHz during 2017–2018 are included in Fig. 2.

2.1.4. GMVA+ALMA 86 GHz

M87 was observed by the Global Millimetre VLBI Array (GMVA) for the first time in concert with the phased Atacama Large Millimetre/submillimetre Array (ALMA) and the Greenland Telescope (GLT) on 14–15 April 2018 (See details in Appendix B.1.4 and Table C.1).

2.1.5. KVN 22, 43, 86, and 129 GHz

M87 was regularly monitored by the Korean VLBI Network (KVN) at 22, 43, 86 and 129 GHz, as part of the Interferometric Monitoring of Gamma-ray Bright Active galactic nuclei (iMOGABA; Lee et al. 2016) program. More information can be found in Appendix B.1.5 and Table C.1. Peak flux density light curves obtained with KVN at 86/129 GHz during 2017–2018 are included in Fig. 2.

2.1.6. Global 43 GHz VLBI

M87 was observed on 1 February 2018 at 43 GHz by a global array of sensitive VLBI antennas in full-track to maximise the imaging sensitivity and resolution. While details of the observation and imaging results will be reported elsewhere (Kim et al., in prep.), we provide a peak flux density of the core and a VLBI-scale total flux density (see Table C.1). We estimate these measurements to be uncertain up to ∼20% considering the absolute flux calibration error and systematic uncertainties in defining the core and the field-of-view (see also Appendix B.1.6).

2.1.7. ALMA 93 GHz and 221 GHz

Observations with phased-ALMA (Matthews et al. 2018; Goddi et al. 2019; Crew et al. 2023) were conducted as part of the 2018 VLBI campaigns with the GMVA at 3 mm/93 GHz (ALMA Band 3) and the EHT at 1.3 mm/221 GHz (ALMA Band 6). The main observational and imaging parameters for the 2018 observations are summarised in Table C.2, while those for 2017 are summarised in M87 MWL2017, as well as included in Fig. 2. For each dataset and corresponding image, the table also reports flux-density values both in the central compact core and the extended jet (see Appendix B.1.7 for more detailed descriptions of data acquisition and reduction).

2.1.8. SMA 200–270 and 345 GHz

The long term 1.1–1.4 mm band (200–270 GHz) and 0.87 mm band (345 GHz) flux density light curve data for M87 shown in Fig. 2, and later in Figs. 3 and 4, were obtained at the Submillimeter Array (SMA) near the summit of Mauna Kea (Hawai’i). 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 millimetre wavelengths (Gurwell et al. 2007). A summary of the measurements made from these data for 2017 and 2018 is shown in Table C.2. The reported core flux for M87 is the vector average of baselines longer than 30 kλ (see Appendix B.1.8 for more details)

Fig. 3. Long-term MWL light curves of M87 from the last two decades taken by the instruments participating in the 2018 observational campaign. Here are listed the observations used in the plot and their reference (in case of already published results): Chandra (0.3–7 keV) (Sun et al. 2018), (2.0–10 keV) (Imazawa et al. 2021); MAGIC, H.E.S.S., VERITAS (2004–2011), Fermi-LAT (2008–2011), VLBA (43 GHz peak) (2006–2011) (Abramowski et al. 2012); H.E.S.S. (before April 2004, 2001–2016, after 2019) (H.E.S.S. Collaboration 2023); VERITAS (2011–2012) (Beilicke & VERITAS Collaboration 2012); VERA (2011–2012), MAGIC (2012–2015) (MAGIC Collaboration 2020); KVN (2012–2016) recalculated from Kim et al. (2018); VLBA, EAVN, VERA (2000–2018) (Cui et al. 2023); All instruments 2017 data (M87 MWL2017). We mark the VHE flare in 2018 with a grey-shaded line in the background. The inset in the third row provides a zoomed-in view of the Chandra measurement of HST-1.

Fig. 4. Multi-wavelength light curves of M87 taken during the observational campaign covering MJD range 58211–58245. The time period containing the 2008 VHE flare is shaded in grey.

2.2. Optical and UV observations

2.2.1. Optical observations by the Kanata telescope

Photometric images of M87 were obtained using the 1.5-m Kanata telescope and the Hiroshima Optical and Near-InfraRed Camera (HONIR), that can perform simultaneous two-band imaging (Akitaya et al. 2014). We used optical R_C (634.9 nm) and near-infrared J-band (1250 nm) filters for the observations. Images with a field of view of 10″×10″ and a spatial resolution of 0.294″ pixel⁻¹ through the band-pass filters were obtained. The central flux of M87 within a radius of 5″ was measured. The measured flux should include three components: radiation from the central region near the BH, radiation from the HST-1 knot (located at a projected distance of ∼70 pc from the core), and thermal radiation from the host galaxy. Thus, the measured flux is an upper limit of the flux of the central nuclear region. Additional details appear in Appendix B.2.1 and Paragraph symbol, cited in the author “M. Santander and G. Puhlhofer”..

2.2.2. UV and optical observations from Swift-UVOT

During the Swift pointings, the UVOT instrument observed M87 in its optical (v, b and u) and UV (w1, m2 and w2) photometric bands (Poole et al. 2008; Breeveld et al. 2010). An aperture of 5″ radius was used for the flux extraction for all UVOT bands. Compared to the values obtained during the EHT 2017 campaign, we observed with UVOT in 2018 an increase of the flux by a factor of ∼3 in optical and a factor of ∼2 in UV bands. However, the high contribution of the host galaxy flux in the V band and partly in the b band contributes to the large uncertainties on the flux. A detailed description of data, analysis and modelling is given in Appendix B.2.2 (for all the results of the Swift-UVOT observations during the EHT 2018 campaign see Table C.5.).

2.2.3. UV observation with HST

We have investigated the HST archive to search for HST UV/optical observations during the 2018 EHT MW campaign. Unfortunately, the HST observation nearest to the campaign was performed more than 3 months later, on 2018 July 30 (MJD 5832.90764). The observations were carried out with the STIS camera (NUV-MAMA detector with F25QTZ aperture) over two HST orbits at the pivot wavelength 2359.7 Å (proposal ID 15358). Since none of the other Optical/UV/NUV instruments used in this campaign allow measurements of the core and HST-1 knot fluxes separately at these wavelengths, we have utilised the HST data to obtain fluxes of the core and HST-1 knot.

As in our previous work M87 MWL2017, we have performed photometric decomposition to exclude the host galaxy flux from the image. The host was approximated with two Sersic functions (Sersic 1968) with the core and jet regions masked out. Unfortunately, due to the small FoV of STIS, no stars were available in the images for constructing a PSF. The HST PSF library does not contain PSFs for the STIS instrument, and the HST team advises against using the TinyTim package¹ to build PSFs for STIS images. Nevertheless, we have performed a decomposition without taking into account the PSF, rather by masking the central part of the host galaxy and performing the decomposition of outer parts of the galaxy and extrapolating the result to the centre.

Figure 5 shows the residual image (after host model subtraction), where the core and the HST-1 knot are clearly visible. To obtain the fluxes of these features, we performed aperture photometry with an aperture radius of $0\stackrel{″}{.}4$ (the same as in M87 MWL2017). Table C.6 gives the fluxes of the core and HST-1 knot before and after host model subtraction. The uncertainties of fluxes before correction for the host galaxy were calculated based on the STIS documentation², while we assume uncertainties of the corrected fluxes of 50%, taking into account the method of the decomposition. However, as can be seen in Table C.6, contributions of the host galaxy are small since the band is blue, while the host galaxy is an old elliptical type. The results indicate that the HST-1 knot is significantly fainter than the core, that most likely was also the case during the EHT MW campaign in 2018.

Fig. 5. Hubble Space Telescope image of M87 at 2359.7 Å with the host galaxy subtracted. The core and HST-1 are indicated; the distance between the features is 1$\stackrel{″}{.}$00±0$\stackrel{″}{.}$05, and the pixel size is 0$\stackrel{″}{.}$0248.

2.2.4. Far-UV observations from AstroSat-UVIT

The UV Imaging Telescope (UVIT) instrument on board the AstroSat mission consists of twin telescopes, each of 38 cm diameter, one dedicated to far-UV (FUV: λ = 1300–1800 Å) and the other to near-UV (NUV: λ = 2000–3000 Å) & Visible (VIS: λ = 3200–5500 Å) channels (Kumar et al. 2012). In this work, the far-UV images from UVIT are used to constrain the jet and the HST-1 knot component in the UV emission. The images created for each orbit of the satellite were aligned and merged to create the final image with an effective exposure time of 14.1 ks. The flux contributions of different known components of M87 are extracted by 2D image decomposition modelling. Figure 6 shows the photometric slice made through the galaxy centre in the jet direction obtained in the BaF2 band (F154W: λ_eff ∼ 1541 Å, with Δλ_eff ∼ 380 Å). In this figure the observed galaxy flux is shown with a thick black line, the total model (a sum of all components) with a red line, the host galaxy model with an orange line, the total jet model, that include the asymmetric Core-Sersic model (Graham & Driver 2005) and a point source for the HST-1 knot, is shown with a green line, and the central source with a blue line. A detailed description of data, analysis and modelling is given in Appendix B.2.3.

Fig. 6. AstroSat-UVIT photometric cut of the decomposition model along the jet direction. The thick black line shows the observed flux along the cut, the blue line is the central point source, the orange line is the host galaxy, the green line is the jet total flux, and the red line is the total model flux. The shaded area shows the masked out region of the jet that has a complex fine structure that cannot be accurately modelled with a simple analytical function and was therefore excluded from the decomposition to reduce the model complexity.

2.3. X-ray observations

2.3.1. Chandra and NuSTAR observations and joint spectral analysis

Observations of M87 with the Advanced CCD Imaging Spectrometer (ACIS) on board the Chandra X-ray Observatory were collected via Director’s Discretionary Time (DDT; PI: Wong) on 22 April 2018 (ObsID 21075 for 9.1 ks starting at UT 00:09) and 24 April 2018 (ObsID 21076 for 9.0 ks starting at UT 13:20). An observation with NuSTAR (Harrison et al. 2013) was also obtained via DDT on 25 April 2018 (ObsID 60466002002 for 21.1 ks beginning at UT 13:56), sufficiently close in time to the second Chandra observation that they can be considered contemporaneous. The data from these observations have been included in previous studies of X-ray flux variability in M87 described in Yang et al. (2019) and Imazawa et al. (2021).

Following the procedure in M87 MWL2017, the two Chandra observations are combined in our spectral analysis such that the reported fluxes represent the average for 22–25 April 2018. The variability of the background-subtracted count rates for each of the spatial components is shown in Fig. 7. As is clear from the Chandra light curve, the core was significantly brighter in both observations than in April 2017. No significant variability is detected within either of the 2018 observations. However, the core flux did increase between the two observations separated by approximately two days, while the other components remained steady (see Fig. 7). Moreover, an increase in the X-ray brightness of M87 relative to 2017 is apparent directly from the NuSTAR spectrum, that exhibits higher S/N out to ≃60 keV than the 2017 data (cf. M87 MWL2017, Fig. 9), despite equal exposure time.

Fig. 7. X-ray light curves of the spatially resolved components of M87 from the Chandra observations in 2018. The core (orange) shows clear variability between the two different days on which observations have been carried out, while the HST-1 knot and the rest of the jet remain fairly steady.

In particular, we detect a significant increase in the flux of the core component relative to 2017 observations. The unabsorbed core flux in the 2–10 keV band averaged between the two Chandra observations in April 2018 was 2.6 × 10⁻¹² erg s⁻¹ cm⁻² that represents an ∼80% increase on a timescale of one year. Assuming a distance of 16.8 Mpc, we find the average 2–10 keV band luminosity of the core to be (8.8 ± 0.4)×10⁴⁰ erg s⁻¹, again a marked increase over (4.8 ± 0.2)×10⁴⁰ erg s⁻¹ observed in 2017. The significantly higher luminosity is not seen in the HST-1 knot or the rest of the jet, that both show ≲10% variability. These results are similar to those obtained by Imazawa et al. (2021). Despite the flux increase, we note that there is no evidence for a change in the shape of the spectrum. Our best fit power-law (PL) photon index for the core is ${\mathrm{\Gamma }}_{\mathrm{core}}=2.{03}_{-0.07}^{+0.12},$ that is statistically consistent with ${\mathrm{\Gamma }}_{\mathrm{core}}=2.{06}_{-0.07}^{+0.10}$ from the 2017 EHT observations.

In addition to the observations taken within the EHT MWL campaign, we also examined data taken with Chandra earlier in 2018: OBSIDs 20488 (4 January 2018) and 20489 (21 March 2018). These observations were processed and analysed in exactly the same way as described above in order to evaluate the core flux using the same set of assumptions. We find that the core flux in the 2–10 keV band was lower in both epochs, ≃2 × 10⁻¹² erg s⁻¹ cm⁻², as shown in Fig. 8.

Fig. 8. Mid-term MWL light curves of M87 taken during the observational campaign covering January 2017 to April 2019. Here are listed the observations used in the plot and their reference (in case of already published results): Chandra (0.3–7 keV) (Sun et al. 2018), (2.0–10 keV) (Imazawa et al. 2021) H.E.S.S. (2019) (H.E.S.S. Collaboration 2023); VLBA, EAVN, VERA (2017–2018) (Cui et al. 2023); All instruments 2017 data (M87 MWL2017). We mark the VHE flare in 2018 with a grey-shaded region in the background.

Although a PL model provides a satisfactory fit to the April 2018 data alone and is suitable given the sensitivity of the spectrum in a typical single observation, the growing archive of NuSTAR observations of M87 permits a deeper look at the shape of the X-ray spectrum of the core. While the shape of the spectrum appears to be stable, there is evidence of deviation from a PL. In particular, our analysis of 14 existing NuSTAR observations (Sheridan et al., in prep.) reveals statistically significant evidence for a break in the PL around E_br = 10 keV: compared to PL fits, models with a break improve the Cash fit statistic by ΔC > 150 with the addition of only two free parameters. Relative to the pure PL model, the spectrum is softer below the break (Γ₁) and harder above the break (Γ₂), though the exact values depend on how the data are combined. For example, when all observations are stacked, we find ${\mathrm{\Gamma }}_1=2.{26}_{-0.05}^{+0.06}$, E_br = 11 ± 1 keV, and ${\mathrm{\Gamma }}_2=1.6_{-0.2}^{+0.1}$. When the spectral shape is tied between observations but the flux is allowed to vary, we find ${\mathrm{\Gamma }}_1=2.{20}_{-0.07}^{+0.08}$, E_br = 12.5 ± 1.2 keV, and Γ₂ = 0.9 ± 0.2. We emphasise that Γ₁, E_br, and the existence of a break are weakly sensitive to the stacking method, while Γ₂ depends more on the details. In our SED modelling (Sect. 4) we therefore consider all three spectral models for the core. Additional details are included in Appendix B.3.1.

2.3.2. Swift-XRT observations and analysis

As a part of the EHT campaign, M87 was observed with Swift-XRT from 18 April (MJD 58226) to 29 April 2018 (MJD 58237). A total of 9 observations were conducted in April, with one additional observation on 18 December 2018 (MJD 58240). To put the 2018 observations in context, we also analyse all Swift-XRT photon counting mode observations for most of the Swift-XRT mission duration, extending back to 2009 and forward to 2021. The details of our data processing and spectral analysis are given in Appendix B.3.2.

Analysing these observations constrains both the long-term and short-term variability of the X-ray emission from M87, even though Swift-XRT cannot spatially resolve the inner jet, including the core and the HST-1 knot. The resulting light curve for the 2018 observations is shown in Fig. 9; we include total unabsorbed flux estimates for M87 with background, as well as net unabsorbed flux estimated by removing the contributions from the background model components. The small-scale components’ flux is found to be consistently between 40% and 60% of the total. However, the background-subtracted fluxes shown in Fig. 9 should only be interpreted as upper limits on the X-ray emission from the core, the HST-1 knot, and the jet.

Fig. 9. Total flux (blue) and unresolved summed component fluxes (red) of M87 in the 2–10 keV band observed with Swift-XRT in 2018.

2.3.3. AstroSat-SXT observations and analysis

AstroSat performed two pointing photon counting (PC) mode observations of M87, the first during 17–18 February 2018 (ObsID: A04_115T02_9000001900; MJD 58166.8) and the second during 21–23 April 2018 (ObsID: A04_115T02_9000002042; MJD 58229.4).

We use a model similar to the one adopted for Chandra and NuSTAR spectral analysis with the addition of a vvapec component to account for the diffuse emission from cluster gas because the SXT samples a larger portion of it. The best-fit normalisation obtained from the SXT spectra is higher than the constraints in M87 MWL2017 by about a factor of three mainly for this reason (additional details appear in Appendix B.3.3).

2.4. Gamma ray observations

2.4.1. Fermi-LAT observations

We performed a dedicated analysis of the Fermi Large Area Telescope (LAT) data (Atwood et al. 2009) of M87 using 12 years of data starting on 4 August 2008, in an energy range between 100 MeV and 1 TeV. We repeated the analysis using a time interval of one month centred on the EHT observation period (i.e. 8 April–8 May 2018). This choice is also motivated by a problem that affected the solar panel on 16 March 2018, that stopped the observations until 8 April 2018. We summarise these analyses here and include additional details in Appendix B.4.1.

The results from 12 years present a significant excess above the noise level (a test statistic of TS = 2015, corresponding to a significance³ of > 40σ with an averaged flux of F = (1.74 ± 0.11)×10⁻⁸ ph cm⁻² s⁻¹. The obtained spectral behaviour (Γ = 2.05 ± 0.03) is in agreement with the 4FGL-DR2 result for this source (Abdollahi et al. 2020).

For the period 8 April – 8 May 2018, the source is detected with a significance of ∼6σ (TS = 38), F = (5.0 ± 1.7)×10⁻⁸ ph cm⁻² s⁻¹, three times higher than the average over 12 years. It also presents a soft spectrum with Γ = 2.40 ± 0.21. Regarding the variability study, using for the light-curve study either one-day or two-day cadence temporal analyses, the source was detected (with more than three-sigma significance) in both of the analyses only in the one-day (and two-day) bin around the VHE flare episode. Namely, on the one-day bin centred on MJD = 58228.0 (20 April 2018), presenting a significance of 3.5σ and flux F = (23 ± 10)×10⁻⁸ ph cm⁻² s⁻¹, as well as in the two-day bin centred on MJD = 58227.5 (19–20 April 2018), with a significance of 3.4σ and flux F = (14 ± 6)×10⁻⁸ ph cm⁻² s⁻¹, respectively. Similarly, on both the one-week time scale and the customised time intervals (see Fig. 10 and Table C.13) the source is significantly detected (with a TS = 20) only in the week MJD = 58224–58231 (18–24 April 2018), as well as in a four-day period of the VHE flare, namely MJD = 58227–58231, presenting a flux of (8.1 ± 2.6)×10⁻⁸ ph cm⁻² s⁻¹ and (10.8 ± 3.8)×10⁻⁸ ph cm⁻² s⁻¹, respectively. The weekly flux represents the highest weekly emission observed for this source so far. The second-highest weekly flux ((7.2 ± 3.1)×10⁻⁸ ph cm⁻² s⁻¹) was registered in the first week of February 2014, during which the only VHE observations were obtained with MAGIC. The MAGIC data did not show a VHE flare at this time. Finally, for the four-day period of the VHE γ-ray flare, we performed a spectral analysis of the LAT data finding a PL index of Γ = 2.18 ± 0.33, compatible with the one measured in the 12-year analysis.

Fig. 10. Fermi-LAT light curve of M87 on customised time intervals during April 2018. The red point refers to the flux estimate, while grey points indicate upper limits. The 1σ upper limit is reported when TS < 10. The vertical blue line (band) indicates the time of the VHE flaring episode, while the horizontal lines represent the averaged flux over 12 years and its uncertainty F12 yr = (1.74 ± 0.11)×10−8 ph cm−2 s−1.

2.4.2. VHE observations: H.E.S.S., MAGIC, VERITAS

H.E.S.S. observations in 2018 were performed with the four 12 m CT1-4 telescopes (see Table C.14). A total of 13.0 hours of live time of good-quality data were obtained. The average zenith angle of the observations was 42.9°. A total statistical significance of 14.8σ was obtained for the 13.0 hours. For the systematic errors we conservatively estimated an uncertainty of 20% on the flux normalisation and 0.1 on the photon spectral index (H.E.S.S. Collaboration 2006b). Detailed information on the H.E.S.S. sensitivity and angular resolution can be found in Parsons & Hinton (2014) and H.E.S.S. Collaboration (2020).

MAGIC observed M87 for 11.3 hours after quality cuts for this campaign. The observations were carried out at a mean zenith angle of 29.6°. The total statistical significance from MAGIC observations during this campaign is 4.0σ. Based on the detailed descriptions in Aleksić et al. (2016) we estimate a systematic uncertainty of ∼30% on the integral flux that consists of 11% uncertainty of the flux normalisation, ±0.15 on the spectral index and 15% uncertainty of the energy scale. The long-term MWL light curves presented later in Figs. 3 and 8 also include the MAGIC data collected in 2019. In 2019 MAGIC observed M87 for a total of 29.1 hours after quality selection. The statistical significance obtained from the 2019 observations is 7.0σ.

VERITAS (Holder et al. 2006) observed M87 during April 2018 for 10.5 hours after quality cuts with an average zenith angle of about 22.9°. The total statistical significance from VERITAS observations during this campaign is 7.8σ. We estimated the systematic uncertainty to be ∼25% and ±0.2 on the absolute flux level and spectral index, respectively (Adams et al. 2022). Further details of the individual analyses can be found in Appendix B.4.2.

The light curves of nightly integrated flux measurements from all VHE instruments above 350 GeV between 10 April 2018 (MJD 58218) and 25 April 2018 (MJD 58233) are shown in Fig. 11 with solid colours. We assumed a PL index of Γ_VHE = 2.3 for the integrated flux calculations. We chose the same analysis threshold of 350 GeV as for our 2017 campaign (EHT MWL Science Working Group 2021). Strong moonlight conditions prevented observations after 25 April 2018. In order to investigate possible flux variability, we fitted the flux measurements from all VHE telescopes with a constant flux model. This yields a χ²/d.o.f. of 42.61/16 corresponding to a probability of p ≈3 × 10⁻⁴. In order to estimate the variability time scale, we fitted the light curve with a Gaussian function:

$$\begin{array}{c}\hfill F(t)=F_0+F_1·exp(-\frac{{(t-t_0)}^2}{2{\sigma }^2}),\end{array}$$(1)

Fig. 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 2018. The solid and semi-transparent colours represent the original data and scaled (with cross-normalisation factors) data, respectively. The black curve and shaded region represent the best Gaussian fit and 1σ error to the original flux points. The grey lines indicate the percentage of the Crab Nebula’s flux level from integrating its SED in the same energy range. The SED of the Crab Nebula is obtained by fitting VERITAS’s six years of observations (Meagher & VERITAS Collaboration 2015). For flux points with significance less than 2σ, upper limits are given in Table C.14.

where F₀ was added as baseline flux. The Gaussian fit provides a better description of the measured data with χ²/d.o.f. of 16.46/13 corresponding to a probability of p ≈ 0.23. The best fit results are F₀ = (1.25 ± 0.46)×10⁻¹² ph cm⁻²s⁻¹, F₁ = (2.83 ± 0.62)×10⁻¹² ph cm⁻²s⁻¹, t₀ = 58229.07 ± 0.39 MJD (21 April 2018), σ = 1.32 ± 0.44 days. Defining the variability time scale as the full width at half maximum (FWHM), we obtain $\mathrm{FWHM}=2\sqrt{2ln2}\sigma \approx 2.35\sigma =3.10\pm 1.03$ days. For comparative purposes, the rising part of light curve was also fitted with an exponential function of the form

$$\begin{array}{c}\hfill F(t)=p_0{e}^{-|t-58226|/\tau },\end{array}$$(2)

where the exponential time scale τ is 1.88 ± 0.27 days for the rising flux portion. Furthermore, a Gaussian function fit to H.E.S.S. data alone yields a χ²/d.o.f. of 0.44/2 corresponding to a probability of p ≈ 0.80, F₀ = (0.77 ± 1.06)×10⁻¹² ph cm⁻²s⁻¹, F₁ = (4.50 ± 0.95)×10⁻¹² ph cm⁻²s⁻¹, t₀ = 58229.15 ± 0.11 MJD, and σ = 1.28 ± 0.34 days.

Light curves with shorter time binning of less than one night were calculated, but no significant intra-night variations were found. Further discussion of possible systematic cross-normalisation factors between IACTs and how the factors affect the fitting results can be found in Appendix B.4.2. To conclude, the statement that a model with a Gaussian peak describes the light curve significantly better than a constant holds even taking potential systematic cross-calibration offsets between the VHE instruments into account. The cross-normalisation factors provide a possible range of fitting parameter values and reduce the tension of the flux difference between VERITAS and H.E.S.S. around MJD 58229 (semi-transparent colours in Fig. 11). However, for this research, there is no compelling reason for a rescaling of the original data, reaffirming the fact that the fit to the flare profile without scaling is acceptable.

Individual SEDs obtained with all VHE instruments are shown in Fig. 12 and results of the spectral fits are given in Table C.15. The photon spectra are described with a simple PL function:

$$\begin{array}{c}\hfill \mathrm{d}N/\mathrm{d}E=f_0·{\left(\frac{E}{E_0}\right)}^{-{\mathrm{\Gamma }}_{\mathrm{VHE}}},\end{array}$$(3)

Fig. 12. Very high-energy SEDs measured with H.E.S.S., MAGIC, and VERITAS during the coordinated MWL campaign in 2018. Error bars are equivalent to the 1σ confidence level. An upper limit is drawn where the statistical significance is less than 2σ (i.e. 95% confidence level). The thick solid lines represent the best spectral fitting results by using a simple PL function. The shaded regions represent the 1σ error region of the fitting results. Historical spectral fits are taken from H.E.S.S. Collaboration (2006a,b), Aleksić et al. (2012), Acciari et al. (2008, 2010), Albert et al. (2008), Aliu et al. (2012), Beilicke & VERITAS Collaboration (2012), MAGIC Collaboration (2020) and the thin black line indicates the SED obtained with the MAGIC measurements presented in M87 MWL2017.

where f₀ is the normalisation at the energy E₀ and Γ_VHE is the photon spectral index. It is important to note that the three VHE SEDs are results of different samplings during the observational campaign. MAGIC observed primarily during nights before and after the VHE flux rise discussed below. About half of the VERITAS data were taken before and the other half during the flare. The H.E.S.S. observations cover the whole time span of the flare and lead to the hardest spectral index amongst the three instruments (Fig. 12). This result points to a harder-when-brighter behaviour in VHE, as previously hinted in H.E.S.S. Collaboration (2006a), Albert et al. (2008), Acciari et al. (2010), Aliu et al. (2012), H.E.S.S. Collaboration (2023). We acknowledge that directly comparing the spectral slopes of H.E.S.S., MAGIC, and VERITAS is also limited by the difference in the energy ranges covered by their respective spectra. Differently than the results from H.E.S.S. Collaboration (2024), the spectrum derived from the H.E.S.S. data in this study is limited to a maximum energy of 6.9 TeV. This limitation ensures that each energy bin in Fig. 12 is detected with a significance of at least 2σ.

3. Results

3.1. Multi-wavelength images in 2018

In Fig. 13, we show a composite of M87 MWL images obtained by radio and X-ray instruments during the 2018 campaign, including the EHT image. The M87 jet is imaged at multiple scales from ∼1 kpc down to a few Schwarzschild radii (defined as r_s = 2GM_BH/c² = 2r_g where M_BH is the mass of M87*, G the gravitational constant, c the speed of light, and r_g is the gravitational radius) summarising a contemporaneous MWL view of M87 during the 2018 EHT campaign.

Fig. 13. Composite of the M87 MWL images at various scales obtained in radio and X-rays during the 2018 campaign. The instrument, observing wavelength, and scale are shown at the top-left side of each image. We note that the colour-scale has been chosen to highlight the observed features for each scale and should not be used for noise levels, dynamic range, or flux density calculation purposes.

To check the structural variations, we compared the data obtained in 2017 and 2018. While the appearance of the kpc-scale jet remained relatively consistent during this period, our high resolution VLBI observations unveiled some morphological evolution in the pc-scale jet between the two years (see Fig. 14). Notably, these changes encompassed a counterclockwise shift in the jet position angle (PA), occurring within a few milli-arcsecondfrom the core (PA ∼ 281°/283° in 2017/2018; see Cui et al. (2023) for details). Furthermore, the data from 2018 reveal a more symmetric brightness ratio of the north-south jet limbs. All of these variations suggest the existence of year-scale structural evolution transverse to the jet. Moreover, regular VLBI monitoring observations reveal a gradual decrease in core flux (by 10–15%) at 22–129 GHz from 2017 to 2018 as shown in Fig. 2. Closer to the SMBH, GMVA 86 GHz observations in concert with ALMA and GLT spatially resolve the radio core, revealing a ring-like structure connecting to an edge-brightened jet (Lu et al. 2023).

Fig. 14. M87 jet structure comparison based on yearly stacked EAVN images at 43 GHz in 2017 and 2018 (see also Cui et al. 2023). The restoring beam size is 0.5 mas circular Gaussian, as indicated in the bottom-left corner. The image has been rotated by −18°. Contour levels are scaled as (1, 2, 4...)×C1st, where C1st is 1.8 mJy beam−1 in 2017 and 1.3 mJy beam−1 in 2018. For M87* mass MBH = 6.5 × 109 M⊙ (Event Horizon Telescope Collaboration 2019a), 1 mas ≈ 125 rs ≈ 0.08 pc, where rs is the Schwarzschild radius.

Contemporaneous observations obtained with the VLBA at 24 and 43 GHz were used for the spectral analysis (Sect. 4). Figure 15 compares maps of the spectral index (α; defined as F_ν ∝ ν^α) between these two frequencies from observations on 5 May 2017 (panel (a)) and 28 April 2018 (panel (b)). The M87 jet has an edge-brightened structure and does not have distinct components at each frequency to anchor an alignment of different images. We model-fitted the jet structure by a number of two-dimensional Gaussian components, aligned images at the position of the 43 GHz core, and shifted 24 GHz images by the value of the core shift. Finally we measured the shifts from the phase referencing observations, that are 0.036 and 0.012 mas in right ascension and declination, according to Hada et al. (2011) and Ro et al. (2023). Even though the amplitude of the core shift may change with time (Plavin et al. 2019), multiple studies of M87 indicate the magnitude of the core-shift variation in 22–43 GHz is within 10 μas and negligible (Hada et al. 2012, 2014; Jiang et al. 2021). The core regions in both 2017 and 2018 show relatively flatter spectral index values (α ≈ 0) than jet regions. However, compared to 2017 observations (Fig. 15a; M87 MWL2017), the value of the spectral index in 2018 steepens faster downstream from the core region to the jet region. It indicates that the radio emission in the jet close to the core becomes relatively stronger at lower frequencies compared to the emission in the core in 2017, that might be related to changes in the physical conditions of the particle accelerator at the core. In the extended jet regions, the spectral index distribution appears to be more complex: it shows flatter values (α > −0.5) at the location of total intensity maxima with steeper values (α < −1) that are localised in between. The M87 jet structure and spectral index distribution based on the 8–15 GHz VLBA observations with the highest-to-date dynamic range (Nikonov et al. 2023) are consistent with the development of Kelvin-Helmholtz instabilities in the jet in regions of 10² − 10⁴ r_g from the central engine. Their propagation can produce higher pressure regions near the jet boundary and in its interior, causing variations of the spectral flattening towards the jet edges, as shown in Fig. 15. The observed variability in flat- and steep-spectrum regions could also be influenced by the jet precession, that modulates the observed jet structure over time. However, for the variability results obtained here, rescaling of the original data is not crucial, reaffirming the fact that such a complicated spectral profile in the edge-brightened jet region could be (at least) partly affected by insufficient interferometric (u,v) coverage (e.g. Pashchenko et al. 2023).

Fig. 15. Spectral index maps obtained between 24 and 43 GHz from VLBA 5 May 2017 (a) and 28 April 2018 (b) observations. The restoring beam size is 0.55 mas ×0.55 mas, that is the circular equivalent-area synthesised beam of VLBA at 24 GHz, drawn as a grey circle in the bottom-right corner. The image has been rotated by −18°. The contours denote the total intensity at 43 GHz and start at 3.5σ rms, increasing in steps of two. The (u, v)-ranges match between all images.

3.2. Change in the mas-scale jet and μas-scale ring asymmetry position angles

As reported in Sect. 3.1, a slight and counterclockwise variation in the PA of the mas-scale jet was detected with 43 GHz VLBI observations in 2017 and 2018. The observed PA evolution from 2017 to 2018 appears to follow the long-term (quasi-)periodic variation (on time scales of ∼10 years) found by the recent VLBI jet monitoring studies (Walker et al. 2018; Cui et al. 2023). Possible scenarios causing such year-scale variations of the jet PA include MHD/HD instabilities (Mizuno et al. 2007; Walker et al. 2018) and Lense-Thirring precession in a tilted accretion disk system (Fragile et al. 2007; Liska et al. 2018; Cui et al. 2023). The observed variation in the north-south limb brightness ratio could also be related to the jet PA variation and/or transient turbulence in the jet flow, which requires further dedicated investigations.

Similarly, the EHT-2018 images (M87 2018 I) reveal a significant shift of the PA of the brightness asymmetry of the ring with respect to that of EHT-2017 images (presenting a counterclockwise variation ∼30° from 2017 to 2018). The shift direction is as expected if the large-scale jet is aligned with the disk rotation axis (M87 2017V). Morphological variations of the asymmetric ring at EHT scales could be related to the week-scale variation caused by magnetorotational instabilities (M87 2017V, M87 2018 I). On the other hand, the observed PA change in the EHT images between 2017 and 2018 may imply the presence of a longer-term (year-scale) variation at this scale as originally suggested by Wielgus & Event Horizon Telescope Collaboration (2020).

Nevertheless, it is worth noting that M87 shows counterclockwise PA variations at both milli-arcsecond and micro-arcsecond scales between the two years, that sheds light on the possible physical connection between the mas-scale jet and μas-scale ring. The recent detection of both the ring-like structure and extended jet in a single image using global VLBI at 86 GHz offers a new potential to explore this challenge by bridging the spatial gap between μas- and mas-scale structures (Lu et al. 2023). Future multi-year comparisons of multi-frequency (230/86/43 GHz) VLBI images will be the key to constraining the origin of variable structures at different scales and their physical connection.

3.3. Multi-wavelength light curve and fractional variability

In Figs. 3, 8, and 4 we show the long-term (2000–2022), mid-term (2017–2019) and short-term (April 2018) MWL light curve behaviour of M87 with the observations collected by the instruments participating in the coordinated MWL campaign in 2018.

Comparing the 2018 EHT-MWL campaign (Fig. 4) to 2017 (M87 MWL2017), the source presents brighter emission at high energies, seen both in the X-ray as well as in the γ-ray energy bands, while in optical and at radio frequencies the core flux remained the same or even a bit dimmer. In particular, Chandra observations during April 2018 revealed a core flux increase of ∼80% with respect to the previous 2017 campaign. While the duration of the VHE γ-ray flaring episode is well constrained (see Sect. 2.4.2 and further discussion in Sect. 3.3.1), the lack of monitoring X-ray observations of the source does not allow us to fully investigate the duration of the X-ray variability, as we further elaborate in Sect. 3.3.2.

The fractional variability (Vaughan et al. 2003) of the VHE γ-ray band in 2017 is compatible with zero within uncertainties. The same applies to the 2018 VHE γ-ray variability. This is caused by the large statistical uncertainties of the VHE data. We note that these large statistical uncertainties are also limiting the calculation of a discrete correlation function (Edelson & Krolik 1988). The fractional variability of the Swift-XRT flux of the small-scale components (namely the core, the HST-1 knot and inner jet) is significantly lower in 2018 than in 2017. For both years these fractional variabilities are slightly higher than in the VHE γ-ray band, although well compatible within statistical uncertainties. Further information on our calculation of the fractional variability are reoperted in Appendix D. The resulting fractional variability values are given in Table D.1.

We combined the spectral slopes and flux normalisations from the three VHE energy spectra obtained in this campaign with the historical measurements to test whether there is a general harder-when-brighter behaviour in the VHE band. The photon indices are plotted against the flux normalisations in Fig. 16, in which it is possible to see that there is no dependence of the spectral index on the flux. However, this comparison is limited by the fact that the VHE spectra shown cover different energy ranges (H.E.S.S. Collaboration 2006a; Acciari et al. 2008; Albert et al. 2008; Acciari et al. 2010; Aliu et al. 2012; Aleksić et al. 2012; Beilicke & VERITAS Collaboration 2012; MAGIC Collaboration 2020; EHT MWL Science Working Group 2021).

Fig. 16. Photon PL index Γ compared to the VHE flux normalisation at 1 TeV. Data from H.E.S.S. Collaboration (2006a), Acciari et al. (2008), Albert et al. (2008), Acciari et al. (2010), Aliu et al. (2012), Aleksić et al. (2012), Beilicke & VERITAS Collaboration (2012), MAGIC Collaboration (2020), EHT MWL Science Working Group (2021). Flux normalisations at different energies have been recalculated at 1 TeV using the respective power-law index. No dependence of the spectral index on the flux was observed.

3.3.1. VHE gamma ray flaring episode

The first strong indication of VHE γ-ray emission from M87 was found in 1998 with the High Energy Gamma Ray Astronomy (HEGRA) telescopes (Aharonian et al. 2003). Since 2004, frequent monitoring (Acciari et al. 2009; Abramowski et al. 2012) of M87 has been conducted by H.E.S.S. (H.E.S.S. Collaboration 2006a), MAGIC (Albert et al. 2008; Aleksić et al. 2012; MAGIC Collaboration 2020), and VERITAS (Acciari et al. 2008, 2010; Aliu et al. 2012; Beilicke & VERITAS Collaboration 2012). These monitoring campaigns have yielded detections of VHE flaring activity from M87 in 2005 (H.E.S.S. Collaboration 2006a), 2008 (Albert et al. 2008), and 2010 (Abramowski et al. 2012). Furthermore, while VERITAS observations carried out in 2011 and 2012 did not show bright flares, they reported a hint of an underlying quiescent increase of the VHE emission evolving on longer time scales than day-scale rapid flares (Beilicke & VERITAS Collaboration 2012).

Despite deep investigations of past flaring episodes of the source as well as of the quiescent state during the 2017 EHT campaign (M87 MWL2017), the origin of the VHE γ-ray emission from M87 is still unclear. A unique opportunity to address this question is provided by the correlated temporal variations between the VHE and lower-energy MWL emission, particularly when the source can be spatially resolved (Rieger & Aharonian 2012; Rieger & Levinson 2018). Fortuitously, as shown in Sect. 2.4.2, the 2018 MWL campaign detected the first fast γ-ray variability observed since 2010, providing a chance to investigate the origin of the γ-ray emission in M87.

During the 2018 MWL campaign the VHE flux increased by a factor ∼3.3 with respect to April 2018 low state (see Sect. 2.4.2, Fig. 11) reaching a peak flux that is still far below the flares observed in the past but a factor ∼4 higher than the quiescent state flux measured during the 2017 campaign (M87 MWL2017). In 2005 the peak flux at energies > 350 GeV reached (1.41 ± 0.39)×10⁻¹¹ ph cm⁻² s⁻¹, almost 20 times the quiescent state flux as measured during the 2017 EHT-MWL campaign (M87 MWL2017). Similarly, in 2008 and in 2010, the observed fluxes were 23 and 37 times the quiescent state flux (see Fig. 3), respectively.

The variability time scales observed during the past VHE flares are similar to the 2018 flare (FWHM approximately three days) presented here, albeit somewhat faster. The 2005 and 2008 flares showed a roughly two-day and one-day variability, respectively (H.E.S.S. Collaboration 2006a; Albert et al. 2008). In 2010 a two-sided exponential fit revealed a rise time of 1.69 ± 0.30 days and a decay time of 0.611 ± 0.080 days (Abramowski et al. 2012). The general light curve shape of the 2018 flare is similar to the flare seen in 2010, with only one isolated peak, while the flares observed in 2005 and 2008 comprise multiple separated high states during a similar time period.

3.3.2. Flux enhancement in the X-ray band

The X-ray flux from the unresolved core as observed with Chandra (within 0.4 arcsec radius), the combined core, HST-1, and jet fluxes from Swift-XRT, and the Swift-UVOT B- and U-band fluxes are all generally higher in April 2018 compared to both 2017 and 2019, as shown in Fig. 8. However, without additional X-ray/UV/optical observations immediately before and/or after the gamma-ray flare period, it is difficult to establish a firm link between the elevated X-ray/UV/optical fluxes and the flare itself. Figure 3 shows that M87’s core X-ray flux can vary to a similar degree even in the absence of a flare, so the elevated X-ray/UV/optical flux near the 2018 flare cannot be definitively linked to the gamma-ray flare itself.

It is, anyway, tempting to relate the observed X-ray flux enhancement to the flaring observed in the VHE band despite the very sparse time coverage of Chandra observations. The difference in core flux between the two Chandra observations during the EHT 2018 campaign is clear but small (≃20%; see Fig. 7) and rising away from the VHE flare that peaked approximately one day and three days before the first and the second Chandra observations, respectively. Two AstroSat-SXT observations, that do not directly resolve the core, suggest a modest increase in flux between February and April 2018.

Denser sampling in time is afforded by our Swift-XRT observations (Fig. 9), providing a measure of variability in the joint flux of the small-scale components (core, the HST-1 knot, and the resolved jet), that cannot be resolved individually. This flux, modelled as a fraction of the total flux observed by Swift-XRT does not show any coherent flux increases comparable to the VHE flare within the ≃11-day campaign window covered with roughly daily cadence. The variance over this period (variability amplitude of 5–10%) is consistent with that of long-term Swift-XRT monitoring. Comparing the Swift-XRT small-scale component fluxes to the Chandra core flux, we find that the Swift-XRT small-scale components fluxes produced by our analysis are typically about double the core fluxes reported in Table 1 of Cheng et al. (2023). Our Swift-XRT small-scale component fluxes and the Chandra core fluxes show a similar trend, presenting a small flux increase between 22 and 24 April 2018.

A model-dependent way to test for an association between the X-ray brightening and the VHE flare would be to search for changes in the X-ray spectrum. For example, a harder X-ray spectrum than observed in non-flaring states (e.g., the 2017 EHT campaign) might be expected in some broadband SED models (discussed in detail in Sect. 4.3). While Cheng et al. (2023) claim to find hardening in Chandra observations from April 2018, we find no evidence for this effect in our joint spectral analysis with NuSTAR (possibly because we account for pileup and use a broader energy range to measure the photon index). In our power law model, we find ${\mathrm{\Gamma }}_{\mathrm{core}}=2.{03}_{-0.07}^{+0.12}$ in 2018, statistically consistent with ${\mathrm{\Gamma }}_{\mathrm{core}}=2.{06}_{-0.07}^{+0.10}$ from the 2017 EHT observations. Thus the only statistically significant change in the X-ray core is the flux increase described above. Finally, we note that this flux enhancement is insensitive to the X-ray spectral model: because the PL break lies above 10 keV, it does not dramatically alter the inferred 2–10 keV core flux.

3.4. Multi-wavelength SED

As discussed in the previous subsections, the source underwent a γ-ray flaring episode during the 2018 EHT-MWL campaign, with a hint of spectral hardening towards the highest TeV energies. While a moderate flux enhancement is also seen at X-ray energies with respect to 2017 EHT-MWL observations (M87 MWL2017), in the radio band, the source presents a comparable, or even lower, emission. These features can be seen in Fig. 17 that shows the 2018 MWL broadband SED (reported also in Table C.1), where the 2017 EHT-MWL campaign SED is shown in the background in grey for comparison. The fact that the source was in a low state at all frequencies while the HST-1 knot emission was subdominant to the core in the X-ray band allowed us to confidently associate the Chandra and NuSTAR total fluxes with the core emission in 2017 data; however, the observed X-ray and VHE variability seen in 2018 leads to uncertainties regarding their localisation.

Fig. 17. Observed broadband SED of M87 contemporaneous with the EHT campaign in April 2018, with the exception of GMVA (carried on February 2018) and HST (carried on July 2018) observations (see Table C.1). The fluxes were measured by various instruments, and they are highlighted with different colours and markers. Filled markers represent flux point estimates, while empty markers indicates ULs. In the X-ray band, the two sets of points represent the core flux from the models with a power law (dark purple) and broken power law (light purple; see Table C.12). Grey points represent the observed broadband SED of the 2017 EHT campaign (EHT MWL Science Working Group 2021). MAGIC observations were performed before the γ-ray flare, VERITAS observations only partially overlap the flare, while all H.E.S.S. observations have been taken during the VHE γ-ray flaring episode (see Table C.1). For the radio/millimetre VLBI and interferometric data, the upper limits on emission size (radius) for several representative frequencies, that are obtained by following the procedures described in M87 MWL2017, are labelled to clarify the constraints used in Sect. 4.

In Table C.1 we present the legacy data set with near-simultaneous broadband SED results collected over a broad range between a frequency of ≃20 GHz and photon energy ≃1 TeV. Similar to M87 MWL2017, we report for each given energy or frequency range the associated spatial scale of the total emission, the time ranges of the observations and fluxes (providing both F_ν and νF_ν). As described in Sect. 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. The angular scale of the emission measured by the different instruments involved in this campaign, ranges between ≃20 μas (EHT) and ≃2° (Fermi-LAT in low-energy γ-rays), namely a factor ≳10⁸ in angular scale.

Besides the EHT observation (for which we list the estimated total flux within a 60 × 60 μas² region as described in Event Horizon Telescope Collaboration 2019d), for other radio observations if the core is not resolved we use the flux peak obtained from the beam size as an upper limit on the core emission. For elliptical beams, we list the average of the axes as a representative angular scale. For the observations carried out in the optical and UV bands, due to the difference between Kanata and Swift-UVOT angular resolution (∼10″and ∼2.5″, respectively), we consider the first as an UL and the second as flux points, similarly to what was done for the Swift-UVOT observation carried out in M87 MWL2017. This difference is related to the larger uncertainties in the subtraction of the host galaxy contribution in the Kanata optical observations. However, we want to highlight that, considering its angular resolution, Swift-UVOT measurements include both the core and the HST-1 knot contribution. At higher energies, the scale typically corresponds to the point spread function (PSF) within a given band or the diameter of the region used for data extraction. For Chandra and NuSTAR, as in EHT MWL Science Working Group (2021), we show the inferred core spectrum from the power-law and broken power-law models described in Sect. 2.3.1. We discuss the implications of the inferred X-ray spectrum on the SED modelling in Sect. 4. Despite the fact that NuSTAR can only separate out the core emission spectroscopically (see Sect. 2.3.1), the reported scale for Chandra and NuSTAR data analysed jointly is related to Chandra’s greater resolving power. Given that Swift-XRT cannot distinguish the core from other components detected by Chandra, we treat its unresolved flux as an upper limit.

Considering the timing relative to the observed γ-ray flare it is important to note that the majority of the observations arranged during our 2018 campaign have been carried out either during or within a few days of the flaring episode, well within the minimum dynamical timescale (5–61 days; Satapathy et al. 2022) for the innermost stable circular orbit of M87. The only exception is represented by global 43 GHz VLBI observations, collected in February 2018, that have been included taking into account the non-variability seen in the radio band on short and long time scales (see Sect. 3.3). Although these observations were not all strictly simultaneous with the γ-ray flare and EHT observations, the non-variability observed for the radio, optical and X-ray bands, suggests that they effectively provide a snapshot of the broadband SED of the M87 core. These results can be used as constraints for any model seeking to explain the source behaviour around the time of the observed γ-ray flare and EHT 2018 image acquisition.

In addition to highlighting the various observatories and instruments with different colours and labels, we also report the upper limits on the size of the emission region for the VLBI measurements. These data are provided to the community in machine-readable form via the EHTC Data Webpage⁴ (or directly via DOI ≪TBD≫) along with associated documentation and files needed for spectral modelling of the X-ray data. A complete list of all supplementary material published along with this paper is provided in Appendix C.

4. Heuristic modelling of the 2018 SED

A variety of theoretical models have been proposed to explain MWL emission from M87 (e.g., Rieger & Aharonian 2012; Prieto et al. 2016, for review). Given the M87 polarisation imaging constraints (Event Horizon Telescope Collaboration 2021a,b), the innermost, ionised regions must be moving along strong field lines. While the accretion flow could contribute in some parts of the SED, the jet radiation is more likely to be dominant in the SED as a whole (e.g., EHT MWL Science Working Group 2021, and references therein). For this work, we choose to model the SED based on single zones, a common approach to many AGN blazar studies; we defer more sophisticated modelling such as structured jet models (e.g., Potter & Cotter 2013; Lucchini et al. 2019) to future works. Models based on single zones enable the baseline properties of the emitting regions to be constrained. In addition, this approach allows for comparisons to other M87 epochs, including the 2017 campaign, as well as to the other AGN observed by the EHT.

The origin of the VHE γ-ray emission from M87 is still unclear. The MWL observation campaign described above (see Sect. 2) offers a unique opportunity to tackle this issue, given the following key points: (i) the duration of the VHE γ-ray flaring episode is well constrained, while we cannot detect similar activity in the other energy bands, and (ii) EHT spatially resolved the photon ring and measured its flux during the flaring episode. With these data, the origin of the VHE γ-ray emission from M87 can in principle be better constrained. We briefly describe a heuristic model to explore this. First, as in M87 MWL2017, we introduce the EHT-oriented model. This model uses the flux density and the size of the emitting region obtained from actual EHT observations and allowed us to constrain the average physical state of the EHT emitting region near the black hole, assuming it is uniform within the region. Within the context of this approach, we find that the high energy γ-ray emission cannot be explained only by the EHT-oriented model. Therefore, the portion of the radiation that cannot be explained by the EHT-oriented model will be explained using a high-energy(HE)-oriented model, whose size is set by the VHE γ-ray flare duration.

Even with the unique 2018 observational data, simple HE-oriented models are degenerate in interpretation. We therefore present two model scenarios (Model A and B), both of which assume a TeV γ-ray flux increase during the flare. Model A is a more straightforward two-zone model that is a combination of the EHT-oriented model and the HE-oriented model. It is characterised by having mildly relativistic bulk flow and an internal emitting particle distribution (assumed leptonic) to explain the observed TeV γ-ray data. In contrast, Model B is a three-zone model with an additional fast-moving component (corresponding to a part of the accelerating jet) that explicitly accounts for the Γ_VHE≈2.0 component of the TeV γ-ray data. In this scenario, the γ-ray emission originates from the fast-moving component, as often proposed in the literature (e.g., Ghisellini et al. 2005; Giannios et al. 2009, and references therein).

Here the shape of the X-ray spectrum also becomes relevant. As discussed in 2.3.1 (see also Sheridan et al., in prep.), a stacked analysis of NuSTAR observations of M87 including the 2018 data suggests a break in the X-ray PL spectrum that (i) appears to be stable from year to year and (ii) bends upward towards the γ-ray band. Extrapolating this broken power-law (BPL) to higher energies results in a substantially larger MeV energy flux compared to the decaying PL model. Because the 2018 data are statistically consistent with both the PL model and the BPL model, in subsequent sections we consider both scenarios when modelling the full SED with Model A (see Table 1, Sect. 4.4, and Fig. 18). We are left with three separate models explored in the sections below: model-A with PL, model-A with BPL, and model-B with PL. Unfortunately, the most notable difference between the PL-based models and BPL-based models is a large MeV bump that is not sampled by our observing campaign (see, e.g., Fig. 18), and so ultimately we are not able to rule out either X-ray model based on SED considerations.

Fig. 18. Modified Model A fitted describing the BPL emission seen in the Chandra–NuSTAR X-ray spectrum.

The SED models include parameters relating to the source bulk physical properties as well as the internal emitting particle properties. The microphysics regarding the acceleration of leptons is phenomenologically included in non-thermal electron distributions (hereafter we refer for ease to these as electron distribution functions; eDFs). For this paper, GR and γ–γ absorption effects are not considered in order to enable comparisons with the models in M87 MWL2017 and other single-zone models. In fact, H.E.S.S. Collaboration (2024) demonstrate that the absorption caused by the Extra-galactic Background Light (EBL) only becomes notable (τ_EBL > 0.2) in M87 for γ-ray energies exceeding 10 TeV. Additionally, the absorption attributable to the infrared light from the galaxy is also found to be negligible (τ_IR < 0.001) for TeV γ-rays. Lastly, the absorption by photons in the vicinity of the SMBH at 10–50 r_g, where r_g = GM/c² is the gravitational radius of the black hole of M87*, might have a more significant impact (Brodatzki et al. 2011; Cui et al. 2012). We have assessed the influence of synchrotron self-absorption (SSA) on the energy flux of TeV photons, utilising the best-fit parameters from our models (see Sect. 4.3). Our analysis indicates that SSA, at most, affects the energy flux by less than an order of magnitude. Consequently, this phenomenon does not significantly affect the interpretation of the findings presented in this study. However, we acknowledge that the exact locations of the individual zones proposed in our models relative to each other are unclear and could potentially impact the overall SED. Due to the complexity involved, this aspect remains unexplored in the current work for the sake of simplicity.

4.1. Observed differences between 2017 and 2018 source emission influencing the SED modelling

Before we begin to explain the SED modelling in detail, it is worth summarising the differences between the 2017 and 2018 source emission that influence the SED modelling. The most notable differences are the flux increase and a hint of a spectral hardening (Γ_VHE≈2.0) observed in the VHE γ-ray energy band during the April 2018 campaign. As reported in Sect. 2.4.2, the FWHM time scale of the VHE flux variation was ∼3 days, providing an important constraint on the size of the γ-ray emitting region in April 2018 that was not available in 2017, albeit degenerate with relativistic beaming effects.

Our interpretation of the X-ray data for the modelling also differs somewhat compared to M87 MWL2017. While some variability is implied by the flux enhancement in the X-ray band in 2018 in comparison with 2017 (see Fig. 8 as well as Sect. 3.3.2), there is no definitive evidence of its association with the VHE γ-ray flare. Therefore, we will examine both of the cases where the X-ray flux is associated with the EHT-oriented model and high-energy (E > 100 MeV) oriented model.

Finally, we treat the optical/IR flux observations as upper limits because we cannot separate out the galaxy contamination due to the inadequate resolution of the Kanata optical telescope. While we choose a value of the SED parameters to keep the model within the range consistent with the Swift-UVOT emission.

4.2. EHT oriented model

The following free parameters characterise the spherical, single-zone emission region: the radius (R), the Doppler factor of the bulk motion (δ), and the averaged magnetic field strength (B). The Doppler factor is defined as δ = 1/Γ_L(1 − β cos θ_view), where Γ_L is the bulk Lorentz factor of the emission region, and we set θ_view = 17° (Event Horizon Telescope Collaboration 2019a). The nonthermal eDF within the emission region is characterised by the electron number density (n_e), the spectral indices of the broken power law (p₁ and p₂), and the minimum/break/maximum Lorentz factors (γ_min, γ_brk, and γ_max), respectively.

Overall, the characteristics of the radio-to-millimetre emission in the observed broadband SED in April 2018 are similar to the SED in April 2017 (M87 MWL2017). Therefore, we followed a similar modelling methodology. First, we set R to be equivalent to the photon-ring size (Event Horizon Telescope Collaboration 2019a), that is, the single-zone effectively encompasses a sphere around the SMBH out to R ≈ 5r_g. Since the bulk speed of the jet has not yet become relativistic, we fix δ ∼ 1. We set γ_min = 1 as it has little influence on the SED, that is, the value of γ_min affects the energy density of nonthermal electrons U_e only if the photon index is softer than 2 (e.g., Kino et al. 2002; Finke et al. 2008). Therefore, we fixed the conservative γ_min for simplicity. Based on observed radio core shift measurements (Hada et al. 2011), the frequencies below ∼230 GHz are optically thick to synchrotron self-absorption (SSA), that constrains the magnetic field strength B> a few G (Event Horizon Telescope Collaboration 2021b). Keeping this limit in mind, one can fit the SED by only adjusting B and the eDF parameters. The best fit model parameters are summarised in Table 1. There are no essential differences between 2017 and 2018 for the EHT-oriented model.

Similar to the results obtained in M87 MWL2017, we find that the entire MWL radiation cannot be explained by the EHT-oriented model alone. In this model U_e and U_B are approximately comparable in the EHT emission region when γ_min = 1. Since we adopt γ_min = 1, U_e provides an upper limit. Although a detailed parameter survey will be left to future studies, U_B > U_e is anticipated for γ_min > 1 in the EHT emission region. This would align more closely with the scenario of a magnetically accelerated jet.

4.3. High-energy oriented model

Here we introduce the HE-oriented model to fit both the observed HE and VHE γ-ray spectrum. The variability time scale indicates that the characteristic size of the VHE γ-ray emitting region is limited to $R_{\mathrm{HE}}\lesssim 8r_g\delta \left(\frac{\mathrm{\Delta }t}{3\mathrm{days}}\right)$. From this, it is clear that there is a degeneracy between the observed variability timescale and δ. Hence, we will also consider cases of δ > 1 below. Similar to the EHT-oriented model, we first focus on the case of δ ∼ 1, because both a jet kinematic analysis using VLBI monitoring observations and an analysis of the brightness ratio of the approaching-jet and counter-jet suggest that the M87 jet is moving at subluminal speeds within the deprojected distance about 10³r_g from the central SMBH (Mertens et al. 2016; Park et al. 2019, and references therein). We took as reference for the Doppler factor a value δ = 1.82, that is consistent with the range allowed by the assumed viewing angle (θ_view = 17°) and able to reproduce the SED.

The 2018 SED VHE active phase further motivates us to explore the case of δ > 1. It is not straightforward to fit both the VHE γ-ray spectrum with Γ_VHE ≈ 2.0 and the Fermi-LAT spectrum with a single emitting zone. We therefore considered two approaches: In Model A, we modified the eDF parameters such that the modelled emission accounts for the tentative hardening of the VHE spectrum in April 2018 (i.e. p₁ > p₂), and in Model B, we took into account a fast moving component⁵.

The best-fit for Model A is shown in Fig. 19. We note that the physical origin is not obvious and further theoretical work is needed to assess its feasibility in the case of M87.

Fig. 19. Model A fitted to the observed contemporaneous broadband SED of M87 from the EHT campaign in April 2018. For each considered region (EHT and HE), the synchrotron (labeled as SYN) and synchrotron self-Compton (labeled as SSC) models as well as the summation of the two models (labeled as Sum) are reported.

In Model B, we retain the loss-dominated approach in the single-zone for the higher energy index (i.e. p₁ < p₂). We examine the phenomenological eDF with p₁ < p₂, that is complementary to the model A, since the radiative cooling rate may dominate the acceleration rate. In this model, we consider a fast moving blob as an additional component of the SED to explain the hard H.E.S.S. γ-ray spectrum in a flaring epoch. It should be mentioned that we need this additional component not due to the softening of the high-energy eDF, but due to the appearance of the decrease of VHE γ-ray emission because of the Klein-Nishina effect. To reproduce the H.E.S.S. γ-ray flux and avoid overshooting the observed X-ray fluxes during the flare, we set a very weak magnetic field strength of B = 4 mG and δ = 2.55 in the high-energy components. This model also looks consistent with the observed SED including the VHE γ-ray flare in Fig. 20.

Fig. 20. Model B fitted to the observed contemporaneous broadband SED of M87 for the EHT campaign in April 2018. For each considered region (EHT, high-energy, and very-high-energy flare), the synchrotron (labeled as SYN) and synchrotron self-Compton (labeled as SSC) models as well as the summation of the SYN and SSC models (solid line) are reported.

The resultant U_e/U_B ≃ 5.8 × 10⁴ shown in Table 1 is very high, requiring an extremely weakly magnetised plasma. The extent to which this could be attributed to magnetic reconnection is not yet clear as current reconnection models tend to suggest that the emission region is likely to be near equipartition (e.g., Petropoulou et al. 2016). Another possibility is that the interaction between the highly magnetised jets and the weakly magnetised wind (e.g., within the jet sheath) can result in the dissipation of the magnetic energy (Kusafuka et al. 2023, see also Komissarov 2012).

The estimated Doppler factor of Model B’s VHE-flare component is δ = 2.55 (the corresponding bulk Lorentz factor is ≲2 if the flaring component is moving along the global jet). The assumed Lorentz factor is consistent with the mildly high bulk Lorentz factor proposed by observations of M87: the gradual acceleration of the jet is suggested by high-cadence VLBI observations (Park et al. 2019), where the estimated Lorentz factor is ∼2 at the distance ≳10³r_g, and even at the HST-1 knot, the X-ray observations also suggest that the Lorentz factor of the jet is high ∼6 (Snios et al. 2019). The required jet power of the VHE flare component of the model B is ∼3.7 × 10⁴² ergs⁻¹. This is compatible with the jet power exhibited in GRMHD simulations.

4.4. X-ray broken power law oriented model

Fig. 18 shows the Model A fit for the case of BPL X-ray radiation. The inverted component on the HE side of the X-ray spectrum can be fitted by increasing the synchrotron flux in the HE-oriented model region. This leads to slightly different best fit values, since the eDF must be fine-tuned to reduce the increase in the synchrotron self-Compton (SSC) emission component corresponding to the increased synchrotron flux. It is clear that model A explains both PL and BPL cases naturally without altering the essential parts of the model. As shown by the figures in Table 1, the difference between the two cases can be attributed to the difference in the physical state of the HE region. When n′_e and B′ in the HE region increase by several times, the low-energy tail part of the emission from the HE region contributes to the X-ray energy band, resulting in the realisation of BPL. Furthermore, it should be emphasised that when BPL is realised in the X-ray energy region, the advantage is that tighter constraints can be placed on the model parameters of the HE emission region through SED fitting. As can be seen from the figure, the high-energy side of the BPL spectrum is a combination of the high-energy bump of the SYN component from the HE region and the low-energy side of the SSC component. This SED fitting can only be accurately reproduced with a delicate balance of model parameters. In other words, it can be said that this imposes stricter constraints on the model parameters. As mentioned in Sect. 4.1, the presence of the BPL spectrum requires the bump to be located in the MeV energy band. However, we face the lack of data at MeV energies that could help provide crucial information to distinguish between the proposed models.

Fig. 21. X-ray spectra of M87 from Chandra and NuSTAR observations in April 2018. Colour-coding for the spatially resolved components is the same as in Fig. 7. The red curves represent the spatially resolved model appropriate for each dataset. The NuSTAR model is the sum of the other spatial components (with cross-normalisation constants to account for instrumental differences between NuSTAR and Chandra as well as between NuSTAR’s two focal plane modules).

4.5. The possible origin of γ-ray emission

One motivation to release this data set is to aid in the development of more complex models that can address the origin of the γ-ray emission, at the same time satisfying the new constraints from VLBI imaging. We discuss here potential scenarios for the generation of HE non-thermal electrons (and positrons) in a compact region with the size of the order of ∼10r_g or a larger region ∼100r_g.

The potential production and acceleration of high-energy electron-positron pairs at the vacuum gap formed at the stagnation surface within an accreting BH magnetosphere, could offer an explanation for the observed gamma-ray flare emission (e.g., Levinson & Rieger 2011; Broderick & Tchekhovskoy 2015; Hirotani 2018; Chen & Yuan 2020; Crinquand et al. 2020; Katsoulakos & Rieger 2020; Kisaka et al. 2022, and references therein). Semi-analytical calculations based on this class of models indicate that the potential VHE emission power from this gap region is sufficient to explain the day-scale VHE variability from M87. While this gap-driven black hole magnetosphere model can explain TeV γ rays, it seems difficult to reproduce GeV γ-ray flux due to a hard photon spectrum and a very high maximum energy of electrons (e.g., Levinson & Rieger 2011; Hirotani 2018). Therefore, other radiation mechanisms seem necessary to reproduce the observed GeV γ-ray emission.

Magnetic reconnection, a commonly occurring process in magnetised plasma, is capable of accelerating particles to high energies (e.g. Zenitani & Hoshino 2001; Lyubarsky & Liverts 2008; Sironi & Spitkovsky 2014; Guo et al. 2014; Werner et al. 2016). Comparison of the linearly polarised, horizon-scale emission of M87 detected by EHT with GRMHD simulations strongly suggests that M87*’s accretion flow is currently in a magnetically arrested disk (MAD) state (Event Horizon Telescope Collaboration 2021a,b). The MAD state (Igumenshchev & Narayan 2002; Narayan et al. 2003, e.g.) occurs when magnetic fields carried in by the accretion flow build up on the horizon, to the point where the magnetic pressure can push back, and truncate, the accretion flow.

The preference for a MAD-like state is also supported by other recent investigations on the properties of magnetic fields and velocity fields downstream of the M87 jet (Chael et al. 2019; Cruz-Osorio et al. 2022; Yuan et al. 2022; Kino et al. 2022; Fromm et al. 2022). Recent ultra-high resolution GRMHD simulations of MAD-state flows have shown that magnetic reconnection can occur in an equatorial current sheet just outside the event horizon as a result of magnetic flux eruptions; this process has been proposed as a possible origin of VHE γ-rays (e.g. Ripperda et al. 2020, 2022; Chashkina et al. 2021; Hakobyan et al. 2023), presenting models that are qualitatively consistent with the M87 MWL campaign observation in 2018. Furthermore, the size of the transient current layer can become of the order of 10 r_g (Ripperda et al. 2022), that is qualitatively consistent with the HE-oriented model of M87 2018 data. Recent studies of Gelles et al. (2022) and Jia et al. (2023) showed that during a reconnection-powered VHE flare, the sub-millimetre/radio core flux can decline due to the temporarily ejected accretion disk, that may be similar to what we observe (see light curves in Figs. 8 and 4).

The flaring VHE γ-rays could also originate via an additional fast moving component (i.e. a structured, multi-zone scenario). There are mainly two types of models for representing the fast moving component: (i) the structured jet with fast spine and slow sheath components⁶ (e.g. Ghisellini et al. 2005) and (ii) a fast moving blob (or mini jet) component (e.g. Lenain et al. 2008; Giannios et al. 2009). The latter is considered in Model B (Sect. 4.3) to account for the VHE γ rays. There are two possible origins of this fast moving component: magnetic reconnection or shock.

Whether magnetic reconnection can account for the implied high U_e/U_B (Table 1) remains to be explored. Current (2d PIC) simulations of relativistic reconnection (e.g., Sironi et al. 2015; Petropoulou et al. 2016) rather suggest U_e ∼ U_B. The shock formed via the interaction between strongly magnetised jets and weakly magnetised winds is another possible origin for the blob component of the Model B. Kusafuka et al. (2023) performed one-dimensional spherically symmetric special-relativistic MHD simulations and found that ∼10% of the magnetic energy can be dissipated. This resulted in the acceleration of a weakly magnetised wind, that may be the emission site of the VHE γ-rays. Since the very high U_e/U_B ∼ 10^4 − 6 is also reported in some blazars (e.g., 1ES 1101-232 and 3C 66A), with simple analysis using single-zone SSC models (Domínguez et al. 2013), it will be important to further explore to what extent a high U_e/U_B can be realised under different conditions.

Although we primarily focus on the leptonic model, hadronic scenarios might also play a non-negligible role in explaining the M87 SED (Reimer et al. 2004; Reynoso et al. 2011; Rieger & Aharonian 2012; Kimura & Toma 2020; Boughelilba et al. 2022; H.E.S.S. Collaboration 2023).

4.6. Observational aid towards unveiling gamma ray origin in M87

While Model A agrees with magnetic reconnection on the point of the variation timescale if the reconnection takes place in the vicinity of the black hole in the MAD (e.g., see Sect. 4.1 in Ripperda et al. 2022), Model B may require shock acceleration because of a large U_e/U_B. X-ray polarisation information can help discriminate these two model scenarios by constraining the non-thermal electron production mechanism. If the flare is caused by magnetic reconnection, the position angle may drastically flip with time, whereas a shock-induced flare may not produce this strong effect. Observations with IXPE or other future X-ray polarimetry missions (such as eXTP; Zhang et al. 2019) could provide important information for unveiling the emission mechanisms in M87.

In addition, observations carried with more sensitive VHE instruments like the upcoming Cherenkov Telescope Array (Cherenkov Telescope Array Consortium 2019; Zanin et al. 2022) are of fundamental interest. Thanks to its better short-timescale sensitivity and spectral resolution, CTA may provide new insights on future observed γ-ray flaring episodes.

Moreover, the sparse sampling and limited data qualities in the currently available EHT observations make it difficult to draw a conclusion on the origin of the long-term PA variation at μas scale. Further EHT images in the subsequent years will be needed to constrain the actual duration of PA variations at EHT scales and their physical origin.

More extended MWL monitoring with higher observation frequency is essential for capturing future MWL flares beyond the VHE one detected from M87 in 2018. In future campaigns, we aim to extend the monitoring period, increasing the chances of identifying potential correlations between the MWL light curves and improving sensitivity to detect up to weeks-long variability. Finally, an enhanced EHT array, as well as potentially future next-generation experiments (e.g., Lico et al. 2023), will allow us to increase the observed field-of-view. This will improve our ability to probe regions near the core (e.g., the inner jet) as potential sites for γ-ray flares.

5. Summary

We have presented the results of the second MWL observational campaign on M87 together with M87* observations performed by the EHT in April 2018. In addition to obtaining the second most extensive quasi-simultaneous broadband spectra of M87 taken yet, we have provided details of the individual observations and light curves, with the aim of investigating the emission variability of M87 and constraining the origin of the γ-ray emission in M87. Similar to M87 MWL2017, the primary outcome of this campaign is represented by the incredible legacy data set, for which all data and some analysis scripts are available to the community via a Cyverse repository (see Data availability section). While we defer detailed modelling to future works by the broader community, we draw preliminary conclusions using a heuristic approach to the SED modelling. As seen in 2017, the HST-1 knot component was subdominant, and the core was dominating the total flux in the radio through the X-ray bands, providing ideal conditions for combining MWL data over a large range of spatial resolutions.

The main results of the 2018 EHT-MWL observational campaign on M87 can be summarised as follows:

• We detected a short (approximately three days) VHE γ-ray flare, the first since 2010. The combined data from all IACTs involved in the campaign, H.E.S.S., MAGIC, and VERITAS, helped fully describe the flare duration and the shape of the light curve and indicated spectral hardening during the flare.

• An HE γ-ray flux increase was also detected by Fermi-LAT during the same days.

• A likely longer-term core flux enhancement was observed in the X-ray band by Chandra, presenting an X-ray spectrum harder than in non-flaring states, as was observed during the EHT 2017 campaign.

• While radio and millimetre core fluxes are compatible with (or even lower than) the emission seen in April 2017, VLBA (centimetre to millimetre) radio observations present a clear change, on an annual basis, of the jet position angle within a few milli-arcsecond from the core. The EHT (millimetre) radio observations also revealed a significant variation in the position angle of the asymmetry of the ring, indicating that this might be related to changes seen at larger scales in the jet’s position angle.

• Although the presence of the flaring episodes allowed us to constrain the size of the VHE γ-ray emitting region in the SED modelling, the location is still uncertain. M87*’s complex, broadband SED cannot be modelled by a one-zone model. Instead, an additional second component (or even a third one as obtained for Model B) is needed to capture its emission.

The results of this broadband MWL campaign of M87 offer a unique opportunity to investigate the origin of the γ-ray emission during a flaring episode. In particular, they highlight the importance of having coordinated MWL observations with well-sampled coverage to fully characterise the spectral variability in M87*, that likely spans several different timescales. The data sets made public as a part of this effort also offer a rich archive for future investigations. As outlined in Sect. 1, these data add to the growing suite of observational constraints that inform GR models as well as phenomenological models, such as those presented in this work. The ∼17% uncertainty on testing GR is currently dominated by theoretical uncertainties, but combining the M87 SED with fits to the EHT image can reduce these errors and improve our understanding of the underlying astrophysics.

The EHT-MWL observational campaigns performed with a more sensitive array in 2021 and 2022 as well as those planned for upcoming years will enable further insights for constraining the physics around M87*, such as the disk-jet connection, and the origin and mechanisms responsible for the emission of γ-ray photons. In addition, the results from the EHT observations of other sources such as Sgr A* (Event Horizon Telescope Collaboration 2022a,b), 3C 279 (Kim et al. 2020), Centaurus A (Janssen et al. 2021), OJ 287, and other targets planned for the upcoming campaigns, together with their contemporaneous MWL coverage, will greatly expand the sample of objects and the impact of the MWL studies undertaken within the EHT-MWL working group. The high resolution EHT millimetre-VLBI images and the broadband SED information will enable deep investigation of the emission mechanisms driving AGN

Data availability

Tables from Appendix C (from C.1 to C.15) are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/692/A140

Data products presented in this paper are available for download through the EHT Collaboration Data Webpage (https://eventhorizontelescope.org/for-astronomers/data), or directly from the CyVerse repository. The repository contains the following data products:

• Broadband spectrum table with frequency, flux density, its uncertainty, and instrument index (format: DAT).

• Sampled posterior distributions of the SED broadband spectral model (A-PL) described in Sect. 4.2 (format: DAT).

• Sampled posterior distributions of the SED broadband spectral model (A-BPL) described in Sect. 4.3 (format: DAT).

• Sampled posterior distributions of the SED broadband spectral model (B) described in Sect. 4.3 (format: DAT).

Acknowledgments

The full acknowledgements are available in Appendix A.

References

Appendix A: Acknowledgements

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

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

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

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

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

G. Principe: project leadership, Fermi-LAT data analysis, MWL variability and SED study, results interpretation, paper writing; J. C. Algaba: MWL images study; C. Arcaro: H.E.S.S. analysis; M. Balokovic: X-ray data analysis and interpretation ; V. Barbosa Martins: H.E.S.S. analysis, results interpretation, and paper writing; S. Chandra: AstroSAT data analysis; Y.-Z. Cui: VERA, EAVN/KaVA data analysis, and radio properties study; F. D’Ammando: Swift-UVOT; A. D. Falcone: Swift-XRT analysis; N. M. Ford: X-ray results interpretation; D. Glawion: H.E.S.S. analysis and paper writing; M. Giroletti: Fermi-LAT data crosscheck and paper writing; C. Goddi: ALMA data analysis; M. A. Gurwell: SMA data analysis; K. Hada: VERA, EAVN/KaVA data analysis, radio results interpretation, and paper writing; D. Haggard: X-ray analysis, MWL results interpretation and paper writing; A. Hahn: MAGIC data analysis, MWL variability study, paper writing, and IACT cross-normalisation analysis; W. Jin: VERITAS data analysis, IACT cross-normalisation analysis, and paper writing; S. Jorstad: Swift-UVOT data analysis; A. Kaur: Swift-XRT data analysis; T. Kawashima: SED Modelling and interpretation; S. Kerby: Swift-XRT data analysis; J.-Y. Kim: KVN data analysis; M. Kino: SED modelling and interpretation; E. V. Kravchenko: VLBA data analysis; S.-S. Lee: KVN data analysis; R.-S. Lu: GMVA data analysis; S. Markoff: MWL results interpretation and paper writing ; D. Mazin: MAGIC data analysis; J. Michail: paper writing; J. Neilsen: Chandra and NuSTAR data analysis and X-ray results interpretation; M. A. Nowak: Swift-XRT analysis; G. Pühlhofer: results interpretation and IACT cross-normalisation analysis; V. Ramakrishnan: paper writing; B. Ripperda: MWL results interpretation and paper writing; M. Santander: VERITAS data analysis; M. Sasada: SED modelling; S. S. Savchenko: Swift-UVOT analysis crosscheck; C. Sheridan: paper writing; all remaining authors: contributed with technical tasks, development, maintenance, and operations of the instruments as well as with reviewing the draft.

This research has made use of data obtained with 12 radio telescopes from the East Asian VLBI Network (EAVN): 2 stations from the Chinese VLBI Network (the Tianma radio telescope operated by Shanghai Astronomical Observatory of Chinese Academy of Sciences (CAS) and the Nanshan radio telescope operated by Xinjiang Astronomical Observatory of CAS), 2 stations from the Japanese VLBI Network (the Hitachi station operated by Ibaraki University and the Kashima station operated by the National Institute of Information and Communications Technology), the Korean VLBI Network operated by the Korea Astronomy and Space Science Institute (KASI), the VLBI Exploration of Radio Astrometry (VERA) operated by the National Astronomical Observatory of Japan (NAOJ), and the Nobeyama 45-meter radio telescope operated by NAOJ. We thank VERA staff members who helped the operation. We are grateful to KVN staff who helped operate the array and to correlate the data for quasi-simultaneous multi-wavelength observations of M87. The KVN is a facility operated by KASI. The KVN operations are supported by KREONET (Korea Research Environment Open NETwork), which is managed and operated by KISTI (Korea Institute of Science and Technology Information)

The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc.

The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Etudes Spatiales in France. This work is performed in part under DOE Contract DE-AC02-76SF00515.

The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l’énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Knut and Alice Wallenberg Foundation, the National Science Centre, Poland grant no. 2016/22/M/ST9/00382, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation.

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22, PID2022-136828NB-C41, PID2022-137810NB-C22, PID2022-138172NB-C41, PID2022-138172NB-C42, PID2022-138172NB-C43, PID2022-139117NB-C41, PID2022-139117NB-C42, PID2022-139117NB-C43, PID2022-139117NB-C44 funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making Europe”; the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-400/18.12.2020 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also been supported by Centros de Excelencia “Severo Ochoa” y Unidades “María de Maeztu” program of the Spanish MCIN/AEI/ 10.13039/501100011033 (CEX2019-000920-S, CEX2019-000918-M, CEX2021-001131-S) and by the CERCA institution and grants 2021SGR00426 and 2021SGR00773 of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2022-10-4595 and the University of Rijeka Project uniri-prirod-18-48; by the Deutsche Forschungsgemeinschaft (SFB1491) and by the Lamarr-Institute for Machine Learning and Artificial Intelligence; by the Polish Ministry Of Education and Science grant No. 2021/WK/08; and by the Brazilian MCTIC, CNPq and FAPERJ.

The Medicina and Noto radio telescopes are funded by the Ministry of University and Research (MUR) and are operated as National Facility by the National Institute for Astrophysics (INAF).

VERITAS is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy’s Office of Science, and resources of the National Energy Research 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.

S.M. is thankful for support from an NWO (Netherlands Organisation for Scientific Research) VICI award, grant Nr. 639.043.513.

D.H. acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant and the Canada Research Chairs program.

G.P. acknowledges support by ICSC - Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union - NextGenerationEU.

J.-C.A. acknowledges support from the Malaysian Fundamental Research Grant Scheme (FRGS) FRGS/1/2019/STG02/UM/02/6.

M.B. acknowledges support from the YCAA Prize Postdoctoral Fellowship and the Black Hole Initiative at Harvard University, which is funded in part by the Gordon and Betty Moore Foundation (grant GBMF8273) and in part by the John Templeton Foundation.

K.H. acknowledges support from JSPS KAKENHI grant Nos. JJP18H03721, JP19H01943, and JP18KK0090.

T.K. acknowledges support from JSPS KAKENHI grant Nos. JP18K13594, JP19H01908, JP19H01906, JP23K03448, JP23H00117, MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Structure and Evolution of680 the Universe Unraveled by Fusion of Simulation and AI;681 Grant Number JPMXP1020230406; Project ID hp230204 , hp240219), by the HPCI System Research Project (Project ID: hp230116, hp240054), and JICFuS. A part of the calculations were carried out on the XC50 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

R.-S.L. is supported by the National Science Fund for Distinguished Young Scholars of China (Grant No. 12325302), the Key Program of the National Natural Science Foundation of China (NSFC, grant No. 11933007), the Key Research Program of Frontier Sciences, CAS (grant No. ZDBS-LY-SLH011), the Shanghai Pilot Program for Basic Research, CAS, Shanghai Branch (JCYJ-SHFY-2021-013) and the Max Planck Partner Group of the MPG and the CAS.

J.N. acknowledges support from SAO award DD7-18089X and NASA award 80NSSC20K0645.

J.P. acknowledges financial support from the Korean National Research Foundation (NRF) via Global PhD Fellowship grant 2014H1A2A1018695 and support through the EACOA Fellowship awarded by the East Asia Core Observatories Association, which consists of the Academia Sinica Institute of Astronomy and Astrophysics, the National Astronomical Observatory of Japan, Center for Astronomical Mega-Science, Chinese Academy of Sciences, and the Korea Astronomy and Space Science Institute.

Y.-Z.C. acknowledges support from the Natural Science Foundation of China (grant 12303021) and the China Postdoctoral Science Foundation (No. 2024T170845).

J.Y.K. is supported for this research by the National Research Foundation of Korea (NRF) grant funded by the Korean government (Ministry of Science and ICT; grant no. 2022R1C1C1005255).

J.M. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2401752.

Appendix B: Detailed analysis descriptions

B.1. Radio observations

B.1.1. VERA 22 GHz

All the data were analysed in with standard VERA data reduction procedures (see Nagai et al. 2013; Hada et al. 2014, for more details). We note that VERA recovers only part of the extended jet emission due to the lack of short baselines, so the total VERA fluxes listed in Table C.1 significantly underestimate the actual total jet fluxes.

B.1.2. EAVN/KaVA 22 and 43 GHz

Especially between March and May of 2017 and 2018 when EHT M87 observations were performed, the observing cadence of EAVN monitoring was intensified. The default array configurations were KVN+VERA+Tianma+Nanshan+Ibaraki at 22 GHz and KVN+VERA+Tianma at 43 GHz respectively, while occasionally a few more stations in East Asia (Sejong, Kashima, Nobeyama) and Italy (Medicina, Sardinia) joined if available. During the period outside March–May of each year, the observations were mostly conducted with the KVN and VERA Array (KaVA; Niinuma et al. 2014; Hada et al. 2017; Park et al. 2019), a core array of EAVN.

Each of the KaVA/EAVN sessions was made in a 5–7-hour continuous run at a data recording rate of 1 Gbps (a total bandwidth of 256 MHz in a single polarisation mode). 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 (Greisen 2003) software package for the initial calibration of visibility amplitude, bandpass and phase calibration. Imaging with the CLEAN algorithm (Högbom 1974) and self-calibration were performed with the Difmap (Shepherd 1997) software.

B.1.3. VLBA 24 and 43 GHz

The total on-source time during the full-track segments is about 1.7 hours at 24 GHz and 6 hours at 43 GHz, while for the short segments it amounts to 30 minutes 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 polarisations at a rate of 2048 Mbps, and correlated with the VLBA software correlator in Socorro (Deller et al. 2011). The initial data reduction was conducted using AIPS, following the standard calibration procedures for VLBI data (Crossley et al. 2012; Walker 2014). 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. Images and polarimetric results from part of these VLBA observations are published in Kravchenko et al. (2020) and Park et al. (2021).

B.1.4. GMVA+ALMA 86 GHz

The participating GMVA stations include the eight VLBA antennas equipped with 3.5 mm receivers, the Effelsberg 100 m telescope, the 100 m Green Bank Telescope, the Metsähovi 14 m telescope, the Onsala 20 m telescope, the IRAM 30 m telescope, and the Yebes 40m telescope (YS). We note that the participation of ALMA and the GLT has improved the north-south resolution by a factor of ∼4 for the observations of M87.

The observations (GMVA project code ML005) were performed for a total of ∼13 h (∼6 h with five telescopes in Europe (Effelsberg, Metstähovi, Onsala, IRAM30m and Yebes), ∼7 h with ALMA, VLBA and GBT, and ∼12 h with the GLT). Two bright calibration sources (3C273 and 3C279) were observed every 20-40 minutes in interleaved scans. The data were correlated with the VLBI correlator at the Max Planck Institute for Radio Astronomy using DiFX (Deller et al. 2011).

During these observations, most stations recorded the signals on a circular polarisation basis. However, the phased ALMA recorded the signals on a linear polarisation basis. Consequently, the visibilities between ALMA and the other stations were correlated on a mixed polarisation basis. After correlation, these visibilities were converted into a pure circular basis using the internal calibration of the ALMA interferometric data with the PolConvert algorithm (Martí-Vidal et al. 2016a). In addition, GLT participated in these observations on a best-effort basis during its commissioning phase. The waveguide phase-shifter of the GLT that was used to convert the on-sky circular-polarisation signals into linear polarisation for detection at the backend was erroneously configured to apply a rotation of 45° (instead of 90°) between the polarisation channels. The instrumental polarisation effects caused by this rotation were removed with a dedicated algorithm after correlation. After this instrumental polarisation calibration, the data were then further calibrated following typical VLBI data reduction procedures in AIPS. Hybrid imaging was performed in the Difmap software with CLEAN and self-calibration. We refer the reader to Lu et al. (2023) for further details of data calibration. Based on the calibration and imaging of the calibrator 3C273, we consider an error budget of 25 % for the flux density estimate of M87. While a dedicated study based on these data is presented in Lu et al. (2023), here we provide a peak flux density of the core and VLBI-scale total flux density (Tables C.1 and C.2).

B.1.5. KVN 22, 43, 86, and 129 GHz

The data were calibrated using a modified version of the KVN pipeline (Hodgson et al. 2016) to include ionospheric delay corrections with total electron content (TEC) maps, parallactic angle corrections, and Earth orientation parameter corrections. Prior to global fringe fitting, the 43, 86, and 129 GHz data were calibrated with the 22 GHz data by means of frequency phase transfer. After the pipeline calibration with AIPS, the inner 75% of each subband (due to flux loss at the bandpass edges) was used for further analysis. Following Kim et al. (2018), the data were averaged by 30 seconds for the 22 GHz and 43 GHz data, and by 10 seconds for the 86 GHz and 129 GHz data. As in M87 MWL2017, we assume a flux density measurement uncertainty of 10% at 43 GHz and 86 GHz, and 30% at 129 GHz. During the observations, a phase instability issue occurred with the 22 GHz receiver at the KVN Yonsei site (KYS) between January and March 2018, and between November 2018 and February 2019. This resulted in a ∼50% flux loss for all KYS baselines. From a sample of compact sources within the iMOGABA program, we derived and applied empirical correction factors to the impacted data. An uncertainty of 15% is assumed for the 22 GHz flux density measurements made during observations where the empirical corrections were applied. For observations that were not impacted by the flux loss, an uncertainty of 10% is adopted. A common 1 mas FWHM circular Gaussian beam is used to convolve the 22, 43, 86, 129 GHz CLEAN maps when determining the peak flux density of the core. This is closest to the natural beams of the KVN at 86 and 129 GHz. In Table C.1, we summarise the data obtained during 2018. In Fig.3, we also add the earlier (2012–2016) data that were originally published in the literature (Kim et al. 2018) but the peak flux densities are recalculated using the procedures applied for 2017–2018 data.

B.1.6. Global 43 GHz VLBI

The participating stations include Effelsberg, Yebes 40 m, Metsahovi, KVN (three stations), 100-m Green Bank Telescope, VLBA (10 stations), and the phased Jansky Very Large Array (JVLA). The data were recorded with a total recording rate of 2 Gbps, corresponding to 256 MHz bandwidth per polarisation. The scheduling, observation, and data reduction followed the standard procedures using AIPS. Similar to other VLBI data reduction and imaging, the CLEAN imaging and self-calibration were made using procedures implemented in Difmap. The resulting FWHM angular resolutions are 0.2 × 0.7 (0.12 × 0.47) mas² with the beam position angle ∼ 0 deg for the natural (uniform) weighting.

B.1.7. ALMA 93 GHz and 221 GHz

The VLBI observations were carried out while the array was in its most compact configuration and only antennas within a radius of 180-m (from the array centre) were used for phasing. The observations were performed in full-polarisation mode to supply the inputs to the polarisation conversion process at the VLBI correlators (Martí-Vidal et al. 2016b).

The setup included four spectral windows (SPWs) of 1875 MHz in each band, two in the lower and two in the upper sideband, centred at 86, 88, 98, and 100 GHz in Band 3, and 213, 215, 227, and 229 GHz in Band 6, respectively. The SPWs were correlated with 240 channels per SPW (corresponding to a spectral resolution of 7.8125 MHz). Details about ALMA operations in VLBI mode 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 package (CASA, version 6.4.1 used; McMullin et al. 2007; CASA Team 2022) task tclean. Only phased antennas were used to produce the final images (with baselines < 360 m), yielding synthesised beam sizes ∼2.5″ in Band 3 and ∼ 1″ in Band 6, respectively. To isolate the core emission from the jet, we compute the sum of the central nine pixels of the model (cleaning component) map, an area of 3 × 3 pixels; the contribution from the jet is accounted for by also summing the clean components along the jet. While the extended emission accounts for less than 20% of the total emission at 1.3 mm, at 3mm it also includes emission from the radio lobes, accounting for more than 60% of the total emission. The absolute amplitude scale calibration has a 5%/10% uncertainty in Band 3/6, respectively. Details about the imaging and flux extraction methods can be found in Goddi et al. (2021).

B.1.8. SMA 200 & 345 GHz

Data from this program are updated regularly and are available at the SMA Observer Center website (SMAOC⁷). Data were primarily obtained in a compact configuration (baselines from 10 to 75 m) though some were obtained at longer baselines up to 210 m. The effective spatial resolution was therefore generally around 3″. We have limited analysis to projected baselines longer than 30 kλ (effective spatial scales smaller than 6.8″). The flux density was obtained by vector averaging of the calibrated visibility data from each observation.

Interferometric observations of M87 were additionally conducted as part of the 2017 and 2018 EHT campaigns, concurrent with the SMA running in phased-array mode operating at 230 GHz. In 2017 the SMA was in compact configuration, while in 2018 the SMA was in extended configuration (baselines 25 m to 210 m). In both years, the interferometer operated in dual-polarisation 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 passband and amplitude calibrated using 3C 279, with flux calibration performed using either Callisto, Ganymede, or Titan. Phase calibration was done through self-calibration of the M87 data itself, assuming a point-source model.

B.2. Optical and UV observations

B.2.1. Optical observations by the Kanata telescope

We adopted standard reduction procedures for the images: bias- or dark-frame subtraction, correction of bad pixels, flattening correction, removal of contamination of CR signals, sky-background subtraction, solving the corresponding coordinates in the World Coordinate System (WCS) using field stars, and stacking of the images aligned with the coordinates of WCS.

We performed an aperture photometry of field stars in the stacked images using Source Extractor software (Bertin & Arnouts 1996). The value of FLUX_AUTO was adopted as the photometric measurements. We estimated zero points of the images in the AB magnitude system using nearby stars whose actual magnitudes have been provided in the references: the Pan-STARRS1 catalogue (Flewelling et al. 2020) of the optical band after conversion of the filter system (Tonry et al. 2012) and the 2MASS catalogue (Skrutskie et al. 2006) of the near-infrared band.

B.2.2. Optical and UV observations from Swift-UVOT

UVOT data in all filters were analysed with the uvotimsum and uvotsource tasks included in the HEASoft package (v6.32) and the 20231208 CALDB-UVOTA release. Source counts were extracted from a circular region of 5 ″radius centred on the source, while background counts were derived from a circular region with 30 ″radius in a source-free region located outside the host galaxy radius (∼ 120 ″). The UVOT magnitudes are corrected for the host galaxy contamination (see Table C.6 ) and the Galactic extinction using an E(B–V) value of 0.020 from Schlafly & Finkbeiner (2011) and the extinction laws from Cardelli et al. (1989) and converted to flux densities using the conversion factors from Breeveld et al. (2010). We have used the same host galaxy model as presented in Table A2 of EHT MWL Science Working Group (2021). However, we have to note that flux densities of the host galaxy listed in Table A2 of EHT MWL Science Working Group (2021) did not have aspect corrections although the correct values of the host galaxy as given in Table C.6 were used in the data reduction of Swift-UVOT data in EHT MWL Science Working Group (2021).

B.2.3. Far-UV observations from AstroSat-UVIT

For the AstroSat-UVIT analysis, the VIS data are used for drift correction of the images in NUV and FUV channels. The FUV and NUV detectors work in photon-counting mode and the VIS detector works in integration mode. Each of these channels is equipped with a number of narrow and wide band filters (for details refer to Tandon et al. 2017, 2020). Out of the two AstroSat pointing observations (i.e. A04_115T02_9000001900 and A04_115T02_9000002042), only the UVIT data of the second one are available. These observations are made in the BaF2 band (F154W: λ_eff ∼ 1541 Å, with Δλ_eff ∼ 380 Å). The level-1 data, downloaded from the AstroSat data archive ‘astrobrowse’, is used to create the image using the CCDLAB software package (Postma & Leahy 2017) incorporating the corrections for the spacecraft drift, geometrical distortions and flat-fielding.

A two-dimensional (2D) image decomposition method was used with IMFIT⁸ (Erwin 2015) tool. The following combination of 2D models is used to create the model image of M87 using UVIT: PointSource + PointSource + Core-Sersic + SersicJet + Flat-Sky. The point sources used here are convolved with the PSF of the UVIT instrument.

The UVIT-BaF2 image, centred at M87 and zooming into a small region around it is shown in Fig. B.1. The two point source models are added to account for the contributions from M87*’s core and the HST-1 knot. The Core-Sersic model (Graham et al. 2004) is used to represent the emission from the host galaxy. In order to take the asymmetric jet into account we created a modified version of the Sersic model by adding an asymmetry term to the originally symmetric Sersic function. But even with the added asymmetry, the shape of the bright spots towards the end of the jet is too complex to develop a good analytical model to describe it. Motivated by the above we masked this irregular complex part of the jet with a rectangular mask. The Flat-Sky model (Mondal et al. 2018) was used to include a component for the sky background. The best fit composite model corresponds to a reduced χ² value of 1.

Fig. B.1. AstroSat UVIT observations of M87. Left: M87 in BaF2 band. Right: Difference image of M87. The patched area is basically the masked region of the jet.

This model is then used to represent the composite radial profile for the host and jet. We then performed aperture photometry of the image using the Python package photutils⁹ with a source region of 5″ centred at the position of M87*. The contributions of other prominent components to this region were estimated using the modelled radial profile for composite emission from the jet, the host galaxy and the sky. The corrected magnitude is then converted to flux using an appropriate flux conversion factor as provided by the instrument team. The corresponding numbers are shown in the Table C.8.

B.3. X-ray observations

B.3.1. Chandra and NuSTAR observations and joint spectral analysis

We keep the data processing and analysis as consistent as possible with our previous work on the Chandra and NuSTAR data taken in parallel with the 2017 EHT campaign described in M87 MWL2017. We only briefly summarise them here. For the Chandra data, we used standard data reduction procedures in CIAO (Fruscione et al. 2006) version 4.9 to extract 0.4–8 keV spectra from the spatially resolved components, the core, the HST-1 knot, and the outer jet, along with the corresponding instrumental response files. The NuSTAR observation was processed following standard procedures outlined in Perri et al. (2017), using NuSTARDAS version 1.8.0 to produce 3–79 keV spectra and instrumental responses with CALDB version 20180312. Source and background regions for Focal Plane Module A and B (FPMA, FPMB) were the same as in M87 MWL2017: 45″ circles that are large enough to include all of the Chandra spatial components mentioned above as well as much of the surrounding diffuse emission from the intracluster medium (ICM).

We proceed to joint spectral modelling of the Chandra and NuSTAR data in order to take advantage of the high spatial resolution of Chandra and the high-energy sensitivity of NuSTAR. Our modelling strategy, described in detail in M87 MWL2017, properly accounts for the effects of pileup in Chandra data and includes consideration of the cross-normalisation between the instruments with respect to the modelled diffuse ICM emission that is differently sampled by Chandra and NuSTAR. For the diffuse component, we adopt a spectral model with two APEC components (Smith et al. 2001) with different temperatures and variable elemental abundances, three additional emission lines, and a cutoff component for the low-mass X-ray binary contribution from within the host galaxy. For the core, the HST-1 knot, and the remainder of the jet, we modelled their continuum emission using PLs. All spectral models include interstellar absorption (Wilms et al. 2000). The absorption column towards the core was previously found to be higher (e.g. Perlman & Wilson 2005), that we allow for by including an extra layer of absorption in the model for this component.

We performed spectral fitting on spectra binned using a custom algorithm that ensures a minimum signal-to-noise ratio within each energy bin (as described in M87 MWL2017). To find the best-fit model parameters we minimised the Cash statistic (Cash 1979) as appropriate for the low-count regime and our binning strategy. Our best fit has a Cash statistic of 1219.023 for 1034 data bins and 37 free parameters. Our multi-component best-fit model is shown along with the data in Fig. 21; red lines represent the model for each dataset, and the NuSTAR model is the sum of all the other spatial components. We derive uncertainties using a Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013) set to run for 60,000 steps with 10 walkers for each of 37 parameters (a total of 22.2 million model evaluations). We discarded the first 185 000 evaluations for calculating chain statistics. Unless otherwise noted, quoted uncertainties represent minimum width intervals containing 90% of the probability function for each spectral parameter.

B.3.2. Swift-XRT observations and analysis

The standard pipeline xrtpipeline version 0.13.5 was used to create level 3 cleaned event files for each Swift-XRT observation. We use a circular source region of radius 35″ for spectral extraction, encompassing all three of the X-ray sources examined individually in the Chandra analysis. For background extraction, we create an annular background region with inner and outer radii of 130″ and 150″, respectively. Using xselect (version 2.4k), we extract spectra for the source and background regions after correcting for pile-up in the source region. If the source region count rate exceeded 0.5 ct s⁻¹ for any observations, we applied a pile-up correction by using an annular region to exclude central pixels affected by pile-up.

Creating ancillary response files for each spectrum with the xrtmkarf task, we grouped the spectral files for each epoch to a minimum of 20 counts per bin using the grppha command to facilitate χ² fitting using xspec. Since the spectral resolution of Chandra is better than Swift-XRT, we adopted the source and background model fitted in the Chandra analysis. All spectral parameters of that model were frozen at the Chandra best-fit values (see Sect. 2.3.1) except the normalisation of the APEC and line emission backgrounds and the joint normalisation of the three power law components for the core, the HST-1 knot, and jet emission. The relative normalisation between the two APEC background components and between the three PL source components were kept frozen at the Chandra values. We use χ² fitting via xspec to conduct this X-ray analysis. In this way, we evaluate the changing total and small-scale components (core, the HST-1 knot and inner jet) fluxes from the source.

B.3.3. AstroSat-SXT observations and analysis

We used sxtpipeline (version 1.4b) along with the standard filtering criteria to generate the cleaned events file. The standard procedure includes time-tagging of events, removing hot pixels, bias corrections, coordinate transformation, filtering (South Atlantic Anomaly passage, Sun-angle, Moon-angle), charge transfer inefficiency corrections, and generation of full-frame products. The processed event lists obtained in this way were then input to xselect to extract spectra, light curves, and images.

We used a circular source region of radius of 15′ centred on the brightest pixel in the SXT image (α = 12:30:50.41, δ = +12:22:32.4), that is consistent with the peak of the flux distribution in Chandra and NuSTAR images. From the PSF size of the SXT we estimate that this region includes over 97% of the total number of detected photons. The background spectrum extracted using deep blank sky observations and distributed by the instrument team¹⁰ is used to model the spectrum. We used the tool sxtARFModule (version 0.3) along with the full-frame ARF (sxt_pc_excl00_v04.arf) distributed by the SXT Payload Operation Center¹¹ to create ARF for the specific source region selection. We utilised the RMF file sxt_pc_mat_g0to12.rmf corresponding to the default 0–12 grade selection.

For spectral fitting we used a model defined with the following expression in xspec: TBabs × (vvapec + vvapec + cutoffpl + TBabs × powerlaw + powerlaw + powerlaw + Gaussian + vvapec). This is similar to the model adopted for Chandra and NuSTAR spectral analysis, except one Gaussian component, that is added in order to fit the broad hump-like residuals around 1.4 keV. An additional vvapec component is included to account for the diffuse emission from cluster gas in the outer regions because the SXT samples a larger portion of it (radius ≃15′) than other focusing X-ray instruments considered in this work. The parameters of the three vvapec components, except the temperature and normalisation, are fixed to the values obtained from joint Chandra and NuSTAR spectral modelling (Sect. 2.3.1). To reduce the number of free parameters, we tied together the temperature and normalisation parameters of the vvapec components with a multiplicative factor. We adopt the sum of these three components with different temperatures to represent the large-scale emission from hot gas in the galaxy and the surrounding cluster. Three PL components are representative of the core, the HST-1 knot, and jet contributions; their PL photon indices are kept fixed based on values determined from the analysis of Chandra and NuSTAR data. We keep the normalisation of the core contribution free in the fitting, while the normalisation of the HST-1 knot and jet components is tied with the core normalisation using the multiplicative factors also employed in Swift-XRT spectral modelling in Sect. 2.3.2, as well as in M87 MWL2017.

As in Sect. 2.3.1, we used the Cash statistic to find the best-fitting composite model. To account for the systematics introduced by the procedures related to ARF generation (such as vignetting model and PSF-based scaling of the effective area, etc.), we included an additional 1% uncertainty to the observed spectra. Constraints on the free parameters of the model are given in Table C.9. Figure B.2 shows the AstroSat-SXT individual PL components to be interpreted as upper limits due to the angular resolution of the instrument.

B.4. γ-ray observations

B.4.1. Fermi-LAT observations

Fermi-LAT (Atwood et al. 2009) is a γ-ray telescope that detects photons by conversion into electron-positron pairs, and it has an operational energy range from 20 MeV to more than 300 GeV. Since its launch in June 2008, it mostly works in survey mode covering the whole sky in about three hours and with a remarkable overall uptime of 99.8% (Ajello et al. 2021).

Regarding the data selection, we used P8R3 (v3) SOURCE class events (Bruel et al. 2018), in the energy range between 100 MeV and 1 TeV, in a region of interest (ROI) of 15° ×15° centred on the position of M87. The low-energy threshold is motivated by the large uncertainties in the arrival directions of the photons below 100 MeV, leading to possible confusion between point-like sources and the Galactic diffuse component (see Principe et al. (2018) for a different analysis implementation to solve this and other issues at low energies with Fermi-LAT).

The analysis (that consists of the following steps: model optimisation, source localisation, spectrum and emission variability studies) was performed with Fermipy¹² (v1.0.1; Wood et al. 2017), that uses the Fermi Tools, version 2-0-18. We created counts maps using 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. 2009). In order to optimise the analysis, we applied different cuts for the data selections at low energies and selected event types with the best PSF¹³. In particular, for energies below 300 MeV, we excluded events with zenith angle larger than 85°, as well as photons from PSF0 and PSF1 event types. Between 300 MeV and 1 GeV we excluded events with zenith angle larger than 95°, as well as photons from the PSF0 event type. Above 1 GeV we used all events with zenith angles less than 105°. The P8R3_Source_V3 instrument response functions (IRFs) were used. The model used to describe the sky includes all point-like and extended sources located at a distance < 15° from the source position and listed in the 4FGL-DR2 (Abdollahi et al. 2020) as well as the Galactic diffuse and isotropic emission. For the latter contributions, we used the same templates¹⁴ adopted to derive the 4FGL catalogue.

Fig. B.2. AstroSat-SXT observations of M87 with SED estimates from spectral modelling. Dashed-dotted and dotted lines are representative of observations made in February and March 2018, respectively. The grey, light-blue and light-red shaded regions around the lines are 1σ uncertainty regime derived using parameter estimations, each representing core, jet and the HST-1 knot sub-components, respectively. The downward arrows are indicative that given the limitations of the instrumental capabilities, these estimations should be treated as upper limits for further investigations.

We performed the spectral analysis, allowing the diffuse background template normalisation and the spectral parameters of the sources to vary within 3° from our targets. For the sources in a radius between 3°–5° and all variable sources only the normalisation was fit, while we fixed to their 4FGL values the parameters of all the remaining sources within the ROI at larger angular distances from our target. We modelled the spectrum of M87 with a PL function.

Finally, for the analysis on the 12 years of LAT data we extracted a light curve using time bins of one month. Similarly, for analysis on a one-month scale centred on EHT observations, the light curve analysis was repeated using time intervals of two days in order to investigate the emission variability on short time scales. We tried also different time intervals with both fixed (1 day and 1 week) or customised duration (discriminating between before, during and after the VHE γ-ray flare). The fluxes in each interval were obtained by leaving only the normalisation free to vary and freezing the other spectral parameters to the best-fit values obtained from the full range analysis.

B.4.2. Very high energy observations: H.E.S.S., MAGIC, VERITAS

H.E.S.S. observations were performed with a wobble mode with an offset of 0.7° from the position of M87, that allows for a simultaneous background estimation. The analysis and reconstruction of the Cherenkov shower images were done with the ImPACT maximum likelihood-based technique (Parsons & Hinton 2014). The separation between γ-rays and hadrons was performed with a boosted decision tree classification method (Ohm et al. 2009). For the signal estimation, the background was estimated using the “ring-background model” whereas the calculation of the light curve and the spectrum was done using the “reflected-background model” (Berge et al. 2007). A source with 1% of the Crab Nebula’s flux can be detected above 100 GeV in ∼10 hours.

The stereoscopic MAGIC observations were conducted in wobble mode with an offset of 0.4°. Image reconstruction and data analysis were performed with the standard MAGIC analysis framework MARS (MAGIC Analysis and Reconstruction Software; Zanin 2013; Aleksić et al. 2016). To estimate the instrument’s PSF in true energy we used the spatial likelihood analysis package SkyPrism from Vovk et al. (2018). The observations on MJD 58232 were taken with a strong moonlight background present and therefore analysed as described in Ahnen et al. (2017). MAGIC’s integral sensitivity is 1% of the Crab Nebula’s flux in ∼26 hours above an analysis energy threshold of 290 GeV.

Data analysis of the observations in 2019 was performed identically to the aforementioned 2018 observations. The reconstructed flux levels in 2019 are compatible with the previously observed quiescent flux (see, e.g. MAGIC Collaboration 2020; EHT MWL Science Working Group 2021). The MAGIC fluxes published in MAGIC Collaboration (2020) use an analysis threshold of 300 GeV for the light curve. Instead of scaling the VHE fluxes between 2012 and 2015 according to the PL spectrum, we repeated the flux calculations with a threshold of 350 GeV, using the original high-level analysis files of the authors of MAGIC Collaboration (2020).

The stereoscopic VERITAS observations were conducted in wobble mode with an offset of 0.5°. The data were analysed using the VERITAS standard analysis procedure (Daniel 2008). The main software package used in this work is called Eventdisplay. An overview of this analysis package can be found in Maier & Holder (2018). In addition, a cross-check analysis was performed using an independent software package called VEGAS (Cogan 2007). VERITAS is designed to measure γ rays with energies from ∼85 GeV and it is able to detect a source with about 1% of the Crab Nebula within 25 hours (Park 2016).

To explore possible systematic flux scale offsets between VERITAS, MAGIC, and H.E.S.S. and investigate how the uncertain flux scales might affect the results in Sect. 2.4.2, we applied two cross-normalisation factors (f₁ and f₂) to the VERITAS and MAGIC fluxes as extra fitting parameters. We note that statistical errors on individual flux points were not altered. We note that this choice is purely arbitrary and that the cross-normalisation factors could be applied to any two of the three IACTs. For an earlier study of systematic offsets between IACTs see Meyer et al. (2010). In our case, the maximum flux scale offsets between the instruments are set to 15%. These offsets are used to constrain the fitting range of f₁ and f₂. The results of the fit of a constant and a Gaussian function (Eq. (1)) are shown in Table B.1.

Appendix C: Tabulated data

This section contains the list of tabulated data from observations described in Sect. 2, as follows:

1. Summary of radio observations and flux densities at centimetre wavelengths in Table C.1 and millimetre wavelengths in Table C.2;

2. Kanata Telescope flux densities at the optical wavelength of 634.9 nm in Table C.3 and near-infrared wavelength of 1250 nm in Table C.4;

3. Optical and ultraviolet flux densities of Swift-UVOT in Table C.5, while Table C.6 contains of the host galaxy;

4. Photometry at the ultraviolet wavelength of 154 nm from AstroSat-UVIT in Table C.8;

5. Optical HST observations in Table C.7;

6. Parameters of the composite spectral model for AstroSat-SXT data in Table C.9;

7. Summary of Swift-XRT observations and fluxes in the 2–10 keV band in Table C.10;

8. Summary of Chandra core fluxes for the observations performed in 2017 and 2018. in Table C.11, while the Chandra+NuSTAR spectral values for stacked observation, normalised on 2018 emission are reported in Table C.12;

9. Summary of Fermi-LAT fluxes for the customised-bin light curve analysis on April 2018. in Table C.13;

10. Summary of VHE γ-ray observations and fluxes above energies of 350 GeV in Table C.14 and parameters of the spectral models fitted to data from H.E.S.S., MAGIC, and VERITAS in Table C.15;

11. Spectral energy distribution for M87 observed during the 2018 EHT campaign in Table C.1;

12. In addition the spectral models (Model-A-PL, Model-A, BPL, Model-B) used to describe the MWL SED (see Sect. 4) are also reported as supplementary material.

Tables C.1 to C.15 are available at the CDS. Furthermore, these data are provided to the community in machine-readable form via the EHTC Data Webpage15 along with associated documentation and files needed for spectral modelling of the X-ray data.

Appendix D: Multi-wavelength light curve variability

We estimate the degree of variability in different wavebands by calculating the fractional variability F_var from Eq. (10) in Vaughan et al. (2003):

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\sqrt{\frac{S^2-\overline{{\sigma }_{\mathrm{err}}^2}}{{\overline{x}}^2}}\end{array}$$(D.1)

where $S^2=\frac{1}{N-1}{\displaystyle \sum _{i=1}^N}{(x_i-\overline{x})}^2$ is the sample variance where $\overline{x}$ is the arithmetic mean of the flux measurements x_i and ${\sigma }_{\mathrm{err}}^2=\frac{1}{N}{\displaystyle \sum _{i=1}^N}{\sigma }_{\mathrm{err},i}^2$ the mean squared measurement uncertainty. The uncertainties of F_var are calculated following Poutanen et al. (2008) as:

$$\begin{array}{c}\hfill \mathrm{\Delta }{F}_{\mathrm{var}}=\sqrt{F_{\mathrm{var}}^2+\mathit{err}\left({\sigma }_{\mathrm{NXS}}^2\right)}-F_{\mathrm{var}}\end{array}$$(D.2)

with

$$\begin{array}{c}\hfill \mathit{err}\left({\sigma }_{\mathrm{NXS}}^2\right)=\sqrt{{\left(\sqrt{\frac{2}{N}}\frac{\overline{{\sigma }_{\mathrm{err}}^2}}{{\overline{x}}^2}\right)}^2+{\left(\sqrt{\frac{\overline{{\sigma }_{\mathrm{err}}^2}}{N}}\frac{2F_{\mathrm{var}}}{{\overline{x}}^2}\right)}^2}.\end{array}$$(D.3)

We present the calculated values of F_var for the 2018 observational campaign in Table D.1 and also show the same values of F_var during the 2017 campaign (M87 MWL2017). As described in Sect. 2.4.2 the observed VHE flare has a time variability scale 𝒪(1 day). Therefore, following the Nyquist-Shannon sampling theorem (Whittaker 1915), flux points closer than 0.5 days are combined using their weighted average value. To be able to compare the fractional variabilities of different bands it is necessary to only select simultaneous data. Therefore, we can only compare the fractional variability of the VHE γ-ray instruments and Swift-XRT listed in Table D.1.

All Tables

All Figures

	Fig. 1. Instrument coverage summary of the 2018 M87 MWL campaign covering MJD range 58198–58270. The grey band indicates the time interval of the VHE γ-ray flare episode.
In the text

	Fig. 2. Radio light curves of the M87 core in 2017 and 2018 at multiple bands. The upper and lower panels are for connected interferometers and VLBI, respectively. The corresponding beam sizes are indicated in Table C.1. KVN data at 22 and 43 GHz are not shown here since KVN captures the data from the shortest baselines of EAVN.
In the text

Fig. 3. Long-term MWL light curves of M87 from the last two decades taken by the instruments participating in the 2018 observational campaign. Here are listed the observations used in the plot and their reference (in case of already published results): Chandra (0.3–7 keV) (Sun et al. 2018), (2.0–10 keV) (Imazawa et al. 2021); MAGIC, H.E.S.S., VERITAS (2004–2011), Fermi-LAT (2008–2011), VLBA (43 GHz peak) (2006–2011) (Abramowski et al. 2012); H.E.S.S. (before April 2004, 2001–2016, after 2019) (H.E.S.S. Collaboration 2023); VERITAS (2011–2012) (Beilicke & VERITAS Collaboration 2012); VERA (2011–2012), MAGIC (2012–2015) (MAGIC Collaboration 2020); KVN (2012–2016) recalculated from Kim et al. (2018); VLBA, EAVN, VERA (2000–2018) (Cui et al. 2023); All instruments 2017 data (M87 MWL2017). We mark the VHE flare in 2018 with a grey-shaded line in the background. The inset in the third row provides a zoomed-in view of the Chandra measurement of HST-1.

In the text

	Fig. 4. Multi-wavelength light curves of M87 taken during the observational campaign covering MJD range 58211–58245. The time period containing the 2008 VHE flare is shaded in grey.
In the text

	Fig. 5. Hubble Space Telescope image of M87 at 2359.7 Å with the host galaxy subtracted. The core and HST-1 are indicated; the distance between the features is 1$\stackrel{″}{.}$00±0$\stackrel{″}{.}$05, and the pixel size is 0$\stackrel{″}{.}$0248.
In the text

Fig. 6. AstroSat-UVIT photometric cut of the decomposition model along the jet direction. The thick black line shows the observed flux along the cut, the blue line is the central point source, the orange line is the host galaxy, the green line is the jet total flux, and the red line is the total model flux. The shaded area shows the masked out region of the jet that has a complex fine structure that cannot be accurately modelled with a simple analytical function and was therefore excluded from the decomposition to reduce the model complexity.

In the text

	Fig. 7. X-ray light curves of the spatially resolved components of M87 from the Chandra observations in 2018. The core (orange) shows clear variability between the two different days on which observations have been carried out, while the HST-1 knot and the rest of the jet remain fairly steady.
In the text

Fig. 8. Mid-term MWL light curves of M87 taken during the observational campaign covering January 2017 to April 2019. Here are listed the observations used in the plot and their reference (in case of already published results): Chandra (0.3–7 keV) (Sun et al. 2018), (2.0–10 keV) (Imazawa et al. 2021) H.E.S.S. (2019) (H.E.S.S. Collaboration 2023); VLBA, EAVN, VERA (2017–2018) (Cui et al. 2023); All instruments 2017 data (M87 MWL2017). We mark the VHE flare in 2018 with a grey-shaded region in the background.

In the text

	Fig. 9. Total flux (blue) and unresolved summed component fluxes (red) of M87 in the 2–10 keV band observed with Swift-XRT in 2018.
In the text

Fig. 10. Fermi-LAT light curve of M87 on customised time intervals during April 2018. The red point refers to the flux estimate, while grey points indicate upper limits. The 1σ upper limit is reported when TS < 10. The vertical blue line (band) indicates the time of the VHE flaring episode, while the horizontal lines represent the averaged flux over 12 years and its uncertainty F12 yr = (1.74 ± 0.11)×10−8 ph cm−2 s−1.

In the text

Fig. 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 2018. The solid and semi-transparent colours represent the original data and scaled (with cross-normalisation factors) data, respectively. The black curve and shaded region represent the best Gaussian fit and 1σ error to the original flux points. The grey lines indicate the percentage of the Crab Nebula’s flux level from integrating its SED in the same energy range. The SED of the Crab Nebula is obtained by fitting VERITAS’s six years of observations (Meagher & VERITAS Collaboration 2015). For flux points with significance less than 2σ, upper limits are given in Table C.14.

In the text

Fig. 12. Very high-energy SEDs measured with H.E.S.S., MAGIC, and VERITAS during the coordinated MWL campaign in 2018. Error bars are equivalent to the 1σ confidence level. An upper limit is drawn where the statistical significance is less than 2σ (i.e. 95% confidence level). The thick solid lines represent the best spectral fitting results by using a simple PL function. The shaded regions represent the 1σ error region of the fitting results. Historical spectral fits are taken from H.E.S.S. Collaboration (2006a,b), Aleksić et al. (2012), Acciari et al. (2008, 2010), Albert et al. (2008), Aliu et al. (2012), Beilicke & VERITAS Collaboration (2012), MAGIC Collaboration (2020) and the thin black line indicates the SED obtained with the MAGIC measurements presented in M87 MWL2017.

In the text

	Fig. 13. Composite of the M87 MWL images at various scales obtained in radio and X-rays during the 2018 campaign. The instrument, observing wavelength, and scale are shown at the top-left side of each image. We note that the colour-scale has been chosen to highlight the observed features for each scale and should not be used for noise levels, dynamic range, or flux density calculation purposes.
In the text

Fig. 14. M87 jet structure comparison based on yearly stacked EAVN images at 43 GHz in 2017 and 2018 (see also Cui et al. 2023). The restoring beam size is 0.5 mas circular Gaussian, as indicated in the bottom-left corner. The image has been rotated by −18°. Contour levels are scaled as (1, 2, 4...)×C1st, where C1st is 1.8 mJy beam−1 in 2017 and 1.3 mJy beam−1 in 2018. For M87* mass MBH = 6.5 × 109 M⊙ (Event Horizon Telescope Collaboration 2019a), 1 mas ≈ 125 rs ≈ 0.08 pc, where rs is the Schwarzschild radius.

In the text

Fig. 15. Spectral index maps obtained between 24 and 43 GHz from VLBA 5 May 2017 (a) and 28 April 2018 (b) observations. The restoring beam size is 0.55 mas ×0.55 mas, that is the circular equivalent-area synthesised beam of VLBA at 24 GHz, drawn as a grey circle in the bottom-right corner. The image has been rotated by −18°. The contours denote the total intensity at 43 GHz and start at 3.5σ rms, increasing in steps of two. The (u, v)-ranges match between all images.

In the text

Fig. 16. Photon PL index Γ compared to the VHE flux normalisation at 1 TeV. Data from H.E.S.S. Collaboration (2006a), Acciari et al. (2008), Albert et al. (2008), Acciari et al. (2010), Aliu et al. (2012), Aleksić et al. (2012), Beilicke & VERITAS Collaboration (2012), MAGIC Collaboration (2020), EHT MWL Science Working Group (2021). Flux normalisations at different energies have been recalculated at 1 TeV using the respective power-law index. No dependence of the spectral index on the flux was observed.

In the text

Fig. 17. Observed broadband SED of M87 contemporaneous with the EHT campaign in April 2018, with the exception of GMVA (carried on February 2018) and HST (carried on July 2018) observations (see Table C.1). The fluxes were measured by various instruments, and they are highlighted with different colours and markers. Filled markers represent flux point estimates, while empty markers indicates ULs. In the X-ray band, the two sets of points represent the core flux from the models with a power law (dark purple) and broken power law (light purple; see Table C.12). Grey points represent the observed broadband SED of the 2017 EHT campaign (EHT MWL Science Working Group 2021). MAGIC observations were performed before the γ-ray flare, VERITAS observations only partially overlap the flare, while all H.E.S.S. observations have been taken during the VHE γ-ray flaring episode (see Table C.1). For the radio/millimetre VLBI and interferometric data, the upper limits on emission size (radius) for several representative frequencies, that are obtained by following the procedures described in M87 MWL2017, are labelled to clarify the constraints used in Sect. 4.

In the text

	Fig. 18. Modified Model A fitted describing the BPL emission seen in the Chandra–NuSTAR X-ray spectrum.
In the text

	Fig. 19. Model A fitted to the observed contemporaneous broadband SED of M87 from the EHT campaign in April 2018. For each considered region (EHT and HE), the synchrotron (labeled as SYN) and synchrotron self-Compton (labeled as SSC) models as well as the summation of the two models (labeled as Sum) are reported.
In the text

	Fig. 20. Model B fitted to the observed contemporaneous broadband SED of M87 for the EHT campaign in April 2018. For each considered region (EHT, high-energy, and very-high-energy flare), the synchrotron (labeled as SYN) and synchrotron self-Compton (labeled as SSC) models as well as the summation of the SYN and SSC models (solid line) are reported.
In the text

Fig. 21. X-ray spectra of M87 from Chandra and NuSTAR observations in April 2018. Colour-coding for the spatially resolved components is the same as in Fig. 7. The red curves represent the spatially resolved model appropriate for each dataset. The NuSTAR model is the sum of the other spatial components (with cross-normalisation constants to account for instrumental differences between NuSTAR and Chandra as well as between NuSTAR’s two focal plane modules).

In the text

	Fig. B.1. AstroSat UVIT observations of M87. Left: M87 in BaF2 band. Right: Difference image of M87. The patched area is basically the masked region of the jet.
In the text

Fig. B.2. AstroSat-SXT observations of M87 with SED estimates from spectral modelling. Dashed-dotted and dotted lines are representative of observations made in February and March 2018, respectively. The grey, light-blue and light-red shaded regions around the lines are 1σ uncertainty regime derived using parameter estimations, each representing core, jet and the HST-1 knot sub-components, respectively. The downward arrows are indicative that given the limitations of the instrumental capabilities, these estimations should be treated as upper limits for further investigations.

In the text