© The Authors 2024

1. Introduction

Black holes are a fundamental prediction of Einstein’s theory of general relativity. A defining feature of black holes is their event horizon, beyond which the escape velocity exceeds the speed of light, resulting in a dark compact object (Schwarzschild 1916). Furthermore, photons approaching non-rotating black holes within a critical impact parameter $R_{\mathrm{c}}=\sqrt{27}GM/c^2$, where G is the gravitational constant, M is the black hole mass, and c is the speed of light, significantly bend in their trajectories by large angles and cannot escape to infinity (e.g., Hilbert 1917; von Laue 1921; Bardeen 1973). The value of R_c changes for spinning black holes, but within ≲4%. The photon

capture and redshift effects produce gravitational signatures in the observed images of astrophysical black holes immersed in background radiation fields, consisting of a characteristic dark center surrounded by ring-like emission (e.g., Luminet 1979), commonly referred to as the black hole shadow. Therefore, observing the shadow, measuring its properties, such as the ring diameter, and comparing it with theoretically predicted shadow morphologies provide unique opportunities to explore curved spacetime and extreme gravitational potential, allowing direct tests of general relativity. This is complementary to tests from gravitational waves, high precision astrometry, and cosmology (e.g., Event Horizon Telescope Collaboration 2022f; Abbott et al. 2016; GRAVITY Collaboration 2018; Ferreira 2019).

From an astrophysical point of view, black holes can be formed in many ways, and the masses of astrophysical black holes vary greatly, from relatively small (≲100 M_⊙; stellar mass) to supermassive scales (∼10⁶ − 10⁹ M_⊙); for example, readers can refer to Greene et al. (2020). Supermassive black holes (SMBHs) are predominantly located at the centers of galaxies, with larger black holes associated with more massive host galaxies (Kormendy & Ho 2013). When surrounding matter accretes onto the central SMBHs, the gravitational potential energy can be released in various forms, making active galactic nuclei (AGNs) shine across the electromagnetic spectrum, depending on the type and rate of mass accretion. Accreting SMBHs can also produce collimated relativistic plasma jets and release enormous amounts of energy into space in the form of strong magnetic fields and relativistic particles. Theoretical studies suggest the ultimate energy source of the jet to be either the accretion flow (Blandford & Payne 1982) or the central black hole itself (Blandford & Znajek 1977), in each case mediated through magnetic fields threading the plasma in the disk or near the horizon. Therefore, high-resolution imaging of accreting SMBHs can also provide unique constraints on modeling the detailed physics of the accretion disk and jet launching in AGNs (e.g., Blandford et al. 2019).

Since the angular diameter of the black hole shadow is d_sh = 2R_cD⁻¹ ≈ 10 GMc⁻²D⁻¹ where D is the distance to the black hole (Bardeen 1973), the imaging capability of the black hole shadow requires a mass-to-distance ratio sufficiently large to provide d_sh resolvable with the existing instruments. In this regard, SMBHs in the center of the radio galaxy M 87 (hereafter M 87^*) and our Galaxy, Sagittarius A^* (Sgr A^*), offer the best opportunities. As for M 87, its distance is measured via multiple methods using primary or secondary distance indicators, and we adopt D = 16.8 ± 0.8 Mpc (Bird et al. 2010; Blakeslee et al. 2009; Cantiello et al. 2018, see Event Horizon Telescope Collaboration 2019f). The central black hole mass is approximately 6 × 10⁹ M_⊙ (Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019f). The combination of D and M yields a d_sh ≈ 40 μas, which is the second largest value on the sky after Sgr A^* (M ≈ 4 × 10⁶ M_⊙, D ≈ 8 kpc, d_sh ≈ 50 μas; see Event Horizon Telescope Collaboration 2022a). M 87 is a representative Fanaroff-Riley I galaxy (Fanaroff & Riley 1974) with a prominent jet extending from the vicinity of the central black hole to kiloparsec scales (Owen et al. 1989; Reid et al. 1989; Junor et al. 1999; Asada & Nakamura 2012; de Gasperin et al. 2012; Hada et al. 2016; Mertens et al. 2016; Kim et al. 2018a; Nakamura et al. 2018; Walker et al. 2018; Park et al. 2019; Goddi et al. 2021; Lu et al. 2023). Multi-wavelength (MWL) studies of the jet from radio to γ-rays show both steady and sporadic bright flaring emission (see EHT MWL Science Working Group 2021). This makes M 87 one of the best candidates for ultra-high resolution imaging of an astrophysical black hole and the jet launching site.

The Event Horizon Telescope (EHT) is a global network of radio telescopes for very long baseline interferometry (VLBI) observations primarily at a wavelength of λ ≈ 1.3 mm or a frequency of ν ≈ 230 GHz (Event Horizon Telescope Collaboration 2019b). In VLBI, we measure complex visibilities or Fourier components of the radio brightness distribution of the sky at spatial frequencies u = (u, v), corresponding to the projected baselines in units of observing wavelength in a plane normal to the direction of the phase reference position (Thompson et al. 2017). These complex visibilities can be separated into amplitude and phase components, and are commonly constructed as closure quantities (Rogers et al. 1974; Readhead et al. 1980; Blackburn et al. 2020), which are immune to station based systematic corruptions. The longest baselines of the EHT result in an array with a theoretical diffraction limited resolution of 1/|u| ≈ 25 μas from the ground, and can directly resolve the black hole shadows in M 87^* and Sgr A^*. The first image of the black hole shadow of M 87^* was obtained by the EHT observations in April 2017 (Event Horizon Telescope Collaboration 2019a,b,c,d,e,f, 2021a,b, 2023, hereafter M 87^* 2017 I–IX). Also, another image of the black hole shadow, that of Sgr A^*, was published (Event Horizon Telescope Collaboration 2022a,b,c,d,e,f, hereafter Sgr A^* 2017 I–VI).

2. Previous EHT results and outline for new results

Below we provide a brief review of the major results on M 87^* from the 2017 EHT observations to introduce the 2018 EHT results, which are the main focus of this paper. For additional details, readers can refer to M 87^∗ 2017 I and Sgr A^∗ 2017 I for an overview of the first black hole shadow imaging of M 87^* and Sgr A^*, respectively.

2.1. Results from previous EHT observations of M 87^*

In April 2017, eight radio telescopes executed the first EHT VLBI observations of M 87^* at a wavelength of 1.3 mm (M 87^∗ 2017 II). Three independent EHT data calibration pipelines were developed and applied to validate the fringe detection and quantify their systematic uncertainties to be used for downstream analysis (M 87^∗ 2017 III). The M 87^* data reveal the presence of two nulls in correlated flux density at ∼3.4 and ∼8.3 Gλ (M 87^∗ 2017 III) and temporal evolution in closure quantities, indicating intrinsic variability of compact structure on a timescale of days or several light-crossing times for a few billion solar-mass black hole (M 87^∗ 2017 IV).

Independent analysis of the calibrated data using imaging and geometric modeling methods produce reconstructions that are well characterized by a bright and asymmetric ring-like structure with a diameter of 42 ± 3 μas, a central flux depression with an intensity ratio of > 10 : 1, and a position angle of the brightness asymetry of ∼180°. These measurements are consistent with expectations for a SMBH of mass M = (6.5 ± 0.7) × 10⁹ M_⊙ (M 87^∗ 2017 IV; M 87^∗ 2017 VI). A small level of interday variability was also detected, which was analyzed in follow-up studies (e.g., Georgiev et al. 2022; Satapathy et al. 2022). Analyses from outside the EHT collaboration featuring both imaging (Carilli & Thyagarajan 2022; Arras et al. 2022) and geometric modeling methods (Lockhart & Gralla 2022) agree that the 2017 EHT data are consistent with a ring-like structure. There is another external analysis of the 2017 data that claims not to support the EHT interpretation (Miyoshi et al. 2022).

In addition, Wielgus et al. (2020) reported results of 1.3 mm VLBI monitoring of the compact core of M 87 with a sparser pre-EHT array. By assuming a ring-like structure, this analysis found a consistent diameter across observations from 2009 to 2017 and strong evidence that the position angle of the brightness asymmetry varies by many degrees from year to year. Linear polarized features were also detected with a ring-like structure (M 87^∗ 2017 VII), along with the first limits on the circular polarization at the horizon scales (M 87^∗ 2017 IX).

2.2. Summary of the new results and paper organization

The EHT array in 2018 was significantly improved since 2017. The recording rate has increased by a factor of two compared to 2017, and the Greenland Telescope (GLT; Inoue et al. 2014; Chen et al. 2023) joined the EHT array for the first time. The addition of this new station results in a better (u, v) coverage compared with the 2017 EHT observations, particularly along the north-south direction. New EHT images of M 87^* can provide further evidence that the event-horizon-scale images from the 2017 observations are consistent with the prediction from general relativity that strongly lensed emission around a black hole should form a persistent ring-like structure. A detection of variability around the ring-like structure can also provide us with a better insight into the dynamics of the underlying plasma and the black hole angular momentum.

The structure of this paper and the main results are as follows: In Sect. 3, we describe the observations and processing of the 2018 EHT campaign data, providing a summary of the data properties and issues compared to that from the 2017 EHT campaign. As in 2017, we see a clear indication of a deep null in the visibility amplitudes around 3.4 Gλ. We also see significant closure phase deviations relative to 2017, indicating non-trivial structural changes.

We explain how we derived pre-imaging constraints on the expected compact flux density and its size in Sect. 4. Combining results from simple visibility modeling and information from lower frequencies, we settled on pre-imaging initial parameters for the compact flux density to be between 0.4 to 1.0 Jy and the size of the source to be ≤100 μas.

In Sect. 5, we describe the imaging procedure and investigate the primary image morphologies. We reconstruct the images using a variety of imaging techniques, ranging from traditional inverse imaging to Bayesian posterior exploration. We compare the recovered image structure across the different imaging methods, frequency bands, calibration pipelines, and observing dates. We recover a clear ring-like structure, with a deep central depression and a diameter of ∼42 μas. We also find that the position angle of the brightest region of the ring is between 200° and 230°, which is different from the position angle measured in 2017. We present the average images of M 87^* across methods on 2018 April 21 compared to 2017 April 11 in Fig. 1.

Inspired by the clear ring-like structure seen in the imaging methods, in Sect. 6 we quantify the support for the ring-like structure by comparing the Bayesian evidence of different geometric models against the data. We find that ring-like models are strongly preferred by the 2018 data. We then check the consistency of important physical quantities such as the diameter and position angle of the emission ring by directly modeling those features in the visibility domain. We find that all direct modeling methods find a position angle around 210°. The direct modeling methods also prefer a slightly higher diameter, around 45 μas. The diameter estimates from imaging and direct modeling live within each other’s error bars, and systematic differences in the median diameters can be attributed to model-dependent resolution effects.

In Sect. 7, we build upon the findings from the 2017 EHT campaign by comparing them against the imaging and direct modeling results from the 2018 data. We explore the persistence of the ring diameter between 2017 and 2018 and discuss the stability of the mass across the two years. We also discuss the possible origin of the shift in the position angle between the two years, which is consistent with our understanding of turbulent material around the black hole. We investigate the compact flux density estimates produced by the Bayesian image reconstruction and modeling methods, and find that the compact flux density is especially difficult to constrain in the 2018 data. The main conclusions of our work are summarized in Sect. 8. A detailed theoretical interpretation will be presented in a later paper.

3. Observations and data processing

In this section, we briefly describe the 2018 EHT observing campaign (Sect. 3.1), the data correlation (Sect. 3.2), the calibration and data reduction (Sect. 3.3), and finally provide some comparisons between the 2017 and 2018 data properties (Sect. 3.4).

3.1. Overview of the 2018 observing campaign

The 2018 EHT observing campaign was scheduled from 2018 April 18 to 29. Based on the expected opacity and weather conditions at each site, observations were triggered on six observing days with a total of nine participating stations, including all the stations which participated in the 2017 EHT campaign (see M 87^∗ 2017 II) with the addition of the GLT. Thus, the 2018 array includes the phased ALMA and the phased Submillimeter Array (SMA), made up of 43 12-m and seven 6-m dishes, respectively, as well as the following seven single-dish stations: Atacama Pathfinder Experiment (APEX), IRAM 30-m telescope (PV), South Pole Telescope (SPT), Large Millimeter Telescope Alfonso Serrano (LMT), Submillimeter Telescope (SMT), James Clerk Maxwell Telescope (JCMT), and the GLT. In Fig. 2, we show a map of the EHT array with the stations that participated in the 2018 campaign.

Four frequency bands centered at sky frequencies of 213.1 (band 1), 215.1 (band 2), 227.1 (band 3), and 229.1 GHz (band 4) were recorded, except for the GLT which observed only in bands 3 and 4. Each band has a bandwidth of 2048 MHz for every station, except for ALMA, which provides an effective bandwidth of 1875 MHz. This represents an improvement in spectral coverage by a factor of two compared to the 2017 observations. Additionally, the LMT had an increased sensitivity, given its increase in effective dish diameter from 32.5 m in 2017 to 50 m in 2018 and a change from a double-sideband (DSB) receiver to a two-single-sideband (2SB) one.

The telescope systems recorded both right-hand and left-hand circular polarization (RCP and LCP), except for JCMT, which observed only RCP throughout the 2018 campaign, and ALMA, which recorded linear polarization. Details on the treatment of different polarization types are given in Sect. 3.2.

Scans on M 87^* were included in four observing days with up to eight participating stations. The scan duration varied between 4 min and 7 min through the observations (although shorter scans were used in two days of the campaign). For calibration purposes, M 87^* scans were interleaved with 3−4 min duration scans on the source 3C 279. The median zenith sky opacities at 1.3 mm throughout the array on each observing day are provided in Table 1. In Fig. 3, we show the scan distribution per participating station in the schedule of the best observing days toward M 87^* and 3C 279.

The observations on April 28 had an overall higher median zenith opacity throughout the array than the other observing days (Table 1) and were notably affected by bad weather at ALMA. This prevented the detection of fringes in most baselines and resulted in many triangles with no closure information. Consequently, observations on April 28 did not pass the required quality checks (including the second phase of the ALMA quality assurance process, i.e., QA2) and were discarded in subsequent analyses.

The operations at the LMT were halted midway through the 2018 EHT campaign owing to local security issues. On the first day of the campaign (April 21), the observations started later than planned due to technical issues, which repeated in the middle of the observing day. At the start of that same day, the first few scans from Hawai‘i stations were lost to bad weather, and the SMA started observations even later due to technical issues. This resulted in less observing time for the LMT and Hawai‘i stations, which explains the differences in the (u, v) coverage in 2018 compared to 2017 (Fig. 4, left panels).

Furthermore, the (u, v) coverage on April 22 is minimal. The observations on April 25 suffered from coherence loss on baselines with APEX, bad weather at PV, and poor phasing efficiencies at the SMA. Hence, April 21 is the best observing day of the 2018 EHT campaign for M 87^*, followed by April 25.

3.2. Correlation

The data from the 2018 observing campaign were correlated at two correlation centers located at the Max-Planck-Institut für Radioastronomie in Bonn, Germany and at the MIT Haystack Observatory in Westford, Massachusetts, USA. The Distributed FX (DiFX, version 2.6.2) correlation package (Deller et al. 2011) was used. The data underwent multiple correlation passes to diagnose and correct data issues, including incorrect polarization labeling on the GLT and LMT. The data used in this paper is based on the fourth revision (Rev4) of the 2018 correlation products, released as part of the 2018 EHT Level 1 (L1) data package that also includes metadata required for absolute flux density calibration (Koay et al. 2023a). The final correlation produced 32 baseband channels, each 58 MHz wide with a spectral resolution of 0.5 MHz and averaged to a 0.4 s accumulation period. After correlation, the recorded RCP and LCP streams are multiplied to create the parallel-hand (RR, LL) and cross-hand (RL, LR) polarization products. Due to the non-matching sampling rates between ALMA and the other stations, a portion of the bandwidths is lost during correlation, resulting in a correlated bandwidth of 1.856 GHz per band (Matthews et al. 2018, M 87^∗ 2017 III). All the EHT stations observe in circular polarizations except ALMA, which observes in linear polarization. The PolConvert v1.7.9 (Martí-Vidal et al. 2016) was then run to convert the mixed polarization data products (XL, XR, YL, YR) to circular polarization basis on baselines to ALMA, described in more detail by Goddi et al. (2019).

3.3. Calibration and post-processing

To maximize fringe amplitudes when averaging the data in time and frequency, residual delays, delay rates, and atmospheric phase turbulence need to be corrected. Additionally, sampling losses and telescope bandpasses need to be calibrated. These signal stabilization (Janssen et al. 2022) data reduction steps are performed with two independent pipelines, EHT-HOPS (Blackburn et al. 2019) and the CASA-based (CASA Team 2022; van Bemmel et al. 2022) rPICARD software (Janssen et al. 2019a). These pipelines are described in more detail in M 87^∗ 2017 III and Sgr A^∗ 2017 II. No significant algorithmic changes were made to the EHT-HOPS pipeline fringe fitting methods compared to the 2017 analysis. For the rPICARD pipeline, we now use a priori determined solutions for the typically ∼7.3 ns single-band and multi-band delay offsets of the phased ALMA instead of solving for these offsets using the ALMA–APEX baseline as was done for the 2017 data.

