Issue 
A&A
Volume 681, January 2024



Article Number  A79  
Number of page(s)  63  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202347932  
Published online  18 January 2024 
The persistent shadow of the supermassive black hole of M 87
I. Observations, calibration, imaging, and analysis^{⋆}
^{1}
Massachusetts Institute of Technology Haystack Observatory, 99 Millstone Road, Westford, MA, 01886, USA
^{2}
National Astronomical Observatory of Japan, 2211 Osawa, Mitaka, Tokyo 1818588, Japan
^{3}
Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
^{4}
Instituto de Astrofísica de AndalucíaCSIC, Glorieta de la Astronomía s/n, 18008 Granada, Spain
^{5}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{6}
Department of Physics, Faculty of Science, Universiti Malaya, 50603 Kuala Lumpur, Malaysia
^{7}
Center for Astrophysics  Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA
^{8}
Department of Physics & Astronomy, The University of Texas at San Antonio, One UTSA Circle, San Antonio, TX 78249, USA
^{9}
Institute of Astronomy and Astrophysics, Academia Sinica, 11F of AstronomyMathematics Building, AS/NTU No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan ROC
^{10}
Departament d’Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50, 46100 Burjassot, València, Spain
^{11}
Observatori Astronòmic, Universitat de València, C. Catedrático José Beltrán 2, 46980 Paterna, València, Spain
^{12}
Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 43992 Onsala, Sweden
^{13}
Steward Observatory and Department of Astronomy, University of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA
^{14}
Yale Center for Astronomy & Astrophysics, Yale University, 52 Hillhouse Avenue, New Haven, CT 06511, USA
^{15}
Astronomy Department, Universidad de Concepción, Casilla 160, Concepción, Chile
^{16}
Department of Physics, University of Illinois, 1110 West Green Street, Urbana, IL 61801, USA
^{17}
Fermi National Accelerator Laboratory, MS209, PO Box 500, Batavia, IL 60510, USA
^{18}
Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA
^{19}
East Asian Observatory, 660 N. A’ohoku Place, Hilo, HI 96720, USA
^{20}
James Clerk Maxwell Telescope (JCMT), 660 N. A’ohoku Place, Hilo, HI 96720, USA
^{21}
California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA
^{22}
Institute of Astronomy and Astrophysics, Academia Sinica, 645 N. A’ohoku Place, Hilo, HI 96720, USA
^{23}
Department of Physics and Astronomy, University of Hawaii at Manoa, 2505 Correa Road, Honolulu, HI 96822, USA
^{24}
Department of Physics, McGill University, 3600 Rue University, Montréal, QC H3A 2T8, Canada
^{25}
Trottier Space Institute at McGill, 3550 Rue University, Montréal, QC H3A 2A7, Canada
^{26}
Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, 38406 SaintMartind’Hères, France
^{27}
Perimeter Institute for Theoretical Physics, 31 Caroline Street North, Waterloo, ON N2L 2Y5, Canada
^{28}
Department of Physics and Astronomy, University of Waterloo, 200 University Avenue West, Waterloo, ON N2L 3G1, Canada
^{29}
Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L 3G1, Canada
^{30}
Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics (IMAPP), Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{31}
Department of Astronomy, University of Massachusetts, Amherst, MA 01003, USA
^{32}
Kavli Institute for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA
^{33}
Department of Physics, University of Chicago, 5720 South Ellis Avenue, Chicago, IL 60637, USA
^{34}
Enrico Fermi Institute, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA
^{35}
Princeton Gravity Initiative, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA
^{36}
Data Science Institute, University of Arizona, 1230 N. Cherry Ave., Tucson, AZ 85721, USA
^{37}
Program in Applied Mathematics, University of Arizona, 617 N. Santa Rita, Tucson, AZ 85721, USA
^{38}
Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, NY 14853, USA
^{39}
Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, PR China
^{40}
Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100101, PR China
^{41}
Korea Astronomy and Space Science Institute, Daedeokdaero 776, Yuseonggu, Daejeon 34055, Republic of Korea
^{42}
Department of Astronomy, Yonsei University, Yonseiro 50, Seodaemungu 03722 Seoul, Republic of Korea
^{43}
Physics Department, Fairfield University, 1073 North Benson Road, Fairfield, CT 06824, USA
^{44}
Department of Astronomy, University of Illinois at UrbanaChampaign, 1002 West Green Street, Urbana, IL 61801, USA
^{45}
Instituto de Astronomía, Universidad Nacional Autónoma de México (UNAM), Apdo Postal 70264, Ciudad de México, Mexico
^{46}
Institut für Theoretische Physik, GoetheUniversität Frankfurt, MaxvonLaueStraße 1, 60438 Frankfurt am Main, Germany
^{47}
Research Center for Intelligent Computing Platforms, Zhejiang Laboratory, Hangzhou 311100, PR China
^{48}
TsungDao Lee Institute, Shanghai Jiao Tong University, Shengrong Road 520, Shanghai 201210, PR China
^{49}
Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, 500 W. 120th Street, New York, NY 10027, USA
^{50}
Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
^{51}
Dipartimento di Fisica “E. Pancini”, Università di Napoli “Federico II”, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, 80126 Napoli, Italy
^{52}
INFN Sez. di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, 80126 Napoli, Italy
^{53}
Wits Centre for Astrophysics, University of the Witwatersrand, 1 Jan Smuts Avenue, Braamfontein, Johannesburg 2050, South Africa
^{54}
Department of Physics, University of Pretoria, Hatfield, Pretoria 0028, South Africa
^{55}
Centre for Radio Astronomy Techniques and Technologies, Department of Physics and Electronics, Rhodes University, Makhanda 6140, South Africa
^{56}
ASTRON, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
^{57}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 Place Jules Janssen, 92195 Meudon, France
^{58}
JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA
^{59}
National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, PR China
^{60}
Las Cumbres Observatory, 6740 Cortona Drive, Suite 102, Goleta, CA 931175575, USA
^{61}
Department of Physics, University of California, Santa Barbara, CA 931069530, USA
^{62}
National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
^{63}
Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, 32D476, 77 Massachusetts Ave., Cambridge, MA 02142, USA
^{64}
Google Research, 355 Main St., Cambridge, MA 02142, USA
^{65}
Institut für Theoretische Physik und Astrophysik, Universität Würzburg, EmilFischerStr. 31, 97074 Würzburg, Germany
^{66}
Department of History of Science, Harvard University, Cambridge, MA 02138, USA
^{67}
Department of Physics, Harvard University, Cambridge, MA 02138, USA
^{68}
NCSA, University of Illinois, 1205 W. Clark St., Urbana, IL 61801, USA
^{69}
Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, R. do Matão, 1226, São Paulo, SP 05508090, Brazil
^{70}
Dipartimento di Fisica, Università degli Studi di Cagliari, SP MonserratoSestu km 0.7, 09042 Monserrato, CA, Italy
^{71}
INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius, CA, Italy
^{72}
INFN, Sezione di Cagliari, 09042 Monserrato, CA, Italy
^{73}
CP3Origins, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark
^{74}
Instituto Nacional de Astrofísica, Óptica y Electrónica, Apartado Postal 51 y 216, 72000 Puebla Pue., Mexico
^{75}
Consejo Nacional de Humanidades, Ciencia y Tecnología, Av. Insurgentes Sur 1582, 03940 Ciudad de México, Mexico
^{76}
Key Laboratory for Research in Galaxies and Cosmology, Chinese Academy of Sciences, Shanghai 200030, PR China
^{77}
Mizusawa VLBI Observatory, National Astronomical Observatory of Japan, 212 Hoshigaoka, Mizusawa, Oshu, Iwate 0230861, Japan
^{78}
Department of Astronomical Science, The Graduate University for Advanced Studies (SOKENDAI), 2211 Osawa, Mitaka, Tokyo 1818588, Japan
^{79}
NOVA Submm Instrumentation Group, Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen, The Netherlands
^{80}
Department of Astronomy, School of Physics, Peking University, Beijing 100871, PR China
^{81}
Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, PR China
^{82}
Department of Astronomy, Graduate School of Science, The University of Tokyo, 731 Hongo, Bunkyoku, Tokyo 1130033, Japan
^{83}
The Institute of Statistical Mathematics, 103 Midoricho, Tachikawa, Tokyo 1908562, Japan
^{84}
Department of Statistical Science, The Graduate University for Advanced Studies (SOKENDAI), 103 Midoricho, Tachikawa, Tokyo 1908562, Japan
^{85}
Kavli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, 515 Kashiwanoha, Kashiwa 2778583, Japan
^{86}
Leiden Observatory, Leiden University, Postbus 2300, 9513 RA Leiden, The Netherlands
^{87}
ASTRAVEO LLC, PO Box 1668, Gloucester, MA 01931, USA
^{88}
Applied Materials Inc., 35 Dory Road, Gloucester, MA 01930, USA
^{89}
Institute for Astrophysical Research, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA
^{90}
University of Science and Technology, Gajeongro 217, Yuseonggu, Daejeon 34113, Republic of Korea
^{91}
Institute for Cosmic Ray Research, The University of Tokyo, 515 Kashiwanoha, Kashiwa, Chiba 2778582, Japan
^{92}
Joint Institute for VLBI ERIC (JIVE), Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
^{93}
Department of Astronomy and Atmospheric Sciences, Kyungpook National University, Daegu 702701, Republic of Korea
^{94}
Department of Physics, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehakro, Yuseonggu, Daejeon 34141, Republic of Korea
^{95}
Kogakuin University of Technology & Engineering, Academic Support Center, 26651 Nakano, Hachioji, Tokyo 1920015, Japan
^{96}
Graduate School of Science and Technology, Niigata University, 8050 Ikarashi 2nocho, Nishiku, Niigata 9502181, Japan
^{97}
Physics Department, National Sun YatSen University, No. 70, LienHai Road, Kaosiung City 80424, Taiwan ROC
^{98}
School of Astronomy and Space Science, Nanjing University, Nanjing 210023, PR China
^{99}
Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Nanjing 210023, PR China
^{100}
INAFIstituto di Radioastronomia, Via P. Gobetti 101, 40129 Bologna, Italy
^{101}
INAFIstituto di Radioastronomia & Italian ALMA Regional Centre, Via P. Gobetti 101, 40129 Bologna, Italy
^{102}
Department of Physics, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd., Taipei 10617, Taiwan ROC
^{103}
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia 58089, Mexico
^{104}
Yunnan Observatories, Chinese Academy of Sciences, 650011 Kunming, Yunnan Province, PR China
^{105}
Center for Astronomical MegaScience, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, PR China
^{106}
Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, 650011 Kunming, PR China
^{107}
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
^{108}
Gravitation and Astroparticle Physics Amsterdam (GRAPPA) Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
^{109}
Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA
^{110}
School of Physics and Astronomy, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240, PR China
^{111}
Institut de Radioastronomie Millimétrique (IRAM), Avenida Divina Pastora 7, Local 20, 18012 Granada, Spain
^{112}
National Institute of Technology, Hachinohe College, 161 Uwanotai, Tamonoki, Hachinohe City, Aomori 0391192, Japan
^{113}
Research Center for Astronomy, Academy of Athens, Soranou Efessiou 4, 115 27 Athens, Greece
^{114}
Department of Physics, Villanova University, 800 Lancaster Avenue, Villanova, PA 19085, USA
^{115}
Physics Department, Washington University, CB 1105, St. Louis, MO 63130, USA
^{116}
School of Physics, Georgia Institute of Technology, 837 State St NW, Atlanta, GA 30332, USA
^{117}
Department of Astronomy and Space Science, Kyung Hee University, 1732, Deogyeongdaero, Giheunggu, Yonginsi, Gyeonggido 17104, Republic of Korea
^{118}
Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
^{119}
Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada
^{120}
Canadian Institute for Advanced Research, 180 Dundas St West, Toronto, ON M5G 1Z8, Canada
^{121}
Department of Physics, National Taiwan Normal University, No. 88, Sec. 4, Tingzhou Rd., Taipei 116, Taiwan ROC
^{122}
Center of Astronomy and Gravitation, National Taiwan Normal University, No. 88, Sec. 4, Tingzhou Road, Taipei 116, Taiwan ROC
^{123}
Finnish Centre for Astronomy with ESO, University of Turku, 20014 Turku, Finland
^{124}
Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland
^{125}
Gemini Observatory/NSF’s NOIRLab, 670 N. A’ohōkū Place, Hilo, HI 96720, USA
^{126}
Frankfurt Institute for Advanced Studies, RuthMoufangStrasse 1, 60438 Frankfurt, Germany
^{127}
School of Mathematics, Trinity College, Dublin 2, Ireland
^{128}
Department of Physics, University of Toronto, 60 St. George Street, Toronto, ON M5S 1A7, Canada
^{129}
Department of Physics, Tokyo Institute of Technology, 2121 Ookayama, Meguroku, Tokyo 1528551, Japan
^{130}
Hiroshima Astrophysical Science Center, Hiroshima University, 131 Kagamiyama, HigashiHiroshima, Hiroshima 7398526, Japan
^{131}
Aalto University Department of Electronics and Nanoengineering, PL 15500, 00076 Aalto, Finland
^{132}
Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK
^{133}
National Biomedical Imaging Center, Peking University, Beijing 100871, PR China
^{134}
College of Future Technology, Peking University, Beijing 100871, PR China
^{135}
Tokyo Electron Technology Solutions Limited, 52 Matsunagane, Iwayado, Esashi, Oshu, Iwate 0231101, Japan
^{136}
Department of Physics and Astronomy, University of Lethbridge, Lethbridge, Alberta T1K 3M4, Canada
^{137}
Netherlands Organisation for Scientific Research (NWO), Postbus 93138, 2509 AC Den Haag, The Netherlands
^{138}
Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 9808578, Japan
^{139}
Astronomical Institute, Tohoku University, Sendai 9808578, Japan
^{140}
Department of Physics and Astronomy, Seoul National University, Gwanakgu, Seoul 08826, Republic of Korea
^{141}
University of New Mexico, Department of Physics and Astronomy, Albuquerque, NM 87131, USA
^{142}
Physics Department, Brandeis University, 415 South Street, Waltham, MA 02453, USA
^{143}
Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Turku, Finland
^{144}
Radboud Excellence Fellow of Radboud University, Nijmegen, The Netherlands
^{145}
School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA
^{146}
School of Physics, Huazhong University of Science and Technology, Wuhan, Hubei 430074, PR China
^{147}
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK
^{148}
School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, PR China
^{149}
Astronomy Department, University of Science and Technology of China, Hefei 230026, PR China
^{150}
Department of Physics and Astronomy, Michigan State University, 567 Wilson Rd, East Lansing, MI 48824, USA
^{151}
University of Vermont, Burlington, VT 05405, USA
^{152}
National ChungShan Institute of Science and Technology, No. 566, Ln. 134, Longyuan Rd., Longtan Dist., Taoyuan City 325, Taiwan ROC
^{153}
Western University, 1151 Richmond Street, London, Ontario N6A 3K7, Canada
^{154}
Department of Physics, Graduate School of Science, Osaka Metropolitan University, 11 Gakuencho, Nakaku, Sakai, Osaka 5998531, Japan
^{155}
Department of Physics, National Cheng Kung University, Tainan 701401, Taiwan ROC
^{156}
Coorporación Mexicana de Investigación en Materiales S.A. de C.V. Ciencia y Tecnología 790, Colonia Saltillo 400, Saltillo, Coahuila CP 25290, Mexico
^{157}
Dipartimento di Fisica, Università di Trieste, 34127 Trieste, Italy
^{158}
INFN Sez. di Trieste, 34127 Trieste, Italy
^{159}
Ground Station for Deep Space Exploration and Telecommunication (GREAT) Project, Institute of Space and Astronautical Science (ISAS), Japan Aerospace Exploration Agency (JAXA), 311 Yoshinodai, Chuoku, Sagamihara, Kanagawa 2525210, Japan
Received:
11
September
2023
Accepted:
9
November
2023
In April 2019, the Event Horizon Telescope (EHT) Collaboration reported the firstever eventhorizonscale images of a black hole, resolving the central compact radio source in the giant elliptical galaxy M 87. These images reveal a ring with a southerly brightness distribution and a diameter of ∼42 μas, consistent with the predicted size and shape of a shadow produced by the gravitationally lensed emission around a supermassive black hole. These results were obtained as part of the April 2017 EHT observation campaign, using a global very long baseline interferometric radio array operating at a wavelength of 1.3 mm. Here, we present results based on the second EHT observing campaign, taking place in April 2018 with an improved array, wider frequency coverage, and increased bandwidth. In particular, the additional baselines provided by the Greenland telescope improved the coverage of the array. Multiyear EHT observations provide independent snapshots of the horizonscale emission, allowing us to confirm the persistence, size, and shape of the black hole shadow, and constrain the intrinsic structural variability of the accretion flow. We have confirmed the presence of an asymmetric ring structure, brighter in the southwest, with a median diameter of 43.3_{−3.1}^{+1.5} μas. The diameter of the 2018 ring is remarkably consistent with the diameter obtained from the previous 2017 observations. On the other hand, the position angle of the brightness asymmetry in 2018 is shifted by about 30° relative to 2017. The perennial persistence of the ring and its diameter robustly support the interpretation that the ring is formed by lensed emission surrounding a Kerr black hole with a mass ∼6.5 × 10^{9} M_{⊙}. The significant change in the ring brightness asymmetry implies a spin axis that is more consistent with the position angle of the largescale jet.
Key words: accretion, accretion disks / black hole physics / gravitation / galaxies: active / galaxies: individual: M 87 / galaxies: jets
The reduced image (FITS file) is available at the CDS via anonymous ftp to https://cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/vizbin/cat/J/A+A/681/A79
© The Authors 2024
Open 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 nonrotating 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 R_{c} changes for spinning black holes, but within ≲4%. The photon
capture and redshift effects produce gravitational signatures in the observed images of astrophysical black holes immersed in background radiation fields, consisting of a characteristic dark center surrounded by ringlike emission (e.g., Luminet 1979), commonly referred to as the black hole shadow. Therefore, observing the shadow, measuring its properties, such as the ring diameter, and comparing it with theoretically predicted shadow morphologies provide unique opportunities to explore curved spacetime and extreme gravitational potential, allowing direct tests of general relativity. This is complementary to tests from gravitational waves, high precision astrometry, and cosmology (e.g., Event Horizon Telescope Collaboration 2022f; Abbott et al. 2016; GRAVITY Collaboration 2018; Ferreira 2019).
From an astrophysical point of view, black holes can be formed in many ways, and the masses of astrophysical black holes vary greatly, from relatively small (≲100 M_{⊙}; stellar mass) to supermassive scales (∼10^{6} − 10^{9} M_{⊙}); for example, readers can refer to Greene et al. (2020). Supermassive black holes (SMBHs) are predominantly located at the centers of galaxies, with larger black holes associated with more massive host galaxies (Kormendy & Ho 2013). When surrounding matter accretes onto the central SMBHs, the gravitational potential energy can be released in various forms, making active galactic nuclei (AGNs) shine across the electromagnetic spectrum, depending on the type and rate of mass accretion. Accreting SMBHs can also produce collimated relativistic plasma jets and release enormous amounts of energy into space in the form of strong magnetic fields and relativistic particles. Theoretical studies suggest the ultimate energy source of the jet to be either the accretion flow (Blandford & Payne 1982) or the central black hole itself (Blandford & Znajek 1977), in each case mediated through magnetic fields threading the plasma in the disk or near the horizon. Therefore, highresolution imaging of accreting SMBHs can also provide unique constraints on modeling the detailed physics of the accretion disk and jet launching in AGNs (e.g., Blandford et al. 2019).
Since the angular diameter of the black hole shadow is d_{sh} = 2R_{c}D^{−1} ≈ 10 GMc^{−2}D^{−1} where D is the distance to the black hole (Bardeen 1973), the imaging capability of the black hole shadow requires a masstodistance ratio sufficiently large to provide d_{sh} resolvable with the existing instruments. In this regard, SMBHs in the center of the radio galaxy M 87 (hereafter M 87^{*}) and our Galaxy, Sagittarius A^{*} (Sgr A^{*}), offer the best opportunities. As for M 87, its distance is measured via multiple methods using primary or secondary distance indicators, and we adopt D = 16.8 ± 0.8 Mpc (Bird et al. 2010; Blakeslee et al. 2009; Cantiello et al. 2018, see Event Horizon Telescope Collaboration 2019f). The central black hole mass is approximately 6 × 10^{9} M_{⊙} (Gebhardt et al. 2011; Event Horizon Telescope Collaboration 2019f). The combination of D and M yields a d_{sh} ≈ 40 μas, which is the second largest value on the sky after Sgr A^{*} (M ≈ 4 × 10^{6} M_{⊙}, D ≈ 8 kpc, d_{sh} ≈ 50 μas; see Event Horizon Telescope Collaboration 2022a). M 87 is a representative FanaroffRiley 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). Multiwavelength (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 ultrahigh 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 lightcrossing times for a few billion solarmass 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 ringlike structure with a diameter of 42 ± 3 μas, a central flux depression with an intensity ratio of > 10 : 1, and a position angle of the brightness asymetry of ∼180°. These measurements are consistent with expectations for a SMBH of mass M = (6.5 ± 0.7) × 10^{9} M_{⊙} (M 87^{∗} 2017 IV; M 87^{∗} 2017 VI). A small level of interday variability was also detected, which was analyzed in followup 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 ringlike 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 preEHT array. By assuming a ringlike 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 ringlike 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 northsouth direction. New EHT images of M 87^{*} can provide further evidence that the eventhorizonscale 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 ringlike structure. A detection of variability around the ringlike 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 nontrivial structural changes.
We explain how we derived preimaging 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 preimaging 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 ringlike 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.
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 ringlike structure seen in the imaging methods, in Sect. 6 we quantify the support for the ringlike structure by comparing the Bayesian evidence of different geometric models against the data. We find that ringlike 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 modeldependent 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 12m and seven 6m dishes, respectively, as well as the following seven singledish stations: Atacama Pathfinder Experiment (APEX), IRAM 30m 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.
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. Colocated 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 dashedlines. 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 doublesideband (DSB) receiver to a twosinglesideband (2SB) one.
The telescope systems recorded both righthand and lefthand 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.
Median zenith sky opacities (1.3 mm) at EHT sites during the 2018 April observations toward M 87^{*}.
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).
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 lowband (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 MaxPlanckInstitut 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 parallelhand (RR, LL) and crosshand (RL, LR) polarization products. Due to the nonmatching 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 postprocessing
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, EHTHOPS (Blackburn et al. 2019) and the CASAbased (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 EHTHOPS 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 singleband and multiband 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, elevationdependent gain curves, and timedependent system temperatures derived from the raw metadata collected for each station, as described in Koay et al. (2023a). A summary of the amplitude calibration parameters of each station and their associated uncertainties is provided in Appendix A. The data are then averaged over time (10 s) and each 1.856 GHz band before being exported by AIPS as UVFITS files. Subsequently, we perform a polarimetric gains ratio calibration (polgainscal) to align the RR and LL phases in the HOPS data. We correct for phase offsets between the two polarization streams and allow the Stokes I visibilities to be obtained^{1}. Network calibration (M 87^{∗} 2017 III) with 10 s solution intervals is then applied to correct the amplitudes of stations with a colocated 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 intrasite baselines and that the flux densities on the intrasite 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 nonclosing systematic uncertainties possibly arising from polarimetric leakages and bandpass effects (see Appendix B). For M 87^{*}, the nonclosing 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 stationbased amplitude gain errors such as those described in Appendices A and E, which are corrected by modeldependent calibration during imaging.
Since some imaging pipelines use timeaveraged 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 nonclosing 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 scanaveraged 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 signaltonoise 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 highbands 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 lowband.
The GLT adds northsouth coverage to the EHT array in bands 3 and 4, probing baseline lengths of ∼7.1 Gλ with the colocated 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 eastwest direction given by PV–Hawai‘i are missing, such that the eastwest 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 networkcalibrated 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 ringlike 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.
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 lowband (gray points). All data shown include a priori stationbased amplitude calibration and network calibration but are prior to any modeldependent selfcalibration. 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 nontrivial 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 daytoday 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.
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 (lowband), gray circles denote data from 2017 April 11 (lowband), 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 networkcalibrated 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 nonBayesian image reconstruction methods (see Sect. 5).
Measurements of the total flux density of the arcsecondscale 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 intrasite baselines (ALMA–APEX and SMA–JCMT) of 0.1−1.1 Mλ and the intersite baselines of 1.3−7.3 Gλ, which are sensitive to the subarcsecond and microarcsecond scale structures, respectively. The longest intrasite baseline is ALMA–APEX, with fringe spacings of 106−131 mas for M 87^{*}. On the other hand, that of the shortest intersite 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 arcsecondscale radio core, constrained by the short intrasite baselines. The total flux density as measured by the ALMAonly interferometric observations was F_{tot} ≈ 1.13 Jy by averaging the values over four bands on April 21 and 22. The second one is the total compact flux density, F_{cpct} ≤ F_{tot}, which is the integrated flux density of the microarcsecondscale structure constrained by intersite baselines.
We estimate the range of F_{cpct} and the size of the compact emission region θ_{cpct}, represented by an equivalent Gaussian fullwidth half maximum (FWHM), following the two independent procedures in M 87^{∗} 2017 IV. One is directly estimated from the 2018 EHT visibility data. With this procedure, we can derive the constraints on F_{cpct} and θ_{cpct} simultaneously. The other is a constraint on the F_{cpct} by utilizing quasisimultaneous VLBI data based on an extrapolation between 1.3 mm and the longer wavelengths. In Appendix F, we derive a series of constraints on F_{cpct} and θ_{cpct} using these procedures.
The EHT constraints give F_{cpct} = 0.30 − 1.13 Jy with a θ_{cpct} between 39 and 98 μas. The source size constraints are roughly consistent with the 2017 constraints. The resulting constraints on F_{cpct} for 2018 are looser than those for 2017 due to the lack of additional constraints using SMT–PV and LMT–PV crossing baseline tracks. Nevertheless, the MWL constraints give a slightly tighter range of F_{cpct} = 0.5 − 0.7 Jy within the central 100 × 100 μas^{2} of the nucleus, which is comparable to the F_{cpct} = 0.56 − 1.03 Jy estimated from the 2017 data. Therefore, we conclude that the total compact flux density of M 87^{*} in 2018 at 1.3 mm is broadly consistent with that of 2017 when considering all these estimates. While these flux density and size constraints in the preimaging 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 preimaging 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 illposed 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 nonunique, 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 selfcalibration 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 CLEANbased 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 ehtimaging (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 (cres90, cres0, cres90, cres180, dblsrc, disk, ring, edisk, point+disk, point+edisk) and a snapshot of a generalrelativistic 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 ringlike and nonring 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 ehtimaging 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 bestfit 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 groundtruth 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 chisquares 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 nonclosing 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.
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.1–5.1.3.
As mentioned above, we also utilized two methods, THEMIS (Broderick et al. 2020a) and Comrade (Tiede 2022), that followed a Bayesian posterior sampling approach to image reconstruction. To quantify the imaging uncertainty, THEMIS and Comrade use Bayesian inference and cast VLBI imaging as a Bayesian inverse problem. Since these Bayesian methods do not require training for selecting model hyperparameters, all synthetic data models are treated as validation data sets. In addition to producing a bestfit 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 pointsource 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 selfcalibration using the same averaging timescale. 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 intrasite 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 selfcalibration 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 selfcalibration, diameter of the CLEAN window, and the powerlaw 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 selfcalibration with the initial imaging result.
We implemented phaseonly selfcalibration 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 phaseonly selfcalibration, 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.
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, lowband). 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. ehtimaging
The Python package ehtimaging (Chael et al. 2018, 2019a,b) is an RMLbased VLBI imaging software capable of producing images by placing different relative weights on the fits to closure quantities and complex visibilities. The results in this paper were produced by performing a parameter survey using the 2017 M 87^{*} imaging pipeline^{2}. 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 preimaging 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 scanaveraging the data. Then the correlated flux densities at intrasite 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 nonclosing errors in the data. As measurements taken by the LMT suffer from large gain uncertainties, we performed an initial amplitudeonly selfcalibration to the LMT data and adopted station based gain corrections for LMT. This selfcalibration 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 preimaging calibration steps, the pipeline proceeded with four iterations of imaging and selfcalibration. The imaging was initialized with a circular Gaussian of FWHM and compact flux density specified by the given parameter combination. The details of the selfcalibration 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 selfcalibration 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 nonzero 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 selfcalibration rounds were performed on the scanaveraged 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 bestfit 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 preimaging calibration step.
5.1.3. SMILI
The imaging pipeline SMILI (Akiyama et al. 2017a,b) is a Python and FORTRAN RMLbased imaging software. Similar to the survey conducted with ehtimaging, 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 scanaveraged the visibilities. Precalibration of the LMT gains and the addition of nonclosing systematic errors were performed as described in the second paragraph of Sect. 5.1.2. We performed four imaging cycles for each selfcalibration stage, selfcalibrating 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 recentered to the image center of mass and blurred with a 20 μas FWHM circular Gaussian. Each imaging cycle performed 1000 iterations of the LBFGSB algorithm (Byrd et al. 1995; Zhu et al. 1997) used for the gradientdescent optimization in SMILI’s image solver. After the three imaging cycles, complex visibilities were selfcalibrated with their output image. In the first two selfcalibration 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 selfcalibration stages. Similar to ehtimaging, SMILI also performed all rounds of selfcalibration on scanaveraged 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 onsky brightness distribution and the timedependent 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 intrasite baselines and the shortest intersite 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 scanaveraged 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 nonparametric modeling and include a rasterized image model similar to the one described in Broderick et al. (2020b). We first scanaveraged 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 intrasite baselines.
Our model for all image reconstructions consisted of a rasterized image model (of dimensions N_{x} × N_{y}), a Gaussian of FWHM = 1 mas, and scanaveraged station gains. The Gaussian was used to model the amplitudes of short intrasite (ALMA–APEX, JCMT–SMA) baselines. In the image model, the grid of raster points was convolved with a thirdorder basisspline (Bspline) 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 Bspline kernel. A detailed description of the image model is given in Appendix L.
The hyperparameters of the model are the FOV and raster size (N_{x} × N_{y}). For April 21 bands 3 and 4, we used a 12 × 12 pixel raster with a 7.5 μas pixel size, which was ∼1/3 of the nominal resolution and FOV of 90 μas. This FOV was chosen by checking the visibilityamplitude 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.
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 f_{g} for the flux density of the largescale 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 loggain 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 unnormalized posterior was formed by taking a product of the prior and the likelihood distributions. To sample from the posterior, we used a twostage strategy. First, to find a reasonable starting location, we adopted the LBFGS 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 NoUTurn 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 EHTlike 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.
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, ehtimaging, 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, ehtimaging, 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, ehtimaging, and SMILI could all produce a nonzero 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 timedependent 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:
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. 
Fig. 9. Visibility amplitudes of band 3 data on April 21 as a function of baseline length compared with corresponding gaincalibrated visibility amplitudes from the five representative image models from each pipeline. The fiducial images are used for DIFMAP, ehtimaging, 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 selfcalibration. ehtimaging scales the intrasite baselines to 0.5 Jy, and SMILI scales the intrasite baselines to 0.6 Jy, and both precalibrate the LMTSMT baselines before fitting. DIFMAP uses 10 s averaged data, and all other methods use scanaveraged data. THEMIS and Comrade apply 1% fractional systematic noise to all baselines, and ehtimaging 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 selfcalibration 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. 
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:
where N_{lnAC} is the number of log closure amplitude data points, is the model closure amplitude, A_{C} 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, ehtimaging, 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 downweighted the ALMA visibilites during selfcalibration, resulting in comparable normalized χ^{2} values to the RML methods.
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, loglikelihoods, 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.
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. ):
where N_{V} is the number of complex visibility data points, N_{model DOF} is the number of model degrees of freedom, N_{gain DOF} is the number of gain degrees of freedom (we multiply the final reduced number by 2 since we fit both the real and imaginary components), 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 ringlike morphology compared to bands 1 and 2.
On the April 25 data, DIFMAP, ehtimaging, 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 ringlike 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 ehtimaging, and 0.6 Jy for SMILI; all of these values fell within the range of allowed compact flux densities derived from the preimaging constraints (Sect. 4). The RML methods only modeled the compact emission (since they rescale the intrasite 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 preimaging 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 intrasite visibilities and verified that the image structure is not significantly impacted. See Sect. 7.3 for a related discussion of the compact flux density estimates from the Bayesian methods, and Appendix M for supplementary synthetic data validation tests at different compact flux densities for the RML and CLEAN methods.
The difference of the flux densities between F_{tot} and F_{cpct} implies there is some emission outside the compact region. However, its structure (e.g., scale and direction) is unconstrained by the data (Sect. 4). The submilliarcsecond 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 submilliarcsecond 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 submilliarcsecond 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 submilliarcsecond 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 smallscale diffuse structure in the 2017 EHT data reconstructions which aligns with the limbbrightened jet observed at 3 mm (Arras et al. 2022; Carilli & Thyagarajan 2022; Broderick et al. 2022a), no constraints on largerscale (> 0.2 mas) jet emission are supported by the data.
5.2.3. Image statistics
The distributions of both the Top Set and posterior images help us understand the image uncertainties associated with each method. Figures 10–12 show the image and visibilitydomain uncertainties associated with the image sets for ehtimaging, THEMIS, and Comrade, respectively.
Fig. 10. Visualization of image statistics calculated using the Top Set images from the ehtimaging 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 ehtimaging 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. 
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. 
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 ehtimaging 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 visibilitydomain 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 ringlike 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 ringlike feature is robustly recovered, finding it has a low fractional standard deviation. Additionally, both THEMIS and Comrade find no evidence for nonring emission. For THEMIS, this can be seen from the fact that the raster is concentrated around the ring. Substantial nonring 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 nonring emission and finds no significant (> 3σ) detection of nonring 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 sparsityinducing 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 visibilitydomain 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 ringlike 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. Negativevalued 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 ehtimaging and is described in detail in M 87^{∗} 2017 IV. The second IDFE tool, Variational ImageDomain 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 mFring 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 ringlike images, VIDA can support a variety of templates, which can be used to evaluate the success of our imaging methods by quantitatively reproducing the nonring 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 ringlike 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 superresolution, it is challenging to combine the width measurements across methods. While some methods trend toward slightly lower widths compared to their 2017 values (DIFMAP, ehtimaging, 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 ringlike structures.
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 (ehtimaging) 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 unshaded region corresponds to the 2018 April 25 measurements. 
Tables 7–10 list the diameter, width, position angle, asymmetry, and fractional central brightness measured with REx and VIDA from all days, bands, a priori calibration, and imaging pipelines. The brightness asymmetry A is ∼0.2 − 0.4, indicating an asymmetric ringlike structure. Moreover, the fractional central brightness f_{C} is ∼0 − 0.5, indicating a clear central depression.
Ring parameters for 2018 April 21 HOPS data.
Ring parameters for the 2018 April 25 HOPS data.
Ring parameters for 2018 April 21 CASA data.
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, ehtimaging, 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 ringlike 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 ringlike 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 imagedomain 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 ringlike features are preferred by the data by comparing the Bayesian evidence of both ringlike and nonring visibilitydomain 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 ringlike 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 mrings (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 stationbased 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 ringlike models provide a much better fit to the data when compared to nonring models. We select the bestperforming crescent and mring models; the THEMIS xsringauss model from M 87^{∗} 2017 VI (11 parameters) and a stretched mring of order 2 (8 parameters). We use these two models as the basis for defining fiducial models for direct feature extraction.
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 ringlike 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 ringlike image feature associated with halfinteger 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 halforbit around the back of the black hole to impact our line of sight, and so on. This feature is consistently seen in analytic, semianalytic, and GRMHD simulations of optically thin accretion flows around black holes, and directly modeling this feature can, in principle, provide strong constraints on spacetime 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 largescale 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 scanaveraged, networkcalibrated 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 xsringauss 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 xsringauss 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 “tophat” 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 xsringauss model construction is in the treatment of the diffuse emission. The xsringauss model used a pair of asymmetric Gaussians to help model any nonring 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 xsringauss, 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 nonring emission, the xsringauss was more agnostic to the precise structure of the ring.
We compare the xsringauss model to the scanaveraged networkcalibrated 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 xsringauss. 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 xsringauss parameters in Appendix N.
6.2.2. Comrade direct modeling method
We also define a fiducial model based on the bestperforming mring model in Fig. 14 for analysis with Comrade. The fiducial model was taken to be an mring 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 mring. The resulting combination was convolved with a circular Gaussian. This construction defines a nineparameter, elliptical, floored mring model, which we refer to hereafter as an “mFring”. An additional pair of elliptical nuisance Gaussians were also included to improve the fit quality. A summary of the mFring 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 mFring and the xsringauss compact flux density. We do not report a compact flux density for the mFring 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 xsringauss, as well as the log closure amplitude and closure phase residuals for the mFring. The gray points represent the scanaveraged data, the colored points represent the maximum likelihood model for the THEMISbased models, and the maximum a posteriori model for the mFring. We find that each model performs very well when compared against the data, with small, unstructured residuals. Appendix N also presents imagedomain representations for each day and all observing bands, demonstrating consistency with the other imaging methods.
Fig. 15. Complex visibility (V) fit comparisons and normalized residuals for the Hybrid Themage and xsringauss, and the logclosure amplitude (lnA_{C}) and closure phase (ψ_{C}) fit comparisons and normalized residuals for the mFring 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 xsringauss, and from the maximum a posteriori model for the mFring. The gray points represent the scanaveraged 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 mFring, and 202° as fitted by the xsringauss. 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 bestfit xsringauss and mFring models produce consistent position angle extraction and thus conclude that the measured position angle from the mFring and the xsringauss model are consistent.
Fig. 16. Fitted features of the xsringauss, mFring, 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 mFring model, while the red band is the equivalent for the xsringauss. 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 xsringauss fits. The hashed region is the 68% posterior probability interval taken from the 2017 M 87^{*} xsringauss fits. 
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 xsringauss and mFring 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 mFring 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 ringlike emission structure around the central SMBH M 87^{*}. This result across various independent methods demonstrates the robustness of our conclusions. The consistency of a ringlike emission structure on a 1year 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.
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 xsringauss, and the mFring 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 ehtimaging, 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 xsringauss and mFring 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 ringlike structure, from the direct modeling methods which do assume a ringlike 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 followup 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 counterclockwise shift in the position angle of the brightness asymmetry in the ring from 2017 to 2018. Finally, we investigate the challenges related to estimating the compact flux density of M 87^{*} on horizon scales with the sparse (u, v) coverage of the EHT.
7.1. Consistency of the ring diameter
The asymmetric ring’s estimated characteristic parameters are provided in Tables 7–10 via image domain feature extraction and Table 11 via the direct modeling results. The parameters are broadly consistent across different calibration schemes and reconstruction methods. There is a tendency for the ring diameters extracted with the image domain feature extraction to be slightly smaller than those from the direct modeling methods. This tendency is also observed in the 2017 results and is related to the effective resolution of each method (M 87^{∗} 2017 VI).
We combine the results from image domain feature extraction and direct modeling to construct a median diameter across all methods (Tables 7–11). We first construct a diameter histogram for each method, normalize the histograms to equally weigh each method, then combine the normalized histograms to produce a single diameter histogram across all imaging and direct modeling methods. We present the resulting median diameter and 68% intervals in the last column of Table 12. This lets us more easily compare the ring diameter of 2018 EHT campaign with that of 2017 (see Table 1 of M 87^{∗} 2017 I). We find the difference in median diameter between these two campaigns is on the order of only ∼1 μas, with error bars that almost completely overlap, indicating a persistent emission structure.
Comparison of the extracted ring parameters for 2017 and 2018.
Since the diameter estimates are modeldependent, 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:
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 electrontoproton 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 recalibrate 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 xsringauss or Table 6 for DIFMAP, ehtimaging, 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, xsringauss, and mFring 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.
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, ehtimaging, SMILI, xsringauss), we use the methodspecific 2017 scaling factor to determine the 2018 d/α_{2017} values. For DIFMAP, ehtimaging, and SMILI we use the scaling factors from Table 6 of M 87^{∗} 2017 VI, and for the xsringauss 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, mFring), 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 modelspecific 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 xsringauss 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 ringlike structure in our images is consistent with the shadow of a supermassive black hole with a mass of ∼6.5 × 10^{9} M_{⊙}. Any discrepancies may disappear in a proper 2018 θ_{g} calibration.
We note that there are at least two new integral field spectroscopic stellar kinematics mass estimates for M 87^{*} recently published. One of these adapts new triaxial orbit models to Keck II telescope observations (Liepold et al. 2023) and favors a mass for M 87^{*} of 5.3 × 10^{9} M_{⊙}. The other estimate uses adaptive optics observations from the VLT instruments MUSE and OASIS (Simon et al. 2024) and favors a mass of 8.7 × 10^{9} M_{⊙}. However, when considering a different stellar mass profile in the inner region (and thus a different masstolight ratio), this method obtains a mass estimate of 5.5 × 10^{9} M_{⊙}. These ∼5 × 10^{9} M_{⊙} mass estimates are within 1.5σ of the 2017 EHT mass estimate. Arguably, there are significant systematic uncertainties in the nonhorizon scale mass estimates of M 87^{*} that require continued investigation.
7.2. Position angle variation
It is worth noting that there is a significant counterclockwise change in the position angle of the brighter region in the ringlike 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 largescale jet is normal to the disk.
The work by Wielgus et al. (2020) shows evidence of yearscale 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 yearscale variation is unambiguously confirmed in the image domain (see Fig. 19). Furthermore, recent longterm monitoring studies of M 87^{*} using longerwavelength VLBI found a systematic position angle oscillation of the parsecscale 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 yearscale morphological variations and a counterclockwise 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 interyear EHT images. Further accumulation of EHT images over forthcoming years, along with parsecscale jet monitoring, is required to understand better the origin of the yearscale variation of the ringlike structure and its possible connection to the largescale jet.
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 milliarcsecondscale 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 nonGaussian, and exhibit large tails. 
In addition to the yearscale 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 timescale 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, highband) 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 onefluid GRMHD simulations are initialized with a weakly magnetized torus of plasma orbiting in the equatorial plane of the black hole, and the rotational axis of the torus is aligned with the black hole axis. After the system evolves to a steady state, we perform radiative transfer to compute the thermal synchrotron emission to generate images with the parameterized electron thermal dynamics (see M 87^{∗} 2017 V, for details). The selected models were prepared by the code iharm (Gammie et al. 2003), and the model images including different black hole spins^{3}, accretion flow types, and assumptions for the ratio between the electron and ion temperature in the simulation (see M 87^{∗} 2017 V, for details).
Fig. 20. Distribution of bestfit position angle (in degrees) of the forward jet by fitting model images with the 2017 (April 6, highband) and 2018 (April 21, band 4) observations. The bestfit 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 milliarcsecondscale 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 electrothermodynamics (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 milliarcsecondscale jet (Walker et al. 2018; Cui et al. 2023), which is ∼288° counterclockwise 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 largescale 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 counterclockwise 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 intrasites 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 longerwavelength and zerobaseline 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 EHTscales.
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 xsringauss 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 xsringauss models in Fig. 21, using the HOPS data. There are two clear modes: a “lowestimate” for compact flux density around 0.6 Jy, and a “highestimate” 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 xsringauss 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.
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 bimodal distributions in the compact flux density, FOV_{y}, and to a lesser degree, in the crescent flux density. This suggests that, at least for this reconstruction, adding the entire raster flux density may not be appropriate for computing the compact flux density; the raster is elongated in one direction and places some flux outside of what we might consider the compact region. This is evident in the April 21 band 1 Hybrid Themage image in Fig. N.2, where the raster extends flux in the northwest direction.
Fig. 22. Triangle plot comparing the compact flux density (), raster FOV in the x and y directions, and crescent flux (I_{X}) parameters for the Hybrid Themage reconstruction of April 21 band 1. The compact flux density and raster FOV parameters exhibit clearly bimodal 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 timedependent 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 ringlike 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 largescale jet seen at longer wavelengths.
Acknowledgments
The Event Horizon Telescope Collaboration thanks the following organizations and programs: the Academia Sinica; the Academy of Finland (projects 274477, 284495, 312496, 315721); the Agencia Nacional de Investigación y Desarrollo (ANID), Chile via NCN19_058 (TITANs), Fondecyt 1221421 and BASAL FB210003; the Alexander von Humboldt Stiftung; an Alfred P. Sloan Research Fellowship; Allegro, the European ALMA Regional Centre node in the Netherlands, the NL astronomy research network NOVA and the astronomy institutes of the University of Amsterdam, Leiden University, and Radboud University; the ALMA North America Development Fund; the Astrophysics and High Energy Physics programme by MCIN (with funding from European Union NextGenerationEU, PRTRC17I1); 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 CASANID folio CAS220010; Chandra DD718089X and TM617006X; the China Scholarship Council; the China Postdoctoral Science Foundation fellowships (2020M671266, 2022M712084); Consejo Nacional de Humanidades, Ciencia y Tecnología (CONAHCYT, Mexico, projects U0004246083, U0004259839, F0003272050, M0037279006, F0003281692, 104497, 275201, 263356); the Colfuturo Scholarship; the Consejería de Economía, Conocimiento, Empresas y Universidad of the Junta de Andalucía (grant P18FR1769), 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émicoUniversidad Nacional Autónoma de México (DGAPAUNAM, 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 MegaScience, 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 PID2019108995GBC22; 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/011838; 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) GrantinAid 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 QYZDJSSWSLH057, QYZDJSSWSYS008, ZDBSLYSLH011); the Leverhulme Trust Early Career Research Fellowship; the MaxPlanckGesellschaft (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 (1032119M001010MY2, 1052112M001025MY3, 1052119M001042, 1062112M001011, 1062119M001013, 1062119M001027, 1062923M001005, 1072119M001017, 1072119M001020, 1072119M001041, 1072119M110005, 1072923M001009, 1082112M001048, 1082112M001051, 1082923M001002, 1092112M001025, 1092124M001005, 1092923M001001, 1102112M001033, 1102124M001007, 1102923M001001, and 1122112M003010MY3); the National Science and Technology Council (NSTC) of Taiwan (1112124M001005, 1122124M001014); 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 HSTHF251431.001A, HSTHF251482.001A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS526555; 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 1112112M001041, NSTC 1112124M001005, NSTC 1122124M001014); the National Science Foundation (NSF, grants AST0096454, AST0352953, AST0521233, AST0705062, AST0905844, AST0922984, AST1126433, AST1140030, DGE1144085, AST1207704, AST1207730, AST1207752, MRI1228509, OPP1248097, AST1310896, AST1440254, AST1555365, AST1614868, AST1615796, AST1715061, AST1716327, OISE1743747, AST1816420, AST1935980, AST2034306, AST2307887); NSF Astronomy and Astrophysics Postdoctoral Fellowship (AST1903847); 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 ScholarshipsDoctoral Program); the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF2015H1A2A1033752, the Korea Research Fellowship Program: NRF2015H1D3A1066561, 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 201700648); 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 PGC2018098915BC21, AYA201680889P, PID2019108995GBC21, PID2020117404GBC21); 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 (JCYJSHFY2021013); 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 (SEV20170709); the Spanish Ministry for Science and Innovation grant CEX2021001131S funded by MCIN/AEI/10.13039/501100011033; the Spinoza Prize SPI 78409; 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. DEAC0500OR22725; the ASTROVIVES FEDER infrastructure, with project code IDIFEDER2021086; 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/011838. APEX is a collaboration between the MaxPlanckInstitut 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 30m telescope on Pico Veleta, Spain is operated by IRAM and supported by CNRS (Centre National de la Recherche Scientifique, France), MPG (MaxPlanckGesellschaft, 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 OPP1852617 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 ACI1548562, and CyVerse, supported by NSF grants DBI0735191, DBI1265383, and DBI1743442. XSEDE Stampede2 resource at TACC was allocated through TGAST170024 and TGAST080026N. 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 LargeScale Community Partnerships allocation AST20023. Frontera is made possible by National Science Foundation award OAC1818253. 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 MarkI 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 opensource license by the Collaboration for Astronomy Signal Processing and Electronics Research (CASPER). The EHT project is grateful to T4Science and Microsemi for their assistance with hydrogen masers. This research has made use of NASA’s Astrophysics Data System. We gratefully acknowledge the support provided by the extended staff of the ALMA, from the inception of the ALMA Phasing Project through the observational campaigns of 2017 and 2018. We would like to thank A. Deller and W. Brisken for EHTspecific support with the use of DiFX. We 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
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [Google Scholar]
 Akiyama, K., Lu, R.S., Fish, V. L., et al. 2015, ApJ, 807, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, AJ, 153, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, ApJ, 838, 1 [NASA ADS] [CrossRef] [Google Scholar]
 An, T., Sohn, B. W., & Imai, H. 2018, Nat. Astron., 2, 118 [Google Scholar]
 Arras, P., Frank, P., Haim, P., et al. 2022, Nat. Astron., 6, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Asada, K., & Nakamura, M. 2012, ApJ, 745, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Bardeen, J. M. 1973, Les Astres Occlus, 215 [Google Scholar]
 Benkevitch, L., Akiyama, K., Lu, R., Doeleman, S., & Fish, V. 2016, ArXiv eprints [arXiv:1609.00055] [Google Scholar]
 Bird, S., Harris, W. E., Blakeslee, J. P., & Flynn, C. 2010, A&A, 524, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blackburn, L., Chan, C.K., Crew, G. B., et al. 2019, ApJ, 882, 23 [Google Scholar]
 Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020, ApJ, 894, 31 [Google Scholar]
 Blakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694, 556 [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., Meier, D., & Readhead, A. 2019, ARA&A, 57, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Bower, G. C., Dexter, J., Markoff, S., et al. 2015, ApJ, 811, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., & Pesce, D. W. 2020, ApJ, 904, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., Gold, R., Karami, M., et al. 2020a, ApJ, 897, 139 [Google Scholar]
 Broderick, A. E., Pesce, D. W., Tiede, P., Pu, H.Y., & Gold, R. 2020b, ApJ, 898, 9 [Google Scholar]
 Broderick, A. E., Pesce, D. W., Gold, R., et al. 2022a, ApJ, 935, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., Tiede, P., Pesce, D. W., & Gold, R. 2022b, ApJ, 927, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Sci. Comput., 16, 1190 [Google Scholar]
 Cantiello, M., Blakeslee, J. P., Ferrarese, L., et al. 2018, ApJ, 856, 126 [Google Scholar]
 Carilli, C. L., & Thyagarajan, N. 2022, ApJ, 924, 125 [NASA ADS] [CrossRef] [Google Scholar]
 CASA Team (Bean, B., et al.) 2022, PASP, 134, 114501 [NASA ADS] [CrossRef] [Google Scholar]
 Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23 [Google Scholar]
 Chael, A., Bouman, K., Johnson, M., et al. 2019a, https://doi.org/10.5281/zenodo.2614016 [Google Scholar]
 Chael, A. A., Bouman, K. L., Johnson, M. D., et al. 2019b, Astrophysics Source Code Library [record ascl:1904.004] [Google Scholar]
 Chael, A., Johnson, M. D., & Lupsasca, A. 2021, ApJ, 918, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, M.T., Asada, K., Matsushita, S., et al. 2023, PASP, 135, 095001 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, I., Jung, T., Zhao, G.Y., et al. 2017, PASJ, 69, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Cho, I., Zhao, G.Y., Kawashima, T., et al. 2022, ApJ, 926, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Clark, B. G. 1980, A&A, 89, 377 [NASA ADS] [Google Scholar]
 Cui, Y.Z., Hada, K., Kino, M., et al. 2021, Res. Astron. Astrophys., 21, 205 [CrossRef] [Google Scholar]
 Cui, Y., Hada, K., Kawashima, T., et al. 2023, Nature, 621, 711 [NASA ADS] [CrossRef] [Google Scholar]
 Cunningham, C. T., & Bardeen, J. M. 1973, ApJ, 183, 237 [Google Scholar]
 de Gasperin, F., Orrú, E., Murgia, M., et al. 2012, A&A, 547, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Deller, A. T., Brisken, W. F., Phillips, C. J., et al. 2011, PASP, 123, 275 [Google Scholar]
 Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517 [NASA ADS] [CrossRef] [Google Scholar]
 Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Science, 338, 355 [NASA ADS] [CrossRef] [Google Scholar]
 EHT MWL Science Working Group (Algaba, J. C., et al.) 2021, ApJ, 911, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Emami, R., Ricarte, A., Wong, G. N., et al. 2023, ApJ, 950, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019a, ApJ, 875, L1 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019b, ApJ, 875, L2 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019c, ApJ, 875, L3 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019d, ApJ, 875, L4 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019e, ApJ, 875, L5 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019f, ApJ, 875, L6 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021a, ApJ, 910, L12 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021b, ApJ, 910, L13 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022a, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022b, ApJ, 930, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022c, ApJ, 930, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022d, ApJ, 930, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022e, ApJ, 930, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022f, ApJ, 930, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2023, ApJ, 957, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Falcke, H., Melia, F., & Agol, E. 2000, ApJ, 528, L13 [Google Scholar]
 Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31P [Google Scholar]
 Ferreira, P. G. 2019, ARA&A, 57, 335 [Google Scholar]
 Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444 [Google Scholar]
 Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119 [Google Scholar]
 Georgiev, B., Pesce, D. W., Broderick, A. E., et al. 2022, ApJ, 930, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Goddi, C., MartíVidal, I., Messias, H., et al. 2019, PASP, 131, 075003 [Google Scholar]
 Goddi, C., MartíVidal, I., Messias, H., et al. 2021, ApJ, 910, L14 [NASA ADS] [CrossRef] [Google Scholar]
 GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 [Google Scholar]
 Greisen, E. W. 2003, Astrophys. Space Sci. Lib., 285, 109 [NASA ADS] [Google Scholar]
 Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131 [Google Scholar]
 Hilbert, D. 1917, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen – Mathematischphysikalische Klasse (Berlin: Weidmannsche Buchhandlung), 53 [Google Scholar]
 Hoffman, M. D., & Gelman, A. 2014, J. Mach. Learn. Res., 15, 1593 [Google Scholar]
 Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
 Hovatta, T., Aller, M. F., Aller, H. D., et al. 2014, AJ, 147, 143 [Google Scholar]
 Inoue, M., AlgabaMarcos, J. C., Asada, K., et al. 2014, Radio Sci., 49, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019a, A&A, 626, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Janssen, M., Blackburn, L., Issaoun, S., et al. 2019b, A Priori Calibration of EHT Stations, Tech. Rep. [Google Scholar]
 Janssen, M., Radcliffe, J. F., & Wagner, J. 2022, Universe, 8, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, M. D., Lupsasca, A., Strominger, A., et al. 2020, Sci. Adv., 6, eaaz1310 [NASA ADS] [CrossRef] [Google Scholar]
 Junor, W., Biretta, J. A., & Livio, M. 1999, Nature, 401, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018a, A&A, 616, A188 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, J.Y., Lee, S.S., Hodgson, J. A., et al. 2018b, A&A, 610, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koay, J. Y., RomeroCañizales, C., Matthews, L. D., et al. 2023a, ArXiv eprints [arXiv:2312.03505] [Google Scholar]
 Koay, J. Y., Asada, K., Matsushita, S., et al. 2023b, ArXiv eprints [arXiv:2312.02759] [Google Scholar]
 Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511 [Google Scholar]
 Lee, S.S., Wajima, K., Algaba, J.C., et al. 2016, ApJS, 227, 8 [Google Scholar]
 Liepold, E. R., Ma, C.P., & Walsh, J. L. 2023, ApJ, 945, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, D. C., & Nocedal, J. 1989, Math. Program., 45, 503 [CrossRef] [Google Scholar]
 Lockhart, W., & Gralla, S. E. 2022, MNRAS, 509, 3643 [Google Scholar]
 Lu, R., Asada, K., Krichbaum, T. P., et al. 2023, Nature, 616, 686 [CrossRef] [Google Scholar]
 Luminet, J.P. 1979, A&A, 75, 228 [NASA ADS] [Google Scholar]
 MartíVidal, I., Roy, A., Conway, J., & Zensus, A. J. 2016, A&A, 587, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018, PASP, 130, 015002 [Google Scholar]
 Mertens, F., Lobanov, A. P., Walker, R. C., & Hardee, P. E. 2016, A&A, 595, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miyoshi, M., Kato, Y., & Makino, J. 2022, ApJ, 933, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Mogensen, P. K., & Riseth, A. N. 2018, J. Open Source Softw., 3, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146 [Google Scholar]
 Owen, F. N., Hardee, P. E., & Cornwell, T. J. 1989, ApJ, 340, 698 [NASA ADS] [CrossRef] [Google Scholar]
 Park, J., Hada, K., Kino, M., et al. 2019, ApJ, 887, 147 [Google Scholar]
 Pordes, R., Petravick, D., Kramer, B., et al. 2007, J. Phys. Conf. Ser., 78, 012057 [NASA ADS] [CrossRef] [Google Scholar]
 Readhead, A. C. S., Walker, R. C., Pearson, T. J., & Cohen, M. H. 1980, Nature, 285, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Biretta, J. A., Junor, W., Muxlow, T. W. B., & Spencer, R. E. 1989, ApJ, 336, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, A. E. E., Hinteregger, H. F., Whitney, A. R., et al. 1974, ApJ, 193, 293 [Google Scholar]
 Rogers, A. E. E., Doeleman, S. S., & Moran, J. M. 1995, AJ, 109, 1391 [NASA ADS] [CrossRef] [Google Scholar]
 Satapathy, K., Psaltis, D., Özel, F., et al. 2022, ApJ, 925, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Schwarzschild, K. 1916, Abh. Konigl. Preuss. Akad. Wissenschaften Jahre, 189 [Google Scholar]
 Sfiligoi, I., Bradley, D. C., Holzman, B., et al. 2009, WRI World Congress on Computer Science and Information Engineering, 428 [Google Scholar]
 Shepherd, M. 2011, Astrophysics Source Code Library [record ascl:1103.001] [Google Scholar]
 Shepherd, M. C. 1997, ASP Conf. Ser., 125, 77 [Google Scholar]
 Simon, D. A., Cappellari, M., & Hartke, J. 2024, MNRAS, 527, 2341 [Google Scholar]
 Skilling, J. 2006, Bayesian Anal., 1, 833 [Google Scholar]
 Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
 Thompson, A. R., Moran, J. M., & Swenson, G. W. 2017, Interferometry and Synthesis in Radio Astronomy, 3rd edn. (Springer) [CrossRef] [Google Scholar]
 Tiede, P. 2022, J. Open Source Softw., 7, 4457 [NASA ADS] [CrossRef] [Google Scholar]
 Tiede, P., Broderick, A. E., & Palumbo, D. C. M. 2022a, ApJ, 925, 122 [CrossRef] [Google Scholar]
 Tiede, P., Johnson, M. D., Pesce, D. W., et al. 2022b, Galaxies, 10, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Tiede, P., Broderick, A. E., Palumbo, D. C. M., & Chael, A. 2022c, ApJ, 940, 182 [NASA ADS] [CrossRef] [Google Scholar]
 van Bemmel, I. M., Kettenis, M., Small, D., et al. 2022, PASP, 134, 114502 [CrossRef] [Google Scholar]
 Vertatschitsch, L., Primiani, R., Young, A., et al. 2015, PASP, 127, 1226 [NASA ADS] [CrossRef] [Google Scholar]
 von Laue, M. 1921, Relativitätstheorie (Braunschweig: Friedrich Vieweg& Sohn) [Google Scholar]
 Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128 [Google Scholar]
 Wielgus, M. 2021, Phys. Rev. D, 104, 124058 [NASA ADS] [CrossRef] [Google Scholar]
 Wielgus, M., Blackburn, L., Issaoun, S., et al. 2019, EHT Data Set Validation and Characterization of Errors, Tech. Rep. 2019CE02, EHT Memo Series [Google Scholar]
 Wielgus, M., Akiyama, K., Blackburn, L., et al. 2020, ApJ, 901, 67 [Google Scholar]
 Xu, K., Ge, H., Tebbutt, W., et al. 2020, Symposium on Advances in Approximate Bayesian Inference, PMLR, 1 [Google Scholar]
 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):
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 g_{E} is the normalized elevationdependent gain curve of the telescope. η_{ph} is the timedependent 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:
where A_{geom} is the geometrical area of the antenna dish in units of m^{2}, and k = 1.38 × 10^{3} is the Boltzmann constant in units of Jy m^{2}/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 suboptimal 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.
EHT station flux density calibration parameters and their uncertainties for the 2018 observing campaign.
Appendix B: Nonclosing systematic error budget
The closure phases (ψ_{C}), log closure amplitudes (lnA_{C}), closure trace phases and closure trace log amplitudes 𝒯, the latter two of which are novel diagnostics described in detail by Broderick & Pesce (2020), are independent of stationbased gain errors, but are sensitive to systematic nonclosing 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 colocated 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 nonclosing systematic errors can lead to deviations from the expected normal distribution, as can be seen for the blue histograms in Fig. B.1.
The nonclosing systematic error, s, can thus be determined by solving for its value while enforcing the condition of:
where σ is the total uncertainty associated with the random variable X, and σ_{th} is the known a priori thermal uncertainty. In this construction, X is a random variable representing the closure products being examined, for instance band 1 − band 2 closure phase difference, RR−LL log closure amplitude difference, etc. mad_{0} is the modified median absolute deviation given by:
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 nonclosing 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.
Nonclosing 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.
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 A_{scan}/A_{2s} for each scan (M 87^{∗} 2017 III). A_{scan} is the debiased amplitude of the coherently averaged visibilities over the entire scan, and A_{2s} is the amplitude derived from data coherently averaged over 2s, then incoherently averaged over the entire scan.
Figure C.1 shows the cumulative distribution of the quantity A_{scan}/A_{2s} for 2342 unique scans (per band, polarization, and baseline) of the M87^{*} April 21 and 25 data with S/N > 7, common to both the HOPS and CASA pipelines. We find that 92.9% of HOPS and 85.6% of CASA scanaveraged data show decoherence losses of less than 10%. When selecting April 21 data only, 97.4% and 87.9% of HOPS and CASA scanaveraged data show better than 90% coherence. For April 25, 86.3% and 82.3% of HOPS and CASA scanaveraged data show better than 90% coherence.
Fig. C.1. Cumulative histograms of the amplitude ratios between coherent averaging for entire scans (A_{scan}), and coherent averaging for 2 s before incoherent averaging over scans (A_{2s}), 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 scanaveraged 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. Followup diagnostics found this issue to be caused by low power levels of the intermediate frequency (IF) signal and too large blockdown 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° peaktopeak 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 postprocessing.
While the above problematic L1 data have either been corrected or excised during postprocessing, 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 modeldependent selfcalibration 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 wellcharacterized and that the amplitude calibration uncertainties are relatively small in comparison, we find the amplitudes on LMT baselines to be undercalibrated 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 undercalibrated, 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 stationbased gains derived from modeldependent calibration, that is selfcalibration 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.
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 eastwest direction, in units of wavelength, demonstrating the undercalibration 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 quasisimultaneous multiwavelength 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 intrasite baseline, SMT–LMT, which is necessary for this analysis. Thus we can derive the three constraints summarized below only for the first two observing days. For reference, the total flux density F_{tot} estimated from the ALMA interferometric data is 1.13 Jy for both April 21 and 22. This value is an average of the total flux values derived for each of the bands 1–4.
Constraint 1. The first constraint is based on the fact that the visibility amplitudes on short baselines can be approximated by a circular Gaussian visibility function,
where V_{0} is the total flux density of the Gaussian source, u is the length of the baseline, and θ is its FWHM in radians. The intermediatetolong 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 finescale source structure:
Here, 𝒱_{i − j} denotes the true visibility on the baseline i − j. Considering the gain uncertainties, the ratio of the true visibility amplitudes is lowerbounded by
where V_{i − j} and Δg_{i} denote the calibrated measured visibility and the gain uncertainty, respectively. We take the median of the visibility amplitude ratio from the collection of constraints derived for each single VLBI scan to derive a robust estimate of the minimum source size.
Constraint 2. The second constraint comes from the curvature of visibility amplitudes between the intrasite baselines and the LMT–SMT baseline. Since the compact flux density should not exceed the total flux density measured with the intrasite baselines, the amplitude drop from the intrasite 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
When we consider the gain uncertainties,
We assumed the uncertainty of the total flux Δg_{tot} = 0.1. Similarly to Constraint 1, the median of the ratio is adopted to mitigate the effects of statistical errors.
Constraint 3. The minimum compact flux density can be derived by the maximum amplitudes of LMT–SMT baseline, since the visibility amplitudes are maximum at zero baseline length and the a priori calibration should not underestimate the station sensitivity. Therefore the compact flux density is constrained as
The equivalent circular Gaussian gives a stronger constraint with the minimum source size (Constraint 1) extrapolated from the LMT–SMT baselines, described by,
These constraints are validated using the various synthetic models described in Appendix G. The upper and lower limits on the source size obtained from Constraints 1 and 2 are along the directions of LMT–SMT and LMT–ALMA baselines.
In summary, Constraints 1−3 lead to the conclusion that 0.30 Jy ≤F_{cpct} ≤ 1.13 Jy on April 21 and 22. If we use Eq. F.7, the minimum compact flux is 0.38 Jy. The compact source size has an equivalent Gaussian FWHM ranging between 39 and 98 μas on the first two observing days.
F.2. Constraints from quasisimultaneous multiwavelength observations
Supplementary, we also make use of quasisimultaneous VLBI data at longer wavelengths obtained through part of a large MWL observing campaign coordinated with the EHT 2018 observations. Our approach assumes that the total compact emission of M87^{*} at ν ≈ 230 GHz (λ ≈ 1.3 mm) F_{tot, 230} originates from the combination of (i) a compact emitter with unknown flux density F_{cpct, 230} and spectral index at 230 GHz, and (ii) a large and diffuse (milliarcsecondscale) jet, which can be characterized by its total flux density F_{jet, 230} = F_{tot, 230} − F_{cpct, 230} and a steep spectral index of α_{jet} ∼ −(0.7 − 1.0), as measured at cm and mmwavelengths (with typical errors of ∼0.3; see Hovatta et al. 2014; Hada et al. 2016). Then, by measuring F_{jet} and F_{tot} at lower frequencies, and extrapolating them to 230 GHz using a powerlaw model, we can estimate the total compact flux density F_{cpct} at 230 GHz: F_{cpct, 230} = F_{tot, 230} − F_{jet, 230}.
We utilize the data obtained with the East Asian VLBI Network (EAVN; An et al. 2018; Cui et al. 2021; Cho et al. 2022) at 22 and 43 GHz, the Korean VLBI Network (KVN; Lee et al. 2016) at 43, 86, and 129 GHz, and the Global Millimeter VLBI Array (GMVA)+ALMA at 86 GHz. All of these observations were conducted between 2018 March 25 and April 21. For details of data reduction and calibration, see Cho et al. (2017, 2021) for EAVN, Lee et al. (2016), Kim et al. (2018b) for KVN, and Lu et al. (2023) for GMVA+ALMA, respectively.
The total flux densities of the nuclear region (i.e., F_{tot} = F_{jet} + F_{cpct}) at 22, 43, and 86 GHz are measured from the EAVN and KVN images by integrating the flux densities over a large window of 500 × 500 μas^{2} centered on the radio core. We obtain total flux densities of F_{tot, 22} ≈ 1.1 Jy, F_{tot, 43} ≈ 0.9 − 1.1 Jy, F_{tot, 86} ≈ 0.9 Jy and F_{tot, 129} ≈ 0.6 − 0.8 Jy at 22, 43, 86, and 129 GHz, respectively. The KVN 129 GHz data suffer from larger uncertainties due to a stronger influence from weather. We fit a powerlaw model to the total flux densities (i.e., F_{tot} ∝ ν^{+αtot}). We find α_{tot} ≈ −0.17 when 129 GHz data are included, while α_{tot} ≈ −0.11 when they are excluded. By extrapolating a powerlaw model with a range α_{tot} ≈ −(0.11 − 0.17) to 230 GHz, we estimate F_{tot, 230} = 0.72 − 0.80 Jy. While this range of values is lower than the ALMA interferometric measurements during the EHT observations (∼1.11−1.18 Jy; Table 3 in Goddi et al. 2021), it only includes the total flux density within the inner 500 × 500 μas^{2}, which is nearly three orders of magnitude below the resolution limit of ALMA.
The compact flux density F_{cpct, 86} at 86 GHz is estimated from a GMVA+ALMA image that has a higher angular resolution than the KVN images (Lu et al. 2023). We obtain F_{cpct, 86} ≈ 0.51 Jy by integrating over a smaller window of 100 × 100 μas^{2} centered on the radio core. We can then derive a jet flux density at 86 GHz to be F_{jet, 86} = F_{tot, 86} − F_{cpct, 86} ≈ 0.39 Jy.
Finally, we compute the expected F_{cpct, 230} by estimating F_{jet, 230} and subtracting it from the previously estimated F_{tot, 230}. As explained above, we assume that F_{jet, 230} follows a power law that can be extrapolated up to 230 GHz (i.e., F_{jet, 230} ∝ ν^{+αjet}). For a range of α_{jet} = −(0.7 − 1.0), we find F_{jet, 230} = 0.15 − 0.20 Jy. This value gives a range of F_{cpct, 230} ≈ 0.5 − 0.7 Jy within the central 100 × 100 μas^{2} region. We emphasize that this estimate is based on lowerfrequency data and assumptions on the jet spectral index. Therefore it should not be regarded as a tight constraint on F_{cpct, 230} but more as a useful reference.
Appendix G: Imaging validation
Here we describe the synthetic data generation process and subsequent reconstructions, as well as the resulting distribution of Top Set images for the real M87^{*} data. We also investigate the properties of the image reconstructions when removing individual baselines.
G.1. Synthetic data
In order to test the efficacy of the imaging pipelines, we generate synthetic data sets using both geometric models and snapshots of GRMHD simulations. All of the geometric models contain 0.6 Jy in compact flux density and 0.5 Jy in a larger extended jet modeled by three Gaussians. We show the ring model as an example in Fig. G.1. The sizes and locations of these three Gaussians are identical to those reported in Table 10 of M 87^{∗} 2017 IV, with the flux densities scaled down to total 0.5 Jy rather than 0.6 Jy. Of the geometric models we generated, four of them were utilized in the procedure to generate the Top Set images. These four geometric models are very similar to those used in M 87^{∗} 2017 IV and are as follows:

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

2.
ring: A thin uniform ring of radius r_{0} = 23 μas blurred by a circular Gaussian beam of FWHM 10 μas.

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

4.
disk: A uniform disk of radius r_{0} = 35 μas blurred by a circular Gaussian beam of FWHM 10 μas.
Groundtruth images of these models are shown in Fig. G.2. Aside from these training data, we also generated seven validation data sets to ensure good performance from the CLEAN and RML method Top Sets. These data, along with reconstructions from each pipeline, are shown in Fig. G.3. The model descriptions are as follows

1.
cres–90, cres0, cres90: Three asymmetric ring models with r_{0} = 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 ehtimaging and all include random thermal noise, stationbased gain and phase errors, and polarimetric leakage. For a more indepth 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, ehtimaging, 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, ehtimaging, 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.
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, ehtimaging, 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 ehtimaging and SMILI show doublering structures although they are a minority in the Top Sets. As discussed above, we see similar minor doublering 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 ehtimaging 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 ehtimaging images reconstructed with fiducial parameters, while the third one presents ehtimaging 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 ehtimaging 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 ehtimaging and in a loss of contrast for Comrade. Removing PV baselines results in an elongation of the ringlike 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 ehtimaging reconstructions without one antenna are affected by a stronger tessellation pattern. However, the morphology of the ringlike structure is more robust to the removal of a single site in 2018 compared to 2017. For example, in 2018, a ringlike 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.
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 groundtruth images. The DIFMAP, ehtimaging, 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. 
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 groundtruth images. The DIFMAP, ehtimaging, 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. 
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 groundtruth structure of each synthetic data model is shown in gray contours. 
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. 
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. 
Fig. G.7. Twenty five randomly selected Top Set images from HOPS band 3 April 21 data for DIFMAP (upper left), ehtimaging (upper right), and SMILI (bottom). 
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 ehtimaging respectively for band 3 on 2018 April 21. The images for the ehtimaging pipeline were reconstructed using the ehtimaging 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 ehtimaging 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 selfcalibration
This Appendix reports the results of the DIFMAP parameter survey, which used an initial phaseonly selfcalibration 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 selfcalibration, this survey requires a careful consideration of the CLEAN window location. This is due to the fact that the selfcalibration 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 selfcalibration 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 selfcalibration. 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 rescaling 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 ringlike morphology, consistent with the results from the survey discussed in Sect. 5.1.1. This result demonstrates that the ringlike structure for M87^{*}, reconstructed by DIFMAP, is stable independently of the geometrical model selected for the first phase of selfcalibration. However, the overall performance of this survey was lower than that of the survey that optimized the geometric models for initial phaseonly selfcalibration, resulting in a smaller number of Top Set images. This can be attributed to the fact that the point source model used as initial phaseonly selfcalibration 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.
Fig. H.1. Fiducial images produced by the DIFMAP pipeline reconstructed through the implementation of an initial phaseonly selfcalibration process using a point source model. The reconstructed images exhibit ringlike structures, which are comparable to the results obtained from the pipeline that utilized closure phases to select the optimal geometric model for the initial phaseonly selfcalibration, 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, lowband) 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.
Fig. I.1. Fiducual images of the 2017 April 11 lowband 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 phaseonly selfcalibration by employing a point source model, whereas the bottom row of images is obtained from the pipeline that identifies the optimal geometric model for selfcalibration. 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 largescale asymmetric Gaussian, and a thin slashed crescent. A more detailed description of the largescale Gaussian and crescent models, as well as details about their technical implementation, are described in Broderick et al. (2020a). A much more detailed discussion of the raster model can be found in Broderick et al. (2020b).
The splined raster model is defined by a rectilinear set of control points which may independently vary in intensity, I_{M, N}. The intensity map is produced using an approximate cubic spline interpolation between the control points. The field of view FOV_{x} and FOV_{y} are permitted to vary, along with the orientation of the grid ϕ and the raster center (x, y). The largescale Gaussian is asymmetric, and is characterized by its total flux I_{G}, a symmetric standard deviation σ_{G}, asymmetry parameter A_{G}, position angle ϕ_{G}, and position (x_{G}, y_{G}). The slashed crescent is based on the xsringauss model also described in Broderick et al. (2020b), which also serves as the core of the THEMIS xsringauss geometric model featured in Sect. 6.2 and in M 87^{∗} 2017 VI.
As a part of the Hybrid Themage model, we eliminate the contributions from the asymmetry and pinned Gaussian parameters by appropriately restricting their priors. We also restrict the priors on the fractional width parameter ψ to range from 0% to 5%. This ensures the ring width is well below that of the raster spacing, to avoid excessive flux trading between the ring and raster components. This leaves the ring flux I_{X}, outer radius R_{out}, linear brightness gradient f_{X}, brightness gradient position angle ϕ_{X}, and position (x_{X}, y_{X}) as searchable parameters.
We use a deterministic evenodd (DEO) swap tempering scheme with a Hamiltonian MonteCarlo sampling kernel based on the Stan package. Chain convergence is assessed based on traditional criteria such as the split , autocorrelation time, and visual inspection of the parameter traces. The parameters and priors for each model component are listed in Table J.1.
THEMIS Raster Model Priors.
Appendix K: THEMIS Raster dimensions
THEMIS raster size (N_{x} × N_{y}) survey for April 21 Band 3.
The best fit raster dimensions N_{x} × N_{y} for the THEMIS Raster models in principle depend on the (u, v) coverage and source complexity. To construct the best model, one should survey the raster dimensions for each unique data set. In order to provide easier comparisons across data sets, we only survey the raster dimensions on our primary science data set, April 21 band 3, and use those best fit raster dimensions for the other data sets. We present the results of this survey in Table K.1, where we survey a 4 × 4, 5 × 5, 6 × 6, and 7 × 7 raster. We rank the models in terms of the difference in the logarithm of the Bayesian Evidence (Δlog(Z)), where positive values are more preferred. The 5 × 5 raster is most preferred by this data, followed by the 4 × 4 raster. The 5 × 5 raster is also the same dimensions as the model used in Broderick et al. (2022a), which applied the Hybrid Themage model to the 2017 data. Thus, the image structures we show in this analysis should be directly comparable to the image structures in the previous paper. We also report the best fit reduced χ^{2} for the surveyed raster sizes, and note that all models provide good fits to the data.
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 rasteronly 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 Bspline kernel of order 3 (see Eqs. L.1 and L.2) to generate flux densities that smoothly vary in all directions. We set the raster points by (l_{i}, m_{j}) given by Eq. L.3, for a given FOV and number of pixels N_{x} and N_{y} for each side. The FOV and number of pixels are the hyperparameters of the model:
where B(x) is a Bspline kernel of order 3 given by,
and raster points are given by,
where s_{x} = FOV/N_{x} and s_{y} = FOV/N_{y}. We restrict the images considered in this paper to having N_{x} = N_{y}, such that the images are square and s_{x} = s_{y}. We add a circular Gaussian with FWHM = mas (G_{1mas}) to the raster to model the amplitudes and station gains of short intrasite (ALMA–APEX, JCMT–SMA) baselines. We form our final image model,
where f is the total flux density of the raster (also described as compact flux) and f_{g} is the fraction of the total flux density that corresponds to the flux density of the Gaussian. We take the Fourier transform of M(I_{ij}, f, f_{g}) to get model Fourier transform 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 f_{g} for the flux density of the largescale Gaussian. For the raster pixels fluxes, I_{i}, we use a symmetric Dirichlet distribution
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 ∑_{i}I_{i} = 1 and 0 ≤ I_{i} ≤ 1. For this work, we set ξ = 1, which is equivalent to a uniform distribution on the simplex.
From the analysis of the crossing baseline tracks and the EHT station flux density calibration parameters (Table A.1), for station gain priors, we used a Normal distribution for the loggain 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 nonimaging constraints on compact flux density permit a wide range of allowable values, and a comparison of the compact flux densities measured by the Bayesian methods suggests that the current data set cannot be used to constrain the compact flux on horizon scales. In this appendix, we test the robustness of the ring reconstructions from the RML and CLEAN methods against different assumptions of the compact flux density.
Firstly, we investigated the dependence of the assumed total compact flux density (F_{cpct}) on the estimated ring parameters of M87^{*} from reconstructed images following the approach adopted in our previous EHT analysis (M 87^{∗} 2017 IV). In this analysis, we only change F_{cpct}, while all other imaging parameters are kept at their fiducial values. We select acceptable images based on two criteria: (1) normalized , (with a 10% systematic uncertainty for DIFMAP), and (2) a corresponding selfcalibration solution with 0.9 ≤ median(1/g_{SMT}) ≤ 1.1, where g_{SMT} is the gain for the SMT station. Figure M.1 illustrates how the estimated ring parameters vary with the assumed total compact flux density. While the ring width for all pipelines and fractional central brightness for RML methods slightly increase with F_{cpct}, the ring diameter and position angle values of the acceptable images do not exhibit significant dependence on F_{cpct}. Therefore, the ring diameter and position angle are robust against different assumptions on the compact total flux density.
Secondly, We tested if our DIFMAP imaging script could distinguish different morphologies for different values of the total compact flux density using synthetic data. We prepared a series of synthetic data with the same morphologies of four geometric models used for Top Set selection in the imaging survey, but with the assumed compact total flux density in range from 0.4 to 1.1 Jy. Here, we used the total flux density, which is the sum of the assumed compact total flux density and flux density of a larger extended jet modeled by three Gaussians, to be 1.1 Jy for all synthetic data. We conduct image reconstruction using the DIFMAP imaging script with the fiducial parameter sets, but adopting the compact flux densities to be the same as that of the assumed values of synthetic data. In Fig. M.2, we show the reconstructed images of geometric models with different values of the total compact flux density. It is clearly demonstrated that our DIFMAP imaging script can distinguish ringlike morphologies from disk or double source morphologies if we properly choose the expected total compact flux density.
Fig. M.1. Measured diameter d, width w, position angle, and fractional central brightness f_{C} of M87^{*}, measured from image reconstructions assuming different total compact flux densities, F_{cpct} for DIFMAP (orange), ehtimaging (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 F_{cpct} that produces images that have (1) normalized , (with 10% systematic uncertainty for DIFMAP), and (2) a corresponding selfcalibration solution with 0.9 ≤ median(1/g_{SMT}) ≤ 1.1. 
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
Fig. N.1. Schematic summary of the xsringauss and mFring model parameters. The xsringauss 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 mFring is an mring 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 mring. 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 mFring in blue. 
Image reconstructions of M87^{*} show strong support for a ringlike structure. We quantify this support through the Bayesian evidence as a robust metric for the goodness of fit.
THEMIS xsringauss model parameters and priors.
mFring parameters and priors.
Fig. N.2. Representative images for each observing band from random posterior samples for the Hybrid Themage, xsringauss, and mFring 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,
which also defines the likelihood of the data for a chosen model and parameter,
the prior probability on model parameters,
and the Bayesian evidence for a chosen model,
The Bayesian evidence allows the relative support between two models to be determined through a Bayes factor,
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,
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 mring of order 2 with a flat emission floor and two elliptical Gaussians (mFring), and an xsringauss 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 mFring and the xsringauss model are their diameters, widths, central brightnesses, and brightness asymmetry position angles. The diameter of the mFring is defined to be the sum of the debiased semimajor and semiminor axis, while that of the xsringauss is taken to be the sum of its inner and outer radii,
Here, the debiased radius () is linked to the defining radius of the mFring (r_{i}), and the FWHM of the blurring Gaussian kernel (σ^{*}) through:
We define the fractional widths of the mFring and the xsringauss respectively to be,
and the brightness asymmetry position angle to be,
The last shared parameter shown between the two models is the fractional central flux depression, or relative central brightness, which is given by,
As described above in the main text, the THEMIS xsringauss 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, xsringauss, and mFring in the visibility domain, we can also produce conventional image domain representations of the models. Figure N.2 shows the onsky representations of the Hybrid Themage, xsringauss, and mFring 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 xsringauss and mFring 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 nonring 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
Median zenith sky opacities (1.3 mm) at EHT sites during the 2018 April observations toward M 87^{*}.
EHT station flux density calibration parameters and their uncertainties for the 2018 observing campaign.
Nonclosing 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.
Complex visibility reduced χ^{2} () for the THEMIS Raster and Hybrid Themage model best fits to the HOPS pipeline data.
All Figures
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 
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. Colocated 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 dashedlines. While the SPT cannot observe M 87^{*}, it observed 3C 279 and was used to calibrate the data. 

In the text 
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 
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 lowband (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 
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 lowband (gray points). All data shown include a priori stationbased amplitude calibration and network calibration but are prior to any modeldependent selfcalibration. 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 
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 (lowband), gray circles denote data from 2017 April 11 (lowband), and blue squares denote data from 2018 April 21 (band 3). Error bars are the ±1σ uncertainties. 

In the text 
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, ehtimaging, 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, ehtimaging, 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 
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 
Fig. 9. Visibility amplitudes of band 3 data on April 21 as a function of baseline length compared with corresponding gaincalibrated visibility amplitudes from the five representative image models from each pipeline. The fiducial images are used for DIFMAP, ehtimaging, 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 selfcalibration. ehtimaging scales the intrasite baselines to 0.5 Jy, and SMILI scales the intrasite baselines to 0.6 Jy, and both precalibrate the LMTSMT baselines before fitting. DIFMAP uses 10 s averaged data, and all other methods use scanaveraged data. THEMIS and Comrade apply 1% fractional systematic noise to all baselines, and ehtimaging 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 selfcalibration 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 
Fig. 10. Visualization of image statistics calculated using the Top Set images from the ehtimaging 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 ehtimaging 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 
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 
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 
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 (ehtimaging) 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 unshaded region corresponds to the 2018 April 25 measurements. 

In the text 
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 ringlike 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 
Fig. 15. Complex visibility (V) fit comparisons and normalized residuals for the Hybrid Themage and xsringauss, and the logclosure amplitude (lnA_{C}) and closure phase (ψ_{C}) fit comparisons and normalized residuals for the mFring 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 xsringauss, and from the maximum a posteriori model for the mFring. The gray points represent the scanaveraged data, with 1% fractional systematic noise added to all baselines. 

In the text 
Fig. 16. Fitted features of the xsringauss, mFring, 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 mFring model, while the red band is the equivalent for the xsringauss. 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 xsringauss fits. The hashed region is the 68% posterior probability interval taken from the 2017 M 87^{*} xsringauss fits. 

In the text 
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 xsringauss, and the mFring 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 ehtimaging, 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 xsringauss and mFring 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 ringlike structure, from the direct modeling methods which do assume a ringlike structure. 

In the text 
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, ehtimaging, SMILI, xsringauss), we use the methodspecific 2017 scaling factor to determine the 2018 d/α_{2017} values. For DIFMAP, ehtimaging, and SMILI we use the scaling factors from Table 6 of M 87^{∗} 2017 VI, and for the xsringauss 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, mFring), 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 modelspecific 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 xsringauss 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 
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 milliarcsecondscale 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 nonGaussian, and exhibit large tails. 

