Press Release
Open Access
Issue
A&A
Volume 681, January 2024
Article Number A79
Number of page(s) 63
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202347932
Published online 18 January 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 , 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 Rc 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 (∼106 − 109M); 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 dsh = 2RcD−1 ≈ 10 GMc−2D−1 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 dsh 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 × 109M (Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019f). The combination of D and M yields a dsh ≈ 40 μas, which is the second largest value on the sky after Sgr A* (M ≈ 4 × 106M, D ≈ 8 kpc, dsh ≈ 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) × 109M (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.

thumbnail Fig. 1.

Representative example images of M 87* from the EHT observations taken on 2017 April 11 and 2018 April 21 (north is up and east is to the left). The 2017 image is generated with the average of fiducial parameter sets from the imaging techniques used in M 87 2017 IV. The 2018 image is created by taking the average of the blurred images generated by the imaging techniques found in Sect. 5. Comparison of the images shows consistency in the diameter across observation epochs, but a shift in position angle of brightness asymmetry. The circle represents a Gaussian blurring kernel with a full width half maximum of 20 μas.

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.

thumbnail Fig. 2.

Map showing the stations that participated in the EHT 2018 campaign (black circles), which differs from the EHT array in 2017 by the addition of the GLT. Co-located sites in Chile and Hawai‘i appear superimposed. The SPT projected location from the back of the map is indicated with a dashed circle, and baselines to this station are represented with dashed-lines. While the SPT cannot observe M 87*, it observed 3C 279 and was used to calibrate the data.

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.

Table 1.

Median zenith sky opacities (1.3 mm) at EHT sites during the 2018 April observations toward M 87*.

thumbnail Fig. 3.

EHT observing schedules for M 87* (blue) and its calibrator 3C 279 (orange) on the 2018 April 21 (left panel) and April 25 (right panel) observing days, which began at the end in UTC of April 20 and 24, respectively. Open rectangles represent scans that were scheduled but not observed owing to weather or technical issues. The filled rectangles mark the scans with detections in the final data set. On these two particular days, scan durations are typically 4 to 5 min each for M 87* and 4 min each for 3C 279, as reflected by the width of each rectangle.

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).

thumbnail Fig. 4.

M 87* (u, v) coverage (colored points) in band 3 (top panels) and band 2 (bottom panels) for observations on 2018 April 21 (left panels) and April 25 (right panels), overlaid on the low-band (u, v) coverage for 2017 April 11 (gray points). The dashed circles show baseline lengths corresponding to fringe spacings of 25 and 50 μas respectively. The (u, v) coverage in band 1 and band 4 is comparable to that in band 2 and band 3, respectively.

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 obtained1. 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.

thumbnail Fig. 5.

Measured correlated flux densities of M 87* as a function of baseline lengths in units of wavelength, for 2018 April 21 (top panels) and April 25 (bottom panels) in band 3, for both HOPS (left panels) and CASA (right panels) outputs. The 2018 data (colored points) are overlaid on the corresponding flux densities of the 2017 April 11 observations in low-band (gray points). All data shown include a priori station-based amplitude calibration and network calibration but are prior to any model-dependent self-calibration. Error bars denote ±1σ from thermal noise. Redundant baselines are shown with different symbols: circles for baselines to ALMA and SMA; diamonds for baselines to APEX and JCMT. The dashed line corresponds to an azimuthally symmetric thin ring model with a 42 μas diameter.

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.

thumbnail Fig. 6.

Variation of closure phases as a function of Greenwich mean sidereal time (GMST) for selected triangles using HOPS calibrated data. Red diamonds denote data from 2017 April 6 (low-band), gray circles denote data from 2017 April 11 (low-band), and blue squares denote data from 2018 April 21 (band 3). Error bars are the ±1σ uncertainties.

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 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 Ftot ≈ 1.13 Jy by averaging the values over four bands on April 21 and 22. The second one is the total compact flux density, Fcpct ≤ Ftot, which is the integrated flux density of the microarcsecond-scale structure constrained by inter-site baselines.

We estimate the range of Fcpct 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 Fcpct and θcpct simultaneously. The other is a constraint on the Fcpct 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 Fcpct and θcpct using these procedures.

The EHT constraints give Fcpct = 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 Fcpct 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 Fcpct = 0.5 − 0.7 Jy within the central 100 × 100 μas2 of the nucleus, which is comparable to the Fcpct = 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 <  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.

Table 2.

Parameter survey results for April 21 band 3 data.

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.15.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 χ2. 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.

Table 3.

Initial DIFMAP geometric models.

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 pipeline2. 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 Nx × Ny), 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 (Nx × Ny). 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.

Table 4.

Comrade raster model hyperparameters.

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 fg 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- (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.

thumbnail Fig. 7.

Representative images recovered from the HOPS and CASA data with all five imaging pipelines for two observing days (April 21 and 25). Each panel shows the fiducial image of the corresponding top set images for the DIFMAP, eht-imaging, and SMILI pipelines, and a random sample from the respective posterior for the THEMIS and Comrade pipelines. We do not have Top Sets for band 1 and band 2 from DIFMAP, eht-imaging, and SMILI pipelines on April 25. The dashed horizontal line in each block separates the DIFMAP and RML methods above from the Bayesian methods below. The circles in the DIFMAP panels represent an effective Gaussian blurring kernel of 20 μas. The solid lines in the THEMIS and Comrade panels represent the size of the blurring kernel used to achieve the same effective resolution as the DIFMAP method.

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 (χ2), normalized by the number of data points to be compared. For example, we define the normalized closure phase as:

thumbnail Fig. 8.

Closure phases plotted as a function of GMST on three selected triangles from the April 21 band 3 observations (black points). The error bars on the data points denote the ±1σ uncertainties. The colored and dashed lines indicate the model closure phase curves from the fiducial images or posterior samples produced by the five imaging pipelines.

thumbnail Fig. 9.

Visibility amplitudes of band 3 data on April 21 as a function of baseline length compared with corresponding gain-calibrated visibility amplitudes from the five representative image models from each pipeline. The fiducial images are used for DIFMAP, eht-imaging, and SMILI, and the maximum likelihood model from the sampling chains are used for THEMIS and maximum a posteriori model for Comrade. The gray points represent the data used by each method after flagging, but before individual self-calibration. eht-imaging scales the intra-site baselines to 0.5 Jy, and SMILI scales the intra-site baselines to 0.6 Jy, and both pre-calibrate the LMT-SMT baselines before fitting. DIFMAP uses 10 s averaged data, and all other methods use scan-averaged data. THEMIS and Comrade apply 1% fractional systematic noise to all baselines, and eht-imaging and SMILI apply 1% and 0% fractional systematic noise respectively to all baselines, seen as minor differences in the gray points. The colored points represent the image model visibilities from each method, with station gains derived from each method’s internal self-calibration procedure applied to the image model visibilities. Below each visibility amplitude figure are the normalized residuals for each image, which is the difference between the gray and the colored points, divided by the uncertainty of the gray data points.

(1)

where Nψ is the number of closure phase data points, 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 as:

(2)

where NlnAC is the number of log closure amplitude data points, is the model closure amplitude, AC is the real data closure amplitude, and 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 χ2 values, which should not be confused for the formal definition of a reduced χ2. In Table 5, we show the normalized χ2 values for the fiducial images and normalized χ2 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 χ2 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 χ2 values to the RML methods.

Table 5.

Fiducual image and Top Set closure quantity normalized χ2 values.

For the Bayesian imaging methods, the fit quality is assessed through comparisons of the complex visibility reduced χ2 for THEMIS, the visibility amplitude and closure phase reduced χ2 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 χ2 values in Table 6.

Table 6.

Reduced χ2 quantities for the THEMIS and Comrade raster models.

We calculate the reduced χ2 by dividing the standard χ2 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 χ2 (Red. ):

(3)

(4)

where NV is the number of complex visibility data points, Nmodel DOF is the number of model degrees of freedom, Ngain 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), is the model complex visibility, V is the data complex visibility, and σV is the complex visibility standard deviation. One can construct similar reduced χ2 quantities for the visibility amplitudes (Red. ) and closure phases (Red. ) 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 χ2 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 χ2 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 Ftot and Fcpct 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 1012 show the image- and visibility-domain uncertainties associated with the image sets for eht-imaging, THEMIS, and Comrade, respectively.

thumbnail Fig. 10.