The stabilized data are then loaded into the AIPS software package (Greisen et al. 2003), where a priori amplitude calibration is applied using the task APCAL to convert the correlation coefficients to a common flux density scale across the very heterogeneous EHT array. These are based on calibration tables in the standard ANTAB format containing information on telescope sensitivity, elevation-dependent gain curves, and time-dependent system temperatures derived from the raw metadata collected for each station, as described in Koay et al. (2023a). A summary of the amplitude calibration parameters of each station and their associated uncertainties is provided in Appendix A. The data are then averaged over time (10 s) and each 1.856 GHz band before being exported by AIPS as UVFITS files. Subsequently, we perform a polarimetric gains ratio calibration (polgainscal) to align the RR and LL phases in the HOPS data. We correct for phase offsets between the two polarization streams and allow the Stokes I visibilities to be obtained¹. Network calibration (M 87^∗ 2017 III) with 10 s solution intervals is then applied to correct the amplitudes of stations with a co-located partner (ALMA and APEX in Chile, SMA and JCMT in Hawai‘i), using redundancies and total flux density measurements from the phased ALMA (EHT MWL Science Working Group et al., in prep.). The network calibration takes advantage of the fact that the source is dominated by an unresolved core component on the intra-site baselines and that the flux densities on the intra-site baselines are comparable to that of the total flux densities of the core component as observed by the phased ALMA. We expect that the difference between the correlated flux density measured by the phased ALMA and the ALMA–APEX baselines should be less than 1%.

At each stage of the calibration, the data products are verified via consistency checks of closure quantities between bands and polarizations, as well as checks of the trivial closure quantities (Fig. B.1, see also M 87^∗ 2017 III and Wielgus et al. 2019, for more details). We also use these closure analyses, and the more novel closure trace analyses (Broderick & Pesce 2020) to roughly quantify the non-closing systematic uncertainties possibly arising from polarimetric leakages and bandpass effects (see Appendix B). For M 87^*, the non-closing systematic errors are estimated to be 1.2° in closure phases and 2.7% in log closure amplitudes, using both the HOPS and CASA data products. Assuming the errors are baseline independent, these translate to 0.7° systematic errors in visibility phases and 1.3% systematic error in visibility amplitudes. These uncertainties do not include station-based amplitude gain errors such as those described in Appendices A and E, which are corrected by model-dependent calibration during imaging.

Since some imaging pipelines use time-averaged data up to the length of a single scan, that is 4 to 5 min, we examine the level of decoherence losses in amplitude when the data are averaged in time. The details are presented in Appendix C. Such amplitude losses arise from imperfectly aligned phases after calibration and can lead to non-closing errors in visibilities. We find that 92.9% of HOPS and 85.6% of CASA M 87^* data have ≥90% coherence (i.e., exhibit decoherence losses of less than 10%) when averaged over the entire scan. When considering April 21 data only, the coherence levels are better, with 97.4% and 87.9% of HOPS and CASA scan-averaged data having better than 90% coherence.

We also check for consistency between the closure phases and log closure amplitudes derived from the HOPS and CASA data products for which signal-to-noise ratio (S/N) > 7 and find that there are < 4% outliers of ≥2σ when σ is equivalent to the thermal noise (σ_th). The fraction of outliers decreases to < 2% when σ includes the systematic uncertainties described above. In Fig. B.1, the fraction of ≥3σ outliers are consistent between HOPS and CASA also indicate that both pipelines produce consistent calibrated data sets. For reference, we provide detailed discussions of data issues in Appendix D, and in particular, highlight gain calibration issues for the GLT and LMT in Appendix E.

3.4. M 87^* 2018 data properties and comparisons with 2017

The 2018 EHT observations were made with an array similar to the one in 2017, except for the addition of the GLT and doubling the bandwidth of the 2017 observations, as mentioned in Sect. 3.1. The central frequency of bands 3 and 4 in 2018 coincides with the low- and high-bands in 2017, respectively. In Fig. 4, we show the (u, v) coverage of the best observing days in 2018 (April 21 and 25) in bands 2 and 3, overlaid on the (u, v) coverage of the best observing day in 2017 (April 11) in the low-band.

The GLT adds north-south coverage to the EHT array in bands 3 and 4, probing baseline lengths of ∼7.1 Gλ with the co-located Chilean sites and baseline lengths of ∼3.8 Gλ and ∼2.8 Gλ with intermediate stations. However, on April 21, the longest baselines in the east-west direction given by PV–Hawai‘i are missing, such that the east-west coverage is worse than that of April 25 and the 2017 observations. The LMT is missing on April 25, and its coverage is significantly reduced on April 21 compared to 2017, which impacts the compact flux density constraints in 2018.

Figure 5 shows the visibility amplitude versus (u, v) distance for 2018 April 21 and 25, overlaid on 2017 April 11. The a priori- and network-calibrated visibility amplitudes (see M 87^∗ 2017 III, for details) display the characteristic secondary peak beyond a deep amplitude minima at ∼3.4 Gλ, also observed in 2017 (M 87^∗ 2017 IV), despite the differences in the (u, v) coverage. A ring-like structure will also produce a second null beyond ∼8.3 Gλ. However, the longest baseline the EHT can probe is 8 Gλ between the Hawai‘i and PV stations, which has only one detection on 2018 April 25. While the S/N at the amplitude minima is lower (≲10) than that at the other regions, they provide very strong upper limits on the amplitudes (≲50 mJy) at these “null” locations. In particular, part of the LMT–Hawai‘i coverage samples similar fringe spacings as the GLT–PV and GLT–SMT baselines (although at almost perpendicular orientations), and part of the Chile–LMT coverage samples similar spacings as the GLT–LMT baseline. The GLT–PV baselines display lower amplitude than the GLT–LMT and GLT–SMT baselines. Thus, in 2018 we also see evidence of source anisotropy around the first null thanks to the GLT, as was observed in 2017 when the GLT was not available, but the corresponding fringe spacings were sampled by the 2017 LMT baselines.

The closure phases in all non-trivial triangles show differences in 2018 compared to 2017, providing evidence for significant changes in structure between these two epochs. In Fig. 6, we compare 2017 and 2018 coherently averaged visibilities on the same triangles (or similar ones) selected in 2017, where day-to-day variability of closure phases was observed (see Fig. 14 in M 87^∗ 2017 III). Here we also include a wide GLT triangle whose closure phases deviate from zero, indicating the presence of resolved asymmetric structure, similar to what has been observed in other triangles in 2017 before the inclusion of the GLT in the EHT array (M 87^∗ 2017 III). Gray circles and red diamonds in Fig. 6 show closure phase measurements from two days in 2017, and blue squares show closure phases measured in 2018. While some triangles show some small deviations in closure phase between days in 2017, the 2018 closure phases are qualitatively different, demonstrating clear indications of different asymmetric structure in the 2018 data.

4. Compact flux density and source size constraints

In this section, we analyze the network-calibrated visibilities of M 87^* to produce estimates of the compact flux density and source size. These are used to generate reasonable synthetic data sets (Sect. 5) and inform our choices of imaging parameters over which to survey for the non-Bayesian image reconstruction methods (see Sect. 5).

Measurements of the total flux density of the arcsecond-scale radio core in M 87 at 1.3 mm are between 1.0 and 2.0 Jy (Bower et al. 2015). Previous VLBI observations of M 87 at 1.3 mm between 2009 and 2017 indicated total compact flux densities at the milliarcsecond and microarcsecond scale fluctuating within a range of 0.5−1.0 Jy (Doeleman et al. 2012; Akiyama et al. 2015; Wielgus et al. 2020). We infer that the discrepancy between the total emission at arcsecond scales and VLBI scales is associated with the unconstrained extended emission outside of the compact emission region, presumably dominated by the jet.

Angular scales in M 87^* are sampled over nearly five orders of magnitude in (u, v) space by the EHT baselines. There is a large gap in the coverage between the intra-site baselines (ALMA–APEX and SMA–JCMT) of 0.1−1.1 Mλ and the inter-site baselines of 1.3−7.3 Gλ, which are sensitive to the sub-arcsecond and microarcsecond scale structures, respectively. The longest intra-site baseline is ALMA–APEX, with fringe spacings of 106−131 mas for M 87^*. On the other hand, that of the shortest inter-site baseline (SMT–LMT) is 139−166 μas. Therefore, the structures on spatial scales of ∼0.2−100 mas are not sampled in the visibility domain, and are thus unconstrained by the current EHT measurements.

The EHT data provide two relevant estimates of the total flux densities at different spatial scales. The first one is the total flux density of the arcsecond-scale radio core, constrained by the short intra-site baselines. The total flux density as measured by the ALMA-only interferometric observations was F_tot ≈ 1.13 Jy by averaging the values over four bands on April 21 and 22. The second one is the total compact flux density, F_cpct ≤ F_tot, which is the integrated flux density of the microarcsecond-scale structure constrained by inter-site baselines.

We estimate the range of F_cpct and the size of the compact emission region θ_cpct, represented by an equivalent Gaussian full-width half maximum (FWHM), following the two independent procedures in M 87^∗ 2017 IV. One is directly estimated from the 2018 EHT visibility data. With this procedure, we can derive the constraints on F_cpct and θ_cpct simultaneously. The other is a constraint on the F_cpct by utilizing quasi-simultaneous VLBI data based on an extrapolation between 1.3 mm and the longer wavelengths. In Appendix F, we derive a series of constraints on F_cpct and θ_cpct using these procedures.

The EHT constraints give F_cpct = 0.30 − 1.13 Jy with a θ_cpct between 39 and 98 μas. The source size constraints are roughly consistent with the 2017 constraints. The resulting constraints on F_cpct for 2018 are looser than those for 2017 due to the lack of additional constraints using SMT–PV and LMT–PV crossing baseline tracks. Nevertheless, the MWL constraints give a slightly tighter range of F_cpct = 0.5 − 0.7 Jy within the central 100 × 100 μas² of the nucleus, which is comparable to the F_cpct = 0.56 − 1.03 Jy estimated from the 2017 data. Therefore, we conclude that the total compact flux density of M 87^* in 2018 at 1.3 mm is broadly consistent with that of 2017 when considering all these estimates. While these flux density and size constraints in the pre-imaging analysis provide a valuable guide for subsequent imaging, analysis, and synthetic data generation, we do not strictly enforce them and leave the compact flux density as a free parameter in the subsequent imaging surveys.

In our pre-imaging constraints obtained above we predict a compact source size less than 100 μas with an unconstrained extended structure larger than ∼0.2 mas. Since there are no visibilities probing structure above ∼0.2 mas, the EHT data cannot be used to constrain jet emission above this scale. Moreover, the emission from the extended jet is significantly lower in surface brightness than the core region, even in 3 mm observations (Hada et al. 2016; Kim et al. 2018a; Lu et al. 2023), so the jet emission is expected to be even weaker at 1.3 mm due to its optically thin (steep) synchrotron spectral nature (EHT MWL Science Working Group 2021). Even though we have improved coverage with the GLT in 2018, the continued lack of short baselines means we do not expect the 2018 EHT data can constrain this very dim jet emission.

5. Imaging and image domain analysis

The EHT’s sparse (u, v) coverage results in an ill-posed inverse problem that prevents the recovery of a unique image from a measured set of visibilities. By using realistic priors, it is nevertheless possible to reconstruct brightness distributions that are consistent with the data. While these reconstructed distributions are still non-unique, they should represent a family of images that are interpretable within our understanding of the physical system. While it is possible to produce reasonable image reconstructions using only closure quantities, all methods featured in this section utilize some data products that are complicated by persistent systematic calibration errors in both the data amplitude and phases. This requires each reconstruction method to undergo a careful self-calibration or station gain reconstruction process in order to extract the most information from the data. As in M 87^∗ 2017 IV, both inverse and forward modeling were used in the image reconstructions. For inverse modeling, a CLEAN-based algorithm (e.g., Högbom 1974; Clark 1980) was employed with DIFMAP software (Shepherd 1997, 2011). The forward modeling techniques used in this study were regularized maximum likelihood (RML) methods represented by eht-imaging (Chael et al. 2018, 2019a,b) and SMILI (Akiyama et al. 2017a,b), and two Bayesian sampling methods to explore the full posterior space with THEMIS (Broderick et al. 2020a,b, 2022a) and Comrade (Tiede 2022).

To explore how the reconstructed images were affected by various imaging and optimization choices, we generated synthetic data sets using ten geometric models (cres-90, cres0, cres90, cres180, dblsrc, disk, ring, edisk, point+disk, point+edisk) and a snapshot of a general-relativistic magnetohydrodynamic (GRMHD) simulation (see Appendix G). The synthetic data were designed to have properties similar to the EHT M 87^* visibility amplitudes, including prominent amplitude nulls, but also reflected a diversity of both ring-like and non-ring source structures. This procedure was similar to previous EHT imaging analyses (M 87^∗ 2017 IV; Sgr A^∗ 2017 III). For the RML and CLEAN imaging methods, four of these ten data sets (cres180, ring, dblsrc, disk) were used as training sets to evaluate the ability of unique combinations of imaging parameters to faithfully reconstruct images close to ground truth. For the Bayesian methods, all eleven data sets were used as validation exercises to confirm the ability of the Bayesian methods to reconstruct diverse image structures.

Section 5.1 presents the image strategies used by each method. Section 5.2 presents the images and discusses consistency between different methods, frequency bands, calibration pipelines, and observation dates. Section 5.3 describes the process by which we measure important image domain quantities such as the diameter and position angle.

5.1. Strategy of the imaging analysis

We conducted parameter surveys with three imaging pipelines; employing CLEAN using the DIFMAP software, and RML methods using eht-imaging and SMILI. We used the four training data sets to explore the impact of different imaging assumptions, such as hyperparameters (weights for both the data and regularizers, see M 87^∗ 2017 IV) and optimization choices, on the resulting image morphology. All images in the surveys were reconstructed from data calibrated as described in Sect. 3.3. From these surveys, we selected a “Top Set” of parameter combinations for each method, which represent the set of best-fit images to the data. Each synthetic data set was blurred with a nominal beam width corresponding to the effective resolution of the longest baseline. We then compared the normalized cross correlation, ρ_NX, of the unblurred ground-truth image against the blurred ground truth to determine a ρ_NX cutoff value. We then found the model ρ_NX value for each set of image parameters against the blurred synthetic data, and a set of model image parameters passed to the next step if the model ρ_NX is above the cutoff value described above. From this pruned set, the Top Set images were then selected by calculating the normalized closure quantity chi-squares of the images for the real M 87^* data, and taking only the images that have a normalized χ² <  2 on data with 0% systematic noise (to account for non-closing errors) for the RML methods and 10% systematic noise for CLEAN (see Sect. 6.3.1 in M 87^∗ 2017 IV and Sect. 5.2.1 in this paper for more details). The closure quantities were averaged over the entire scan before comparison to ensure sufficient S/N (Rogers et al. 1995; Blackburn et al. 2020). The distribution of images in the resulting Top Sets on the real M 87^* data is a proxy for the uncertainties due to different imaging strategies and assumptions. For each method we also defined a “fiducial” image by converting the model ρ_NX of each Top Set image to an effective blurring width, averaged this width across all synthetic data sets, and found the set of parameters that produced the minimum effective blurring width. A summary of the Top Set parameters for each method is shown in Table 2.

Validation of the Top Set parameter combinations was then performed by imaging the remaining six geometric models as well as a GRMHD snapshot image and ensuring that the resulting images closely match their ground truths (see Appendix G.1 and Fig. G.3). We also found very few artifacts in the image reconstructions across the Top Set, as shown by stacking Top Set images for each synthetic data set as shown in Fig. G.4. The detailed imaging strategies of the three pipelines are summarized in Sects. 5.1.1–5.1.3.

As mentioned above, we also utilized two methods, THEMIS (Broderick et al. 2020a) and Comrade (Tiede 2022), that followed a Bayesian posterior sampling approach to image reconstruction. To quantify the imaging uncertainty, THEMIS and Comrade use Bayesian inference and cast VLBI imaging as a Bayesian inverse problem. Since these Bayesian methods do not require training for selecting model hyperparameters, all synthetic data models are treated as validation data sets. In addition to producing a best-fit image to the relevant visibility data, these approaches produce an image family that captures the image uncertainty arising from the measurement uncertainty of the data given the assumptions of the model. This permits us to quantitatively assess the significance of image structures, array calibration quantities, and the relationship between these features. More detailed descriptions of these methods are given in Sects. 5.1.4 and 5.1.5.

5.1.1. DIFMAP

The CLEAN algorithm (e.g., Högbom 1974; Clark 1980) has been widely used for image reconstruction in radio interferometric imaging. This algorithm uses an inverse modeling approach to derive a sparse reconstruction on the image domain using the interferometer’s point-source response (i.e., the dirty beam). The CLEAN algorithm used in this paper assumes that the sky brightness distribution is a collection of point sources. The imaging process involves generating point sources (CLEAN components) at the location of the brightness peak in the dirty image (defined as the 2D Fourier transform of the measured visibilities) and iteratively subtracting them until a stopping criterion is reached. The final image is obtained by convolving the CLEAN components with a Gaussian restoring beam and adding the residual image to represent the residual noise. Restrictions, known as CLEAN boxes or CLEAN windows, can be placed on the area of the map in which the CLEAN components are searched and are especially important for data with sparse (u, v) coverage.