In the text 
Fig. 20. Distribution of bestfit position angle (in degrees) of the forward jet by fitting model images with the 2017 (April 6, highband) and 2018 (April 21, band 4) observations. The bestfit 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 milliarcsecondscale 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 electrothermodynamics (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 
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 
Fig. 22. Triangle plot comparing the compact flux density (), raster FOV in the x and y directions, and crescent flux (I_{X}) parameters for the Hybrid Themage reconstruction of April 21 band 1. The compact flux density and raster FOV parameters exhibit clearly bimodal distributions in this parameter space. The color regions represent the 99%, 90%, and 50% quantile regions. 

In the text 
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 
Fig. C.1. Cumulative histograms of the amplitude ratios between coherent averaging for entire scans (A_{scan}), and coherent averaging for 2 s before incoherent averaging over scans (A_{2s}), 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 
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 eastwest direction, in units of wavelength, demonstrating the undercalibration of amplitudes on LMT baselines. HOPS data are shown in this figure, but are consistent with that from CASA. 

In the text 
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 
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 groundtruth images. The DIFMAP, ehtimaging, 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 
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 groundtruth images. The DIFMAP, ehtimaging, 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 
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 groundtruth structure of each synthetic data model is shown in gray contours. 

In the text 
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 
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 
Fig. G.7. Twenty five randomly selected Top Set images from HOPS band 3 April 21 data for DIFMAP (upper left), ehtimaging (upper right), and SMILI (bottom). 

In the text 
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 ehtimaging respectively for band 3 on 2018 April 21. The images for the ehtimaging pipeline were reconstructed using the ehtimaging 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 ehtimaging 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 
Fig. H.1. Fiducial images produced by the DIFMAP pipeline reconstructed through the implementation of an initial phaseonly selfcalibration process using a point source model. The reconstructed images exhibit ringlike structures, which are comparable to the results obtained from the pipeline that utilized closure phases to select the optimal geometric model for the initial phaseonly selfcalibration, as demonstrated in Fig. 7. The circle represents a Gaussian blurring kernel with a FWHM of 20 μas. 

In the text 
Fig. I.1. Fiducual images of the 2017 April 11 lowband 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 phaseonly selfcalibration by employing a point source model, whereas the bottom row of images is obtained from the pipeline that identifies the optimal geometric model for selfcalibration. The circles represent a Gaussian blurring kernel with a FWHM of 20 μas. 

In the text 
Fig. M.1. Measured diameter d, width w, position angle, and fractional central brightness f_{C} of M87^{*}, measured from image reconstructions assuming different total compact flux densities, F_{cpct} for DIFMAP (orange), ehtimaging (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 F_{cpct} that produces images that have (1) normalized , (with 10% systematic uncertainty for DIFMAP), and (2) a corresponding selfcalibration solution with 0.9 ≤ median(1/g_{SMT}) ≤ 1.1. 

In the text 
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 
Fig. N.1. Schematic summary of the xsringauss and mFring model parameters. The xsringauss 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 mFring is an mring 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 mring. 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 mFring in blue. 

In the text 
Fig. N.2. Representative images for each observing band from random posterior samples for the Hybrid Themage, xsringauss, and mFring models, for the April 21 and April 25 HOPS data. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.