Visualization of image statistics calculated using the Top Set images from the eht-imaging pipeline for observations taken on April 21 band 3. We emphasize that these images do not represent the posterior probability space for the reconstructions. Each image reconstructed using eht-imaging is the maximum a posteriori (MAP) image for a given parameter set. Thus, the statistics shown represent uncertainties that arise from different choices of regularizer weights, not from an exploration of posterior space. The top row shows top statistics in the image domain while the bottom row shows the visibility domain. Overlaid on the visibility domain panels is the (u, v) coverage for the April 21 observation. From left to right, we present the mean image; the standard deviation; the normalized standard deviation, calculated by rescaling each image to the flux of the mean image; and the fractional standard deviation, calculated by dividing the standard deviation by the mean. The fractional standard deviation panel has been clipped to a maximum value of 1. Portions of the image exhibit large fractional standard deviations due to pixel values very close to zero in the mean image. In the top row, image contours are drawn at 10%, 20%, 40%, and 80% of the peak values from the mean image. In the bottom row, the gray contours represent 0.1%, 1%, and 10% of the peak while the black contours represent 10 and 100 mJy (left three panels) and 0.1 (right most panel). The complex visibilities are calculated by taking a Fourier transform of the images and then calculating the mean and standard deviation. The absolute value of the mean and standard deviation of the complex visibilities are used to calculate visibility amplitudes.

thumbnail Fig. 11.

Visualization of image uncertainties using images from the posterior of THEMIS pipeline for observations taken with band 3 on April 21. The contour lines shown are drawn as described in Fig. 10.

thumbnail Fig. 12.

Visualization of image uncertainties using images from the posterior of Comrade pipeline for observations taken with band 3 on April 21. The contour lines shown are the same as described in Fig. 10.

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.

thumbnail Fig. 13.

Ring characteristics from all bands, observational days, and imaging pipelines coming form HOPS and CASA data (colored points). Median values and 68% confidence intervals are shown for diameter and width. Circular mean and standard deviation values are shown for the position angle. Circle and triangle markers correspond to REx and VIDA respectively. Red (eht-imaging) and blue (DIFMAP) dashed lines and shaded regions show the ring parameters and 68% confidence intervals from the 2017 April 11 measurements. The vertical gray shaded region corresponds to the 2018 April 21 measurements, and the vertical un-shaded region corresponds to the 2018 April 25 measurements.

Tables 710 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 fC is ∼0 − 0.5, indicating a clear central depression.

Table 7.

Ring parameters for 2018 April 21 HOPS data.

Table 8.

Ring parameters for the 2018 April 25 HOPS data.

Table 9.

Ring parameters for 2018 April 21 CASA data.

Table 10.

Ring parameters for 2018 April 25 CASA data.

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.

thumbnail Fig. 14.

Comparison of Δlog(Z) for a series of geometric models ordered by model complexity. The Bayesian evidence for each model is evaluated with data generated from closure amplitude and phases. The number of parameters needed to define each model is given in parentheses. Circle markers denote ring-like models and hourglass markers denote other models. colors denote models with similar construction. The right panel shows a zoom in of the gray shaded region of the left panel.

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.

thumbnail Fig. 15.

Complex visibility (V) fit comparisons and normalized residuals for the Hybrid Themage and xs-ringauss, and the log-closure amplitude (lnAC) and closure phase (ψC) fit comparisons and normalized residuals for the mF-ring model. Filled points correspond to the real components of the complex visibility and open points correspond to the imaginary components. The colored points come from the maximum likelihood model for the Hybrid Themage and xs-ringauss, and from the maximum a posteriori model for the mF-ring. The gray points represent the scan-averaged data, with 1% fractional systematic noise added to all baselines.

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 ( 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.

thumbnail Fig. 16.

Fitted features of the xs-ringauss, mF-ring, and Hybrid Themage fiducial models to the 2018 M 87* HOPS data. We include fits from bands 1 through 4 on April 21 and April 25. Each point shows the median value of the posterior distribution and the error bars indicate the 68% posterior probability range centered around the median. The blue line and band represent the median and 68% confidence interval for the posterior generated by combining all bands and all days for the mF-ring model, while the red band is the equivalent for the xs-ringauss. The pink lines and bars represent the statistical mean and standard deviations of the Hybrid Themage method. The black line indicates the median over all days and bands of the 2017 M 87* analysis from xs-ringauss fits. The hashed region is the 68% posterior probability interval taken from the 2017 M 87* xs-ringauss fits.

Table 11.

Summary of the direct modeling parameters.

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.

thumbnail Fig. 17.

Representative band 3 images of M 87* on April 21 from each imaging pipeline (top row). Fiducial images are shown for DIFMAP and the RML methods. The DIFMAP image is restored with a circular 20 μas beam, as shown by the circle in the lower right corner. For THEMIS, Comrade, Hybrid Themage, the xs-ringauss, and the mF-ring a random sample drawn from the posterior is shown. Representative band 3 images of M 87* on April 21 from each imaging pipeline after blurring them with a circular Gaussian beam are shown on the bottom row – the FWHM of each beam is shown with the horizontal bar in the bottom right of each image. The eht-imaging, SMILI, THEMIS, Comrade, and Hybrid Themage images have been restored with 16.9, 16.1, 19.5, 18.8, and 14.2 μas FWHM Gaussian beams, respectively, to match the resolution of the DIFMAP reconstruction, whereas the xs-ringauss and mF-ring images were restored with a 20 μas FWHM Gaussian beam. The vertical dashed line separates the DIFMAP and RML methods from the Bayesian methods, and the solid vertical line separates the imaging methods which do not assume a ring-like structure, from the direct modeling methods which do assume a ring-like structure.

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 710 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 711). 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.

Table 12.

Comparison of the extracted ring parameters for 2017 and 2018.

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:

(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/α2017, 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/α2017 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.

thumbnail Fig. 18.

Comparisons of d/α2017, which serves as a proxy for θg, using multiple methods for 2017 (squares) and 2018 April 21 band 3 (diamonds). For the 2018 methods which have a corresponding θg calibration from the 2017 analysis (DIFMAP, eht-imaging, SMILI, xs-ringauss), we use the method-specific 2017 scaling factor to determine the 2018 d/α2017 values. For DIFMAP, eht-imaging, and SMILI we use the scaling factors from Table 6 of M 87 2017 VI, and for the xs-ringauss we use the scaling factors from Table 4 of M 87 2017 VI. For the 2018 methods which do not have a 2017 θg calibration (THEMIS, Comrade, mF-ring), we use α = 11.0, coming from the median α across all 2017 methods (M 87 2017 I). The error bars for the 2018 points are representative of the 68% confidence intervals of the model-specific diameter estimates. The two image domain feature extraction methods are shown with red points for REx and blue points for VIDA. For the 2017 points, we show the measured θg (black squares), the uncertainty due to differences in the 2017 observational details (σobs, solid black error bars), and the uncertainty due to the diversity of the 2017 GRMHD library (σthy, gray dashed error bars). The gray horizontal line and shaded region represent the 2017 θg value and σthy uncertainty for the xs-ringauss model. While the uncertainty in diameter for some of the 2018 methods is larger than that for the same methods in 2017, the uncertainty related to the different GRMHD models used as the calibration set dominates and spans all methods.

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 × 109M. 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 × 109M. 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 × 109M. 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 × 109M. These ∼5 × 109M 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.

thumbnail Fig. 19.

Comparison of the brightness position angle measured in the EHT observations during 2009−2018. The orange shadow covering (288 ± 10)° indicates the observational position angle range of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023). The red points for 2018 are the median values after combining the posteriors of all bands on April 21 and April 25 (see Table 11). The blue points for April 6 and April 11 in 2017 are adopted from M 87 2017 IV, M 87 2017 VI and Wielgus et al. (2020, see the Table 4). The gray points for 2009−2013 are adopted from Table 3 in Wielgus et al. (2020), and represent the 68% confidence intervals. The posterior shapes for the 2009−2013 are non-Gaussian, and exhibit large tails.

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 spins3, accretion flow types, and assumptions for the ratio between the electron and ion temperature in the simulation (see M 87 2017 V, for details).

thumbnail Fig. 20.

Distribution of best-fit position angle (in degrees) of the forward jet by fitting model images with the 2017 (April 6, high-band) and 2018 (April 21, band 4) observations. The best-fit 10% of images (solid lines) among all (∼18 000) images (dashed lines) are also shown. For reference, the vertical line shows the position angle ∼288° of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023), with the orange shadow area from (288 − 10)° to (288 + 10)°. While the fitted model images include different black hole spins, accretion types, and different electro-thermodynamics (M 87 2017 V; M 87 2017 VIII), the black hole spin vector is pointing away from Earth in all images. The 2018 EHT results are consistent with a black hole spin vector pointing away from Earth.

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.

thumbnail Fig. 21.

Model compact flux density estimates from the Bayesian Image reconstruction methods for the April 21 HOPS data. The error bars represent the 68% confidence intervals around the median (square points).

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, FOVy, 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.

thumbnail Fig. 22.