We conducted imaging with CLEAN using DIFMAP (Shepherd 1997, 2011) based on the imaging pipelines utilized for the 2017 EHT images of M 87^* (M 87^∗ 2017 IV) and of Sgr A^* (Sgr A^∗ 2017 III). We perform the DIFMAP analysis using data averaged to 10 s, and perform self-calibration using the same averaging time-scale. The presented pipeline has a few minor modifications compared to the 2017 M 87^* pipeline. We omitted setting the estimated flux density expected from a baseline of zero length since there are actual data from intra-site baselines at very short lengths. Based on the gain error budget presented in Appendix E, we set acceptable limits on amplitude gain correction factors during the self-calibration process, which were within the range of 0.5−2.0, instead of the 0.83−1.2 range used in the 2017 pipeline. We found large gain correction factors for the GLT, LMT, and PV at least for several scans. With the addition of the GLT as a new station in the array, the (u, v) coverage of the 2018 EHT data was improved, especially on intermediate baselines. This helped to suppress the side lobes in the synthesized beam, improving the CLEAN image reconstruction. The pipeline surveyed five parameters: the total assumed compact flux density, cleaning stopping condition, relative weight correction factor for ALMA in self-calibration, diameter of the CLEAN window, and the power-law scaling of the (u, v) density weighting function. The same parameter ranges used for the 2017 M 87^* imaging were used, with a compact flux density of 0.4 Jy to maintain consistency with RML imaging methods. As explained in M 87^∗ 2017 IV, the amplitudes measured from LMT suffer from relatively poor a priori calibration and thus require self-calibration with the initial imaging result.

We implemented phase-only self-calibration with geometric models to mitigate the impact of a priori calibration uncertainties in EHT measurements during the imaging process. Two strategies were employed for selecting the initial model in the parameter survey using synthetic and real data sets. The first strategy involved employing a point source with a flux density of 1 Jy positioned at the origin of the dirty image, similar to the approach used in the 2017 M 87^* imaging (M 87^∗ 2017 IV). The second strategy aimed to choose a better geometric model for the initial phase-only self-calibration, following the methodology of the 2017 Sgr A^* imaging (Sgr A^∗ 2017 III). For the latter strategy, synthetic visibilities were generated for various geometric models, including a Gaussian with 15 μas FWHM (representing an unresolved symmetric model), a uniform disk with sizes ranging from 56 to 84 μas in steps of 4 μas, and a uniform ring with sizes ranging from 36 to 68 μas (also in steps of 4 μas, without width). The best model was selected based on the closure phase normalized χ². Unlike the Sgr A^* imaging, we did not incorporate further CLEAN reconstruction after determining the best initial model because it resulted in instability in the best initial model for the data presented in this paper. The initial models chosen for each data set can be found in Table 3. Our analysis revealed that selecting the geometric model based on closure phase fitting outperformed the point source model strategy for synthetic data reconstruction, yielding more Top Set images. This outcome was anticipated, considering the more complex source geometries in synthetic and real M 87^* data compared to a simple point source. Even though the point source model strategy produced fewer Top Set images, the images themselves look similar to the Top Set images produced using the closure phase fitting strategy. We present the results from closure phase fitting in the main text and include the point source modeling strategy in Appendix H.

The DIFMAP pipelines presented in this paper have been updated compared to those used in previous EHT imaging analyses (M 87^∗ 2017 IV; Sgr A^∗ 2017 III). We performed a backward compatibility test for these pipelines on the 2017 M 87^* data (April 11, low-band). We found that the resulting images agree with those presented in M 87^∗ 2017 IV. We present the fiducial images of the M 87^* data and the synthetic data in Appendix I.

5.1.2. eht-imaging

The Python package eht-imaging (Chael et al. 2018, 2019a,b) is an RML-based VLBI imaging software capable of producing images by placing different relative weights on the fits to closure quantities and complex visibilities. The results in this paper were produced by performing a parameter survey using the 2017 M 87^* imaging pipeline². The surveyed parameters consist of the total assumed compact flux density, the fractional systematic error on the measured visibilities, the FWHM of the circular Gaussian used as the initial and prior image, and weights for four of the regularizers. See Table 2 for a summary of the surveyed parameter space and Appendix A in M 87^∗ 2017 IV for detailed definitions of each regularizer. The range of surveyed compact flux density values were chosen based on the pre-imaging constraints outlined in Sect. 4, while the range of surveyed regularizer weights were chosen based on experience from values surveyed in M 87^∗ 2017 IV. All images were reconstructed with a 128 μas FOV and a 64 × 64 pixel grid.

The imaging pipeline starts by loading and coherently scan-averaging the data. Then the correlated flux densities at intra-site baselines were rescaled by the given compact flux density to remove the contributions from unresolved extended emission outside the FOV. We added an additional fractional systematic error term to the visibilities’ error budget to account for unknown amounts of non-closing errors in the data. As measurements taken by the LMT suffer from large gain uncertainties, we performed an initial amplitude-only self-calibration to the LMT data and adopted station based gain corrections for LMT. This self-calibration is performed with a circular Gaussian geometric model with FWHM of 60 μas and flux density of 0.6 Jy, chosen to fall in the center of the compact flux density limits derived in Sect. 4. Lastly, the visibility amplitudes were inverse tapered with a 5 μas FWHM circular Gaussian to enforce an angular resolution limit on the final reconstructed image.

After these pre-imaging calibration steps, the pipeline proceeded with four iterations of imaging and self-calibration. The imaging was initialized with a circular Gaussian of FWHM and compact flux density specified by the given parameter combination. The details of the self-calibration and the relative weights placed on fits to the various data products were modified between each iteration to reflect progressing amounts of confidence in the gain and phase solutions. The first two rounds of self-calibration were performed only on the phases while the last two rounds were performed on amplitudes and phases. For the relative data weights, we began the first round of imaging by placing unity weight on the closure quantities, a fifth of that on the visibility amplitudes, and no weight on complex visibilities. As we progressed through iterations, we removed weight on the visibility amplitudes and allowed non-zero weight on complex visibilities. The ratio between weights placed on close quantities compared to complex visibilities decreased in later iterations as we converged on a phase solution. All self-calibration rounds were performed on the scan-averaged time intervals.

Each iteration involves several attempts at producing an image to prevent the imaging function from getting stuck in a local minimum. Each attempt utilized the previous best-fit image blurred to the nominal array resolution as the initial image. At the end of all four iterations of imaging, we ensured consistency with the original data and limited the angular resolution of reconstructed features by convolving the final image with the same 5 μas Gaussian used for inverse tapering in the pre-imaging calibration step.

5.1.3. SMILI

The imaging pipeline SMILI (Akiyama et al. 2017a,b) is a Python and FORTRAN RML-based imaging software. Similar to the survey conducted with eht-imaging, the SMILI parameter survey allows variation of the following parameters: the total assumed compact flux density, the FWHM of the circular Gaussian used for the prior image, the fractional systematic error on the measured visibilities, and the weights of three regularizers. See Table 2 and Appendix A in M 87^∗ 2017 IV for details of the SMILI parameter ranges and descriptions of the regularizers, respectively. All images are reconstructed with a 128 μas FOV and a 64 × 64 pixel grid.

Before imaging, the script coherently scan-averaged the visibilities. Pre-calibration of the LMT gains and the addition of non-closing systematic errors were performed as described in the second paragraph of Sect. 5.1.2. We performed four imaging cycles for each self-calibration stage, self-calibrating only to the final reconstructed images in each of the three cycles. Reconstructions at each stage were initialized with a circular Gaussian of FWHM = 20 μas containing the compact flux density of the given parameter combination. Fractional uncertainties of 50%, 30%, and 5% were added in quadrature to the error of visibility amplitudes on LMT, GLT, and baselines to other antennas, respectively, to account for residual errors in the amplitude calibration. After the first imaging attempt in each stage, subsequent initializations used the previously obtained image re-centered to the image center of mass and blurred with a 20 μas FWHM circular Gaussian. Each imaging cycle performed 1000 iterations of the L-BFGS-B algorithm (Byrd et al. 1995; Zhu et al. 1997) used for the gradient-descent optimization in SMILI’s image solver. After the three imaging cycles, complex visibilities were self-calibrated with their output image. In the first two self-calibration stages, the imaging step used closure data (closure amplitudes and phases) and visibility amplitudes. The imaging used closure data and complex visibilities for the final two self-calibration stages. Similar to eht-imaging, SMILI also performed all rounds of self-calibration on scan-averaged timescales.

5.1.4. THEMIS

We employ THEMIS, a generic modular parameter estimation framework specifically developed to compare parameterized models with VLBI data produced by the EHT (Broderick et al. 2020a). The basic THEMIS image model aims to reproduce the on-sky brightness distribution and the time-dependent complex station gains. For reproducing generic brightness distributions, THEMIS utilizes an adaptive splined raster model defined by a set of brightness control points arranged on an adjustable rectilinear grid (Broderick et al. 2020b), which we call a “Themage” (THEMIS image, see Appendix J for details). We used a raster resolution of 5 × 5, which was chosen after a survey of other resolutions ranked by the Bayesian evidence and described in more detail in Appendix K. Since there was a significant difference in the flux density between the intra-site baselines and the shortest inter-site baselines, we also included a large scale asymmetric Gaussian component to approximate the flux density contribution from unresolved emission outside of the FOV.

The THEMIS image model’s raster and Gaussian components are constructed in the visibility domain and directly fit against the scan averaged full complex visibilities generated from the standard 10 s averaged network calibrated data described in Sect. 3.3, to which we added an additional 1% fractional systematic error in quadrature to the thermal uncertainties. When fitting complex visibilities, thermal errors are strictly Gaussian and improve the smoothness of the likelihood surface. When we reconstructed the scan-averaged station gains, we imposed Gaussian priors on the logarithmic gain amplitudes and flat priors on the gain phases. Informed by the analysis of crossing baseline tracks and the EHT station flux density calibration parameters (Table A.1), for the April 21 data, we imposed prior widths of 10% for all stations except LMT and GLT, which used 30% and 100% prior widths respectively. We used the same gain priors for the April 25 data except for permitting a 100% gain prior for PV, since it was noted that PV exhibited large amplitude fluctuations after UTC 02:00 due to poor weather (see Appendix D).

5.1.5. Comrade

Comrade is a Bayesian differentiable modular modeling framework written in the Julia programming language for use with VLBI (Tiede 2022). In this work, we applied non-parametric modeling and include a rasterized image model similar to the one described in Broderick et al. (2020b). We first scan-averaged the standard 10 s averaged network calibrated data described in Sect. 3.3. We added an additional 1% fractional systematic error in quadrature to the visibility thermal error and then extracted the closure phases and visibility amplitudes. We removed closure phases with S/N <  3 and any closure phases with intra-site baselines.

Our model for all image reconstructions consisted of a rasterized image model (of dimensions N_x × N_y), a Gaussian of FWHM = 1 mas, and scan-averaged station gains. The Gaussian was used to model the amplitudes of short intra-site (ALMA–APEX, JCMT–SMA) baselines. In the image model, the grid of raster points was convolved with a third-order basis-spline (B-spline) kernel (Eq. (L.2)) to generate flux densities that smoothly varied in all directions. Comrade’s raster model is similar to THEMIS in Sect. 5.1.4, except that we used the B-spline kernel. A detailed description of the image model is given in Appendix L.

The hyperparameters of the model are the FOV and raster size (N_x × N_y). For April 21 bands 3 and 4, we used a 12 × 12 pixel raster with a 7.5 μas pixel size, which was ∼1/3 of the nominal resolution and FOV of 90 μas. This FOV was chosen by checking the visibility-amplitude residuals to incorporate all the flux density within the FOV. If the FOV is too small (< 65 μas) then the residuals are large. If the FOV is too large (> 90 μas), we do not have information on those scales from the visibilities and we get poor reconstructions for the synthetic data sets. For bands 1 and 2, we shrunk the FOV to 75 μas due to the lack of intermediate baselines from GLT; maintaining the same pixel size resulted in a 10 × 10 raster. The hyperparameters are the same for M 87^* and the synthetic data. The hyperparameters were changed depending on (u, v) coverage for different bands and days, and the specific details are given in Table 4.

We formed the visibility amplitude likelihood and closure phase likelihood as described in Appendix F of Sgr A^∗ 2017 IV. We used a uniform distribution 𝒰(0.0, 1.5) Jy for the prior on the total flux density f and 𝒰(0.0, 1.0) for the fraction of the total flux density f_g for the flux density of the large-scale Gaussian. For the raster pixel fluxes, we used a symmetric Dirichlet distribution on the simplex (see Eq. (L.5)). For the station gain priors, we used a normal distribution for the log-gain amplitudes for each station. The prior width of station gain priors was similar to THEMIS in Sect. 5.1.4. A detailed description of the prior distributions is given in Appendix L.

Finally, the un-normalized posterior was formed by taking a product of the prior and the likelihood distributions. To sample from the posterior, we used a two-stage strategy. First, to find a reasonable starting location, we adopted the L-BFGS optimizer (Liu & Nocedal 1989; Mogensen & Riseth 2018), running it 20 times, and initializing it from random draws from the prior. We then selected the location that optimized the log joint probability density at our starting location for the Hamiltonian Monte Carlo No-U-Turn Sampler (Hoffman & Gelman 2014; Xu et al. 2020). We ran the sampler for 12 000 steps; the first 10 000 are adaptation steps. To check for Markov chain Monte Carlo (MCMC) convergence, we measured the split-$\widehat{R}$ (the average variance of draws in one chain compared to the variance across all chains) and the effective sample size. After sampling, we calculated statistics of the posterior, such as the mean and standard deviation of the images and station gains. Since the phase center for EHT-like data is unconstrained, we used the image centroid to align the images.

5.2. Presentation of images

5.2.1. Comparison between methods, bands, and days

Figure 7 shows representative images of M 87^* produced with each of the five imaging methods using data from both calibration pipelines. Images are shown for all four bands on April 21 and April 25. For the DIFMAP and RML methods, we present fiducial images. For THEMIS and Comrade, we display a random sample from the posterior.

We first discuss images from the band 3 data on April 21, which represent (together with band 4) the best (u, v) coverage and the most stable imaging results. On April 21, DIFMAP, eht-imaging, and SMILI could all produce a non-zero number of Top Set images for all bands. A visual inspection shows that all images display similar characteristics, including diameter, a central flux depression, and a brightness asymmetry in roughly similar positions (see Appendix G.2). Apparent differences in detailed structure between methods can be attributed to differences in the effective resolutions of the imaging pipelines. For example, a 20 μas deconvolution beam was used for DIFMAP imaging, so DIFMAP images tend to have a larger ring width and weaker central depression. In addition to good agreement in structure among the image reconstructions, we also find good agreement in the reconstruction of the time-dependent station gains for both the synthetic and real data (see Appendix G.1).

Figures 8 and 9 show the data visibility amplitudes and closure phases with the corresponding model visibility amplitudes and closure phases computed from each April 21 band 3 image. These figures demonstrate that the images produced by all pipelines fit the data well. When comparing the model closure quantities of a Top Set image against the real data closure quantities, we used the mean squared standardized residual (χ²), normalized by the number of data points to be compared. For example, we define the normalized closure phase ${\chi }_{\mathrm{CP}}^2$ as:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{CP}}^2=\frac{1}{N_{\psi }}{\displaystyle \sum \frac{|e^{i{\widehat{\psi }}_{\mathrm{C}}}-e^{i{\psi }_{\mathrm{C}}}{|}^2}{{\sigma }_{{\psi }_{\mathrm{C}}}^2}},\end{array}$$(1)

where N_ψ is the number of closure phase data points, ${\widehat{\psi }}_{\mathrm{C}}$ is the model closure phase, ψ_C is the real data closure phase, and σ_ψC is the closure phase standard deviation. Similarly, we can construct the normalized log closure amplitude ${\chi }_{log\mathrm{CA}}^2$ as:

$$\begin{array}{c}\hfill {\chi }_{log\mathrm{CA}}^2=\frac{1}{N_{ln{A}_{\mathrm{C}}}}{\displaystyle \sum \frac{|ln{\widehat{A}}_{\mathrm{C}}-ln{A}_{\mathrm{C}}{|}^2}{{\sigma }_{ln{A}_{\mathrm{C}}}^2}},\end{array}$$(2)

where N_lnAC is the number of log closure amplitude data points, ${\widehat{A}}_{\mathrm{C}}$ is the model closure amplitude, A_C is the real data closure amplitude, and ${\sigma }_{\ln A_{\text{C}}}^2$ is the log closure amplitude standard deviation. Readers can refer to Sect. 2.1 of M 87^∗ 2017 IV for additional details about constructing normalized closure χ² values, which should not be confused for the formal definition of a reduced χ². In Table 5, we show the normalized χ² values for the fiducial images and normalized χ² statistics across the Top Sets for DIFMAP, eht-imaging, and SMILI. The fiducial images from the RML methods are consistent with the data roughly within the thermal noise, and the normalized χ² values have little scatter across the Top Sets. We added 10% systematic uncertainty to the real data in the evaluation process for DIFMAP because the image generation process did not use closure quantities (unlike the RML methods), and down-weighted the ALMA visibilites during self-calibration, resulting in comparable normalized χ² values to the RML methods.

For the Bayesian imaging methods, the fit quality is assessed through comparisons of the complex visibility reduced χ² for THEMIS, the visibility amplitude and closure phase reduced χ² for Comrade, log-likelihoods, the logarithm of the Bayesian evidence (log(Z)), and distribution of the residuals. For this fit quality assessment, the maximum likelihood model from the sampling chains were used for THEMIS and maximum a posteriori model for Comrade. We present the relevant reduced χ² values in Table 6.

We calculate the reduced χ² by dividing the standard χ² by the sum of the number of independent model and gain parameters subtracted from the number of data points, or for the case of the complex visibility reduced χ² (Red. ${\chi }_{\mathrm{V}}^2$):

$$\begin{array}{cc}& \mathrm{Red}.N_{\mathrm{V}}=N_{\mathrm{V}}-N_{\mathrm{model}\mathrm{DOF}}-N_{\mathrm{gain}\mathrm{DOF}},\hfill \end{array}$$(3)

$$\begin{array}{cc}& \mathrm{Red}.{\chi }_{\mathrm{V}}^2=\frac{1}{2(\mathrm{red}.N_{\mathrm{V}})}{\displaystyle \sum \frac{|\widehat{V}{-V|}^2}{{\sigma }_{\mathrm{V}}^2}},\hfill \end{array}$$(4)

where N_V is the number of complex visibility data points, N_model DOF is the number of model degrees of freedom, N_gain DOF is the number of gain degrees of freedom (we multiply the final reduced number by 2 since we fit both the real and imaginary components), $\widehat{V}$ is the model complex visibility, V is the data complex visibility, and σ_V is the complex visibility standard deviation. One can construct similar reduced χ² quantities for the visibility amplitudes (Red. ${\chi }_{\mathrm{A}}^2$) and closure phases (Red. ${\chi }_{\mathrm{CP}}^2$) using the appropriate model and data products, and properly counting the number of data points, model, and gain degrees of freedom to construct the reduced denominator.

The THEMIS and Comrade raster reconstructions produce low reduced χ² across all bands and both days. Since the Comrade raster used a substantially larger raster grid than THEMIS, we expected the Comrade fits to produce reduced χ² well below unity. In the THEMIS fits, we find that band 1 reconstructions feature systematically worse performance. This may be indicative of additional systematic errors in those data sets.

Details of the observations contribute to the differences between images from different bands and days. The improved (u, v) coverage in bands 3 and 4, given by the participation of the GLT, allows for improved reconstructions of M 87^* images. The GLT is especially important in probing the null point near 3.4 Gλ. This is proven by the increased number of Top Set images for the trained methods and the cleaner reconstruction of the ring-like morphology compared to bands 1 and 2.

On the April 25 data, DIFMAP, eht-imaging, and SMILI struggled to produce a significant number of Top Set images – none of the methods could produce Top Set images for bands 1 and 2. The Bayesian methods also struggle to produce a ring-like morphology for data taken on this day. This performance issue is mainly due to a lack of data from LMT on this day. The LMT–SMT baseline provided the only probe of the visibility structure around 1 Gλ; the lack of this baseline hampered imaging, especially the capacity to select a preferred compact flux density.

5.2.2. Compact flux density for RML and CLEAN

The M 87^* fiducial images were reconstructed with a total compact flux of 0.4 Jy for DIFMAP, 0.5 Jy for eht-imaging, and 0.6 Jy for SMILI; all of these values fell within the range of allowed compact flux densities derived from the pre-imaging constraints (Sect. 4). The RML methods only modeled the compact emission (since they rescale the intra-site baselines), so it is necessary to confirm that the image structure would not significantly change within the range of compact flux densities allowed by the pre-imaging constraints. Though each pipeline surveyed over several different assumptions on the total compact flux density (see Table 2 for a summary), we found no strong preference for a particular value. As a complementary check we generated images without intra-site visibilities and verified that the image structure is not significantly impacted. See Sect. 7.3 for a related discussion of the compact flux density estimates from the Bayesian methods, and Appendix M for supplementary synthetic data validation tests at different compact flux densities for the RML and CLEAN methods.

The difference of the flux densities between F_tot and F_cpct implies there is some emission outside the compact region. However, its structure (e.g., scale and direction) is unconstrained by the data (Sect. 4). The sub-milliarcsecond scale jet is not detected by the 2017 and 2018 EHT data, but is clearly seen in the 3 mm images (Hada et al. 2016; Kim et al. 2018a; Lu et al. 2023). This sub-milliarcsecond scale jet is presumably resolved out due to lack of visibility and closure information from short enough (u, v) spacing (e.g., ≲1 Gλ). The contribution from the sub-milliarcsecond scale jet becomes less important at (u, v) distance ≥1 Gλ in 3 mm data (Lu et al. 2023). Assuming the jet emission is optically thin synchrotron radiation (EHT MWL Science Working Group 2021), the sub-milliarcsecond jet can be less bright at higher frequencies, so we expect that it would be challenging for the EHT to detect without much better (u, v) coverage at short baselines.

As we did in M 87^∗ 2017 IV, the geometric models in the synthetic data sets included emission from a simple extended jet modeled out to ∼2 mas. None of the image reconstructions on the synthetic data recover this extended component. Similarly, no method was able to recover diffuse or extended emission displaced from the ring seen in the GMRHD synthetic data. While several studies suggest the presence of a small-scale diffuse structure in the 2017 EHT data reconstructions which aligns with the limb-brightened jet observed at 3 mm (Arras et al. 2022; Carilli & Thyagarajan 2022; Broderick et al. 2022a), no constraints on larger-scale (> 0.2 mas) jet emission are supported by the data.

5.2.3. Image statistics

The distributions of both the Top Set and posterior images help us understand the image uncertainties associated with each method. Figures 10–12 show the image- and visibility-domain uncertainties associated with the image sets for eht-imaging, THEMIS, and Comrade, respectively.

The uncertainties shown for eht-imaging are a proxy for the uncertainties in choosing the regularizer weights and hyperparameters. Similar to the corresponding figure in M 87^∗ 2017 IV, we find that the high image uncertainties correspond to locations of high brightness temperature and visibility-domain uncertainties primarily due to gaps in (u, v) coverage. We see in the image domain figure of normalized standard deviation that small concentrations of uncertainties exist at various locations along the ring. However, they are less pronounced than the “knots” in 2017. The very small amount of uncertainty in the central flux depression indicates the robustness of the ring-like feature, as it occurs in nearly every Top Set image. The normalized standard deviation visibility domain image shows a similar concentration of uncertainty around 2 Gλ. The fractional standard deviation image shows a sharp increase in uncertainty at the boundary between (u, v) distances probed by the EHT versus those outside the maximum probed distance. This indicates that the RML methods are not assigning high confidence to image features smaller than the minimum scale probed by the observations.

Comparatively, the uncertainties derived from THEMIS and Comrade stem from statistical uncertainty of the image posterior space given the model assumptions. Both THEMIS and Comrade find that the ring-like feature is robustly recovered, finding it has a low fractional standard deviation. Additionally, both THEMIS and Comrade find no evidence for non-ring emission. For THEMIS, this can be seen from the fact that the raster is concentrated around the ring. Substantial non-ring emission would cause the THEMIS raster to expand to cover the region, as seen in Fig. G.3. Comrade’s larger assumed FOV allows for a direct quantitative estimate of the significance of non-ring emission and finds no significant (> 3σ) detection of non-ring emission. Moving to the visibility domain, we see that the uncertainties reported by THEMIS and Comrade are qualitatively different and generally find that THEMIS’ uncertainty is much smaller than Comrade’s. This disparity reflects the differences in each method’s model specification. For short (u, v) distances, the lower standard deviation for THEMIS is due to its sparsity-inducing prior and fitted raster dimension, which tends to produce smaller image FOVs than Comrade, reducing the uncertainty on scales ≲3 Gλ. Additionally, since THEMIS uses a relatively large pixel size (≳12 μas), it is unable to produce significant amplitudes on scales ≳10 Gλ, which drastically decreases its measured visibility-domain uncertainty compared to Comrade.

5.3. Image domain feature extraction (IDFE)

While the imaging methods described above produce excellent fits to the visibility data and produce a consistent ring-like structure, the specific choice of regularizers and weights, masks, and parameter combinations used to produce these image reconstructions are generally agnostic to the image morphology and the hyperparameters may not directly relate to physical quantities. To produce measurements of ring diameter, width, flux asymmetry, depth of the central depression, etc., we pass the reconstructed images through feature extraction algorithms similar to the analyses featured in M 87^∗ 2017 IV and Sgr A^∗ 2017 III. The entire Top Set is used for the CLEAN and RML methods to obtain these measurements, while 500 random samples drawn from the respective posteriors are used for THEMIS and Comrade. Negative-valued pixels in THEMIS images are replaced with zeros in line with the previous Sgr A^* analysis (Sgr A^∗ 2017 III; Sgr A^∗ 2017 IV).

We use two image domain feature extraction tools, REx and VIDA. The Ring Extractor (REx) is available as part of eht-imaging and is described in detail in M 87^∗ 2017 IV. The second IDFE tool, Variational Image-Domain Analysis (VIDA, Tiede et al. 2022a), uses parameterized templates to approximate an image and adjusts the parameter values until a specified cost function, in the form of a probability divergence that provides a distance metric between the image and template, is minimized. We use an updated version of the SymCosineRingwFloor template used in Sgr A^∗ 2017 IV to match the mF-ring model described later in Sect. 6.2. We set m = 2 (the azimuthal brightness mode) for consistency with the geometric modeling analyses (see Sect. 6.1). While REx is intended to work only on ring-like images, VIDA can support a variety of templates, which can be used to evaluate the success of our imaging methods by quantitatively reproducing the non-ring synthetic data sets such as the disk and the double source.

Figure 13 shows the ring parameters (diameter d, width w, and position angle) measured with REx and VIDA for all bands, days, and imaging pipelines for HOPS and CASA data reconstructions. On April 21, all bands and pipelines show ring-like structures with a roughly consistent diameter of ∼42 μas as expected from the image morphology discussed in Sect. 5.2. Position angle values are around 215° and are consistent among all bands and pipelines for April 21. Width measurements have relatively large differences among the imaging pipelines, especially between DIFMAP and the other methods. This is primarily due to the beam convolution effect in CLEAN imaging. Since measurements of the ring width are dependent on an individual method’s ability to leverage super-resolution, it is challenging to combine the width measurements across methods. While some methods trend toward slightly lower widths compared to their 2017 values (DIFMAP, eht-imaging, THEMIS), others are more consistent with the 2017 values (SMILI, Comrade). While the results from two independent IDFE methods, REx and VIDA, generally agree, there do exist some discrepancies. These discrepancies between the IDFE methods primarily arise due to the difference in the functional responses to an image with double rings or otherwise corrupted ring-like structures.

Tables 7–10 list the diameter, width, position angle, asymmetry, and fractional central brightness measured with REx and VIDA from all days, bands, a priori calibration, and imaging pipelines. The brightness asymmetry A is ∼0.2 − 0.4, indicating an asymmetric ring-like structure. Moreover, the fractional central brightness f_C is ∼0 − 0.5, indicating a clear central depression.

Comparing the extracted features between April 21 and April 25 is challenging due to the lack of Top Set images for DIFMAP, eht-imaging, and SMILI on April 25 bands 1 and 2 as well as the sparser (u, v) coverage on April 25 leading to higher uncertainty. For the available measurements, the features look consistent within their errors between days. Major discrepancies in the measurements, especially in the ring width, arise from the broken ring structures, as shown in Fig. 7. Ring parameters between HOPS and CASA agree well with each other, and we obtain a diameter of ∼42 μas and position angle of ∼215°. All parameter measurements except for the position angle are consistent between 2017 and 2018. The position angle has changed from ∼180° in 2017 to ∼215° in 2018. The interpretation of this change is discussed in Sect. 7.

In summary, the 2018 EHT data is well described by a ring-like crescent structure. Where the data quality and coverage are good, we find a consistent diameter compared to the 2017 images and a slightly shifted position angle. Where the (u, v) coverage is poor, such as in bands 1 and 2 on April 25, we expect our reconstruction methods to struggle, and indeed we have much larger uncertainties in the image structure for those data sets. Now that we have determined the presence of a ring-like structure in the 2018 data, in the next section we can compare and verify the measured image domain parameters by directly modeling those features with an assumed ring in the reconstructions.

6. Direct modeling of physical parameters

In addition to the image-domain feature extraction methods, we also attempted to model the physical ring parameters from the visibility domain information directly. We use a procedure similar to one described in M 87^∗ 2017 VI where we construct a family of geometric models which closely match the observed M 87^* visibilities. This procedure can be seen as another form of image reconstruction with stronger assumptions on the basic image structure. First, in Sect. 6.1, we assess and quantify the extent to which ring-like features are preferred by the data by comparing the Bayesian evidence of both ring-like and non-ring visibility-domain geometric models. We then describe our fiducial models for direct physical parameter extraction and analyze their inferred ring features in Sect. 6.2. Finally, in Sect. 6.3, we check for consistency in the ring parameters across models and compare the results to M 87^∗ 2017 VI.

6.1. Geometric modeling evidence exploration

We estimated the preference for ring-like models by comparing geometric models of varying degrees of complexity to the 2018 M 87^* EHT data. This is very similar to the procedures found in M 87^∗ 2017 VI and Sgr A^∗ 2017 IV to evaluate the performance of different models against the EHT data. In this instance, we used different numbers and combinations of Gaussians, disks, rings, crescents, and m-rings (see Johnson et al. 2020 and Sect. 4.3 in Sgr A^∗ 2017 IV) as our model library. We evaluate each model based on the difference in the log Bayesian evidence, Δlog(Z).

Our Bayesian evidence analysis was performed on data generated from closure quantities constructed with band 3 observations made on April 21 after applying an S/N >  3 cut. This means that these reconstructions are not susceptible to station-based corruption effects. The Bayesian evidences and posteriors for each model were evaluated using Comrade with a Nested Sampling (NS; Skilling 2006) posterior reconstruction scheme through dynesty (Speagle 2020).

The resulting Bayesian evidence of each model is shown in Fig. 14, where it is evident that ring-like models provide a much better fit to the data when compared to non-ring models. We select the best-performing crescent and m-ring models; the THEMIS xs-ringauss model from M 87^∗ 2017 VI (11 parameters) and a stretched m-ring of order 2 (8 parameters). We use these two models as the basis for defining fiducial models for direct feature extraction.

6.2. Direct geometric modeling

6.2.1. THEMIS direct modeling methods

We expect that if the object we see at the center of M 87 is a black hole (or sufficiently similar object), and the material surrounding it is optically thin, it should have a ring-like image feature associated with half-integer orbits (n) of photons around the central gravitational object. Thus, the n = 0 “ring” corresponds to direct emission from around the black hole, the n = 1 ring corresponds to photons that complete a half-orbit around the back of the black hole to impact our line of sight, and so on. This feature is consistently seen in analytic, semi-analytic, and GRMHD simulations of optically thin accretion flows around black holes, and directly modeling this feature can, in principle, provide strong constraints on space-time properties (Johnson et al. 2020; Wielgus 2021; Broderick et al. 2022b). Motivated by this strong physical prediction, we proceed to fit a thin geometric crescent model with an asymmetric azimuthal brightness distribution combined with the standard splined raster and large-scale Gaussian model described above, collectively referred to as a “Hybrid Themage”. A more detailed explanation of the model parameters and prior distributions is given in Appendix J.

There have been important discussions about how well the Hybrid Themage model recovers the properties of the n = 1 photon ring and by how much the recovered parameters are biased by contributions from the n = 0 ring (Tiede et al. 2022b). For this reason, we cannot assert that this thin geometric component represents the n = 1 photon ring. Any model with an explicit ring component is likely to fare well when applied to M 87^* EHT data with sufficient coverage. In a future paper, we intend to conduct a series of detailed synthetic data tests to better understand the behavior of the thin ring component.

For the Hybrid Themage model, we fit the same scan-averaged, network-calibrated data we use in the standard Themage raster fitting, along with the same priors on the station gain parameters. We used the same sampler and adaptive tempering scheme to fit the Hybrid Themage models we also used to fit the standard Themage models. This type of image model has consistently produced good fits when applied to M 87^* data, both in previous papers analyzing the 2017 data (Lockhart & Gralla 2022; Broderick et al. 2022a) and in the current data set as seen in Appendix K.

As mentioned above, we also repeat the same procedure done in M 87^∗ 2017 VI to directly model physical parameters in the visibility data. We used the THEMIS xs-ringauss model to provide direct continuity with the 2017 EHT analysis, which was also the best performing crescent model featured in Sect. 6.1. A detailed description of the construction of the xs-ringauss model is described in Sect. 5.1 of M 87^∗ 2017 VI, but we also briefly summarize the main components here.

The primary ring was constructed by subtracting a uniform disk offset from a larger uniform disk, creating a “top-hat” crescent, then “slashed” to create a linear brightness gradient. To quantify the depth of the central flux depression, we added a circular Gaussian to act as a “floor” in the crescent center. We also added an additional asymmetric Gaussian pinned to the inner edge of the crescent at the point where the width is largest, inspired by the model from Benkevitch et al. (2016). The overall “image” was then smoothed with a Gaussian kernel, and two additional “nuisance” asymmetric Gaussian components were introduced to help model diffuse emission near the ring. A cartoon describing each model component is shown in Fig. N.1.

The main difference between the Hybrid Themage and the xs-ringauss model construction is in the treatment of the diffuse emission. The xs-ringauss model used a pair of asymmetric Gaussians to help model any non-ring emission, and the Hybrid Themage used the adaptive raster components. The ring model in the Hybrid Themage was constructed using the same model object in the THEMIS code as the xs-ringauss, but with restricted priors to suppress the “pinned” Gaussian and enforce a small ring width. While the Hybrid Themage utilized a more flexible method to model the non-ring emission, the xs-ringauss was more agnostic to the precise structure of the ring.

We compare the xs-ringauss model to the scan-averaged network-calibrated data, similar to how we fit the raster Themage models. We utilized the same parallel tempered Hamiltonian MCMC sampler to explore the likelihood space for the xs-ringauss. We also fit the complex station gains in the same way as the raster models, using the same priors for the gain amplitudes and phases. We describe the detailed priors on the xs-ringauss parameters in Appendix N.