Triangle plot comparing the compact flux density (), raster FOV in the x and y directions, and crescent flux (IX) parameters for the Hybrid Themage reconstruction of April 21 band 1. The compact flux density and raster FOV parameters exhibit clearly bi-modal distributions in this parameter space. The color regions represent the 99%, 90%, and 50% quantile regions.

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 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.


1

In the CASA data, a priori corrections for the field rotation angle variations are applied upstream and the RCP-LCP phase offsets are solved as part of the instrumental calibration steps.

3

The dimensionless spin parameter a* and inclination angle i between the observer’s line of sight and the spin axis of the accretion flow are: a* = ( − 0.94, −0.5) for i = 17°, and a* = (0, 0.5, 0.94) for i = 163°.

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

  1. Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [Google Scholar]
  2. Akiyama, K., Lu, R.-S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
  3. Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, AJ, 153, 159 [NASA ADS] [CrossRef] [Google Scholar]
  4. Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, ApJ, 838, 1 [NASA ADS] [CrossRef] [Google Scholar]
  5. An, T., Sohn, B. W., & Imai, H. 2018, Nat. Astron., 2, 118 [Google Scholar]
  6. Arras, P., Frank, P., Haim, P., et al. 2022, Nat. Astron., 6, 259 [NASA ADS] [CrossRef] [Google Scholar]
  7. Asada, K., & Nakamura, M. 2012, ApJ, 745, L28 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bardeen, J. M. 1973, Les Astres Occlus, 215 [Google Scholar]
  9. Benkevitch, L., Akiyama, K., Lu, R., Doeleman, S., & Fish, V. 2016, ArXiv e-prints [arXiv:1609.00055] [Google Scholar]
  10. Bird, S., Harris, W. E., Blakeslee, J. P., & Flynn, C. 2010, A&A, 524, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Blackburn, L., Chan, C.-K., Crew, G. B., et al. 2019, ApJ, 882, 23 [Google Scholar]
  12. Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ, 894, 31 [Google Scholar]
  13. Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694, 556 [Google Scholar]
  14. Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
  15. Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
  16. Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bower, G. C., Dexter, J., Markoff, S., et al. 2015, ApJ, 811, L6 [NASA ADS] [CrossRef] [Google Scholar]
  18. Broderick, A. E., & Pesce, D. W. 2020, ApJ, 904, 126 [NASA ADS] [CrossRef] [Google Scholar]
  19. Broderick, A. E., Gold, R., Karami, M., et al. 2020a, ApJ, 897, 139 [Google Scholar]
  20. Broderick, A. E., Pesce, D. W., Tiede, P., Pu, H.-Y., & Gold, R. 2020b, ApJ, 898, 9 [Google Scholar]
  21. Broderick, A. E., Pesce, D. W., Gold, R., et al. 2022a, ApJ, 935, 61 [NASA ADS] [CrossRef] [Google Scholar]
  22. Broderick, A. E., Tiede, P., Pesce, D. W., & Gold, R. 2022b, ApJ, 927, 6 [NASA ADS] [CrossRef] [Google Scholar]
  23. Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Sci. Comput., 16, 1190 [Google Scholar]
  24. Cantiello, M., Blakeslee, J. P., Ferrarese, L., et al. 2018, ApJ, 856, 126 [Google Scholar]
  25. Carilli, C. L., & Thyagarajan, N. 2022, ApJ, 924, 125 [NASA ADS] [CrossRef] [Google Scholar]
  26. CASA Team (Bean, B., et al.) 2022, PASP, 134, 114501 [NASA ADS] [CrossRef] [Google Scholar]
  27. Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23 [Google Scholar]
  28. Chael, A., Bouman, K., Johnson, M., et al. 2019a, https://doi.org/10.5281/zenodo.2614016 [Google Scholar]
  29. Chael, A. A., Bouman, K. L., Johnson, M. D., et al. 2019b, Astrophysics Source Code Library [record ascl:1904.004] [Google Scholar]
  30. Chael, A., Johnson, M. D., & Lupsasca, A. 2021, ApJ, 918, 6 [NASA ADS] [CrossRef] [Google Scholar]
  31. Chen, M.-T., Asada, K., Matsushita, S., et al. 2023, PASP, 135, 095001 [NASA ADS] [CrossRef] [Google Scholar]
  32. Cho, I., Jung, T., Zhao, G.-Y., et al. 2017, PASJ, 69, 87 [NASA ADS] [CrossRef] [Google Scholar]
  33. Cho, I., Zhao, G.-Y., Kawashima, T., et al. 2022, ApJ, 926, 108 [NASA ADS] [CrossRef] [Google Scholar]
  34. Clark, B. G. 1980, A&A, 89, 377 [NASA ADS] [Google Scholar]
  35. Cui, Y.-Z., Hada, K., Kino, M., et al. 2021, Res. Astron. Astrophys., 21, 205 [CrossRef] [Google Scholar]
  36. Cui, Y., Hada, K., Kawashima, T., et al. 2023, Nature, 621, 711 [NASA ADS] [CrossRef] [Google Scholar]
  37. Cunningham, C. T., & Bardeen, J. M. 1973, ApJ, 183, 237 [Google Scholar]
  38. de Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275 [Google Scholar]
  40. Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
  41. Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Science, 338, 355 [NASA ADS] [CrossRef] [Google Scholar]
  42. EHT MWL Science Working Group (Algaba, J. C., et al.) 2021, ApJ, 911, L11 [NASA ADS] [CrossRef] [Google Scholar]
  43. Emami, R., Ricarte, A., Wong, G. N., et al. 2023, ApJ, 950, 38 [NASA ADS] [CrossRef] [Google Scholar]
  44. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [Google Scholar]
  45. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L2 [Google Scholar]
  46. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019c, ApJ, 875, L3 [Google Scholar]
  47. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019d, ApJ, 875, L4 [Google Scholar]
  48. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019e, ApJ, 875, L5 [Google Scholar]
  49. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019f, ApJ, 875, L6 [Google Scholar]
  50. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021a, ApJ, 910, L12 [Google Scholar]
  51. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021b, ApJ, 910, L13 [Google Scholar]
  52. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022a, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
  53. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022b, ApJ, 930, L13 [NASA ADS] [CrossRef] [Google Scholar]
  54. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022c, ApJ, 930, L14 [NASA ADS] [CrossRef] [Google Scholar]
  55. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022d, ApJ, 930, L15 [NASA ADS] [CrossRef] [Google Scholar]
  56. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022e, ApJ, 930, L16 [NASA ADS] [CrossRef] [Google Scholar]
  57. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022f, ApJ, 930, L17 [NASA ADS] [CrossRef] [Google Scholar]
  58. Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2023, ApJ, 957, L20 [NASA ADS] [CrossRef] [Google Scholar]
  59. Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13 [Google Scholar]
  60. Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P [Google Scholar]
  61. Ferreira, P. G. 2019, ARA&A, 57, 335 [Google Scholar]
  62. Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444 [Google Scholar]
  63. Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119 [Google Scholar]
  64. Georgiev, B., Pesce, D. W., Broderick, A. E., et al. 2022, ApJ, 930, L20 [NASA ADS] [CrossRef] [Google Scholar]
  65. Goddi, C., Martí-Vidal, I., Messias, H., et al. 2019, PASP, 131, 075003 [Google Scholar]
  66. Goddi, C., Martí-Vidal, I., Messias, H., et al. 2021, ApJ, 910, L14 [NASA ADS] [CrossRef] [Google Scholar]
  67. GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
  69. Greisen, E. W. 2003, Astrophys. Space Sci. Lib., 285, 109 [NASA ADS] [Google Scholar]
  70. Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131 [Google Scholar]
  71. Hilbert, D. 1917, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen – Mathematisch-physikalische Klasse (Berlin: Weidmannsche Buchhandlung), 53 [Google Scholar]
  72. Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593 [Google Scholar]
  73. Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
  74. Hovatta, T., Aller, M. F., Aller, H. D., et al. 2014, AJ, 147, 143 [Google Scholar]
  75. Inoue, M., Algaba-Marcos, J. C., Asada, K., et al. 2014, Radio Sci., 49, 564 [NASA ADS] [CrossRef] [Google Scholar]
  76. Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019a, A&A, 626, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Janssen, M., Blackburn, L., Issaoun, S., et al. 2019b, A Priori Calibration of EHT Stations, Tech. Rep. [Google Scholar]
  78. Janssen, M., Radcliffe, J. F., & Wagner, J. 2022, Universe, 8, 527 [NASA ADS] [CrossRef] [Google Scholar]
  79. Johnson, M. D., Lupsasca, A., Strominger, A., et al. 2020, Sci. Adv., 6, eaaz1310 [NASA ADS] [CrossRef] [Google Scholar]
  80. Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891 [NASA ADS] [CrossRef] [Google Scholar]
  81. Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018a, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Kim, J.-Y., Lee, S.-S., Hodgson, J. A., et al. 2018b, A&A, 610, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Koay, J. Y., Romero-Cañizales, C., Matthews, L. D., et al. 2023a, ArXiv e-prints [arXiv:2312.03505] [Google Scholar]
  84. Koay, J. Y., Asada, K., Matsushita, S., et al. 2023b, ArXiv e-prints [arXiv:2312.02759] [Google Scholar]
  85. Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
  86. Lee, S.-S., Wajima, K., Algaba, J.-C., et al. 2016, ApJS, 227, 8 [Google Scholar]
  87. Liepold, E. R., Ma, C.-P., & Walsh, J. L. 2023, ApJ, 945, L35 [NASA ADS] [CrossRef] [Google Scholar]
  88. Liu, D. C., & Nocedal, J. 1989, Math. Program., 45, 503 [CrossRef] [Google Scholar]
  89. Lockhart, W., & Gralla, S. E. 2022, MNRAS, 509, 3643 [Google Scholar]
  90. Lu, R., Asada, K., Krichbaum, T. P., et al. 2023, Nature, 616, 686 [CrossRef] [Google Scholar]
  91. Luminet, J.-P. 1979, A&A, 75, 228 [NASA ADS] [Google Scholar]
  92. Martí-Vidal, I., Roy, A., Conway, J., & Zensus, A. J. 2016, A&A, 587, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002 [Google Scholar]
  94. Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, A&A, 595, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Miyoshi, M., Kato, Y., & Makino, J. 2022, ApJ, 933, 36 [NASA ADS] [CrossRef] [Google Scholar]
  96. Mogensen, P. K., & Riseth, A. N. 2018, J. Open Source Softw., 3, 615 [NASA ADS] [CrossRef] [Google Scholar]
  97. Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146 [Google Scholar]
  98. Owen, F. N., Hardee, P. E., & Cornwell, T. J. 1989, ApJ, 340, 698 [NASA ADS] [CrossRef] [Google Scholar]
  99. Park, J., Hada, K., Kino, M., et al. 2019, ApJ, 887, 147 [Google Scholar]
  100. Pordes, R., Petravick, D., Kramer, B., et al. 2007, J. Phys. Conf. Ser., 78, 012057 [NASA ADS] [CrossRef] [Google Scholar]
  101. Readhead, A. C. S., Walker, R. C., Pearson, T. J., & Cohen, M. H. 1980, Nature, 285, 137 [NASA ADS] [CrossRef] [Google Scholar]
  102. Reid, M. J., Biretta, J. A., Junor, W., Muxlow, T. W. B., & Spencer, R. E. 1989, ApJ, 336, 112 [NASA ADS] [CrossRef] [Google Scholar]
  103. Rogers, A. E. E., Hinteregger, H. F., Whitney, A. R., et al. 1974, ApJ, 193, 293 [Google Scholar]
  104. Rogers, A. E. E., Doeleman, S. S., & Moran, J. M. 1995, AJ, 109, 1391 [NASA ADS] [CrossRef] [Google Scholar]
  105. Satapathy, K., Psaltis, D., Özel, F., et al. 2022, ApJ, 925, 13 [NASA ADS] [CrossRef] [Google Scholar]
  106. Schwarzschild, K. 1916, Abh. Konigl. Preuss. Akad. Wissenschaften Jahre, 189 [Google Scholar]
  107. Sfiligoi, I., Bradley, D. C., Holzman, B., et al. 2009, WRI World Congress on Computer Science and Information Engineering, 428 [Google Scholar]
  108. Shepherd, M. 2011, Astrophysics Source Code Library [record ascl:1103.001] [Google Scholar]
  109. Shepherd, M. C. 1997, ASP Conf. Ser., 125, 77 [Google Scholar]
  110. Simon, D. A., Cappellari, M., & Hartke, J. 2024, MNRAS, 527, 2341 [Google Scholar]
  111. Skilling, J. 2006, Bayesian Anal., 1, 833 [Google Scholar]
  112. Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
  113. Thompson, A. R., Moran, J. M., & Swenson, G. W. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd edn. (Springer) [CrossRef] [Google Scholar]
  114. Tiede, P. 2022, J. Open Source Softw., 7, 4457 [NASA ADS] [CrossRef] [Google Scholar]
  115. Tiede, P., Broderick, A. E., & Palumbo, D. C. M. 2022a, ApJ, 925, 122 [CrossRef] [Google Scholar]
  116. Tiede, P., Johnson, M. D., Pesce, D. W., et al. 2022b, Galaxies, 10, 111 [NASA ADS] [CrossRef] [Google Scholar]
  117. Tiede, P., Broderick, A. E., Palumbo, D. C. M., & Chael, A. 2022c, ApJ, 940, 182 [NASA ADS] [CrossRef] [Google Scholar]
  118. van Bemmel, I. M., Kettenis, M., Small, D., et al. 2022, PASP, 134, 114502 [CrossRef] [Google Scholar]
  119. Vertatschitsch, L., Primiani, R., Young, A., et al. 2015, PASP, 127, 1226 [NASA ADS] [CrossRef] [Google Scholar]
  120. von Laue, M. 1921, Relativitätstheorie (Braunschweig: Friedrich Vieweg& Sohn) [Google Scholar]
  121. Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [Google Scholar]
  122. Wielgus, M. 2021, Phys. Rev. D, 104, 124058 [NASA ADS] [CrossRef] [Google Scholar]
  123. Wielgus, M., Blackburn, L., Issaoun, S., et al. 2019, EHT Data Set Validation and Characterization of Errors, Tech. Rep. 2019-CE-02, EHT Memo Series [Google Scholar]
  124. Wielgus, M., Akiyama, K., Blackburn, L., et al. 2020, ApJ, 901, 67 [Google Scholar]
  125. Xu, K., Ge, H., Tebbutt, W., et al. 2020, Symposium on Advances in Approximate Bayesian Inference, PMLR, 1 [Google Scholar]
  126. Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Trans. Math. Softw., 23, 550 [Google Scholar]

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):