6.2.2. Comrade direct modeling method

We also define a fiducial model based on the best-performing m-ring model in Fig. 14 for analysis with Comrade. The fiducial model was taken to be an m-ring of order 2 that was allowed to have an asymmetric stretch. We included a flat emission floor in this model and constrained it to fill the interior of the m-ring. The resulting combination was convolved with a circular Gaussian. This construction defines a nine-parameter, elliptical, floored m-ring model, which we refer to hereafter as an “mF-ring”. An additional pair of elliptical nuisance Gaussians were also included to improve the fit quality. A summary of the mF-ring and its defining parameters is given in Fig. N.1, with additional details in Appendix N. We have also fit the fractional stretch of the mF-ring and the xs-ringauss compact flux density. We do not report a compact flux density for the mF-ring since we fit this model to closure quantities which are insensitive to the total flux.

6.3. Direct modeling results

Similar to the imaging methods in Sect. 5, each of these direct modeling methods fits the data well. We compare the model performance against the data in Fig. 15, where we show the complex visibility residuals for the Hybrid Themage and xs-ringauss, as well as the log closure amplitude and closure phase residuals for the mF-ring. The gray points represent the scan-averaged data, the colored points represent the maximum likelihood model for the THEMIS-based models, and the maximum a posteriori model for the mF-ring. We find that each model performs very well when compared against the data, with small, unstructured residuals. Appendix N also presents image-domain representations for each day and all observing bands, demonstrating consistency with the other imaging methods.

Table 11 gives a summary of the physical model parameters, and Fig. 16 shows the model parameters and their comparison with the same model parameters from 2017. The most noteworthy change from 2017 is that of the brightness asymmetry position angle ($\widehat{\varphi }$ in Table 11), which has shifted from 160° to 212° as fitted by the mF-ring, and 202° as fitted by the xs-ringauss. The Hybrid Themage finds a position angle of about 215°. The slight discrepancy between the three models is likely due to our models definition of position angle, which is ambivalent to any flux contribution from the nuisance Gaussians and raster components. We note that image domain fits of the best-fit xs-ringauss and mF-ring models produce consistent position angle extraction and thus conclude that the measured position angle from the mF-ring and the xs-ringauss model are consistent.

We also measure a slight increase in diameter from 43 μas as measured by the direct modeling methods in 2017 (M 87^∗ 2017 VI) to a median of 44.6 μas and 44.9 μas for the xs-ringauss and mF-ring models, respectively. The Hybrid Themage has a mean diameter even higher, at 45.5 μas, but the known inverse correlation between diameter and width for thin rings (M 87^∗ 2017 VI) can partially explain this. It is expected that most of the observed 1.3 mm emission originates from the inner portion of the accretion flow–in a region near the vicinity of the photon sphere (M 87^∗ 2017 V; Dexter et al. 2012; Emami et al. 2023)–though this emission morphology produces images that are slightly larger than the predicted size of the black hole shadow. Moreover, this emission is subjected to vary in size resulting from accretion flow physics (Tiede et al. 2022b). Section 5.3 of M 87^∗ 2017 VI gives an estimate of the theoretical scaling and uncertainties associated with inferring physical features from the fitted parameters of the geometric models.

The fractional stretch of the mF-ring allows for a measurement of image ellipticity which could originate from many possible features such as the shadow (Falcke et al. 2000; Cunningham & Bardeen 1973) or inner shadow (Chael et al. 2021). However, the ability to extract this feature from data can be complicated by gaps in the EHT (u, v) coverage. It thus may not be directly related to any intrinsic image geometry (Tiede et al. 2022c). The depth of the central brightness depression and the fraction width measurements in 2018 are consistent with their respective measurements from 2017.

7. Discussion

In this work we utilize the EHT observations from April 2018 to obtain a new image of M 87^*, building upon the findings from the 2017 EHT campaign (M 87^∗ 2017 I; M 87^∗ 2017 II; M 87^∗ 2017 III; M 87^∗ 2017 IV; M 87^∗ 2017 V; M 87^∗ 2017 VI; M 87^∗ 2017 VII; M 87^∗ 2017 VIII). As shown in Fig. 17, representative images from band 3 April 21 data produced by both the imaging pipelines and direct modeling methods reveal the existence of a ring-like emission structure around the central SMBH M 87^*. This result across various independent methods demonstrates the robustness of our conclusions. The consistency of a ring-like emission structure on a 1-year time scale supports the interpretation of this structure as the black hole shadow of M 87^* (M 87^∗ 2017 I; M 87^∗ 2017 VI), consistent with the predictions from general relativity.

In this section, we discuss the main implications from the new images in 2018. We start by discussing the significance of the new diameter estimate and compare it to the previous diameter and mass estimate from the analysis of the 2017 data. We do not produce a new mass estimate, since this requires a dedicated calibration exercise using GRMHD simulations (M 87^∗ 2017 VI; Sgr A^∗ 2017 IV; Sgr A^∗ 2017 VI), which will be featured in a follow-up paper. Instead we use the scaling factor from the 2017 analysis to construct a proxy for the gravitational diameter (M 87^∗ 2017 VI), allowing us to compare the mass estimate from 2017 with the diameter measurements from 2018. Then we discuss the significance of the counter-clockwise shift in the position angle of the brightness asymmetry in the ring from 2017 to 2018. Finally, we investigate the challenges related to estimating the compact flux density of M 87^* on horizon scales with the sparse (u, v) coverage of the EHT.

7.1. Consistency of the ring diameter

The asymmetric ring’s estimated characteristic parameters are provided in Tables 7–10 via image domain feature extraction and Table 11 via the direct modeling results. The parameters are broadly consistent across different calibration schemes and reconstruction methods. There is a tendency for the ring diameters extracted with the image domain feature extraction to be slightly smaller than those from the direct modeling methods. This tendency is also observed in the 2017 results and is related to the effective resolution of each method (M 87^∗ 2017 VI).

We combine the results from image domain feature extraction and direct modeling to construct a median diameter across all methods (Tables 7–11). We first construct a diameter histogram for each method, normalize the histograms to equally weigh each method, then combine the normalized histograms to produce a single diameter histogram across all imaging and direct modeling methods. We present the resulting median diameter and 68% intervals in the last column of Table 12. This lets us more easily compare the ring diameter of 2018 EHT campaign with that of 2017 (see Table 1 of M 87^∗ 2017 I). We find the difference in median diameter between these two campaigns is on the order of only ∼1 μas, with error bars that almost completely overlap, indicating a persistent emission structure.

Since the diameter estimates are model-dependent, we can instead convert the diameter to a common physical scale. For black holes, a natural quantity is the angular gravitational radius θ_g. We can construct this quantity by using the diameter values (d) obtained via image reconstruction and direct modeling along with a scaling factor (α) related to GRMHD simulations:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{g}}=\frac{d}{\alpha }·\end{array}$$(5)

Since we know the true angular gravitational radius in the GRMHD simulations, we can produce image reconstructions of GRMHD synthetic data and directly determine the scaling factor between the reconstructed diameter and the true θ_g. This scaling factor is a function of the GRMHD black hole spin, inclination, and temperature ratio for the electron-to-proton coupling in both the weakly and strongly magnetized regions. The scaling factor also depends on the particular image reconstruction model, and so in principle it should be unique for each method.

In our previous 2017 analysis, we derived α values between 10.7 and 11.5 with associated errors of ∼10% across all methods, from imaging to direct modeling (M 87^∗ 2017 I; M 87^∗ 2017 VI). Because the α value has some dependence on the observational properties such as coverage and S/N, and because we use some methods in 2018 that were not used in 2017, it is necessary to re-calibrate the 2018 methods to determine α values for a proper 2018 θ_g estimate. This would require image and model reconstructions of a large number of GRMHD synthetic data sets. This is a significant undertaking, and since this paper is primarily focused on presenting the first 2018 images, it is beyond the scope of this work to produce a 2018 θ_g calibration.

Nevertheless, it is still valuable to produce some comparison with the 2017 results, and to that end we construct a proxy for θ_g, d/α₂₀₁₇, which combines the diameters we measure in 2018 with the α values from 2017. For the 2018 methods which have a 2017 θ_g calibration, we will use the corresponding 2017 α values from M 87^∗ 2017 VI, either in Table 4 for the xs-ringauss or Table 6 for DIFMAP, eht-imaging, and SMILI. For the 2018 methods which do not have a corresponding 2017 θ_g calibration, we will adopt α = 11.0, which is the median value across all 2017 methods (see Table 1 in M 87^∗ 2017 I).

We show a comparison between the 2018 April 21 band 3 d/α₂₀₁₇ and 2017 θ_g values in Fig. 18 for each method. We show the 2018 points with colored diamonds and the 2017 points with black squares. Red (REx) and blue VIDA points correspond to the diameters estimated via image domain feature extraction, and the pink, purple, and dark blue points correspond to the Hybrid Themage, xs-ringauss, and mF-ring direct modeling methods. The uncertainty in the 2018 points is directly related to the 68% confidence intervals of the April 21 band 3 diameter estimates (see Tables 7 and 11). For the 2017 points, the uncertainty related to the 2017 observational systematics are slightly larger than the measurement uncertainty (see Sect. 7 of M 87^∗ 2017 VI), and is represented by the solid black error bars. The uncertainty related to the diversity of GMRHD models used in the calibration set dominates over all other uncertainties, and is shown by the dashed gray error bars.

While the diameter uncertainties for some of the 2018 reconstructions are larger than the observational uncertainties in 2017, all methods are broadly consistent with each other, and all live within in the GRMHD uncertainty. This supports our previous interpretation that the ring-like structure in our images is consistent with the shadow of a supermassive black hole with a mass of ∼6.5 × 10⁹ M_⊙. Any discrepancies may disappear in a proper 2018 θ_g calibration.

We note that there are at least two new integral field spectroscopic stellar kinematics mass estimates for M 87^* recently published. One of these adapts new triaxial orbit models to Keck II telescope observations (Liepold et al. 2023) and favors a mass for M 87^* of 5.3 × 10⁹ M_⊙. The other estimate uses adaptive optics observations from the VLT instruments MUSE and OASIS (Simon et al. 2024) and favors a mass of 8.7 × 10⁹ M_⊙. However, when considering a different stellar mass profile in the inner region (and thus a different mass-to-light ratio), this method obtains a mass estimate of 5.5 × 10⁹ M_⊙. These ∼5 × 10⁹ M_⊙ mass estimates are within 1.5σ of the 2017 EHT mass estimate. Arguably, there are significant systematic uncertainties in the non-horizon scale mass estimates of M 87^* that require continued investigation.

7.2. Position angle variation

It is worth noting that there is a significant counter-clockwise change in the position angle of the brighter region in the ring-like structure from 2017 to 2018 (∼30°). The shift direction is consistent with the prediction reported in M 87^∗ 2017 V if the alignment of the large-scale jet is normal to the disk.

The work by Wielgus et al. (2020) shows evidence of year-scale position angle variation of the ring’s peak brightness, based on simple modeling of prototype EHT data. The present study is the first case where such a year-scale variation is unambiguously confirmed in the image domain (see Fig. 19). Furthermore, recent long-term monitoring studies of M 87^* using longer-wavelength VLBI found a systematic position angle oscillation of the parsec-scale jet (Walker et al. 2018; Cui et al. 2023), which could be caused by flow instabilities (Walker et al. 2018) or precession of the central compact source (Cui et al. 2023). M 87^* exhibits year-scale morphological variations and a counter-clockwise shift from 2017 to 2018 at μas scales. Nevertheless, whether the time variation we see in M 87^* and its jet are physically linked with each other or not, there could be some biases, given the large spatial gap between them and also the lack of inter-year EHT images. Further accumulation of EHT images over forthcoming years, along with parsec-scale jet monitoring, is required to understand better the origin of the year-scale variation of the ring-like structure and its possible connection to the large-scale jet.

In addition to the year-scale position angle variation of the brightest spot in the asymmetric ring (Wielgus et al. 2020), GRMHD simulations for the accretion environment around M 87^* also show that the position angle of the brightest location of the asymmetric ring can vary due to the turbulent, magnetized accretion environment, with a time-scale much smaller than the observational cadence between 2017 and 2018 EHT observations (see also Fig. 6 of M 87^∗ 2017 V). As will be discussed further in a forthcoming paper, it is possible to apply the 2017 and 2018 EHT observation results as independent constraints for GRMHD models of the black hole in M 87.

Figure 20 presents the position angle distribution of the fitted jet directions, assuming the black hole spin vector is pointing away from Earth. To this end, from the image library applied in M 87^∗ 2017 V and M 87^∗ 2017 VIII, we select a group of model images which have the black hole spin axis pointing away from Earth, then fit all the selected images to the 2017 (April 6, high-band) and 2018 (April 21, band 4) EHT data by applying the snapshot model (SSM) scoring procedure introduced in M 87^∗ 2017 V and M 87^∗ 2017 VI. For a summary of the model image library preparation in M 87^∗ 2017 V and M 87^∗ 2017 VIII), the one-fluid GRMHD simulations are initialized with a weakly magnetized torus of plasma orbiting in the equatorial plane of the black hole, and the rotational axis of the torus is aligned with the black hole axis. After the system evolves to a steady state, we perform radiative transfer to compute the thermal synchrotron emission to generate images with the parameterized electron thermal dynamics (see M 87^∗ 2017 V, for details). The selected models were prepared by the code iharm (Gammie et al. 2003), and the model images including different black hole spins³, accretion flow types, and assumptions for the ratio between the electron and ion temperature in the simulation (see M 87^∗ 2017 V, for details).

Constrained by the 2018 EHT observation of M 87^*, the distribution of the forward jet direction in the model images is consistent with the observed milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023), which is ∼288° counter-clockwise from the north. This supports the interpretation that the black hole spin vector points away from Earth, which was also strongly favored by the first M 87^* EHT results (M 87^∗ 2017 V). Additionally, with the model images considered, the position angle distribution of the forward jet direction in the 2018 observations has its peak closer to ∼288° (vertical dashed line in Fig. 20). Assuming the large-scale jet is aligned perpendicular to the accretion disk (M 87^∗ 2017 V), with many yearly observations to accumulate a statistically significant number of independent images, we expect to see the position angle of the brightest region to be more similar to that seen in 2018 than in 2017. Since the timescale for position angle shifts in GRMHD simulations is small compared to the yearly cadence of EHT observations, we do not necessarily associate this counter-clockwise shift in the brightness position angle with a global rotation of the accretion flow, but instead consider the 2018 observations to be a snapshot of the most common orientation of the accretion flow. Additionally, we will also directly score a new GRMHD library against the 2018 data using similar procedures to the ones found in M 87^∗ 2017 V; M 87^∗ 2017 VIII; Sgr A^∗ 2017 V, and leverage the 2017 and 2018 data to score the GRMHD simulations across multiple years.

7.3. Compact flux density estimates from Bayesian imaging methods

The EHT has very few baselines that probe the intermediate region in (u, v) space between the intra-sites and about 2 Gλ. The only baseline that probe this region of (u, v) space is LMT–SMT. Even though the LMT has a very sensitive 50 m dish, it typically experiences relatively large pointing errors, increasing the recovered station gain amplitudes. While the LMT seems to have much better pointing in 2018 than in 2017, we have about half the scans in 2018 compared to 2017 and no LMT participation on April 25. This limited data means it is much more challenging to constrain the compact flux density in M 87^* during the 2018 EHT observation campaign. The longer-wavelength and zero-baseline analysis described in Sect. 4 are consistent with the findings from the Bayesian methods: there is a large amount of uncertainty in what we should expect the compact flux density to be at EHT-scales.

To illustrate this, we investigate the compact flux density estimates produced by the various Bayesian methods that fit the complex visibilities (THEMIS) or the visibility amplitudes (Comrade). For the THEMIS raster and Hybrid Themage models, we quantify the compact flux density by summing the flux density from the raster and the raster plus crescent, respectively. For Comrade, the compact flux density is a model parameter which sets the total flux density of the raster as described in Sect. 5.1.5. The THEMIS xs-ringauss model constructs the compact flux density estimate by summing the flux density from the crescent and nuisance Gaussians.

We summarize the compact flux density estimates across all observing bands on April 21 from the THEMIS raster, Hybrid Themage, Comrade raster, and THEMIS xs-ringauss models in Fig. 21, using the HOPS data. There are two clear modes: a “low-estimate” for compact flux density around 0.6 Jy, and a “high-estimate” for compact flux density around 1.0 Jy. The THEMIS raster model is the only model that consistently produces compact flux density estimates below 0.7 Jy for all bands. In comparison, the THEMIS xs-ringauss and Comrade raster produce compact flux density estimates that are greater than 0.9 Jy. The Comrade raster prefers larger compact flux density because of the larger FOV and the choice of prior distribution. We found that the compact flux density strongly depends on the Comrade prior assumptions because of the lack of intermediate baselines. The Hybrid Themage model prefers lower compact flux densities in band 3 and band 4, but a higher compact flux density in band 1, with uncertainty that spans 0.5 Jy. The error bars in Fig. 21 are taken from the 68% confidence intervals and that, generally, the methods that produce low compact flux density estimates do not overlap with the high compact flux density estimates.

There is no consensus among the Bayesian methods about what the compact flux density should be in the 2018 EHT data. The chain parameters for the Hybrid Themage reconstruction of the April 21 band 1 data give us clues about how we might achieve such high values for compact flux density. In Fig. 22, we show the total compact flux density, the raster FOV, and the crescent flux density. Recall that for the Hybrid Themage, we construct the compact flux density by adding the crescent flux density and the flux density from all the raster points. We see clear bi-modal distributions in the compact flux density, FOV_y, and to a lesser degree, in the crescent flux density. This suggests that, at least for this reconstruction, adding the entire raster flux density may not be appropriate for computing the compact flux density; the raster is elongated in one direction and places some flux outside of what we might consider the compact region. This is evident in the April 21 band 1 Hybrid Themage image in Fig. N.2, where the raster extends flux in the northwest direction.

The most crucial parameters of interest, such as the diameter and the position angle, do not seem to depend strongly on the recovered compact flux density. This is consistent with the findings from the RML and CLEAN methods, which are not able to independently constrain the ratio between the total flux density and the compact flux density. The measured diameters and position angles are entirely consistent between models that produce different values of compact flux density (see Appendix M).

8. Conclusion

We present the results of the 2018 EHT observations of M 87^* at 1.3 mm in order to further investigate the nature of the black hole shadow, and the year to year variation in the image structure. During the 2018 observations, the GLT participated for the first time as part of the VLBI array and provided additional baselines to improve the (u, v) coverage. Following the example of the 2017 data analysis, we use multiple independent calibration, imaging, and modeling methods to analyze and interpret the 2018 EHT data. We analyze data on two independent epochs (April 21 and April 25) across four independent frequency bands. The 2018 data show a clear null in the visibility amplitudes at 3.4 Gλ, and exhibit significant variation in the time-dependent closure phases compared to 2017. The differences in the closure phases indicate a significant change in the asymmetric structure from 2017 to 2018.

All of the imaging and modeling methods indicate the presence of a ring-like structure, and the extracted ring characteristics are consistent across all bands and days. Similar to the 2017 data, the current EHT data does not provide strong constraints on extended or diffuse emission outside the ring. All methods on April 21 data suggest a ring diameter of ∼42 μas for the image domain analyses (Sect. 5.3) and ∼45 μas for the direct modeling methods (Sect. 6.3). When combining all methods with equal weights, the median diameter of the ring in 2018 is $43.3_{-3.1}^{+1.5}\mathrm{\mu }$as, in good agreement with the diameter measured in the analysis of the 2017 data.

The lack of variability in the ring diameter between 2017 and 2018 is consistent with the predictions from general relativity for strongly lensed emission around a black hole. In contrast with the 2017 data, we detect a significant shift in the position angle of the ring brightness asymmetry from ∼180° as measured in 2017 to ∼210° in the 2018 data. This shift in position angle is consistent with the expected variability from GRMHD simulations. In particular, when converting the brightness asymmetry position angle to a nominal black hole spin direction, the 2018 image is more consistent with the orientation of the large-scale jet seen at longer wavelengths.

---

Acknowledgments

The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19_058 (TITANs), Fondecyt 1221421 and BASAL FB210003; the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the ALMA North America Development Fund; the Astrophysics and High Energy Physics programme by MCIN (with funding from European Union NextGenerationEU, PRTR-C17I1); the Black Hole Initiative, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation (although the opinions expressed in this work are those of the author(s) 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; the Fondo CAS-ANID folio CAS220010; 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 Horizon ERC Grants 2021 programme under grant agreement No. 101040021; the Generalitat Valenciana (grants APOSTD/2018/177 and ASFAE/2022/018) and GenT Program (project CIDEGENT/2018/021); MICINN Research Project PID2019-108995GB-C22; 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 FAPESP (Fundação de Amparo á Pesquisa do Estado de São Paulo) under grant 2021/01183-8; 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 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 National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-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, OISE-1743747, AST-1816420, AST-1935980, AST-2034306, 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); 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); and the YCAA Prize Postdoctoral Fellowship. We thank the staff at the participating observatories, correlation centers, and institutions for their enthusiastic support. This paper makes use of the following ALMA data: ADS/JAO.ALMA#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 US 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.

References

Appendix A: Flux density calibration parameters and uncertainties

The visibility amplitudes need to be converted from correlation coefficients to units of flux density, considering the different telescope sensitivities across the array. This is achieved by multiplying the normalized visibility amplitudes with the geometric mean of the derived system equivalent flux densities (SEFDs) of the two stations of each respective baseline. The SEFD is given by (e.g., M 87^∗ 2017 III):

$$\begin{array}{c}\hfill \mathrm{SEFD}=\frac{1}{{\eta }_{\mathrm{ph}}}\frac{T_{\mathrm{sys}}^{\ast }}{\mathrm{DPFU}\times g_E},\end{array}$$(A.1)

where $T_{\mathrm{sys}}^{\ast }$ is the effective system noise temperature, corrected for atmospheric attenuation. The degrees per flux density unit, DPFU, is a conversion factor from units of temperature (K) to units of flux density (Jy), and g_E is the normalized elevation-dependent gain curve of the telescope. η_ph is the time-dependent phasing efficiency for a phased array, and has a value of unity for single dish telescopes. The DPFU can be estimated from the antenna aperture efficiency (η_A) as:

$$\begin{array}{c}\hfill \mathrm{DPFU}=\frac{{\eta }_{\mathrm{A}}{A}_{\mathrm{geom}}}{2k},\end{array}$$(A.2)

where A_geom is the geometrical area of the antenna dish in units of m², and k = 1.38 × 10³ is the Boltzmann constant in units of Jy m²/K.

A summary of the station DPFUs and gain curve parameters for the EHT array in 2018 is provided in Table A.1, together with their estimated uncertainties. Details on their derivation for each station are described by Koay et al. (2023a). We expect the DPFU uncertainties to be the dominant source of systematic errors in the flux density calibration of EHT data. $T_{\mathrm{sys}}^{\ast }$ values are not expected to vary much within the typical scan duration of a few minutes, so their uncertainties are expected to be no larger than a few percent. The GLT is an exception, due to sub-optimal performance in 2018 when it was still only partially commissioned, where $T_{\mathrm{sys}}^{\ast }$ uncertainties are also significant (Koay et al. 2023b).

We note in particular that the quoted errors for the GLT and LMT in Table A.1 are very likely to be underestimated, due to insufficient a priori data collected by the stations, as described in Appendix D. We therefore obtain an independent set of constraints on the amplitude calibration uncertainties for these two stations in Appendix E.

Appendix B: Non-closing systematic error budget

The closure phases (ψ_C), log closure amplitudes (lnA_C), closure trace phases and closure trace log amplitudes 𝒯, the latter two of which are novel diagnostics described in detail by Broderick & Pesce (2020), are independent of station-based gain errors, but are sensitive to systematic non-closing systematic errors. Figure B.1 shows the normalized distributions of the differences in the closure quantities between bands (band 1 − band 2, band 3 − band 4), polarizations (RR−LL), as well as the distributions of the trivial closure quantities derived from triangles or quadrangles with at least two co-located stations.

In an ideal scenario where polarimetric leakages are absent, the differences in closure quantities between bands and polarizations, as well as the trivial closure quantities, is expected to follow a normal distribution with a mean value of zero. The presence of non-closing systematic errors can lead to deviations from the expected normal distribution, as can be seen for the blue histograms in Fig. B.1.

The non-closing systematic error, s, can thus be determined by solving for its value while enforcing the condition of:

$$\begin{array}{c}\hfill {\mathrm{mad}}_0\left(\frac{X}{\sigma }\right)={\mathrm{mad}}_0\left(\frac{X}{\sqrt{{\sigma }_{\mathrm{th}}^2+s^2}}\right)=1,\end{array}$$(B.1)

where σ is the total uncertainty associated with the random variable X, and σ_th is the known a priori thermal uncertainty. In this construction, X is a random variable representing the closure products being examined, for instance band 1 − band 2 closure phase difference, RR−LL log closure amplitude difference, etc. mad₀ is the modified median absolute deviation given by:

$$\begin{array}{c}\hfill {\mathrm{mad}}_0(Y)=1.4826\mathrm{med}(|Y|),\end{array}$$(B.2)

where med denotes the median, and the factor 1.4826 scales the result so that it acts as a robust estimator of standard deviation for a normally distributed random variable Y with zero mean.

The non-closing systematic errors estimated for M87^* and 3C 279 using the various tests are presented in Table B.1. To avoid biases arising from low S/N, only closure quantities above the S/N > 6 threshold are used.

Appendix C: Temporal coherence

We quantify the level of decoherence loss of the data calibrated using the HOPS and CASA pipelines using the ratio A_scan/A_2s for each scan (M 87^∗ 2017 III). A_scan is the debiased amplitude of the coherently averaged visibilities over the entire scan, and A_2s is the amplitude derived from data coherently averaged over 2s, then incoherently averaged over the entire scan.

Figure C.1 shows the cumulative distribution of the quantity A_scan/A_2s for 2342 unique scans (per band, polarization, and baseline) of the M87^* April 21 and 25 data with S/N >  7, common to both the HOPS and CASA pipelines. We find that 92.9% of HOPS and 85.6% of CASA scan-averaged data show decoherence losses of less than 10%. When selecting April 21 data only, 97.4% and 87.9% of HOPS and CASA scan-averaged data show better than 90% coherence. For April 25, 86.3% and 82.3% of HOPS and CASA scan-averaged data show better than 90% coherence.

The level of coherence in the 2018 data is slightly lower than that of 2017 M 87^∗ 2017 III, due to (i) overall poorer weather at some sites, including Chile, Hawai‘i, and at the SMT (especially on April 25), as well as (ii) the poor sensitivity levels of the GLT, as described in Appendix E. Overall, CASA data have lower coherence levels than HOPS data, mostly for baselines with lower S/N such as those including the GLT. For example, when we exclude GLT baselines, the percentage of CASA scan-averaged data showing decoherence losses of less than 10% increases by a few percent, and becomes comparable to that of HOPS.

Appendix D: Data issues

In this appendix, we provide details on data issues that were identified in the 2018 L1 data package pertaining to M87^* (and 3C 279), and steps taken to address these issues. Our analysis indicated that the flux densities on APEX baselines are underestimated after a priori amplitude calibration, when compared to the flux densities of ALMA baselines to other stations. We also found 25% differences between RCP and LCP amplitudes on all APEX baselines in band 1, and a smaller but still significant 13% difference between the RCP and LCP amplitudes in band 4. Follow-up diagnostics found this issue to be caused by low power levels of the intermediate frequency (IF) signal and too large block-down converter attenuation settings for the ROACH2 Digital Back End (R2DBE; Vertatschitsch et al. 2015) associated with RCP in bands 1 and 4. These issues are corrected at the network calibration stage. There are also large amplitude drops on APEX baselines after UTC 01:00 on April 25, found to be caused by 180° peak-to-peak phase modulations whose origin is not yet determined to date. We therefore flagged and excised all the affected APEX data on April 25. The APEX baselines also show large amplitude jumps in several scans, typically at the beginning and end of a scan, which we also excised during post-processing.

While the above problematic L1 data have either been corrected or excised during post-processing, there remain a number of issues in the data used in our analyses. The GLT was still under commissioning and participated on a best effort basis in 2018. Due to astigmatism and a low ∼22% aperture efficiency, the visibilities on GLT baselines have lower S/N compared to that of the other stations. The telescope sensitivity was also not well characterized in 2018, so large uncertainties in the amplitude calibration are expected (Koay et al. 2023b). Due to the low telescope sensitivity, pointing sources could not be tracked down to elevations of < 10°, resulting in a poorly constrained pointing model at those low elevations and additional amplitude losses for sources observed at those elevations. This does not significantly affect M87^* which was observed at higher elevations, but severely compromises the quality of GLT data for 3C 279 which was observed between 6° to 8°.

The abrupt halt of the LMT observations in the middle of the 2018 campaign meant that there were no antenna sensitivity measurements recorded at the time. The amplitude calibration applied uses sensitivity measurements obtained in 2020, based on the assumption that the antenna characteristics have not changed in the intervening period. Also there were focusing and pointing issues being corrected on the fly while observing M87^* and 3C 279 on April 21, possibly leading to larger amplitude calibration errors. No elevation dependent gain curves were applied owing to a lack of aperture efficiency measurements at low elevations, where large gain errors of 20% (or possibly higher) are expected. This is not expected to significantly affect the M87^* and 3C 279 observations, since these targets were observed at high elevations of 50° to 77°.

The SMA phasing efficiencies on April 25 are poor (mostly below 50% between UTC 03:00 and 05:30) and have a large scatter, resulting in a large scatter in amplitudes on SMA baselines. The PV station also experienced poor weather on April 25 after UTC 02:00, leading to a large scatter in amplitudes and phases on associated baselines. There are large drops in amplitudes on all ALMA baselines on April 21 between UTC 04:12 to 04:17, attributed to phasing issues on ALMA. These are partly corrected during a priori amplitude calibration and network calibration, but a residual 40% scatter or drop in amplitude remains on the corresponding scan, which can be fixed during model-dependent self-calibration downstream.

Appendix E: Additional LMT and GLT gain constraints from overlapping visibilities

As mentioned in Appendix D, the GLT and LMT station sensitivities are not well constrained (see also Appendix A), and their residual gain uncertainties are expected to be larger than the amplitude calibration uncertainties provided in Table A.1 in Appendix A. To obtain additional a priori constraints on the amplitude calibration uncertainties for baselines to these two stations, we examine and compare amplitudes on the baselines that overlap with that of the LMT and GLT in (u, v) space.

The ALMA–LMT and ALMA–SPT baselines overlap for the calibrator source 3C 279, (Fig. E.1 left panel). Assuming that the SPT sensitivity is well-characterized and that the amplitude calibration uncertainties are relatively small in comparison, we find the amplitudes on LMT baselines to be under-calibrated by ∼34% in bands 1 and 2, ∼23% in band 3 (Fig. E.1, right panel), and ∼16% in band 4.

For the GLT, there are no overlapping baselines for M87^*. Amplitudes on GLT baselines for 3C 279 cannot be used to constrain its amplitude uncertainties for M87^* due to the large pointing offsets (see Appendix D). The closest baselines in (u, v) space where the amplitudes can be compared for M87^* are between ALMA–LMT and GLT–PV, where we find the amplitudes on the GLT–PV baselines are lower than ALMA–LMT amplitudes by up to ∼50% in bands 3 and 4. Considering that the LMT itself is under-calibrated, the amplitude errors for GLT are likely to be larger than 50%.

These a priori constraints are very crude, zeroth order estimates, but provide useful consistency checks with downstream station-based gains derived from model-dependent calibration, that is self-calibration during the imaging analysis. They are also important to better characterize the station uncertainties used in the generation of the synthetic data in Appendix G, so that they are reasonably consistent with that of the actual M87^* data.

Appendix F: Derivation of constraints on the total compact flux density and source size

We estimate constraints on the total compact flux density by using the EHT data themselves, and by comparing the flux density from quasi-simultaneous multi-wavelength observations. These two methods are similar to those used in M 87^∗ 2017 IV and Sgr A^∗ 2017 II. We provide the details below.

F.1. Constraints from EHT observations

Here, we derive constraints on the total compact flux density and the source size. We use the procedure described in Appendix B.1 of M 87^∗ 2017 IV, which was updated in Sgr A^∗ 2017 II to incorporate our uncertainties for the antenna gains using the DPFU uncertainties. This procedure included four constraints for the previous M87^* analysis. However, the 2018 data do not contain crossing points in the (u, v) space for the SMT–PV and LMT–PV baselines; therefore we cannot apply the 4th constraint. Moreover, only the observations on April 21 and 22 contain the short intra-site baseline, SMT–LMT, which is necessary for this analysis. Thus we can derive the three constraints summarized below only for the first two observing days. For reference, the total flux density F_tot estimated from the ALMA interferometric data is 1.13 Jy for both April 21 and 22. This value is an average of the total flux values derived for each of the bands 1–4.

Constraint 1. The first constraint is based on the fact that the visibility amplitudes on short baselines can be approximated by a circular Gaussian visibility function,

$$\begin{array}{c}\hfill V_{\mathrm{G}}(\mathbf{u};V_0,\theta )=V_{0}exp\left(-\frac{{(\pi \theta |\mathbf{u}|)}^2}{4ln2}\right),\end{array}$$(F.1)

where V₀ is the total flux density of the Gaussian source, |u| is the length of the baseline, and θ is its FWHM in radians. The intermediate-to-long baselines tend to measure larger correlated flux density than what is expected for an equivalent Gaussian source.

We can deduce that the measured amplitude ratio of the ALMA–LMT over SMT–LMT baselines will be larger than the corresponding ratio from a circular Gaussian source model. Consequently, the FWHM size of a circular Gaussian determined by the amplitude ratio between SMT–LMT and ALMA–LMT baselines provides an estimate of the minimum compact source size θ_cpct that is not significantly affected by the intrinsic fine-scale source structure:

$$\begin{array}{c}\hfill {\theta }_{\mathrm{cpct}}\gtrsim \sqrt{\frac{4ln2ln\left(\frac{|{\mathcal{V}}_{\mathrm{SMT}-\mathrm{LMT}}|}{|{\mathcal{V}}_{\mathrm{ALMA}-\mathrm{LMT}}|}\right)}{{\pi }^2{(|\mathbf{u}|}_{\mathrm{ALMA}-\mathrm{LMT}}^2-{|\mathbf{u}|}_{\mathrm{SMT}-\mathrm{LMT}}^2)}}.\end{array}$$(F.2)

Here, 𝒱_i − j denotes the true visibility on the baseline i − j. Considering the gain uncertainties, the ratio of the true visibility amplitudes is lower-bounded by