(A.1)

where 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 gE 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:

(A.2)

where Ageom is the geometrical area of the antenna dish in units of m2, and k = 1.38 × 103 is the Boltzmann constant in units of Jy m2/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. 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 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.

Table A.1.

EHT station flux density calibration parameters and their uncertainties for the 2018 observing campaign.

Appendix B: Non-closing systematic error budget

The closure phases (ψC), log closure amplitudes (lnAC), 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:

(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. mad0 is the modified median absolute deviation given by:

(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.

Table B.1.

Non-closing systematic uncertainties, s (and in units of thermal noise, s/σth), for M87* and its calibrator 3C 279 estimated using various statistical tests on both CASA and HOPS products.

thumbnail Fig. B.1.

Diagnostic plots showing the normalized distributions of various closure quantities of M87* and 3C 279 data from both the CASA (top four rows) and HOPS (bottom four rows) reduction pipelines. For each block corresponding to each pipeline, the first two rows show band 1 − band 2 (first column), band 3 − band 4 (second column), RR−LL (third column), and trivial (fourth column) closure quantities. The bottom two rows show band 1 − band 2 (first column), band 3 − band 4 (second column), and trivial (third column) closure trace quantities. Only data from April 21 and 25 have been used. The distributions prior to (blue) and after (red) accounting for the estimated systematic uncertainties, s, are shown. The values of s for each source and reduction pipeline are given in Table B.1. In the top left corner of each distribution, the number of > 3σ outliers are given considering thermal noise only, followed by the number of outliers considering thermal plus systematic noise for σ in parenthesis. These numbers are followed by the total number of data points after a slash.

Appendix C: Temporal coherence

We quantify the level of decoherence loss of the data calibrated using the HOPS and CASA pipelines using the ratio Ascan/A2s for each scan (M 87 2017 III). Ascan is the debiased amplitude of the coherently averaged visibilities over the entire scan, and A2s 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 Ascan/A2s 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.

thumbnail Fig. C.1.

Cumulative histograms of the amplitude ratios between coherent averaging for entire scans (Ascan), and coherent averaging for 2 s before incoherent averaging over scans (A2s), for both HOPS (blue) and CASA (orange) M87* data on both April 21 and April 25. The gray histogram shows the corresponding ratios from the HOPS data on the same days with no atmospheric phase corrections applied. For each pipeline, the fraction of data with coherence above 90% (i.e., with decoherence losses of no more than 10%) when averaged over the full scan is indicated.

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.

thumbnail Fig. E.1.

Left: (u, v) coverage for 3C 279 observations in band 3 on 2018 April 21 (colored points), showing the overlap between ALMA–LMT and ALMA–SPT baselines (highlighted with a box). Right: band 3 flux densities on ALMA–LMT (blue) and ALMA–SPT (orange) as a function of the projected baseline length in the east-west direction, in units of wavelength, demonstrating the under-calibration of amplitudes on LMT baselines. HOPS data are shown in this figure, but are consistent with that from CASA.

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 Ftot 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,

(F.1)

where V0 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:

(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

(F.3)

where Vi − j and Δgi 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

(F.4)

When we consider the gain uncertainties,

(F.5)

We assumed the uncertainty of the total flux Δgtot = 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

(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,

(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 ≤Fcpct ≤ 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) Ftot, 230 originates from the combination of (i) a compact emitter with unknown flux density Fcpct, 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 Fjet, 230 = Ftot, 230 − Fcpct, 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 Fjet and Ftot at lower frequencies, and extrapolating them to 230 GHz using a power-law model, we can estimate the total compact flux density Fcpct at 230 GHz: Fcpct, 230 = Ftot, 230 − Fjet, 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., Ftot = Fjet + Fcpct) 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 μas2 centered on the radio core. We obtain total flux densities of Ftot, 22 ≈ 1.1 Jy, Ftot, 43 ≈ 0.9 − 1.1 Jy, Ftot, 86 ≈ 0.9 Jy and Ftot, 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., Ftot ∝ ν+α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 Ftot, 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 μas2, which is nearly three orders of magnitude below the resolution limit of ALMA.

The compact flux density Fcpct, 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 Fcpct, 86 ≈ 0.51 Jy by integrating over a smaller window of 100 × 100 μas2 centered on the radio core. We can then derive a jet flux density at 86 GHz to be Fjet, 86 = Ftot, 86 − Fcpct, 86 ≈ 0.39 Jy.

Finally, we compute the expected Fcpct, 230 by estimating Fjet, 230 and subtracting it from the previously estimated Ftot, 230. As explained above, we assume that Fjet, 230 follows a power law that can be extrapolated up to 230 GHz (i.e., Fjet, 230 ∝ ν+αjet). For a range of αjet = −(0.7 − 1.0), we find Fjet, 230 = 0.15 − 0.20 Jy. This value gives a range of Fcpct, 230 ≈ 0.5 − 0.7 Jy within the central 100 × 100 μas2 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 Fcpct, 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 r0 = 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 r0 = 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 r0 = 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 r0 = 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.

thumbnail Fig. G.1.

Example geometric ring model used to generate the synthetic data. The left panel exhibits the morphological characteristics of the ring on a linear scale, encompassing a FOV of 130 μas. The right panel depicts the logarithmic scale representation of the extended jet model (FOV = 2900 μas).

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.

thumbnail Fig. G.2.

Four training geometric models as imaged by each imaging method. The first row shows the visibility amplitudes of the model (orange) compared to the visibility amplitudes measured for M87* (blue). The second row shows the ground-truth images. The DIFMAP, eht-imaging, and SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS and Comrade rows show a random draw from the posterior. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 μas. The white lines in the THEMIS and Comrade panels represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution.

thumbnail Fig. G.3.

Seven geometric validation models plus one GRMHD snapshot as imaged by each method. The first row shows the visibility amplitudes of the model (orange) compared to the visibility amplitudes measured for M87* (blue). The second row shows the ground-truth images. The DIFMAP, eht-imaging, and SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS and Comrade rows show a random draw from the posterior. The THEMIS reconstruction of the GRMHD synthetic data uses a Hybrid Themage model. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 μas. The white lines in the THEMIS and Comrade panels represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution.

thumbnail Fig. G.4.

Stacked images with a dynamic range cutoff used to investigate the structural deviations in the Top Set (from SMILI as a representative). First, each Top Set image is normalized with its peak intensity, and the pixels with a flux larger than a certain threshold (e.g., 0.1) are converted to unity (unless 0). Then all the images across the Top Set are stacked and normalized by the maximum number of Top Set images so that the fraction of reliable features over the entire Top Set is emphasized. The ground-truth structure of each synthetic data model is shown in gray contours.

thumbnail Fig. G.5.

Multiplicative gain correction factors at each station as a function of time. The ground truth gains (gray circles) are used to generate the synthetic data with the cres180 model for April 21, band 3. Model gains from each pipeline are derived from the fiducial images or posterior samples of the cres180 synthetic data reconstructions.

thumbnail Fig. G.6.

Multiplicative gain correction factors at each station as a function of time from the fiducual images or posterior samples of the image reconstructions for the M87* HOPS band 3 data.

thumbnail Fig. G.7.

Twenty five randomly selected Top Set images from HOPS band 3 April 21 data for DIFMAP (upper left), eht-imaging (upper right), and SMILI (bottom).

thumbnail Fig. G.8.

Example reconstructions of M87* for 2018 April 21 and 2017 April 11 after omitting visibilities to each geographical site. The images presented show reconstructions that exclude all baselines to the indicated site (i.e., mimicking an observation without that site). Top and middle rows show the reconstructions for Comrade and eht-imaging respectively for band 3 on 2018 April 21. The images for the eht-imaging pipeline were reconstructed using the eht-imaging fiducial parameters (see Sect. 5.2). The images shown for Comrade are the mean images from the posterior. The bottom row shows the reconstructions for 2017 April 11 using the 2018 eht-imaging pipeline but 2017 fiducial parameters (M 87 2017 IV). The ellipse in each panel shows the corresponding synthesized beam with uniform weighting, but the image is not convolved with this beam.

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 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 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.

thumbnail Fig. H.1.

Fiducial images produced by the DIFMAP pipeline reconstructed through the implementation of an initial phase-only self-calibration process using a point source model. The reconstructed images exhibit ring-like structures, which are comparable to the results obtained from the pipeline that utilized closure phases to select the optimal geometric model for the initial phase-only self-calibration, as demonstrated in Fig. 7. The circle represents a Gaussian blurring kernel with a FWHM of 20 μas.

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.

thumbnail Fig. I.1.

Fiducual images of the 2017 April 11 low-band data obtained from the 2018 DIFMAP pipelines discussed in Sect. 5.1.1, as a backwards compatibility check. The top row of images is obtained using the pipeline that conducts initial phase-only self-calibration by employing a point source model, whereas the bottom row of images is obtained from the pipeline that identifies the optimal geometric model for self-calibration. The circles represent a Gaussian blurring kernel with a FWHM of 20 μas.

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, IM, N. The intensity map is produced using an approximate cubic spline interpolation between the control points. The field of view FOVx and FOVy 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 IG, a symmetric standard deviation σG, asymmetry parameter AG, position angle ϕG, and position (xG, yG). 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 IX, outer radius Rout, linear brightness gradient fX, brightness gradient position angle ϕX, and position (xX, yX) 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 , auto-correlation time, and visual inspection of the parameter traces. The parameters and priors for each model component are listed in Table J.1.

Table J.1.

THEMIS Raster Model Priors.

Appendix K: THEMIS Raster dimensions

Table K.1.

THEMIS raster size (Nx × Ny) survey for April 21 Band 3.

The best fit raster dimensions Nx × Ny 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 χ2 for the surveyed raster sizes, and note that all models provide good fits to the data.

Table K.2.

Complex visibility reduced χ2 () for the THEMIS Raster and Hybrid Themage model best fits to the HOPS pipeline 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 χ2 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 (li, mj) given by Eq. L.3, for a given FOV and number of pixels Nx and Ny for each side. The FOV and number of pixels are the hyperparameters of the model:

(L.1)

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

(L.2)

and raster points are given by,

(L.3)

where sx = FOV/Nx and sy = FOV/Ny. We restrict the images considered in this paper to having Nx = Ny, such that the images are square and sx = sy. We add a circular Gaussian with FWHM = mas (G1mas) 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,

(L.4)

where f is the total flux density of the raster (also described as compact flux) and fg is the fraction of the total flux density that corresponds to the flux density of the Gaussian. We take the Fourier transform of M(Iij, f, fg) to get model Fourier transform for each baseline ab. The Fourier amplitudes are then corrupted to model the individual station gain amplitudes |g| to construct the visibility amplitudes,

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 fg for the flux density of the large-scale Gaussian. For the raster pixels fluxes, Ii, we use a symmetric Dirichlet distribution

(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 ∑iIi = 1 and 0 ≤ Ii ≤ 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 (Fcpct) 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 Fcpct, while all other imaging parameters are kept at their fiducial values. We select acceptable images based on two criteria: (1) normalized ,   (with a 10% systematic uncertainty for DIFMAP), and (2) a corresponding self-calibration solution with 0.9 ≤ median(1/|gSMT|) ≤ 1.1, where gSMT 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 Fcpct, the ring diameter and position angle values of the acceptable images do not exhibit significant dependence on Fcpct. 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.

thumbnail Fig. M.1.

Measured diameter d, width w, position angle, and fractional central brightness fC of M87*, measured from image reconstructions assuming different total compact flux densities, Fcpct for DIFMAP (orange), eht-imaging (blue), and SMILI (green) (see Sect. 5.3). All measurements are made using REx. All other imaging parameters were set to the fiducial parameters of the corresponding pipeline. DIFMAP values were measured after restoring with a 20 μas FWHM Gaussian beam. The solid lines indicate the measured value, and the shaded regions give the ±1σ uncertainty of the REx measurement. The darker colored regions correspond to values of Fcpct that produces images that have (1) normalized ,   (with 10% systematic uncertainty for DIFMAP), and (2) a corresponding self-calibration solution with 0.9 ≤ median(1/|gSMT|) ≤ 1.1.

thumbnail Fig. M.2.

Reconstructed DIFMAP images of the four geometric models with different values of the total compact flux density. We use the fiducial parameters from the 2018 Top Set, but with different assumed total compact flux densities. The top label above each column indicates the assumed total compact flux density in units of Jy.

Appendix N: Geometric modeling details

thumbnail Fig. N.1.

Schematic summary of the xs-ringauss and mF-ring model parameters. The xs-ringauss is defined to be a disk with a circular hole removed, that is allowed to be offset from the center of the disk and rotated. A circular Guassian is pinned to the center of the disk as an emission floor, and an additional elliptical Gaussian is pinned to the inner edge of the offset hole. The mF-ring is an m-ring of order 2 that has an elliptical stretch which is allowed to be rotated. A flat emission floor is included to match the interior of the m-ring. Both models have an additional blur defined on them with a Gaussian blurring kernel. The intensity profiles of the fiducial models along two axes are indicated in red and green. An angular intensity profile is also shown for the mF-ring in blue.

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.

Table N.1.

THEMIS xs-ringauss model parameters and priors.

Table N.2.

mF-ring parameters and priors.

thumbnail Fig. N.2.

Representative images for each observing band from random posterior samples for the Hybrid Themage, xs-ringauss, and mF-ring models, for the April 21 and April 25 HOPS data.

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,

(N.1)

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

(N.2)

the prior probability on model parameters,

(N.3)

and the Bayesian evidence for a chosen model,

(N.4)

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

(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,

(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,

(N.7)

Here, the debiased radius () is linked to the defining radius of the mF-ring (ri), and the FWHM of the blurring Gaussian kernel (σ*) through:

(N.8)

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

(N.9)

and the brightness asymmetry position angle to be,

(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,

(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

Table 1.

Median zenith sky opacities (1.3 mm) at EHT sites during the 2018 April observations toward M 87*.

Table 2.

Parameter survey results for April 21 band 3 data.

Table 3.

Initial DIFMAP geometric models.

Table 4.

Comrade raster model hyperparameters.

Table 5.

Fiducual image and Top Set closure quantity normalized χ2 values.

Table 6.

Reduced χ2 quantities for the THEMIS and Comrade raster models.

Table 7.

Ring parameters for 2018 April 21 HOPS data.

Table 8.

Ring parameters for the 2018 April 25 HOPS data.

Table 9.

Ring parameters for 2018 April 21 CASA data.

Table 10.

Ring parameters for 2018 April 25 CASA data.

Table 11.

Summary of the direct modeling parameters.

Table 12.

Comparison of the extracted ring parameters for 2017 and 2018.

Table A.1.

EHT station flux density calibration parameters and their uncertainties for the 2018 observing campaign.

Table B.1.

Non-closing systematic uncertainties, s (and in units of thermal noise, s/σth), for M87* and its calibrator 3C 279 estimated using various statistical tests on both CASA and HOPS products.

Table J.1.

THEMIS Raster Model Priors.

Table K.1.

THEMIS raster size (Nx × Ny) survey for April 21 Band 3.

Table K.2.

Complex visibility reduced χ2 () for the THEMIS Raster and Hybrid Themage model best fits to the HOPS pipeline data.

Table N.1.

THEMIS xs-ringauss model parameters and priors.

Table N.2.

mF-ring parameters and priors.

All Figures

thumbnail Fig. 1.

Representative example images of M 87* from the EHT observations taken on 2017 April 11 and 2018 April 21 (north is up and east is to the left). The 2017 image is generated with the average of fiducial parameter sets from the imaging techniques used in M 87 2017 IV. The 2018 image is created by taking the average of the blurred images generated by the imaging techniques found in Sect. 5. Comparison of the images shows consistency in the diameter across observation epochs, but a shift in position angle of brightness asymmetry. The circle represents a Gaussian blurring kernel with a full width half maximum of 20 μas.

In the text
thumbnail Fig. 2.

Map showing the stations that participated in the EHT 2018 campaign (black circles), which differs from the EHT array in 2017 by the addition of the GLT. Co-located sites in Chile and Hawai‘i appear superimposed. The SPT projected location from the back of the map is indicated with a dashed circle, and baselines to this station are represented with dashed-lines. While the SPT cannot observe M 87*, it observed 3C 279 and was used to calibrate the data.

In the text
thumbnail Fig. 3.

EHT observing schedules for M 87* (blue) and its calibrator 3C 279 (orange) on the 2018 April 21 (left panel) and April 25 (right panel) observing days, which began at the end in UTC of April 20 and 24, respectively. Open rectangles represent scans that were scheduled but not observed owing to weather or technical issues. The filled rectangles mark the scans with detections in the final data set. On these two particular days, scan durations are typically 4 to 5 min each for M 87* and 4 min each for 3C 279, as reflected by the width of each rectangle.

In the text
thumbnail Fig. 4.

M 87* (u, v) coverage (colored points) in band 3 (top panels) and band 2 (bottom panels) for observations on 2018 April 21 (left panels) and April 25 (right panels), overlaid on the low-band (u, v) coverage for 2017 April 11 (gray points). The dashed circles show baseline lengths corresponding to fringe spacings of 25 and 50 μas respectively. The (u, v) coverage in band 1 and band 4 is comparable to that in band 2 and band 3, respectively.

In the text
thumbnail Fig. 5.

Measured correlated flux densities of M 87* as a function of baseline lengths in units of wavelength, for 2018 April 21 (top panels) and April 25 (bottom panels) in band 3, for both HOPS (left panels) and CASA (right panels) outputs. The 2018 data (colored points) are overlaid on the corresponding flux densities of the 2017 April 11 observations in low-band (gray points). All data shown include a priori station-based amplitude calibration and network calibration but are prior to any model-dependent self-calibration. Error bars denote ±1σ from thermal noise. Redundant baselines are shown with different symbols: circles for baselines to ALMA and SMA; diamonds for baselines to APEX and JCMT. The dashed line corresponds to an azimuthally symmetric thin ring model with a 42 μas diameter.

In the text
thumbnail Fig. 6.

Variation of closure phases as a function of Greenwich mean sidereal time (GMST) for selected triangles using HOPS calibrated data. Red diamonds denote data from 2017 April 6 (low-band), gray circles denote data from 2017 April 11 (low-band), and blue squares denote data from 2018 April 21 (band 3). Error bars are the ±1σ uncertainties.

In the text
thumbnail Fig. 7.

Representative images recovered from the HOPS and CASA data with all five imaging pipelines for two observing days (April 21 and 25). Each panel shows the fiducial image of the corresponding top set images for the DIFMAP, eht-imaging, and SMILI pipelines, and a random sample from the respective posterior for the THEMIS and Comrade pipelines. We do not have Top Sets for band 1 and band 2 from DIFMAP, eht-imaging, and SMILI pipelines on April 25. The dashed horizontal line in each block separates the DIFMAP and RML methods above from the Bayesian methods below. The circles in the DIFMAP panels represent an effective Gaussian blurring kernel of 20 μas. The solid lines in the THEMIS and Comrade panels represent the size of the blurring kernel used to achieve the same effective resolution as the DIFMAP method.

In the text
thumbnail Fig. 8.

Closure phases plotted as a function of GMST on three selected triangles from the April 21 band 3 observations (black points). The error bars on the data points denote the ±1σ uncertainties. The colored and dashed lines indicate the model closure phase curves from the fiducial images or posterior samples produced by the five imaging pipelines.

In the text
thumbnail Fig. 9.

Visibility amplitudes of band 3 data on April 21 as a function of baseline length compared with corresponding gain-calibrated visibility amplitudes from the five representative image models from each pipeline. The fiducial images are used for DIFMAP, eht-imaging, and SMILI, and the maximum likelihood model from the sampling chains are used for THEMIS and maximum a posteriori model for Comrade. The gray points represent the data used by each method after flagging, but before individual self-calibration. eht-imaging scales the intra-site baselines to 0.5 Jy, and SMILI scales the intra-site baselines to 0.6 Jy, and both pre-calibrate the LMT-SMT baselines before fitting. DIFMAP uses 10 s averaged data, and all other methods use scan-averaged data. THEMIS and Comrade apply 1% fractional systematic noise to all baselines, and eht-imaging and SMILI apply 1% and 0% fractional systematic noise respectively to all baselines, seen as minor differences in the gray points. The colored points represent the image model visibilities from each method, with station gains derived from each method’s internal self-calibration procedure applied to the image model visibilities. Below each visibility amplitude figure are the normalized residuals for each image, which is the difference between the gray and the colored points, divided by the uncertainty of the gray data points.

In the text
thumbnail Fig. 10.

Visualization of image statistics calculated using the Top Set images from the eht-imaging pipeline for observations taken on April 21 band 3. We emphasize that these images do not represent the posterior probability space for the reconstructions. Each image reconstructed using eht-imaging is the maximum a posteriori (MAP) image for a given parameter set. Thus, the statistics shown represent uncertainties that arise from different choices of regularizer weights, not from an exploration of posterior space. The top row shows top statistics in the image domain while the bottom row shows the visibility domain. Overlaid on the visibility domain panels is the (u, v) coverage for the April 21 observation. From left to right, we present the mean image; the standard deviation; the normalized standard deviation, calculated by rescaling each image to the flux of the mean image; and the fractional standard deviation, calculated by dividing the standard deviation by the mean. The fractional standard deviation panel has been clipped to a maximum value of 1. Portions of the image exhibit large fractional standard deviations due to pixel values very close to zero in the mean image. In the top row, image contours are drawn at 10%, 20%, 40%, and 80% of the peak values from the mean image. In the bottom row, the gray contours represent 0.1%, 1%, and 10% of the peak while the black contours represent 10 and 100 mJy (left three panels) and 0.1 (right most panel). The complex visibilities are calculated by taking a Fourier transform of the images and then calculating the mean and standard deviation. The absolute value of the mean and standard deviation of the complex visibilities are used to calculate visibility amplitudes.

In the text
thumbnail Fig. 11.

Visualization of image uncertainties using images from the posterior of THEMIS pipeline for observations taken with band 3 on April 21. The contour lines shown are drawn as described in Fig. 10.

In the text
thumbnail Fig. 12.

Visualization of image uncertainties using images from the posterior of Comrade pipeline for observations taken with band 3 on April 21. The contour lines shown are the same as described in Fig. 10.

In the text
thumbnail Fig. 13.

Ring characteristics from all bands, observational days, and imaging pipelines coming form HOPS and CASA data (colored points). Median values and 68% confidence intervals are shown for diameter and width. Circular mean and standard deviation values are shown for the position angle. Circle and triangle markers correspond to REx and VIDA respectively. Red (eht-imaging) and blue (DIFMAP) dashed lines and shaded regions show the ring parameters and 68% confidence intervals from the 2017 April 11 measurements. The vertical gray shaded region corresponds to the 2018 April 21 measurements, and the vertical un-shaded region corresponds to the 2018 April 25 measurements.

In the text
thumbnail Fig. 14.

Comparison of Δlog(Z) for a series of geometric models ordered by model complexity. The Bayesian evidence for each model is evaluated with data generated from closure amplitude and phases. The number of parameters needed to define each model is given in parentheses. Circle markers denote ring-like models and hourglass markers denote other models. colors denote models with similar construction. The right panel shows a zoom in of the gray shaded region of the left panel.

In the text
thumbnail Fig. 15.

Complex visibility (V) fit comparisons and normalized residuals for the Hybrid Themage and xs-ringauss, and the log-closure amplitude (lnAC) and closure phase (ψC) fit comparisons and normalized residuals for the mF-ring model. Filled points correspond to the real components of the complex visibility and open points correspond to the imaginary components. The colored points come from the maximum likelihood model for the Hybrid Themage and xs-ringauss, and from the maximum a posteriori model for the mF-ring. The gray points represent the scan-averaged data, with 1% fractional systematic noise added to all baselines.

In the text
thumbnail Fig. 16.

Fitted features of the xs-ringauss, mF-ring, and Hybrid Themage fiducial models to the 2018 M 87* HOPS data. We include fits from bands 1 through 4 on April 21 and April 25. Each point shows the median value of the posterior distribution and the error bars indicate the 68% posterior probability range centered around the median. The blue line and band represent the median and 68% confidence interval for the posterior generated by combining all bands and all days for the mF-ring model, while the red band is the equivalent for the xs-ringauss. The pink lines and bars represent the statistical mean and standard deviations of the Hybrid Themage method. The black line indicates the median over all days and bands of the 2017 M 87* analysis from xs-ringauss fits. The hashed region is the 68% posterior probability interval taken from the 2017 M 87* xs-ringauss fits.

In the text
thumbnail Fig. 17.

Representative band 3 images of M 87* on April 21 from each imaging pipeline (top row). Fiducial images are shown for DIFMAP and the RML methods. The DIFMAP image is restored with a circular 20 μas beam, as shown by the circle in the lower right corner. For THEMIS, Comrade, Hybrid Themage, the xs-ringauss, and the mF-ring a random sample drawn from the posterior is shown. Representative band 3 images of M 87* on April 21 from each imaging pipeline after blurring them with a circular Gaussian beam are shown on the bottom row – the FWHM of each beam is shown with the horizontal bar in the bottom right of each image. The eht-imaging, SMILI, THEMIS, Comrade, and Hybrid Themage images have been restored with 16.9, 16.1, 19.5, 18.8, and 14.2 μas FWHM Gaussian beams, respectively, to match the resolution of the DIFMAP reconstruction, whereas the xs-ringauss and mF-ring images were restored with a 20 μas FWHM Gaussian beam. The vertical dashed line separates the DIFMAP and RML methods from the Bayesian methods, and the solid vertical line separates the imaging methods which do not assume a ring-like structure, from the direct modeling methods which do assume a ring-like structure.

In the text
thumbnail Fig. 18.

Comparisons of d/α2017, which serves as a proxy for θg, using multiple methods for 2017 (squares) and 2018 April 21 band 3 (diamonds). For the 2018 methods which have a corresponding θg calibration from the 2017 analysis (DIFMAP, eht-imaging, SMILI, xs-ringauss), we use the method-specific 2017 scaling factor to determine the 2018 d/α2017 values. For DIFMAP, eht-imaging, and SMILI we use the scaling factors from Table 6 of M 87 2017 VI, and for the xs-ringauss we use the scaling factors from Table 4 of M 87 2017 VI. For the 2018 methods which do not have a 2017 θg calibration (THEMIS, Comrade, mF-ring), we use α = 11.0, coming from the median α across all 2017 methods (M 87 2017 I). The error bars for the 2018 points are representative of the 68% confidence intervals of the model-specific diameter estimates. The two image domain feature extraction methods are shown with red points for REx and blue points for VIDA. For the 2017 points, we show the measured θg (black squares), the uncertainty due to differences in the 2017 observational details (σobs, solid black error bars), and the uncertainty due to the diversity of the 2017 GRMHD library (σthy, gray dashed error bars). The gray horizontal line and shaded region represent the 2017 θg value and σthy uncertainty for the xs-ringauss model. While the uncertainty in diameter for some of the 2018 methods is larger than that for the same methods in 2017, the uncertainty related to the different GRMHD models used as the calibration set dominates and spans all methods.

In the text
thumbnail Fig. 19.

Comparison of the brightness position angle measured in the EHT observations during 2009−2018. The orange shadow covering (288 ± 10)° indicates the observational position angle range of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023). The red points for 2018 are the median values after combining the posteriors of all bands on April 21 and April 25 (see Table 11). The blue points for April 6 and April 11 in 2017 are adopted from M 87 2017 IV, M 87 2017 VI and Wielgus et al. (2020, see the Table 4). The gray points for 2009−2013 are adopted from Table 3 in Wielgus et al. (2020), and represent the 68% confidence intervals. The posterior shapes for the 2009−2013 are non-Gaussian, and exhibit large tails.

In the text
thumbnail Fig. 20.

Distribution of best-fit position angle (in degrees) of the forward jet by fitting model images with the 2017 (April 6, high-band) and 2018 (April 21, band 4) observations. The best-fit 10% of images (solid lines) among all (∼18 000) images (dashed lines) are also shown. For reference, the vertical line shows the position angle ∼288° of the milliarcsecond-scale jet (Walker et al. 2018; Cui et al. 2023), with the orange shadow area from (288 − 10)° to (288 + 10)°. While the fitted model images include different black hole spins, accretion types, and different electro-thermodynamics (M 87 2017 V; M 87 2017 VIII), the black hole spin vector is pointing away from Earth in all images. The 2018 EHT results are consistent with a black hole spin vector pointing away from Earth.

In the text
thumbnail Fig. 21.

Model compact flux density estimates from the Bayesian Image reconstruction methods for the April 21 HOPS data. The error bars represent the 68% confidence intervals around the median (square points).

In the text
thumbnail Fig. 22.

Triangle plot comparing the compact flux density (), raster FOV in the x and y directions, and crescent flux (IX) parameters for the Hybrid Themage reconstruction of April 21 band 1. The compact flux density and raster FOV parameters exhibit clearly bi-modal distributions in this parameter space. The color regions represent the 99%, 90%, and 50% quantile regions.

In the text
thumbnail Fig. B.1.

Diagnostic plots showing the normalized distributions of various closure quantities of M87* and 3C 279 data from both the CASA (top four rows) and HOPS (bottom four rows) reduction pipelines. For each block corresponding to each pipeline, the first two rows show band 1 − band 2 (first column), band 3 − band 4 (second column), RR−LL (third column), and trivial (fourth column) closure quantities. The bottom two rows show band 1 − band 2 (first column), band 3 − band 4 (second column), and trivial (third column) closure trace quantities. Only data from April 21 and 25 have been used. The distributions prior to (blue) and after (red) accounting for the estimated systematic uncertainties, s, are shown. The values of s for each source and reduction pipeline are given in Table B.1. In the top left corner of each distribution, the number of > 3σ outliers are given considering thermal noise only, followed by the number of outliers considering thermal plus systematic noise for σ in parenthesis. These numbers are followed by the total number of data points after a slash.

In the text
thumbnail Fig. C.1.

Cumulative histograms of the amplitude ratios between coherent averaging for entire scans (Ascan), and coherent averaging for 2 s before incoherent averaging over scans (A2s), for both HOPS (blue) and CASA (orange) M87* data on both April 21 and April 25. The gray histogram shows the corresponding ratios from the HOPS data on the same days with no atmospheric phase corrections applied. For each pipeline, the fraction of data with coherence above 90% (i.e., with decoherence losses of no more than 10%) when averaged over the full scan is indicated.

In the text
thumbnail Fig. E.1.

Left: (u, v) coverage for 3C 279 observations in band 3 on 2018 April 21 (colored points), showing the overlap between ALMA–LMT and ALMA–SPT baselines (highlighted with a box). Right: band 3 flux densities on ALMA–LMT (blue) and ALMA–SPT (orange) as a function of the projected baseline length in the east-west direction, in units of wavelength, demonstrating the under-calibration of amplitudes on LMT baselines. HOPS data are shown in this figure, but are consistent with that from CASA.

In the text
thumbnail Fig. G.1.

Example geometric ring model used to generate the synthetic data. The left panel exhibits the morphological characteristics of the ring on a linear scale, encompassing a FOV of 130 μas. The right panel depicts the logarithmic scale representation of the extended jet model (FOV = 2900 μas).

In the text
thumbnail Fig. G.2.

Four training geometric models as imaged by each imaging method. The first row shows the visibility amplitudes of the model (orange) compared to the visibility amplitudes measured for M87* (blue). The second row shows the ground-truth images. The DIFMAP, eht-imaging, and SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS and Comrade rows show a random draw from the posterior. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 μas. The white lines in the THEMIS and Comrade panels represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution.

In the text
thumbnail Fig. G.3.

Seven geometric validation models plus one GRMHD snapshot as imaged by each method. The first row shows the visibility amplitudes of the model (orange) compared to the visibility amplitudes measured for M87* (blue). The second row shows the ground-truth images. The DIFMAP, eht-imaging, and SMILI rows show a fiducial image made from the same parameter sets as the images shown in Fig. 17. The THEMIS and Comrade rows show a random draw from the posterior. The THEMIS reconstruction of the GRMHD synthetic data uses a Hybrid Themage model. The circle in the DIFMAP panels represent a Gaussian blurring kernel of 20 μas. The white lines in the THEMIS and Comrade panels represent the size of FWHM for the Gaussian blurring kernel needed to match the DIFMAP effective resolution.

In the text
thumbnail Fig. G.4.

Stacked images with a dynamic range cutoff used to investigate the structural deviations in the Top Set (from SMILI as a representative). First, each Top Set image is normalized with its peak intensity, and the pixels with a flux larger than a certain threshold (e.g., 0.1) are converted to unity (unless 0). Then all the images across the Top Set are stacked and normalized by the maximum number of Top Set images so that the fraction of reliable features over the entire Top Set is emphasized. The ground-truth structure of each synthetic data model is shown in gray contours.

In the text
thumbnail Fig. G.5.

Multiplicative gain correction factors at each station as a function of time. The ground truth gains (gray circles) are used to generate the synthetic data with the cres180 model for April 21, band 3. Model gains from each pipeline are derived from the fiducial images or posterior samples of the cres180 synthetic data reconstructions.

In the text
thumbnail Fig. G.6.

Multiplicative gain correction factors at each station as a function of time from the fiducual images or posterior samples of the image reconstructions for the M87* HOPS band 3 data.

In the text
thumbnail Fig. G.7.

Twenty five randomly selected Top Set images from HOPS band 3 April 21 data for DIFMAP (upper left), eht-imaging (upper right), and SMILI (bottom).

In the text
thumbnail Fig. G.8.

Example reconstructions of M87* for 2018 April 21 and 2017 April 11 after omitting visibilities to each geographical site. The images presented show reconstructions that exclude all baselines to the indicated site (i.e., mimicking an observation without that site). Top and middle rows show the reconstructions for Comrade and eht-imaging respectively for band 3 on 2018 April 21. The images for the eht-imaging pipeline were reconstructed using the eht-imaging fiducial parameters (see Sect. 5.2). The images shown for Comrade are the mean images from the posterior. The bottom row shows the reconstructions for 2017 April 11 using the 2018 eht-imaging pipeline but 2017 fiducial parameters (M 87 2017 IV). The ellipse in each panel shows the corresponding synthesized beam with uniform weighting, but the image is not convolved with this beam.

In the text
thumbnail Fig. H.1.

Fiducial images produced by the DIFMAP pipeline reconstructed through the implementation of an initial phase-only self-calibration process using a point source model. The reconstructed images exhibit ring-like structures, which are comparable to the results obtained from the pipeline that utilized closure phases to select the optimal geometric model for the initial phase-only self-calibration, as demonstrated in Fig. 7. The circle represents a Gaussian blurring kernel with a FWHM of 20 μas.

In the text
thumbnail Fig. I.1.

Fiducual images of the 2017 April 11 low-band data obtained from the 2018 DIFMAP pipelines discussed in Sect. 5.1.1, as a backwards compatibility check. The top row of images is obtained using the pipeline that conducts initial phase-only self-calibration by employing a point source model, whereas the bottom row of images is obtained from the pipeline that identifies the optimal geometric model for self-calibration. The circles represent a Gaussian blurring kernel with a FWHM of 20 μas.

In the text
thumbnail Fig. M.1.

Measured diameter d, width w, position angle, and fractional central brightness fC of M87*, measured from image reconstructions assuming different total compact flux densities, Fcpct for DIFMAP (orange), eht-imaging (blue), and SMILI (green) (see Sect. 5.3). All measurements are made using REx. All other imaging parameters were set to the fiducial parameters of the corresponding pipeline. DIFMAP values were measured after restoring with a 20 μas FWHM Gaussian beam. The solid lines indicate the measured value, and the shaded regions give the ±1σ uncertainty of the REx measurement. The darker colored regions correspond to values of Fcpct that produces images that have (1) normalized ,   (with 10% systematic uncertainty for DIFMAP), and (2) a corresponding self-calibration solution with 0.9 ≤ median(1/|gSMT|) ≤ 1.1.

In the text
thumbnail Fig. M.2.

Reconstructed DIFMAP images of the four geometric models with different values of the total compact flux density. We use the fiducial parameters from the 2018 Top Set, but with different assumed total compact flux densities. The top label above each column indicates the assumed total compact flux density in units of Jy.

In the text
thumbnail Fig. N.1.

Schematic summary of the xs-ringauss and mF-ring model parameters. The xs-ringauss is defined to be a disk with a circular hole removed, that is allowed to be offset from the center of the disk and rotated. A circular Guassian is pinned to the center of the disk as an emission floor, and an additional elliptical Gaussian is pinned to the inner edge of the offset hole. The mF-ring is an m-ring of order 2 that has an elliptical stretch which is allowed to be rotated. A flat emission floor is included to match the interior of the m-ring. Both models have an additional blur defined on them with a Gaussian blurring kernel. The intensity profiles of the fiducial models along two axes are indicated in red and green. An angular intensity profile is also shown for the mF-ring in blue.

In the text
thumbnail Fig. N.2.

Representative images for each observing band from random posterior samples for the Hybrid Themage, xs-ringauss, and mF-ring models, for the April 21 and April 25 HOPS data.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.