$$\begin{array}{c}\hfill \frac{|{\mathcal{V}}_{\mathrm{SMT}-\mathrm{LMT}}|}{|{\mathcal{V}}_{\mathrm{ALMA}-\mathrm{LMT}}|}\ge (1-\sqrt{\mathrm{\Delta }{g}_{\mathrm{ALMA}}^2+\mathrm{\Delta }{g}_{\mathrm{SMT}}^2})\frac{|V_{\mathrm{SMT}-\mathrm{LMT}}|}{|V_{\mathrm{ALMA}-\mathrm{LMT}}|},\end{array}$$(F.3)

where V_i − j and Δg_i denote the calibrated measured visibility and the gain uncertainty, respectively. We take the median of the visibility amplitude ratio from the collection of constraints derived for each single VLBI scan to derive a robust estimate of the minimum source size.

Constraint 2. The second constraint comes from the curvature of visibility amplitudes between the intra-site baselines and the LMT–SMT baseline. Since the compact flux density should not exceed the total flux density measured with the intra-site baselines, the amplitude drop from the intra-site to LMT–SMT baselines gives the maximum limit of that from the compact flux density to LMT–SMT baseline. Therefore, it gives the maximum limit of the source FWHM size with an equivalent circular Gaussian as

$$\begin{array}{c}\hfill {\theta }_{\mathrm{cpct}}\le \frac{2\sqrt{ln2}}{\pi |\mathbf{u}|}\sqrt{ln\frac{F_{\mathrm{tot}}}{|{\mathcal{V}}_{\mathrm{LMT}-\mathrm{SMT}}|}}.\end{array}$$(F.4)

When we consider the gain uncertainties,

$$\begin{array}{c}\hfill \frac{1}{|{\mathcal{V}}_{\mathrm{LMT}-\mathrm{SMT}}|}\le \frac{1+\sqrt{\mathrm{\Delta }{g}_{\mathrm{LMT}}^2+\mathrm{\Delta }{g}_{\mathrm{SMT}}^2+\mathrm{\Delta }{g}_{\mathrm{tot}}^2}}{|V_{\mathrm{LMT}-\mathrm{SMT}}|}.\end{array}$$(F.5)

We assumed the uncertainty of the total flux Δg_tot = 0.1. Similarly to Constraint 1, the median of the ratio is adopted to mitigate the effects of statistical errors.

Constraint 3. The minimum compact flux density can be derived by the maximum amplitudes of LMT–SMT baseline, since the visibility amplitudes are maximum at zero baseline length and the a priori calibration should not underestimate the station sensitivity. Therefore the compact flux density is constrained as

$$\begin{array}{c}\hfill F_{\mathrm{cpct}}\ge (1-\sqrt{\mathrm{\Delta }{g}_{\mathrm{LMT}}^2+\mathrm{\Delta }{g}_{\mathrm{SMT}}^2})|V_{\mathrm{LMT}-\mathrm{SMT}}|.\end{array}$$(F.6)

The equivalent circular Gaussian gives a stronger constraint with the minimum source size (Constraint 1) extrapolated from the LMT–SMT baselines, described by,

$$\begin{array}{c}\hfill F_{\mathrm{cpct}}\ge |{\mathcal{V}}_{\mathrm{LMT}-\mathrm{SMT}}|exp\left(\frac{(\pi |\mathbf{u}|{\theta }_{\mathrm{cpct}}{)}^2}{4ln2}\right).\end{array}$$(F.7)

These constraints are validated using the various synthetic models described in Appendix G. The upper and lower limits on the source size obtained from Constraints 1 and 2 are along the directions of LMT–SMT and LMT–ALMA baselines.

In summary, Constraints 1−3 lead to the conclusion that 0.30 Jy ≤F_cpct ≤ 1.13 Jy on April 21 and 22. If we use Eq. F.7, the minimum compact flux is 0.38 Jy. The compact source size has an equivalent Gaussian FWHM ranging between 39 and 98 μas on the first two observing days.

F.2. Constraints from quasi-simultaneous multi-wavelength observations

Supplementary, we also make use of quasi-simultaneous VLBI data at longer wavelengths obtained through part of a large MWL observing campaign coordinated with the EHT 2018 observations. Our approach assumes that the total compact emission of M87^* at ν ≈ 230 GHz (λ ≈ 1.3 mm) F_tot, 230 originates from the combination of (i) a compact emitter with unknown flux density F_cpct, 230 and spectral index at 230 GHz, and (ii) a large and diffuse (milliarcsecond-scale) jet, which can be characterized by its total flux density F_jet, 230 = F_tot, 230 − F_cpct, 230 and a steep spectral index of α_jet ∼ −(0.7 − 1.0), as measured at cm- and mm-wavelengths (with typical errors of ∼0.3; see Hovatta et al. 2014; Hada et al. 2016). Then, by measuring F_jet and F_tot at lower frequencies, and extrapolating them to 230 GHz using a power-law model, we can estimate the total compact flux density F_cpct at 230 GHz: F_cpct, 230 = F_tot, 230 − F_jet, 230.

We utilize the data obtained with the East Asian VLBI Network (EAVN; An et al. 2018; Cui et al. 2021; Cho et al. 2022) at 22 and 43 GHz, the Korean VLBI Network (KVN; Lee et al. 2016) at 43, 86, and 129 GHz, and the Global Millimeter VLBI Array (GMVA)+ALMA at 86 GHz. All of these observations were conducted between 2018 March 25 and April 21. For details of data reduction and calibration, see Cho et al. (2017, 2021) for EAVN, Lee et al. (2016), Kim et al. (2018b) for KVN, and Lu et al. (2023) for GMVA+ALMA, respectively.

The total flux densities of the nuclear region (i.e., F_tot = F_jet + F_cpct) at 22, 43, and 86 GHz are measured from the EAVN and KVN images by integrating the flux densities over a large window of 500 × 500 μas² centered on the radio core. We obtain total flux densities of F_tot, 22 ≈ 1.1 Jy, F_tot, 43 ≈ 0.9 − 1.1 Jy, F_tot, 86 ≈ 0.9 Jy and F_tot, 129 ≈ 0.6 − 0.8 Jy at 22, 43, 86, and 129 GHz, respectively. The KVN 129 GHz data suffer from larger uncertainties due to a stronger influence from weather. We fit a power-law model to the total flux densities (i.e., F_tot ∝ ν^+α_tot). We find α_tot ≈ −0.17 when 129 GHz data are included, while α_tot ≈ −0.11 when they are excluded. By extrapolating a power-law model with a range α_tot ≈ −(0.11 − 0.17) to 230 GHz, we estimate F_tot, 230 = 0.72 − 0.80 Jy. While this range of values is lower than the ALMA interferometric measurements during the EHT observations (∼1.11−1.18 Jy; Table 3 in Goddi et al. 2021), it only includes the total flux density within the inner 500 × 500 μas², which is nearly three orders of magnitude below the resolution limit of ALMA.

The compact flux density F_cpct, 86 at 86 GHz is estimated from a GMVA+ALMA image that has a higher angular resolution than the KVN images (Lu et al. 2023). We obtain F_cpct, 86 ≈ 0.51 Jy by integrating over a smaller window of 100 × 100 μas² centered on the radio core. We can then derive a jet flux density at 86 GHz to be F_jet, 86 = F_tot, 86 − F_cpct, 86 ≈ 0.39 Jy.

Finally, we compute the expected F_cpct, 230 by estimating F_jet, 230 and subtracting it from the previously estimated F_tot, 230. As explained above, we assume that F_jet, 230 follows a power law that can be extrapolated up to 230 GHz (i.e., F_jet, 230 ∝ ν^+α_jet). For a range of α_jet = −(0.7 − 1.0), we find F_jet, 230 = 0.15 − 0.20 Jy. This value gives a range of F_cpct, 230 ≈ 0.5 − 0.7 Jy within the central 100 × 100 μas² region. We emphasize that this estimate is based on lower-frequency data and assumptions on the jet spectral index. Therefore it should not be regarded as a tight constraint on F_cpct, 230 but more as a useful reference.

Appendix G: Imaging validation

Here we describe the synthetic data generation process and subsequent reconstructions, as well as the resulting distribution of Top Set images for the real M87^* data. We also investigate the properties of the image reconstructions when removing individual baselines.

G.1. Synthetic data

In order to test the efficacy of the imaging pipelines, we generate synthetic data sets using both geometric models and snapshots of GRMHD simulations. All of the geometric models contain 0.6 Jy in compact flux density and 0.5 Jy in a larger extended jet modeled by three Gaussians. We show the ring model as an example in Fig. G.1. The sizes and locations of these three Gaussians are identical to those reported in Table 10 of M 87^∗ 2017 IV, with the flux densities scaled down to total 0.5 Jy rather than 0.6 Jy. Of the geometric models we generated, four of them were utilized in the procedure to generate the Top Set images. These four geometric models are very similar to those used in M 87^∗ 2017 IV and are as follows:

• 1.

  cres180: An asymmetric ring model with r₀ = 23 μas, a brightness position angle oriented south, and blurred by a circular Gaussian beam of FWHM 10 μas.

• 2.

  ring: A thin uniform ring of radius r₀ = 23 μas blurred by a circular Gaussian beam of FWHM 10 μas.

• 3.

  dblsrc: Two circular Gaussian components each with FWHM of 20 μas. One is located at the origin with a flux density of 0.27 Jy, while the second is positioned at Δ_R.A. = 30 μas and Δ_decl. = −12 μas with a flux density of 0.33 Jy.

• 4.

  disk: A uniform disk of radius r₀ = 35 μas blurred by a circular Gaussian beam of FWHM 10 μas.

Ground-truth images of these models are shown in Fig. G.2. Aside from these training data, we also generated seven validation data sets to ensure good performance from the CLEAN and RML method Top Sets. These data, along with reconstructions from each pipeline, are shown in Fig. G.3. The model descriptions are as follows

• 1.

  cres–90, cres0, cres90: Three asymmetric ring models with r₀ = 23 μas, and brightness position angles oriented east, north, and west, blurred by a circular Gaussian beam of FWHM = 10 μas.

• 2.

  edisk: A uniform elliptical disk model with a major axis of 66 μas, a minor to major axis ratio of 0.65, a major axis position angle of 60°, and a blurring of 10 μas.

• 3.

  point+disk: A point source plus symmetric disk model containing a 10 μas point source centered in 100 μas diameter disk. The point source to disk flux ratio is chosen to be 0.192.

• 4.

  point+edisk: A point source plus elliptical disk model containing a 10 μas point source centered in an ellipse of major axis of 96 μas, a minor to major axis ratio of 0.8, a major axis position angle of 60°, and a blurring of 10 μas. The point source to disk flux ratio is 0.16.

• 5.

  GRMHD: A snapshot from a GRMHD simulation of M87^*.

The synthetic visibilities were generated using eht-imaging and all include random thermal noise, station-based gain and phase errors, and polarimetric leakage. For a more in-depth explanation of synthetic visibility generation, we refer readers to Appendix C.2 of M 87^∗ 2017 IV.

Figures G.2 and G.3 show the image reconstructions for four training and seven validation data sets including the GRMHD model from each imaging pipeline, respectively. The image reconstructed with the fiducial parameters is displayed for DIFMAP, eht-imaging, and SMILI. A single image randomly selected from the posterior is displayed for THEMIS and Comrade. The THEMIS reconstruction of the GRMHD synthetic data uses a Hybrid Themage model. This figure clearly shows that a very wide range of morphology types can be recovered by all imaging pipelines. These simulations are based on the same (u, v) coverage of the 2018 EHT observation, and thus have a common point spread function (PSF). This ensures that the recovered structures by imaging pipelines are not artificially caused by the (u, v) coverage and/or the PSF, but rather come from the real structure imprinted in the data.

Based on ρ_NX values and ring parameters, we again confirmed that most of the Top Set parameters can reconstruct the correct morphology. For instance, 95, 96, 97% of cres–90 April 21 band 3 images reconstructed with the Top Set parameters from DIFMAP, eht-imaging, SMILI passed the ρ_NX criteria imposed on the training data sets though some images have deviated structure as seen in the training data sets. Figure G.4 shows the significance of the deviated structure in the Top Set of synthetic data reconstructions with SMILI. Above a certain dynamic range threshold (e.g., 10% of peak intensity), artificial structures which deviate from the ground truth are shown but its fraction in the Top Set is minor, especially the repeated central structure. In addition, the fiducial images recover the station gains as presented in Fig. G.5. This shows the multiplicative gains at each telescope, in comparison with true gains from the cres180 model on April 21 band 3. In a same manner, the station gains are derived for M87^* which are consistent across different imaging approaches, as presented in Fig. G.6. Certain imaging pipelines may perform poorly at certain epochs and bands. In case of point+disk April 21 band 3, only 40% of DIFMAP images reconstructed with the Top Set parameters passed the ρ_NX criteria.

G.2. The distribution of M87^* Top Set images

As shown in the main text, all Top Set M87^* images have consistent ring features. However, there are slight differences among the Top Set reconstructions. Figure G.7 shows 25 randomly selected Top Set reconstructions from HOPS band 3 April 21 data for DIFMAP, eht-imaging, and SMILI. We see the ∼40 μas ring for all images, and most images have a similar brightness feature, but the ring width seems to have relatively large uncertainty depending on the image assumptions of each imaging method. Also, a few images from eht-imaging and SMILI show double-ring structures although they are a minority in the Top Sets. As discussed above, we see similar minor double-ring structures in ring and crescent synthetic data reconstructions, so the second ring is unlikely to be real and we anticipate this does not affect our conclusion. These trends are also seen on the other observational dates and bands. Image variations among Top Sets are relatively larger, especially for data with poor coverage.

G.3. Image dependence on sites

In addition to analyzing the dependence of M87^* images on frequency bands and changes in coverage from different epochs, it is necessary to understand how much the reconstructed image relies on data from single antenna sites. For this purpose, we imaged M87^* excluding all the data from baselines connected to one site, and repeated the imaging for all six sites. The test was performed on April 21 band 3 data, using the eht-imaging and Comrade pipelines as representative for the RML and Bayesian methods. The first row of Fig. G.8 shows the mean image reconstructed by Comrade, the second row shows eht-imaging images reconstructed with fiducial parameters, while the third one presents eht-imaging fiducial images of the 2017 April 11 data, reconstructed with the 2018 pipeline. We see that in most cases, the removal of one antenna in 2018 affects the quality of the image reconstruction, both for the Comrade and the eht-imaging pipelines.

The strongest image degradation happens when removing the Chile sites, which are fundamental to reconstruct the basic image morphology, but also imaging without the SMT or the Hawai‘i sites results in the appearance of a tessellation pattern for eht-imaging and in a loss of contrast for Comrade. Removing PV baselines results in an elongation of the ring-like structure for both pipelines, while the absence of either the GLT or the LMT stations only introduces minor artifacts, without changing the overall image morphology.

Compared to 2017, the 2018 eht-imaging reconstructions without one antenna are affected by a stronger tessellation pattern. However, the morphology of the ring-like structure is more robust to the removal of a single site in 2018 compared to 2017. For example, in 2018, a ring-like structure is still distinguishable even without the Chile site. Also the removal of the LMT antenna does not affect the orientation of the ring, and the removal of PV does not strongly affect the morphology.

Appendix H: Imaging with DIFMAP using a point source model for initial self-calibration

This Appendix reports the results of the DIFMAP parameter survey, which used an initial phase-only self-calibration process assuming a point source starting model. Unlike the parameter survey discussed in the main text (Sect. 5.1.1), which seeks to determine the best geometric model for initial self-calibration, this survey requires a careful consideration of the CLEAN window location. This is due to the fact that the self-calibration process results in the brightest part of the source being placed at the map’s coordinate origin, requiring a shift in the position of the window center’s position for the subsequent CLEAN and self-calibration procedures.

To ensure optimal performance, the ideal window center location should be specified for each source model and not treated as a global free parameter. In this analysis, we employed an automated optimal mask search procedure during the parameter survey using a point source model for initial self-calibration. The procedure involved a systematic search of the parameter space by shifting the window center from the map origin in increments of 5 μas along the right ascension and declination directions. The shift is limited to a maximum value equal to the size of the CLEAN window. To maintain consistency, the five imaging parameters presented in Sect. 5.1.1 were kept constant at their fiducial parameters which were obtained from the CLEAN imaging of the 2017 EHT M87^* data (M 87^∗ 2017 IV). These parameters included a compact flux density of 0.5 Jy, stopping CLEANing when the required compact flux density was reached, an ALMA weight re-scaling factor of 0.1, a window diameter of 60 μas, and a (u, v) weight exponent of −1. The procedure computed the ${\chi }_{\mathrm{CP}}^2$ value for each image in the survey, as described in Sect. 5.2, and determined the optimal window position based on the image with the lowest ${\chi }_{\mathrm{CP}}^2$ value. This approach enabled the determination of the optimal window position for each source model and window diameter.

Next, we conducted a parameter survey by varying the five imaging parameters (Sect. 5.1.1), assuming the optimal window position obtained during the automatic mask search procedure. We followed the image selection criteria outlined in Sect. 5.2 to select both Top Set and fiducial images from the survey. The fiducial images resulting from this survey are displayed in Fig. H.1. Our findings indicate that the survey was able to accurately reconstruct the majority of the synthetic data models. For the real M87^* data, the Top Set images displayed a ring-like morphology, consistent with the results from the survey discussed in Sect. 5.1.1. This result demonstrates that the ring-like structure for M87^*, reconstructed by DIFMAP, is stable independently of the geometrical model selected for the first phase of self-calibration. However, the overall performance of this survey was lower than that of the survey that optimized the geometric models for initial phase-only self-calibration, resulting in a smaller number of Top Set images. This can be attributed to the fact that the point source model used as initial phase-only self-calibration model in this survey is overly simplistic for many geometric models and actual M87^* source structure. Despite encountering challenges with the uniform disk model, resulting in a lower number of Top Set images obtained from the survey, successful reconstruction of both the ring and crescent synthetic models was achieved. These models were found to closely resemble the true source structure.

Appendix I: Testing the backward compatibility of the DIFMAP pipelines on the 2017 data

The DIFMAP pipelines presented in Sect. 5.1.1 are updated as compared to the DIFMAP pipelines used in previous EHT imaging analyses (M 87^∗ 2017 IV; Sgr A^∗ 2017 III). We conducted a backward compatibility test of these pipelines on the 2017 EHT M87^* data, including the real data from the HOPS pipeline (April 11, low-band) and the synthetic data sets for the training. Our analysis revealed that the resulting images align well with those previously reported in (M 87^∗ 2017 IV). Our fiducial images of the M87^* data and synthetic data can be found in Fig. I.1.

Appendix J: THEMIS imaging model details and priors

The THEMIS image reconstructions are composed of three different primary model components, the splined raster model, a large-scale asymmetric Gaussian, and a thin slashed crescent. A more detailed description of the large-scale Gaussian and crescent models, as well as details about their technical implementation, are described in Broderick et al. (2020a). A much more detailed discussion of the raster model can be found in Broderick et al. (2020b).

The splined raster model is defined by a rectilinear set of control points which may independently vary in intensity, I_M, N. The intensity map is produced using an approximate cubic spline interpolation between the control points. The field of view FOV_x and FOV_y are permitted to vary, along with the orientation of the grid ϕ and the raster center (x, y). The large-scale Gaussian is asymmetric, and is characterized by its total flux I_G, a symmetric standard deviation σ_G, asymmetry parameter A_G, position angle ϕ_G, and position (x_G, y_G). The slashed crescent is based on the xs-ringauss model also described in Broderick et al. (2020b), which also serves as the core of the THEMIS xs-ringauss geometric model featured in Sect. 6.2 and in M 87^∗ 2017 VI.

As a part of the Hybrid Themage model, we eliminate the contributions from the asymmetry and pinned Gaussian parameters by appropriately restricting their priors. We also restrict the priors on the fractional width parameter ψ to range from 0% to 5%. This ensures the ring width is well below that of the raster spacing, to avoid excessive flux trading between the ring and raster components. This leaves the ring flux I_X, outer radius R_out, linear brightness gradient f_X, brightness gradient position angle ϕ_X, and position (x_X, y_X) as searchable parameters.

We use a deterministic even-odd (DEO) swap tempering scheme with a Hamiltonian Monte-Carlo sampling kernel based on the Stan package. Chain convergence is assessed based on traditional criteria such as the split $\widehat{R}$, auto-correlation time, and visual inspection of the parameter traces. The parameters and priors for each model component are listed in Table J.1.

Appendix K: THEMIS Raster dimensions

The best fit raster dimensions N_x × N_y for the THEMIS Raster models in principle depend on the (u, v) coverage and source complexity. To construct the best model, one should survey the raster dimensions for each unique data set. In order to provide easier comparisons across data sets, we only survey the raster dimensions on our primary science data set, April 21 band 3, and use those best fit raster dimensions for the other data sets. We present the results of this survey in Table K.1, where we survey a 4 × 4, 5 × 5, 6 × 6, and 7 × 7 raster. We rank the models in terms of the difference in the logarithm of the Bayesian Evidence (Δlog(Z)), where positive values are more preferred. The 5 × 5 raster is most preferred by this data, followed by the 4 × 4 raster. The 5 × 5 raster is also the same dimensions as the model used in Broderick et al. (2022a), which applied the Hybrid Themage model to the 2017 data. Thus, the image structures we show in this analysis should be directly comparable to the image structures in the previous paper. We also report the best fit reduced χ² for the surveyed raster sizes, and note that all models provide good fits to the data.

We can also compare the fit quality between the standard THEMIS raster and Hybrid Themage models by similarly comparing the Δlog(Z). In Table K.2 we show the reduced χ² and difference in Bayesian evidence between the Hybrid Themage and raster-only model for both days and all bands. Here we see that the Hybrid Themage model is generally preferred in nearly all data sets, except for the April 21 band 1 data.

Appendix L: Comrade image model and prior

For all Comrade image reconstructions, we construct our image model in the following steps. We begin with a rasterized image model convolved with a B-spline kernel of order 3 (see Eqs. L.1 and L.2) to generate flux densities that smoothly vary in all directions. We set the raster points by (l_i, m_j) given by Eq. L.3, for a given FOV and number of pixels N_x and N_y for each side. The FOV and number of pixels are the hyperparameters of the model:

$$\begin{array}{c}\hfill I(l,m)={\displaystyle \sum _{i,j}}B\left(\frac{l-l_i}{s_x}\right)B\left(\frac{m-m_j}{s_y}\right)I_{\mathit{ij}},\end{array}$$(L.1)

where B(x) is a B-spline kernel of order 3 given by,

$$\begin{array}{c}\hfill B(x)=\{\begin{array}{cc}{x}^2/2\hfill & 0\le x<1\hfill \\ (-2x^2+6x-3)/2\hfill & 1\le x<2\hfill \\ {(3-x)}^2/2\hfill & 2\le x<3\hfill \end{array}\end{array}$$(L.2)

and raster points are given by,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill l_i& =-\mathrm{FOV}/2+s_x/2+s_x(i-1)\hfill \\ \hfill m_j& =-\mathrm{FOV}/2+s_y/2+s_y(j-1)\hfill \end{array}\end{array}$$(L.3)

where s_x = FOV/N_x and s_y = FOV/N_y. We restrict the images considered in this paper to having N_x = N_y, such that the images are square and s_x = s_y. We add a circular Gaussian with FWHM = $2\sqrt{2ln}\sigma =1$ mas (G_1mas) to the raster to model the amplitudes and station gains of short intra-site (ALMA–APEX, JCMT–SMA) baselines. We form our final image model,

$$\begin{array}{c}\hfill M(I_{\mathit{ij}},f,f_g)=f(1-f_g)\times I(I_{\mathit{ij}})+f_{g}f\times G_{1\mathrm{mas}}\end{array}$$(L.4)

where f is the total flux density of the raster (also described as compact flux) and f_g is the fraction of the total flux density that corresponds to the flux density of the Gaussian. We take the Fourier transform of M(I_ij, f, f_g) to get model Fourier transform ${\stackrel{\sim }{I}}_{\mathit{ab}}$ for each baseline ab. The Fourier amplitudes are then corrupted to model the individual station gain amplitudes |g| to construct the visibility amplitudes, $|{\widehat{\mathcal{V}}}_{\mathit{ab}}|=|{\mathit{g}}_a||{\mathit{g}}_b||{\stackrel{\sim }{I}}_{\mathit{ab}}|$

We form the visibility amplitude likelihood and closure phase likelihood as described in Appendix F of Sgr A^∗ 2017 IV. Next, we form our prior distribution by using a uniform distribution 𝒰(0.0, 1.5) Jy for the prior on the total flux density f and 𝒰(0.0, 1.0) for the fraction of the total flux density f_g for the flux density of the large-scale Gaussian. For the raster pixels fluxes, I_i, we use a symmetric Dirichlet distribution

$$\begin{array}{c}\hfill p(\mathbf{I}|\xi )=\frac{\mathrm{\Gamma }(K\xi )}{\mathrm{\Gamma }{(\xi )}^K}{\displaystyle \prod _{i=1}^K}{I}_i^{\xi -1},\end{array}$$(L.5)

where I is the flattened vector of pixel fluxes, K is the total number of pixels in the image, Γ is gamma distribution, and ξ is the concentration parameter that controls the sparsity and smoothness of the image. The Dirichlet support is on the simplex, meaning that ∑_iI_i = 1 and 0 ≤ I_i ≤ 1. For this work, we set ξ = 1, which is equivalent to a uniform distribution on the simplex.

From the analysis of the crossing baseline tracks and the EHT station flux density calibration parameters (Table A.1), for station gain priors, we used a Normal distribution for the log-gain amplitudes for each station. For the April 21 data, we imposed prior widths of 10% or all stations, that is a 𝒩(0.0, 0.1) prior, except 30% and 100% prior widths for LMT and GLT, respectively. We used the same gain priors for the April 25 data except for PV which used a 100% prior width, since it was found that PV exhibited large amplitude fluctuations after UTC 02:00 due to poor weather as mentioned in Appendix D.

Appendix M: Impact of compact flux density on image reconstruction

The non-imaging constraints on compact flux density permit a wide range of allowable values, and a comparison of the compact flux densities measured by the Bayesian methods suggests that the current data set cannot be used to constrain the compact flux on horizon scales. In this appendix, we test the robustness of the ring reconstructions from the RML and CLEAN methods against different assumptions of the compact flux density.

Firstly, we investigated the dependence of the assumed total compact flux density (F_cpct) on the estimated ring parameters of M87^* from reconstructed images following the approach adopted in our previous EHT analysis (M 87^∗ 2017 IV). In this analysis, we only change F_cpct, while all other imaging parameters are kept at their fiducial values. We select acceptable images based on two criteria: (1) normalized ${\chi }_{\mathrm{CP}}^2$,  ${\chi }_{\log \mathrm{CA}}^2<2$ (with a 10% systematic uncertainty for DIFMAP), and (2) a corresponding self-calibration solution with 0.9 ≤ median(1/|g_SMT|) ≤ 1.1, where g_SMT is the gain for the SMT station. Figure M.1 illustrates how the estimated ring parameters vary with the assumed total compact flux density. While the ring width for all pipelines and fractional central brightness for RML methods slightly increase with F_cpct, the ring diameter and position angle values of the acceptable images do not exhibit significant dependence on F_cpct. Therefore, the ring diameter and position angle are robust against different assumptions on the compact total flux density.

Secondly, We tested if our DIFMAP imaging script could distinguish different morphologies for different values of the total compact flux density using synthetic data. We prepared a series of synthetic data with the same morphologies of four geometric models used for Top Set selection in the imaging survey, but with the assumed compact total flux density in range from 0.4 to 1.1 Jy. Here, we used the total flux density, which is the sum of the assumed compact total flux density and flux density of a larger extended jet modeled by three Gaussians, to be 1.1 Jy for all synthetic data. We conduct image reconstruction using the DIFMAP imaging script with the fiducial parameter sets, but adopting the compact flux densities to be the same as that of the assumed values of synthetic data. In Fig. M.2, we show the reconstructed images of geometric models with different values of the total compact flux density. It is clearly demonstrated that our DIFMAP imaging script can distinguish ring-like morphologies from disk or double source morphologies if we properly choose the expected total compact flux density.

Appendix N: Geometric modeling details

Image reconstructions of M87^* show strong support for a ring-like structure. We quantify this support through the Bayesian evidence as a robust metric for the goodness of fit.

The Bayesian approach is reliant on making an appropriate estimate of the posterior distribution, P(Θ|M, D), for some parameters Θ of a chosen model, M, that is conditioned on some data, D. The posterior is defined by the relationship,

$$\begin{array}{c}\hfill P(\mathrm{\Theta }|\mathbf{M},\mathbf{D})=\frac{P(\mathbf{D}|\mathbf{M},\mathrm{\Theta })P(\mathrm{\Theta }|\mathbf{M})}{P(\mathbf{D}|\mathbf{M})}\equiv \frac{\mathcal{L}(\mathrm{\Theta }|\mathbf{M})\pi (\mathrm{\Theta }|\mathbf{M})}{Z(\mathbf{M})},\end{array}$$(N.1)

which also defines the likelihood of the data for a chosen model and parameter,

$$\begin{array}{c}\hfill \mathcal{L}(\mathrm{\Theta }|\mathbf{M})\equiv P(\mathbf{D}|\mathbf{M},\mathrm{\Theta })\end{array}$$(N.2)

the prior probability on model parameters,

$$\begin{array}{c}\hfill \pi (\mathrm{\Theta }|\mathbf{M})=P(\mathrm{\Theta }|\mathbf{M})\end{array}$$(N.3)

and the Bayesian evidence for a chosen model,

$$\begin{array}{c}\hfill Z(\mathbf{M})\equiv P(\mathbf{D}|\mathbf{M})={\displaystyle \int \mathcal{L}}(\mathrm{\Theta }|\mathbf{M})\pi (\mathrm{\Theta }|\mathbf{M})d\mathrm{\Theta }.\end{array}$$(N.4)

The Bayesian evidence allows the relative support between two models to be determined through a Bayes factor,

$$\begin{array}{c}\hfill R({\mathbf{M}}_1,{\mathbf{M}}_2)\equiv \frac{Z({\mathbf{M}}_1)}{Z({\mathbf{M}}_1)}\frac{\pi ({\mathbf{M}}_2)}{\pi ({\mathbf{M}}_2)}\end{array}$$(N.5)

where π(M) is the probability prior defined on the set of models. We use the Bayesian evidence with a constant prior to define the goodness of fit for a model M ∈ {M} as,

$$\begin{array}{c}\hfill \mathrm{\Delta }log(Z)=logZ(\mathbf{M})-\mathrm{maximum}[logZ(\{\mathbf{M}\})],\end{array}$$(N.6)

where the argument of Δlog(Z) is the Bayes factor between M and the best performing model in the set.

The best performing models are selected with this metric to construct fiducial models for feature comparison; namely, a stretched m-ring of order 2 with a flat emission floor and two elliptical Gaussians (mF-ring), and an xs-ringauss with two elliptical Gaussians an a large scale circular Gaussian. The schematic summary of defining parameters of these models is shown in Fig. N.1. See Sect. 6.2 for details on these models.

The shared geometric parameters of interest between both the mF-ring and the xs-ringauss model are their diameters, widths, central brightnesses, and brightness asymmetry position angles. The diameter of the mF-ring is defined to be the sum of the debiased semi-major and semi-minor axis, while that of the xs-ringauss is taken to be the sum of its inner and outer radii,

$$\begin{array}{c}\hfill \widehat{d}=\{\begin{array}{cc}{\widehat{r}}_0+{\widehat{r}}_1,\hfill & \mathrm{m}F\text{-}ring\hfill \\ R_{\mathrm{out}}+R_{\mathrm{in}},\hfill & \mathrm{xs}\text{-}\mathrm{ringauss}.\hfill \end{array}\end{array}$$(N.7)

Here, the debiased radius (${\widehat{r}}_i$) is linked to the defining radius of the mF-ring (r_i), and the FWHM of the blurring Gaussian kernel (σ^*) through:

$$\begin{array}{c}\hfill {\widehat{r}}_i=r_i-\frac{{\sigma }^{\ast 2}}{16ln(2)r_i}.\end{array}$$(N.8)

We define the fractional widths of the mF-ring and the xs-ringauss respectively to be,

$$\begin{array}{c}\hfill f_w=\{\begin{array}{cc}\frac{{\sigma }^{\ast }}{\widehat{d}},\hfill & \mathrm{m}F\text{-}ring\hfill \\ [7pt]\frac{R_{\mathrm{out}}-R_{\mathrm{in}}+{\sigma }^{\ast }}{\widehat{d}},\hfill & \mathrm{xs}\text{-}\mathrm{ringauss}.\hfill \end{array}\end{array}$$(N.9)

and the brightness asymmetry position angle to be,

$$\begin{array}{c}\hfill \widehat{\varphi }=\{\begin{array}{cc}{\beta }_1+{\varphi }_{\mathrm{m}},\hfill & \mathrm{m}F\text{-}ring\hfill \\ {\varphi }_{\mathrm{X}},\hfill & \mathrm{xs}\text{-}\mathrm{ringauss}.\hfill \end{array}\end{array}$$(N.10)

The last shared parameter shown between the two models is the fractional central flux depression, or relative central brightness, which is given by,

$$\begin{array}{c}\hfill {\widehat{f}}_{\mathrm{C}}=\{\begin{array}{cc}\frac{\mathrm{brightness}\mathrm{of}\mathrm{emission}\mathrm{floor}}{\mathrm{mean}\mathrm{brightness}\mathrm{of}\mathrm{m}\text{-}\mathrm{ring}},\hfill & \mathrm{m}F\text{-}ring\hfill \\ [7pt]\frac{\mathrm{brightness}\mathrm{of}\mathrm{central}\mathrm{Gaussian}}{\mathrm{mean}\mathrm{brightness}\mathrm{of}\mathrm{crescent}},\hfill & \mathrm{xs}\text{-}\mathrm{ringauss}.\hfill \end{array}\end{array}$$(N.11)

As described above in the main text, the THEMIS xs-ringauss model used to fit the data from the 2018 observations is the same model used to fit the data from the 2017 observations. We list the model parameters and their priors in Tables N.1 and N.2.

Even though we construct and fit the Hybrid Themage, xs-ringauss, and mF-ring in the visibility domain, we can also produce conventional image domain representations of the models. Figure N.2 shows the on-sky representations of the Hybrid Themage, xs-ringauss, and mF-ring models. We find that these models reproduce many of the same features seen in the more agnostic models shown in Sect. 5, such as a consistent diameter and position angle of the brightest part of the ring in the southwest. The xs-ringauss and mF-ring models use a pair of asymmetric Gaussians as nuisance components, designed to help fit the residual data that cannot be captured by the ring component. The raster component in the Hybrid Themage serves a similar purpose. We note that for April 21 band 3 and band 4, the non-ring components settle in the west and southwest part of the image, but these components occupy more varied locations in the band 1, band 2, and April 25 reconstructions.

All Tables

All Figures