Monitoring the Morphology of M87* in 2009-2017 with the Event Horizon Telescope
Maciek Wielgus, Kazunori Akiyama, Lindy Blackburn, Chi-kwan Chan, Jason Dexter, Sheperd S. Doeleman, Vincent L. Fish, Sara Issaoun, Michael D. Johnson, Thomas P. Krichbaum, Ru-Sen Lu, Dominic W. Pesce, George N. Wong, Geoffrey C. Bower, Avery E. Broderick, Andrew Chael, Koushik Chatterjee, Charles F. Gammie, Boris Georgiev, Kazuhiro Hada, Laurent Loinard, Sera Markoff, Daniel P. Marrone, Richard Plambeck, Jonathan Weintroub, Matthew Dexter, David H. E. MacMahon, Melvyn Wright, Event Horizon Telescope Collaboration
DDraft version September 25, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Monitoring the Morphology of M87 ∗ in 2009–2017 with the Event Horizon Telescope Maciek Wielgus,
1, 2
Kazunori Akiyama,
3, 4, 5, 1
Lindy Blackburn,
1, 2
Chi-kwan Chan,
6, 7
Jason Dexter, Sheperd S. Doeleman,
1, 2
Vincent L. Fish, Sara Issaoun, Michael D. Johnson,
1, 2
Thomas P. Krichbaum, Ru-Sen Lu ( 路 如 森 ),
11, 10
Dominic W. Pesce,
1, 2
George N. Wong,
12, 13
Geoffrey C. Bower, Avery E. Broderick,
15, 16, 17
Andrew Chael,
18, 19
Koushik Chatterjee, Charles F. Gammie,
12, 21
Boris Georgiev,
16, 17
Kazuhiro Hada,
22, 23
Laurent Loinard,
24, 25
Sera Markoff,
20, 26
Daniel P. Marrone, Richard Plambeck, Jonathan Weintroub,
1, 2
Matthew Dexter, David H. E. MacMahon, Melvyn Wright, Antxon Alberdi, Walter Alef, Keiichi Asada, Rebecca Azulay,
30, 31, 10
Anne-Kathrin Baczko, David Ball, Mislav Baloković,
1, 2
Enrico Barausse,
32, 33
John Barrett, Dan Bintley, Wilfred Boland, Katherine L. Bouman,
1, 2, 36
Michael Bremer, Christiaan D. Brinkerink, Roger Brissenden,
1, 2
Silke Britzen, Dominique Broguiere, Thomas Bronzwaer, Do-Young Byun,
38, 39
John E. Carlstrom,
40, 41, 42, 43
Shami Chatterjee, Ming-Tang Chen, Yongjun Chen ( 陈 永 军 ),
11, 45
Ilje Cho,
38, 39
Pierre Christian,
6, 2
John E. Conway, James M. Cordes, Geoffrey B. Crew, Yuzhu Cui,
22, 23
Jordy Davelaar, Mariafelicia De Laurentis,
47, 48, 49
Roger Deane,
50, 51
Jessica Dempsey, Gregory Desvignes, Sergio A. Dzib, Ralph P. Eatough, Heino Falcke, Ed Fomalont, Raquel Fraga-Encinas, Per Friberg, Christian M. Fromm, Peter Galison,
1, 53, 54
Roberto García, Olivier Gentaz, Ciriaco Goddi,
9, 55
Roman Gold,
56, 15
José L. Gómez, Arturo I. Gómez-Ruiz,
57, 58
Minfeng Gu ( 顾 敏 峰 ),
11, 59
Mark Gurwell, Michael H. Hecht, Ronald Hesper, Luis C. Ho ( 何 子 山 ),
61, 62
Paul Ho, Mareki Honma,
22, 23
Chih-Wei L. Huang, Lei Huang ( 黄 磊 ),
11, 59
David H. Hughes, Makoto Inoue, David J. James, Buell T. Jannuzi, Michael Janssen, Britton Jeter,
16, 17
Wu Jiang ( 江 悟 ), Alejandra Jimenez-Rosales, Svetlana Jorstad,
66, 67
Taehyun Jung,
38, 39
Mansour Karami,
15, 16
Ramesh Karuppusamy, Tomohisa Kawashima, Garrett K. Keating, Mark Kettenis, Jae-Young Kim, Junhan Kim,
6, 36
Jongsoo Kim, Motoki Kino,
5, 69
Jun Yi Koay, Patrick M. Koch, Shoko Koyama, Michael Kramer, Carsten Kramer, Cheng-Yu Kuo, Tod R. Lauer, Sang-Sung Lee, Yan-Rong Li ( 李 彦 荣 ), Zhiyuan Li ( 李 志 远 ),
73, 74
Michael Lindqvist, Rocco Lico, Kuo Liu, Elisabetta Liuzzo, Wen-Ping Lo,
29, 76
Andrei P. Lobanov, Colin Lonsdale, Nicholas R. MacDonald, Jirong Mao ( 毛 基 荣 ),
77, 78, 79
Nicola Marchili,
75, 10
Alan P. Marscher, Iván Martí-Vidal,
30, 31
Satoki Matsushita, Lynn D. Matthews, Lia Medeiros,
80, 6
Karl M. Menten, Yosuke Mizuno,
48, 81
Izumi Mizuno, James M. Moran,
1, 2
Kotaro Moriyama,
4, 22
Monika Moscibrodzka, Cornelia Müller,
10, 9
Gibwa Musoke,
20, 9
Hiroshi Nagai,
5, 23
Neil M. Nagar, Masanori Nakamura, Ramesh Narayan,
1, 2
Gopal Narayanan, Iniyan Natarajan, Antonios Nathanail, Roberto Neri, Chunchong Ni,
16, 17
Aristeidis Noutsos, Hiroki Okino,
22, 84
Héctor Olivares, Gisela N. Ortiz-León, Tomoaki Oyama, Feryal Özel, Daniel C. M. Palumbo,
1, 2
Jongho Park, Nimesh Patel, Ue-Li Pen,
15, 85, 86, 87
Vincent Piétu, Aleksandar PopStefanija, Oliver Porth,
20, 48
Ben Prather, Jorge A. Preciado-López, Dimitrios Psaltis, Hung-Yi Pu,
15, 88, 29
Venkatessh Ramakrishnan, Ramprasad Rao, Mark G. Rawlings, Alexander W. Raymond,
1, 2
Luciano Rezzolla,
48, 89
Bart Ripperda,
90, 91
Freek Roelofs, Alan Rogers, Eduardo Ros, Mel Rose, Arash Roshanineshat, Helge Rottmann, Alan L. Roy, Chet Ruszczyk, Benjamin R. Ryan,
13, 92
Kazi L. J. Rygl, Salvador Sánchez, David Sánchez-Arguelles,
63, 58
Mahito Sasada,
22, 94
Tuomas Savolainen,
95, 96, 10
F. Peter Schloerb, Karl-Friedrich Schuster, Lijing Shao,
10, 62
Zhiqiang Shen ( 沈 志强 ),
11, 45
Des Small, Bong Won Sohn,
38, 39, 97
Jason SooHoo, Fumie Tazaki, Paul Tiede,
16, 17
Remo P. J. Tilanus,
9, 55, 98, 6
Michael Titus, Kenji Toma,
99, 100
Pablo Torne,
10, 93
Tyler Trent, Efthalia Traianou, Sascha Trippe,
Shuichiro Tsuda, Ilse van Bemmel, Huib Jan van Langevelde,
68, 102
Daniel R. van Rossum, Jan Wagner, John Wardle,
Derek Ward-Thompson,
Norbert Wex, Robert Wharton, Qingwen Wu ( 吴 庆 文 ), Doosoo Yoon, André Young, Ken Young, [email protected] a r X i v : . [ a s t r o - ph . H E ] S e p Ziri Younsi,
Feng Yuan ( 袁 峰 ),
11, 59, 107
Ye-Fei Yuan ( 袁 业 飞 ), J. Anton Zensus, Guangyao Zhao, Shan-Shan Zhao,
9, 73 and Ziyan Zhu Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA National Radio Astronomy Observatory, 520 Edgemont Rd, Charlottesville, VA 22903, USA Massachusetts Institute of Technology Haystack Observatory, 99 Millstone Road, Westford, MA 01886, USA National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan Steward Observatory and Department of Astronomy, University of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA Data Science Institute, University of Arizona, 1230 N. Cherry Ave., Tucson, AZ 85721, USA JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics (IMAPP), Radboud University, P.O. Box9010, 6500 GL Nijmegen, The Netherlands Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People’s Republic of China Department of Physics, University of Illinois, 1110 West Green St, Urbana, IL 61801, USA CCS-2, Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA Institute of Astronomy and Astrophysics, Academia Sinica, 645 N. A’ohoku Place, Hilo, HI 96720, USA Perimeter Institute for Theoretical Physics, 31 Caroline Street North, Waterloo, ON, N2L 2Y5, Canada Department of Physics and Astronomy, University of Waterloo, 200 University Avenue West, Waterloo, ON, N2L 3G1, Canada Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L 3G1 Canada Princeton Center for Theoretical Science, Jadwin Hall, Princeton University, Princeton, NJ 08544, USA NASA Hubble Fellowship Program, Einstein Fellow Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 West Green Street, Urbana, IL 61801, USA Mizusawa VLBI Observatory, National Astronomical Observatory of Japan, 2-12 Hoshigaoka, Mizusawa, Oshu, Iwate 023-0861, Japan Department of Astronomical Science, The Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo181-8588, Japan Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México, Morelia 58089, México Instituto de Astronomía, Universidad Nacional Autónoma de México, CdMx 04510, México Gravitation Astroparticle Physics Amsterdam (GRAPPA) Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam,The Netherlands Radio Astronomy Laboratory, University of California, Berkeley, CA 94720, USA Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain Institute of Astronomy and Astrophysics, Academia Sinica, 11F of Astronomy-Mathematics Building, AS/NTU No. 1, Sec. 4,Roosevelt Rd, Taipei 10617, Taiwan, R.O.C. Departament d’Astronomia i Astrofísica, Universitat de València, C. Dr. Moliner 50, E-46100 Burjassot, València, Spain Observatori Astronòmic, Universitat de València, C. Catedrático José Beltrán 2, E-46980 Paterna, València, Spain SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy East Asian Observatory, 660 N. A’ohoku Place, Hilo, HI 96720, USA Nederlandse Onderzoekschool voor Astronomie (NOVA), PO Box 9513, 2300 RA Leiden, The Netherlands California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA Institut de Radioastronomie Millimétrique, 300 rue de la Piscine, F-38406 Saint Martin d’Hères, France Korea Astronomy and Space Science Institute, Daedeok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea University of Science and Technology, Gajeong-ro 217, Yuseong-gu, Daejeon 34113, Republic of Korea Kavli Institute for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA Department of Physics, University of Chicago, 5720 South Ellis Avenue, Chicago, IL 60637, USA Enrico Fermi Institute, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, NY 14853, USA Key Laboratory of Radio Astronomy, Chinese Academy of Sciences, Nanjing 210008, People’s Republic of China Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, SE-43992 Onsala,Sweden Dipartimento di Fisica “E. Pancini”, Universitá di Napoli “Federico II”, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia,I-80126, Napoli, Italy Institut für Theoretische Physik, Goethe-Universität Frankfurt, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany INFN Sez. di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT Department of Physics, University of Pretoria, Lynnwood Road, Hatfield, Pretoria 0083, South Africa Centre for Radio Astronomy Techniques and Technologies, Department of Physics and Electronics, Rhodes University, Grahamstown6140, South Africa LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195Meudon, France Department of History of Science, Harvard University, Cambridge, MA 02138, USA Department of Physics, Harvard University, Cambridge, MA 02138, USA Leiden Observatory—Allegro, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands CP3-Origins, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark Instituto Nacional de Astrofísica, Óptica y Electrónica, Luis Enrique Erro 1, Tonantzintla, Puebla, C.P. 72840, México Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582, 03940, Ciudad de México, México Key Laboratory for Research in Galaxies and Cosmology, Chinese Academy of Sciences, Shanghai 200030, People’s Republic of China NOVA Sub-mm Instrumentation Group, Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD Groningen,The Netherlands Department of Astronomy, School of Physics, Peking University, Beijing 100871, People’s Republic of China Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People’s Republic of China Instituto Nacional de Astrofísica, Óptica y Electrónica. Apartado Postal 51 y 216, 72000. Puebla Pue., México ASTRAVEO LLC, PO Box 1668, MA 01931 Max-Planck-Institut für Extraterrestrische Physik, Giessenbachstr. 1, D-85748 Garching, Germany Institute for Astrophysical Research, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA Astronomical Institute, St. Petersburg University, Universitetskij pr., 28, Petrodvorets,198504 St.Petersburg, Russia Joint Institute for VLBI ERIC (JIVE), Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands Kogakuin University of Technology & Engineering, Academic Support Center, 2665-1 Nakano, Hachioji, Tokyo 192-0015, Japan Physics Department, National Sun Yat-Sen University, No. 70, Lien-Hai Rd, Kaosiung City 80424, Taiwan, R.O.C National Optical Astronomy Observatory, 950 North Cherry Ave., Tucson, AZ 85719, USA Key Laboratory for Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan Road,Shijingshan District, Beijing, People’s Republic of China School of Astronomy and Space Science, Nanjing University, Nanjing 210023, People’s Republic of China Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Nanjing 210023, People’s Republic of China Italian ALMA Regional Centre, INAF-Istituto di Radioastronomia, Via P. Gobetti 101, I-40129 Bologna, Italy Department of Physics, National Taiwan University, No.1, Sect.4, Roosevelt Rd., Taipei 10617, Taiwan, R.O.C. Yunnan Observatories, Chinese Academy of Sciences, 650011 Kunming, Yunnan Province, People’s Republic of China Center for Astronomical Mega-Science, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing, 100012, People’sRepublic of China Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, 650011 Kunming, People’s Republicof China School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA Tsung-Dao Lee Institute, Shanghai Jiao Tong University, Shanghai, 200240, People’s Republic of China Astronomy Department, Universidad de Concepción, Casilla 160-C, Concepción, Chile Department of Astronomy, University of Massachusetts, 01003, Amherst, MA, USA Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada Canadian Institute for Advanced Research, 180 Dundas St West, Toronto, ON M5G 1Z8, Canada Department of Physics, National Taiwan Normal University, No. 88, Sec.4, Tingzhou Rd., Taipei 116, Taiwan, R.O.C. School of Mathematics, Trinity College, Dublin 2, Ireland Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos, NM, 87545, USA Instituto de Radioastronomía Milimétrica, IRAM, Avenida Divina Pastora 7, Local 20, E-18012, Granada, Spain Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland Department of Astronomy, Yonsei University, Yonsei-ro 50, Seodaemun-gu, 03722 Seoul, Republic of Korea Netherlands Organisation for Scientific Research (NWO), Postbus 93138, 2509 AC Den Haag, The Netherlands Frontier Research Institute for Interdisciplinary Sciences, Tohoku University, Sendai 980-8578, Japan
Astronomical Institute, Tohoku University, Sendai 980-8578, Japan
Department of Physics and Astronomy, Seoul National University, Gwanak-gu, Seoul 08826, Republic of Korea
Leiden Observatory, Leiden University, Postbus 2300, 9513 RA Leiden, The Netherlands
Physics Department, Brandeis University, 415 South Street, Waltham, MA 02453, USA
Jeremiah Horrocks Institute, University of Central Lancashire, Preston PR1 2HE, UK
School of Physics, Huazhong University of Science and Technology, Wuhan, Hubei, 430074, People’s Republic of China
Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK
School of Astronomy and Space Sciences, University of Chinese Academy of Sciences, No. 19A Yuquan Road, Beijing 100049,People’s Republic of China
Astronomy Department, University of Science and Technology of China, Hefei 230026, People’s Republic of China
AbstractThe Event Horizon Telescope (EHT) has recently delivered the first resolved images of M87 ∗ , thesupermassive black hole in the center of the M87 galaxy. These images were produced using 230 GHzobservations performed in 2017 April. Additional observations are required to investigate the persis-tence of the primary image feature – a ring with azimuthal brightness asymmetry – and to quantifythe image variability on event horizon scales. To address this need, we analyze M87 ∗ data collectedwith prototype EHT arrays in 2009, 2011, 2012, and 2013. While these observations do not containenough information to produce images, they are sufficient to constrain simple geometric models. Wedevelop a modeling approach based on the framework utilized for the 2017 EHT data analysis andvalidate our procedures using synthetic data. Applying the same approach to the observational datasets, we find the M87 ∗ morphology in 2009–2017 to be consistent with a persistent asymmetric ringof ∼ µ as diameter. The position angle of the peak intensity varies in time. In particular, we finda significant difference between the position angle measured in 2013 and 2017. These variations arein broad agreement with predictions of a subset of general relativistic magnetohydrodynamic simula-tions. We show that quantifying the variability across multiple observational epochs has the potentialto constrain the physical properties of the source, such as the accretion state or the black hole spin. Keywords: black holes – accretion, accretion disks – galaxies: active – galaxies: individual: M87 –Galaxy: center – techniques: interferometric INTRODUCTIONThe compact radio source in the center of the M87galaxy, hereafter M87 ∗ , has been observed at 1.3 mil-limeter wavelength (230 GHz frequency) using very longbaseline interferometry (VLBI) since 2009. These obser-vations, performed by early configurations of the EventHorizon Telescope (EHT, Doeleman et al. 2009) ar-ray, measured the size of the compact emission to be ∼ µ as, with large systematic uncertainties relatedto the limited baseline coverage (Doeleman et al. 2012;Akiyama et al. 2015). The addition of new sites andsensitivity improvements leading up to the April 2017observations yielded the first resolved images of thesource (Event Horizon Telescope Collaboration et al.2019a,b,c,d,e,f, hereafter EHTC I-VI). These images re-vealed an asymmetric ring (a crescent) with a diameter d = 42 ± µ as and a position angle of the bright side φ B between ◦ and ◦ east of north (counterclockwisefrom north/up as seen on the sky, EHTC VI), see theleft panel of Figure 1. The apparent size and appearanceof the observed ring agree with theoretical expectationsfor a . × M (cid:12) black hole driving a magnetized ac-cretion inflow/outflow system, inefficiently radiating viasynchrotron emission (Yuan & Narayan 2014, EHTC V).Trajectories of the emitted photons are subject to strongdeflection in the vicinity of the event horizon, resultingin a lensed ring-like feature seen by a distant observer – the anticipated shadow of a black hole (Bardeen 1973;Luminet 1979; Falcke et al. 2000; Broderick & Loeb2009).General relativistic magnetohydrodynamic(GRMHD) simulations of relativistic plasma in theaccretion flow and jet-launching region close to theblack hole (EHTC V, Porth et al. 2019) predict that theM87 ∗ source structure will exhibit a prominent asym-metric ring throughout multiple years of observations,with a mean diameter d primarily determined by theblack hole mass-to-distance ratio and a position angle φ B primarily determined by the orientation of the blackhole spin axis. The detailed appearance of M87 ∗ mayalso be influenced by many poorly constrained effects,such as the black hole spin magnitude, magnetic fieldstructure in the accretion flow (Narayan et al. 2012,EHTC V), the electron heating mechanism (e.g., Moś-cibrodzka et al. 2016; Chael et al. 2018a), nonthermalelectrons (e.g., Davelaar et al. 2019), and misalignmentbetween the jet and the black hole spin (White et al.2020; Chatterjee et al. 2020). Moreover, turbulence inthe accretion flow, perhaps driven by the magnetorota-tional instability (Balbus & Hawley 1991), is expectedto cause stochastic variability in the image with correla-tion timescales of up to a few weeks ( ∼ dynamical timefor M87 ∗ ). The model uncertainties and expected time-dependent variability of these theoretical predictions onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT GRMHD φ jet φ spin φ B , exp Observation φ jet φ B , exp µ as Figure 1.
Left panel: one of the images of M87 ∗ obtainedin EHTC IV (see Section 4.2 for details). A 42 µ as circle isplotted with a dashed line for reference. The observed po-sition angle of the approaching jet φ jet is 288 ◦ east of north(Walker et al. 2018). Under the assumed physical interpre-tation of the ring, we expect to find the bright side of thecrescent on average approximately 90 ◦ clockwise from φ jet (EHTC V). We assume a convention φ B , exp = 198 ◦ , indi-cated with a blue dashed line. Right panel: a random snap-shot (note that this is not a fit to the EHT image) froma GRMHD simulation adopting the expected properties ofM87 ∗ (Section 4.1). The spin vector of the black hole ispartially directed into the page, counteraligned with the ap-proaching jet (and aligned with the deboosted receding jet);its projection onto the observer’s screen is located at theposition angle of φ spin = φ jet − ◦ . strongly motivate the need for additional observationsof M87 ∗ , especially on timescales long enough to yielduncorrelated snapshots of the turbulent flow.To this end, we analyze archival EHT observations ofM87 ∗ from observing campaigns in 2009, 2011, 2012,and 2013. While these observations do not have enoughbaseline coverage to form images (EHTC IV), they aresufficient to constrain simple geometrical models, follow-ing procedures similar to those presented in EHTC VI.We employ asymmetric ring models that are motivatedby both results obtained with the mature 2017 array andthe expectation from GRMHD simulations that the ringfeature is persistent.We begin, in Section 2, by summarizing the details ofthese archival observations with the “proto-EHT” arrays.In Section 3, we describe our procedure for fitting simplegeometrical models to these observations. In Section 4,we test this procedure using synthetic proto-EHT obser-vations of GRMHD snapshots and of the EHT imagesof M87 ∗ . We then use the same procedure to fit modelsto the archival observations of M87 ∗ in Section 5. Wediscuss the implications of these results for our theo-retical understanding of M87 ∗ in Section 6, and brieflysummarize our findings in Section 7. OBSERVATIONS AND DATA Our analysis covers five separate 1.3 mm VLBI observ-ing campaigns conducted in 2009, 2011, 2012, 2013, and2017. The M87 ∗ data from 2011 and 2013 have not beenpublished previously. For all campaigns except 2012,M87 ∗ was observed on multiple nights. For the proto-EHT data sets (2009–2013) we simultaneously utilize theentire data set from each year, fitting to data from mul-tiple days with a single source model, when available.This is motivated by the M87 ∗ dynamical timescale ar-gument, little visibility amplitude variation reported byEHTC III on a one-week timescale, as well as by the lim-ited amount of available data and lack of evidence forinterday variability in the proto-EHT data sets. We useincoherent averaging to estimate visibility amplitudes oneach scan ( ∼ few minutes of continuous observation) andbispectral averaging to estimate closure phases (Rogerset al. 1995; Johnson et al. 2015; Fish et al. 2016). Thefrequency setup in 2009–2013 consisted of two 480 MHzbands, centered at 229.089 and 229.601 GHz. Wheneverboth bands or both parallel-hand polarization compo-nents were available, we incoherently averaged all simul-taneous visibility amplitudes. The data sets are sum-marized in Table 1, where the number of detections fornonredundant baselines of different projected baselinelengths is given, with the corresponding ( u, v ) -coverageshown in Figure 2. Redundant baselines yield indepen-dent observations of the same visibility. In Table 1 wealso indicate the number of available nonredundant clo-sure phases (CPs, not counting redundant and intrasitebaselines, minimal set, see Blackburn et al. 2020). Asis the case for non-phase-referenced VLBI observations(Thompson et al. 2017), we do not have access to abso-lute visibility phases. All visibility amplitudes observedin 2009–2013 are presented in Figure 3.A more detailed summary of the observational setupof the proto-EHT array in 2009–2013 and the associateddata reduction procedures can be found in Fish et al.(2016). All data sets discussed in this paper are publiclyavailable . 2.1. Prior to 2013, the proto-EHT array included tele-scopes at three geographical locations: (1) the Com-bined Array for Research in Millimeter-wave Astronomy(CARMA, CA) in Cedar Flat, California, (2) the Sub-millimeter Telescope (SMT, AZ) on Mt. Graham in Ari-zona, and (3) the Submillimeter Array (SMA, SM), theJames Clerk Maxwell Telescope (JCMT, JC), and theCaltech Submillimeter Observatory (CSO, CS) on Mau-nakea in Hawai’i. These arrays were strongly east-westoriented, and the longest projected baselines, betweenSMT and Hawai’i, reached about 3.5 G λ , correspondingto the instrument resolution (maximum fringe spacing)of ∼ µ as. https://eventhorizontelescope.org/for-astronomers/data Table 1.
M87 ∗ data sets analyzed in this paper.Detections on Nonredundant BaselinesYear Telescopes Dates Baselines a Zero Short Medium b Long c Total CPs < . G λ < G λ < . G λ > . G λ > . G λ d d a theoretically available / with detections / nonredundant, nonzero with detections, b all / SMT-Hawai’i, c all / SMT-Chile, d single-day data set The 2011 observations of M87 ∗ have not been pub-lished but follow the data reduction procedures de-scribed in Lu et al. (2013). The 2009 and 2012 observa-tions and data processing of M87 ∗ have been publishedin Doeleman et al. (2012) and Akiyama et al. (2015),respectively. However, our analysis uses a modified pro-cessing of the 2012 data because the original processingerroneously applied the same correction for atmosphericopacity at the SMT twice. The SMT calibration pro-cedures have been updated since then (Issaoun et al.2017).Each observation included multiple subarrays ofCARMA as well as simultaneous measurements ofthe total source flux density with CARMA acting asa connected-element interferometer; these propertiesthen allow the CARMA amplitude gains to be “net-work calibrated” (Fish et al. 2011; Johnson et al. 2015,EHTC III). Of these three observing campaigns, only2012 provides closure phase information for M87 ∗ , andall closure phases measured on the single, narrow tri-angle SMT–SMA–CARMA were consistent with zero towithin σ (Akiyama et al. 2015), see Figure 4.2.2. The 2013 observing epoch did not include the CSO,but added the Atacama Pathfinder Experiment fa-cility (APEX, AP) in the Atacama Desert in Chile.This additional site brought for the first time thelong ( ≈ − λ ) baselines CARMA–APEX and SMT–APEX, that are roughly orthogonal to the CARMA–Hawai’i and SMT–Hawai’i baselines, see Figure 2. Theaddition of APEX increased the instrument resolution(maximum fringe spacing) to ∼ µ as. While the 2013 An opacity correction raises visibility amplitudes on SMT base-lines by ∼
10% in nominal conditions; our visibility amplitudeson SMT baselines are, thus, slightly lower than those reportedby Akiyama et al. (2015). However, the calibration error doesnot change the primary conclusions of Akiyama et al. (2015). observations of Sgr A ∗ were presented in several publi-cations (Johnson et al. 2015; Fish et al. 2016; Lu et al.2018), the M87 ∗ observations obtained during the 2013campaign have not been published previously.The proto-EHT array observed M87 ∗ on March 21st,22nd, 23rd, and 26th 2013. CARMA–APEX detectionswere found on March 22nd (11 detections) and 23rd(7 detections) with a single SMT–APEX detection onMarch 23rd. March 23rd (MJD 36374) was the onlyday with detections on baselines to each of the fourgeographical sites. No detections between Hawai’i andAPEX were found, and there were no simultaneous de-tections over a closed triangle that would allow for themeasurement of closure phase.2.3. In 2017, the EHT observed M87 ∗ with five geo-graphical sites (EHTC I; EHTC II), without CSO andCARMA, but with the addition of the Large Millime-ter Telescope Alfonso Serrano (LMT, LM) on the Vol-cán Sierra Negra in Mexico, the IRAM 30-m telescope(PV) on Pico Veleta in Spain, and the phased-up At-acama Large Millimeter/submillimeter Array (ALMA,AA, Matthews et al. 2018; Goddi et al. 2019). The ex-pansion of the array resulted in significant improvementsin ( u, v ) -coverage, shown with gray lines in Figure 2, andinstrument resolution raised to ∼ µ as. In addition tohardware setup developments (EHTC II), the recordedbandwidth was increased from 2 × × ∗ image by con- onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT − u (G λ ) − v ( G λ ) − u (G λ ) − − u (G λ ) − − u (G λ ) − intrasiteCARMA-SMT CARMA-Hawai’iSMT-Hawai’i APEX-CARMAAPEX-SMT EHT 2017 Figure 2. ( u, v ) -coverage of the M87 ∗ observations performed in 2009–2013 with various proto-EHT arrays. Gray lines indicatedetections obtained during the 2017 observations with a mature EHT array, including several new sites, but without the baselinesto CARMA. Dashed circles correspond to angular scales of 50 µ as (inner) and 25 µ as (outer). straining the set of physical (EHTC V) and geometric(EHTC VI) models representing the source morphology.2.4. M87 ∗ data properties VLBI observations sample the Fourier transform ofthe intensity distribution on the sky I ( x, y ) via thevan Cittert–Zernike theorem (van Cittert 1934; Zernike1938) V ( u, v ) = (cid:90) (cid:90) I ( x, y ) e − π i ( xu + vy ) d x d y , (1)where the measured Fourier coefficients V ( u, v ) are re-ferred to as “visibilities” (Thompson et al. 2017). Whenan array of N telescopes observes a source, N ( N − / independent visibility measurements are obtained, pro-vided detections on all baselines are found. Certainproperties of the geometry described by I ( x, y ) can beinferred directly from inspecting the visibility data.In the top panel of Figure 3 we summarize allthe M87 ∗ detections obtained during 2009–2013 ob-servations as a function of projected baseline length √ u + v . Dashed lines represent (cid:98) R , the analyticFourier transform of an infinitely thin ring with a to-tal intensity I and a diameter d , (cid:98) R (cid:16)(cid:112) u + v (cid:17) = I J (cid:16) πd (cid:112) u + v (cid:17) , (2)where J is a zeroth-order Bessel function of the firstkind. This simple analytic model predicts the visibilitynull located at b ≈ . (cid:18) d µ as (cid:19) − G λ , (3)and a wide plateau around the first maximum, locatedat b ≈ . (cid:18) d µ as (cid:19) − G λ , (4) recovering about 40% of the flux density seen on shortbaselines. In Figure 3 we use d = 45 . µ as and show (cid:98) R curves corresponding to I = λ ) in all data sets and the flux densityrecovery on long baselines to APEX in 2013, are roughlyconsistent with that of a simple ring model. Moreover,all detections on baselines to APEX have a similar fluxdensity of ∼ . Jy, while the projected baseline lengthvaries between 5.2 and 6.1 G λ . In the analytic thinring model framework, this can be readily understood,because baselines to APEX sample the wide plateauaround the maximum located at b , Equation 4.The gray dots in Figure 3 correspond to the sourcemodel constructed based on the 2017 EHT observations– the mean of the four images reconstructed for April5th, 6th, 10th and 11th 2017 with the eht-imaging pipeline (Chael et al. 2016, EHTC IV). In the 2017model, east–west baselines, such as SMT–Hawai’i, probea deep visibility null located around b (Equation 3),where sampled amplitudes drop below 0.01 Jy. North–south baselines do not show a similar feature, which indi-cates source asymmetry. Irrespective of the orientation,visibility amplitudes flatten out around b . Gray dotswith black envelopes represent the 2017 source modelsampled at the ( u, v ) -coordinates of the past observa-tions, for which all medium-length baselines were ori-ented in the east–west direction.One can immediately notice interesting discrepancies.The visibility amplitudes measured on long baselines toAPEX (projected baseline length ∼ b ) in 2013 wereabout a factor of 2 larger than the corresponding 2017source model predictions. At the same time, the fluxdensity on the short CARMA–SMT baseline is con-sistent between the 2013 measurements and the 2017model predictions. This shows that the image on the Baseline (G λ ) − − V i s i b ili t y A m p li t ud e ( J y ) L o n g b a s e li n e s M e d i u m b a s e li n e s Sh o rt b a s e li n e s I n tr a s i t e b a s e li n e s b = . G λ b = . G λ Modified Julian Date . . . . V i s i b ili t y A m p li t ud e ( J y ) I n tr a s i t e Sh o rt Figure 3.
Top:
Visibility amplitudes of M87 ∗ detections in 2009–2013 as a function of projected baseline length √ u + v . Thesource model derived from the EHT 2017 observations is shown with gray dots. Gray dots with black borders show the predictedvisibility amplitudes of the source model at the baselines of the prior observations in 2009–2013. Dashed black lines correspondto the family of Fourier transforms of a symmetric, infinitely thin ring of diameter d = 45 . µ as. Bottom: total arcsecond-scale flux density (on intrasite baselines, network-calibrated), and compact emission flux density from the short CARMA–SMTbaseline. In the case of the short baselines in 2017, predictions of the source model are given. onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT ( u, v ) -distances of 3.2-3.5 G λ on Hawai’i–USA baselines,record flux density above 0.1 Jy. At the same time, the2017 model predicts that these baselines sample a vis-ibility null region around b , with a flux density lowerby an order of magnitude. However, the compact flux(on short baselines) did not change by more than a fac-tor of two, remaining between 0.5 and 1.0 Jy throughoutthe 2009–2017 observations, see the bottom panel of Fig-ure 3. This suggests that the null location in the past (ifpresent) was different than that observed in 2017, whichmay correspond to a fluctuation of the crescent positionangle or a changing degree of source symmetry.Apart from the visibility amplitude data, a limitednumber of closure phases from the narrow triangle SMT–SMA–CARMA has been obtained from the 2012 dataset (Akiyama et al. 2015). All of these closure phasesare measured to be consistent with zero, which suggestsa high degree of east–west symmetry in the geometry ofthe source observed in 2012. While the closure phaseson this triangle were not observed in 2017 (CARMAwas not part of the 2017 array), we can numerically re-sample the 2017 images ( eht-imaging reconstructions,EHTC IV) to verify the consistency. In Figure 4 we showthe closure phases obtained in 2012, averaged betweenbands, the two CARMA subarrays shown separately.Near-zero closure phases observed in 2012 are roughlyconsistent with at least some models from 2017. Unfor-tunately, technical difficulties that occurred during the2012 campaign precluded obtaining measurements be-tween UTC 7.5 and 10.5, where nonzero closure phasesare predicted by all 2017 models.Altogether, we see strong suggestions that the 2009–2013 data sets describe a similar geometry to the 2017results, but there are also substantial hints that the de-tailed properties of the source structure evolved betweenobservations. These differences can be quantified withgeometric modeling of the source morphology. MODELING APPROACHThe sparse nature of the pre-2017 data sets precludesreconstructing images in the manner employed for the2017 data (EHTC IV). However, the earlier data arestill capable of providing interesting constraints on morestrictly parameterized classes of models. Figure 5 showsthe 2013 data set overplotted with a best-fit ring model (in blue; 5 degrees of freedom) and asymmetric Gaus-sian model (in red; 4 degrees of freedom). Both modelsattain similar fit qualities, as determined by Bayesianand Akaike information criteria (see, e.g., Liddle 2007). This is the maximum likelihood estimator for the slashed thickring model (RT), as discussed in Section 3.1. − − C l o s u r e P h a s e ( d e g ) model 17/04/05model 17/04/06model 17/04/10model 17/04/11 model RTmodel RGdata 2012 C1data 2012 C2 Figure 4.
Consistency of the closure phases on the SMT–SMA–CARMA triangle between the values observed in 2012(Akiyama et al. 2015) and numerically resampled sourcemodels constructed based on the 2017 observations. Thepredictions of the asymmetric ring models RT and RG fittedto 2012 observations are also given, see Section 3 and Sec-tion 5. Data corresponding to two CARMA subarrays, C1and C2, are shown separately.
In the absence of prior information, we would be unableto confidently select a preferred model. However, therobust image morphology reconstructed from the 2017data provides a natural and strong prior for selectingan appropriate parameterization. The “generalized cres-cent” (GC) geometric models considered in EHTC VIyielded fit qualities comparable to those of image recon-structions for the 2017 data, and in this paper, we applytwo variants of the GC model to the pre-2017 data sets.Owing to data sparsity, we restrict the parameter spaceof the models to a subset of that considered in EHTC VIcontaining only a handful of parameters of interest.Throughout this paper, we use perceptually uniformcolor maps from the ehtplot library to display the im-ages. In some of the figures we present models blurred toa resolution of 15 µ as, adopted in this paper as the effec-tive resolution of the EHT. The EHT instrument reso-lution measured as the maximum fringe spacing in 2017is about 25 µ as, however, for the image reconstructionmethods employed in EHTC IV a moderate effect of su-perresolution can be expected (Honma et al. 2014; Chaelet al. 2016). The effect may be much more prominentfor the geometric models, which are not fundamentallylimited by the resolution. https://github.com/liamedeiros/ehtplot λ )10 − V i s i b ili t y A m p li t ud e ( J y ) Gaussianringlarge ring µ as Figure 5.
Comparison of maximum likelihood (ML) asym-metric Gaussian and asymmetric ring models fitted to the2013 data. Data are shown as points with error bars corre-sponding to the thermal errors. The shaded regions coverall amplitudes for a given model. The red and blue linesrepresent models evaluated at the ( u, v ) -coordinates of theobservations. Both models offer a very similar fit quality.ML estimators are shown as inset figures. The model ofa ring with roughly double the diameter (dashed curve) fitsthe intermediate baselines, but is excluded by long-baselineamplitudes. Model specification
Given the image morphology inferred from the 2017data, the primary parameters of interest we would liketo constrain using the earlier data sets are the size ofthe source, the orientation of any asymmetry, and thepresence or absence of a central flux depression. Theanalyses presented in this paper use two simple ring-like models – similar to those presented in Kamruddin& Dexter (2013) and Benkevitch et al. (2016) – both ofwhich are subsets of the GC models from EHTC VI.The first model we consider is a concentric “slashed”ring, where the ring intensity is modulated by a lineargradient, hereafter denoted as RT. In this model, theflux is contained within a circular annulus with the innerand outer radii R in and R out , respectively. The modelis described by five parameters:1. the mean diameter of the ring d = R in + R out ,2. the position angle of the bright side of the ring ≤ φ B < π ,3. the fractional thickness of the ring < ψ = 1 − R in /R out < ,4. the total intensity < I < Jy, and5. β , an intensity gradient (“slash”) across the ringin the direction given by φ B , corresponding to the ratio between the dimmest and brightest points onthe ring, < β < . A ring of uniform brightnesshas β = 1 , while a ring with vanishing flux at thedimmest part has β = 0 .This model reduces to a slashed disk for ψ → . Theassumed definition of mean diameter is consistent withthe one used in EHTC VI, allowing for direct compar-isons. Except where otherwise specified, we use this firstmodel for the analyses discussed in this paper.As a check against model-specific biases, we considera second model consisting of an infinitesimally thinslashed ring, blurred with a Gaussian kernel (EHTC IV).The equivalent five parameters for this model are:1. the mean diameter of the ring d = 2 R in = 2 R out ,2. the position angle of the bright side of the ring ≤ φ B < π ,3. the width of the Gaussian blurring kernel < σ < µ as,4. the total intensity < I < Jy, and5. the slash < β < .This second model, hereafter referred to as RG, reducesto a circular Gaussian for d (cid:28) σ .Both the RT and RG models provide a measure of thesource diameter ( d ), the orientation of the brightnessasymmetry ( φ B ), and the presence of a central flux de-pression. We quantify the latter property using the fol-lowing general measure of relative ring thickness (fromEHTC VI) f w = R out − R in + 2 σ √ d , (5)where R out = R in for the RG model and σ = 0 for theRT model.All data sets except 2009 contain observations fromintrasite baselines (“zero baselines”), see Table 1. ForM87 ∗ , these baselines are sensitive to the flux from theextended jet emission on ∼ arcsecond scales (EHTC IV,see also the bottom panel of Figure 3) and do not di-rectly inform us about the compact source structure onscales of ∼ tens of microarcseconds. However, the intr-asite baselines still provide useful constraints on stationgain parameters (see Section 3.2), and hence, we do notflag them. Rather, we parameterize this large-scale fluxusing a large symmetric Gaussian component consist-ing of two parameters, flux and size. This componentis entirely resolved out on intersite baselines and thushas no direct impact on the compact source geometry.Ultimately, the models that we use have 5 geometricparameters for the 2009 data set and 5+2=7 geometricparameters in all other cases. onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT Fitting procedure and priors
We perform the parameter estimation for this pa-per using
Themis , an analysis framework developed byBroderick et al. (2020) for the specific requirements ofEHT data analysis.
Themis operates within a Bayesianformalism, employing a differential evolution Markovchain Monte Carlo (MCMC) sampler to explore the pos-terior space. Prior to model fitting, the data productsare prepared in a manner similar to that described inEHTC VI. Descriptions of the likelihood constructionsfor different classes of data products are given in Brod-erick et al. (2020).One important difference between the 2017 and pre-2017 data sets is that the latter contain almost exclu-sively visibility amplitude information, rather than hav-ing access to the robust closure quantities in both phaseand amplitude (Thompson et al. 2017; Blackburn et al.2020) that aided interpretation of the 2017 data. In ad-dition to thermal noise, visibility amplitudes suffer fromuncertainties in the absolute flux calibration, includingpotential systematic effects such as losses related to tele-scope pointing imperfections (EHTC III). These uncer-tainties are parameterized within
Themis using station-based amplitude gain factors g i , representing the scalingbetween the geometric model amplitudes | ¯ V ij | and thegains-corrected model amplitudes | ˆ V ij | , | ˆ V ij | = (1 + g i )(1 + g j ) | ¯ V ij | . (6)Model amplitudes | ˆ V ij | are then compared with the mea-sured visibility amplitudes | V ij | . Within Themis , thenumber of amplitude gain parameters N g is equal tothe number of (station, scan) pairs, i.e., the gains areassumed to be constant across a single scan but uncor-related from one scan to another. By explicitly mod-eling station-based gains, we correctly account for theotherwise covariant algebraic structure of the visibil-ity calibration errors (Blackburn et al. 2020). At eachMCMC step, Themis marginalizes over the gain am-plitude parameters (subject to Gaussian priors) usinga quadratic expansion of the log-likelihood around itsmaximum given the current parameter vector; see Brod-erick et al. (2020) for details. For the analysis of the2009–2013 data sets presented in this paper, we haveadopted rather conservative 15% amplitude gain un-certainties for each station, represented by symmetricGaussian priors with a mean value of 0.0 and standarddeviation of 0.15. The width of these priors reflects ourconfidence in the flux density calibration rather than thestatistical variation in the visibility data.The RT model is parameterized within
Themis interms of R out , φ B , ψ , I , and β . Uniform priors are usedfor each of these parameters, with ranges of [0 , µ asfor R out , [0 , π ] for φ B , [0 , for ψ , [0 , Jy for I , and [0 , for β . We achieve the “infinitesimally thin” ringof the RG model within Themis by imposing a strictprior on ψ of [10 − , − ] , and the prior on σ is uniform in the range [0 , µ as. Because d and f w are derivedparameters, we do not impose their priors directly butrather infer them from appropriate transformations ofthe priors on R out , ψ , and σ . The effective prior on d = R out (2 − ψ ) is given by π ( d ) = ln(2)200 , µ as ≤ d < µ as ln (cid:0) µ as d (cid:1) , µ as ≤ d < µ as , otherwise , (7)which is uniform within the range [0 , µ as. For theRT model, the effective prior on f w = ψ/ (2 − ψ ) is givenby π RT ( f w ) = f w ) , ≤ f w ≤ , otherwise , (8)which is not uniform but rather increases toward smallervalues. For the RG model, the effective prior on f w = σ (cid:112) /R out is given by π RG ( f w ) = α , ≤ f w ≤ α α f , f w > α , otherwise , (9)where α = (cid:112) / ≈ . for our specified priorson σ and R out ; this prior is uniform within the range [0 , α ] . 3.3. Degeneracies and limitations
Modeling tests revealed the presence of a large-diameter secondary ring mode in the posterior distri-butions for the 2009–2012 data sets, corresponding tothe dashed green line in Figure 5. This mode is ex-cluded by the detections on long baselines (APEX base-lines in 2013, multiple baselines in 2017) and detectionson medium-length ( ∼ . G λ ) baselines (LMT–SMT in2017). Excising this secondary mode, as we do for theposteriors presented in Figure 6, effectively limits the di-ameter d to be less than ∼ µ as. In all cases, the priorrange is sufficient to capture the entire volume of theprimary posterior mode, corresponding to an emissionregion of radius ∼ µ as. We have verified numericallythat this procedure produces the same results as restrict-ing the priors on R out to [0 , µ as for the analysis ofthe 2009–2012 data sets.As a consequence of the Fourier symmetry of a realdomain input signal, we have V ∗ ( u, v ) = V ( − u, − v ) .Hence, visibility amplitude data alone cannot break thedegeneracy between the orientation of φ B and φ (cid:48) B = φ B + 180 ◦ , and effectively, we only constrain the axis ofthe crescent asymmetry. This is how the reported uncer-tainties should be interpreted. Having that in mind, forthe 2009, 2011, and 2013, consisting exclusively of the2visibility amplitude data, we choose the reported φ B us-ing the prior information about the position angle of thejet φ jet to select the φ B such that φ jet − ◦ < φ B < φ jet ,where φ jet = 288 ◦ (Walker et al. 2018). In other words,between the orientations φ B and φ (cid:48) B , we choose theone that is closer to the expected bright side position φ B , exp = 198 ◦ . This is motivated by the theoretical in-terpretation of the asymmetric ring feature (EHTC V).In the case of the 2012 data set, for which a very lim-ited number of closure phases is available, we report theorientation φ B of the maximum likelihood (ML) esti-mators, noting the bimodal character of the posteriordistributions and the aforementioned ◦ degeneracy.These caveats do not apply to the 2017 data set, forwhich substantial closure phase information is availableand breaks the degeneracy.It is important to recognize that the parameters ofa geometric model have no direct relation to the phys-ical parameters of the source, unlike direct fitting us-ing GRMHD simulation snapshots (Dexter et al. 2010;Kim et al. 2016; Fromm et al. 2019, EHTC V) or ray-traced geometric source models (Broderick & Loeb 2009;Broderick et al. 2016; Vincent et al. 2020) to the data.The crescent model is a phenomenological descriptionof the source morphology in the observer’s plane. Ifphysical parameters (such as black hole mass) are to beextracted, additional calibration, in general affected bythe details of the assumed theoretical model and the ( u, v ) -coverage, needs to be performed (EHTC VI). MODELING SYNTHETIC DATAIn order to verify whether the ( u, v ) -coverage andsignal-to-noise ratio of the 2009–2013 observations aresufficient to constrain source geometric properties withsimple asymmetric ring models, we have designed testsusing synthetic VLBI observations. The synthetic obser-vations are generated with the eht-imaging software(Chael et al. 2016, 2018b) by sampling four emissionmodels (GRMHD1, GRMHD2, MODEL1, MODEL2)with the ( u, v ) -coverage and thermal error budget re-ported for past observations. Additionally, corruptionfrom time-dependent station-based gain errors has beenfolded into the synthetic observations. The ground-truth images that we use correspond to ray-traced snap-shots of a GRMHD simulation and published images ofM87 ∗ (EHTC IV), reconstructed based on the 2017 ob-servations. 4.1. GRMHD snapshots
For the first two synthetic data tests, we use a randomsnapshot from a GRMHD simulation of a low magneticflux standard and normal evolution (SANE) accretiondisk (Narayan et al. 2012; Sądowski et al. 2013) arounda black hole with spin a ∗ ≡ Jc/GM = 0 . , shown inFigure 1 (second panel), and in Figure 7 (first panel).The GRMHD simulation was performed with the iharm code (Gammie et al. 2003), and the ray tracing was done with ipole (Mościbrodzka & Gammie 2018). FollowingMościbrodzka et al. (2016) and EHTC V, we assumea thermal electron energy distribution function and re-late the local ratio of ion ( T i ) and electron temperature( T e ) to the plasma parameter β p representing the ratioof gas to magnetic pressure T i T e = R high β β + R low
11 + β , (10)with R high = 40 and R low = 1 for the considered snap-shot. The prescription given by Equation 10 parame-terizes complex plasma microphysics, allowing us to ef-ficiently survey different models of electron heating, re-sulting in a different geometry of the radiating region.As an example, for the SANE models with large R high the emission originates predominantly in the stronglymagnetized jet base region, while for a small R high diskemission dominates (EHTC V).The image considered here is a higher resolution ver-sion (1280 × . × M (cid:12) black hole at a distance D = 16.9 Mpc. This choice results in an M/D ratio of 3.80 µ as and an observed black hole shadow that isnearly circular with an angular diameter not substan-tially different from the Schwarzschild case, which is √ M/D = 39 . µ as (Bardeen 1973). For reference,the dashed circles plotted in Figure 1 and Figure 7 havea diameter of 42.0 µ as. These parameters were chosen tobe consistent with the ones inferred from the EHT 2017observations (EHTC I). The camera is oriented with aninclination angle of ◦ . The viewing angle was chosento agree with the expected inclination of the M87 ∗ jet(Walker et al. 2018). The choices of spin a ∗ , electrontemperature parameter R high , and the SANE accretionstate are arbitrary. The choice of R low follows the as-sumptions made in EHTC V. We also assume that theaccretion disk plane is perpendicular to the black holejet (the disk is not tilted). The first image, GRMHD1,has been rotated in such a way that the projection ofthe simulated black hole spin axis counteraligns withthe observed position angle of the approaching M87 ∗ jet, φ jet = 288 ◦ (Walker et al. 2018; Kim et al. 2018).The GRMHD2 test corresponds to the same snapshot,but rotated counterclockwise by 90 ◦ to φ jet = 18 ◦ ,hence displaying a brightness asymmetry in the east-west rather than in the north-south direction, see thesecond panel of Figure 7. Because the ( u, v ) -coverage in2009–2013 was highly anisotropic, a dependence of thefidelity of the results on the image orientation may beexpected. 4.2. M87 ∗ images Hereafter, we use the natural units in which G = c = 1 . onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT
35 40 45 50 550 . . . . P r o b a b ili t y d e n s i t y
100 200 3000 . . . GRMHD1 φ jet = 288 ◦
35 40 45 50 55d ( µ as)0 . . . . P r o b a b ili t y d e n s i t y
200 300 400 φ B (deg)0 . . . GRMHD2 φ jet = 18 ◦ . . . . P r o b a b ili t y d e n s i t y . . . . MODEL117/04/06
35 40 45 50 55d ( µ as)0 . . . . P r o b a b ili t y d e n s i t y
100 200 300 φ B (deg)0 . . . . MODEL217/04/11
Figure 6.
Top two rows: marginalized distributions of the mean diameter d and brightness maximum position angle φ B forthe RT model fits to GRMHD simulation snapshots GRMHD1 (first row) and GRMHD2 (second row). The 2017 posteriorsare contained within the dark gray bands. The dashed vertical line in the left panels denotes diameter of √ M/D . Thevertical dashed lines in the right panels denote the convention angle φ B , exp , φ B , exp − ◦ , and the approaching jet position angle φ jet = φ B , exp + 90 ◦ . The range of ( φ jet − ◦ , φ jet ) is highlighted. Two bottom rows: similar as above, but for the RT modelfitted to MODEL1 and MODEL2. Lightly shaded areas correspond to values reported in EHTC I, diameter d = 42 ± µ as andposition angle ◦ < φ B < ◦ . ∗ generated based on the 2017 EHT observations,published in EHTC IV. We consider two days of the 2017observations with good coverage and reported structuralsource differences (EHTC III; EHTC IV, Arras et al.2020) - 2017 April 6th (MODEL1) and 2017 April 11th(MODEL2), see the first row of Figure 7. MODEL2was also shown in the first panel of Figure 1. Whilethese models are constructed based on the observationaldata, we resample them numerically to obtain syntheticdata sets considered in this section. Synthetic closurephases on the SMT–SMA–CARMA triangle computedfrom these models were shown in Figure 4. Note thatthere is a subtle difference between resampling a modelconstructed based on the 2017 data with a numericalmodel of the 2017 array and direct modeling of the ac-tual 2017 data, considered in Section 5. The sampled im-ages were generated utilizing the eht-imaging pipelinethrough the procedure outlined in EHTC IV, with a res-olution of 64 ×
64 pixels. This test can be viewed as anattempt to evaluate what the outcome of the modelingefforts would have been had the 2017 EHT observationsbeen carried out with one of the proto-EHT 2009–2013arrays rather than with the mature 2017 array.4.3.
Results for the synthetic data sets
Figure 7 shows a summary of the maximum likeli-hood (ML) RT model fits to the synthetic data sets;each column shows the fits for a single ground-truthimage, and each row shows the fits for a single arrayconfiguration. Though our simple ring models cannotfully reproduce the properties of the abundant and highsignal-to-noise data sampled with the 2017 array (i.e.,fits to these data sets are characterized by poor reduced- χ values of χ n ∼ ), they nevertheless recover diame-ter and orientation values that are reasonably consistentwith those reported in EHTC VI for the 2017 observa-tions, including the counterclockwise shift of the bright-ness position angle between 2017 April 6th (MODEL1)and 11th (MODEL2). We note that the underfitting ofthe data set sampled with the 2017 array results in anartificial narrowing of the parameter posteriors (the fullposterior is captured in EHTC VI by considering a morecomplicated GC model). ML estimators for 2009–2013data sets, on the other hand, typically fit the data muchcloser than the thermal error budget. Because we modeltime-dependent station gains in a small array (often onlytwo to three telescopes observing at the same time in2009–2013 with missing detections on some baselines),the number of model parameters may be formally largerthan the number of data points, and this complicates ourestimation of the number of effective degrees of freedom(see also Section 5.3). As a consequence, we cannot gen- erally utilize a χ n goodness-of-fit statistic as was donein EHTC IV and EHTC VI. The relevant parameter estimates and uncertaintiesfrom the RT model fits are listed in Table 2. In Figure 6we show the marginalized posteriors for the diameterand position angle parameters. For each synthetic dataset we also indicate the image domain position angle η ,defined as η = Arg (cid:18) (cid:80) k I k e iφ k (cid:80) k I k (cid:19) , (11)where I k is the intensity and φ k is the position angleof the k th pixel in the image. A similar image domainposition angle estimator was considered in EHTC IV.We notice that the image- and model-based estimatorsmay occasionally display significant differences (e.g.,GRMHD1). However, they are both sensitive to globalproperties of the brightness distribution, unlike someother estimators that could be considered, such as, e.g.,the location of the brightest pixel. For the diameter d estimates reported in Table 2 we list both the me-dian and ML values, with 68% and 95% confidence in-tervals, respectively. For the orientation angle φ B welist the ML values with 68% confidence intervals, andfor the fractional thickness f w we list the 95th distribu-tion percentile. Values of φ B contained in parenthesesindicate that the 68% confidence interval exceeds 100 ◦ ,in which case we have concluded that the orientationis effectively unconstrained. We find that the diameteris well constrained in general, with the GRMHD datasets recovering a typical value of ∼ µ as and 95% con-fidence intervals that never exceed ± µ as from thisvalue; the analogous measurement for the MODEL datasets is ± µ as. Biases related to the array orientationcan be seen - particularly with the 2013 coverage, theGRMHD2 test estimates an appreciably larger diameterthan GRMHD1, inconsistent within the 68% confidenceinterval.The orientation φ B is poorly constrained, with pos-terior distributions that depend strongly on the detailsof the ( u, v ) -coverage. Nevertheless, the 2009–2013 MLestimates provide orientations of the axis of asymme-try that are consistent within ± ◦ with the results ob-tained using the 2017 synthetic coverage. We note thatin 3 out of 4 synthetic data sets, the limited number ofclosure phases provided by the simulated data sets withthe 2012 coverage is enough to correctly break the degen-eracy in the position angle φ B , discussed in Section 3.3.For the synthetic GRMHD data sets, the preference forthe correct brightness position angle is very strong (seeFigure 6). For the MODEL data sets, the effect of clo-sure phases is much less prominent, the distributionsremain bimodal, and in the case of the MODEL1 data See, e.g., Andrae et al. (2010) for further comments about theproblems with the χ n metric and counting the degrees of freedom. onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT GRMHD1 I N P U T µ as GRMHD2 MODEL1 MODEL2
Figure 7.
ML estimators corresponding to the fits to four synthetic images, shown in the first row (no blurring). Estimatorswere obtained through synthetic VLBI observations with the ( u, v ) -coverage and uncertainties identical to those of the realobservations performed in 2009–2017. The thick ring model (RT) was used, and the presented images of ML estimators wereblurred to µ as resolution. Blue dashed lines indicate the convention for the expected position angle of the bright component φ B , exp . The gray bar represents the ML estimate of φ B . For the 2009, 2011, and 2013 data sets, the orientation is determinedassuming that | φ B , exp − φ B | < ◦ . The dashed circles correspond to a diameter of 42 µ as. Table 2.
Parameter estimates from fitting the RT model tothe synthetic data sets.Coverage d ( µ as) φ B (deg) f w Estimator Median ML ML At MostConfidence 68% 95% 68% 95%GRMHD1 2009 . +5 . − . . +4 . − . +56 − . +5 . − . . +5 . − . +55 − η a = 170 ◦ . +4 . − . . +4 . − . +25 − . +2 . − . . +2 . − . (246) b . +0 . − . . +0 . − . +1 − . +3 . − . . +11 . − . (302) . +2 . − . . +2 . − . (265) η a = 260 ◦ . +3 . − . . +10 . − . +6 − . +1 . − . . +2 . − . +38 − b . +0 . − . . +0 . − . +1 − . +3 . − . . +4 . − . +10 − . +2 . − . . +2 . − . +65 − η a = 155 ◦ . +3 . − . . +9 . − . +38 − . +2 . − . . +1 . − . +36 − b . +0 . − . . +0 . − . +1 − . +3 . − . . +11 . − . +30 − . +2 . − . . +2 . − . +72 − η a = 172 ◦ . +3 . − . . +9 . − . (179) . +1 . − . . +1 . − . +2 − b . +0 . − . . +0 . − . +1 − a η calculated with Equation 11, b using the ( u, v ) -coverageof 2017 April 6th set, the ML estimator points at the wrong orientation,suggesting brightness located in the north.We also consider the fractional thickness f w of thering, as defined in Equation 5. The fractional thicknessprovides a measure of whether the data support the pres-ence of a central flux depression, a signature feature ofthe black hole shadow, or if it is consistent with a disk-like morphology (i.e., f w ≈ for the RT model). Wefind that f w is less well constrained than the diame-ter d , consistently with the conclusions of EHTC VI. Insome cases the ML estimator corresponds to a limit ofa disk-like source morphology without a central depres-sion (see Figure 7). Only the 2013 and 2017 syntheticdata sets allow us to confidently establish the presenceof a central flux depression, with posterior distributionsexcluding f w > . for all synthetic data sets (see Ta-ble 2). For the 2009–2012 coverage synthetic data sets, f w is not constrained sufficiently well to permit similarstatements. We find that the RG model produces re- sults that are typically consistent with those of the RTmodel (see Appendix A, Figure 16). MODELING REAL DATAEncouraged by the results of the tests on syntheticdata sets, we performed the same analysis on the 2009–2013 proto-EHT M87 ∗ observations. We also present theanalysis of the 2017 observations with the RT and RGmodels. In the latter case only lower band data (2 GHzbandwidth centered at 227 GHz) were used.5.1. Source geometry estimators
In the first row of Figure 8 we show the ML estima-tors obtained by fitting the RT model to each observa-tional data set; in Figure 9 we show the same for theRG model. For the 2009, 2011, and 2013 data sets,which only constrain the axis of the crescent asymmetry,the orientation of the brightness peak was selected witha prior derived from the approaching jet orientation onthe sky, φ B ∈ ( φ jet − ◦ , φ jet ) , Section 3.3. The 2012data set, for which some closure phases are available,indicates a weak preference toward the brightness posi-tion located in the north rather than in the south (Fig-ure 10, and Figures 17-18 in the Appendix B). However,the posterior distribution remains bimodal and 47% ofits volume remains consistent with the jet orientationbased prior. Hence, the distinction is not very signif-icant – it is entirely dependent on the sign of closurephases shown in Figure 4, which are all consistent withzero to within 2 σ . Closure phases predicted by the MLestimators for the RT and RG models are indicated inFigure 4. The 2017 data sets are consistent with the ori-entation imposed by the jet position prior. The secondrows of Figures 8 and 9 show the ML estimators blurredto a resolution of µ as. In the third rows of Figures8 and 9 we present “mean images” for each data set,obtained by sampling 2 × sets of model parametersfrom the MCMC chains and averaging the correspond-ing images. The mean images highlight structure thatis “typical” of a random draw from the posterior distri-bution, though we note that a mean image itself doesnot necessarily provide a good fit to the data. Becauseof the rotational degeneracy, the orientation is alwaysassumed to be the one closer to the orientation givenby the ML estimate of φ B for the construction of theseimages. 5.2. Estimated parameters
The marginalized posteriors for the mean diameter d and position angle φ B for the observational data setsare shown in Figure 10 for both the RT and RG models,and tabulated values of the relevant estimates for the RTmodel are given in Table 3. The posterior distributionsfor the 2009–2012 data sets have complex shapes, notall parameters are well constrained, and ML estimatorsdo not necessarily coincide with the marginalized poste-riors maxima of the individual model parameters, which onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT M L µ as M L + B L U R M E AN Figure 8.
First row:
ML estimators obtained from fitting the RT model to the 2009–2017 observations. The position angle φ B isindicated with a bar. For the 2009, 2011, and 2013 data sets the orientation is determined assuming that φ jet − ◦ < φ B < φ jet ,where φ jet = 288 ◦ . Position angle 68% confidence intervals are shown for the 2009–2013 data sets. Second row:
RT modelsfrom the first row blurred to a 15 µ as resolution, indicated with a beam circle in the bottom-right corner of the first panel. Thedashed circle of 42 µ as diameter is plotted for reference. Third row: mean of the 2 × images drawn from the posterior of theRT model fits. M L M L + B L U R M E AN Figure 9.
Same as Figure 8, but for the RG model. The 2009–2012 posterior distributions contain a Gaussian mode, manifestingas a bright ring interior in the mean images. This is related to the lower spatial resolution of the 2009–2012 observations. . . . P r o b a b ili t y d e n s i t y A p r . . . . A p r A p r RT
30 35 40 45 50d ( µ as)0 . . . P r o b a b ili t y d e n s i t y A p r
100 200 300 φ B (deg)0 . . . . A p r A p r RG Figure 10.
Top: marginalized distributions of mean diameter d and brightness maximum position angle φ B for the RT modelfits to 2009–2017 observational data sets. Lightly shaded areas correspond to values reported in EHTC I, diameter d = 42 ± µ asand position angle ◦ < φ B < ◦ . The 2017 posteriors are contained within the dark gray bands. The dashed vertical linein the left panels denotes the diameter of √ M/D . The vertical dashed lines in the right panels denote the convention angle φ B , exp = 198 ◦ , φ B , exp − ◦ , and the approaching jet position angle φ jet = 288 ◦ . The range of ( φ jet − ◦ , φ jet ) is highlighted. Bottom: same, but for the RG model. can be seen in the corner plots (Figures 17-18 in the Ap-pendix B). The behavior of the posterior distributionsis much improved for the 2013 data set, and becomes ex-emplary in the case of the 2017 data sets (Figures 19-20in the Appendix B).
Table 3.
Parameter estimates from fitting the RT model tothe observational data sets. d ( µ as) φ B (deg) f w Estimator Median ML ML At MostConfidence 68% 95% 68% 95%2009 . +5 . − . . +3 . − . +68 − . +2 . − . . +2 . − . (130) a . +2 . − . . +2 . − . +42 − . +2 . − . . +3 . − . +17 − b . +0 . − . . +0 . − . +1 − c . +0 . − . . +0 . − . +1 − a secondary mode present at φ B − ◦ (see the text), b c Similar to the case of the synthetic data sets, we findthat the diameter d is well constrained; the RT model95% confidence intervals across all observational datasets always fall within ± µ as from d = 40 µ as. The2013 proto-EHT observations provide meaningful con-straints on φ B , indicating that the source asymmetry in2013 was in the east-west direction, rather than in thenorth-south direction, as in the case of the 2017 dataset. The 2009 and 2011 data sets do not constrain theorientation well.All ML estimators and mean images from the RTmodel fits show a clear shadow feature, indicating thata disk-like, filled-in structure is disfavored by all ofthe observations (however, for 2009–2012 it cannot beexcluded with high confidence based on the relativethickness parameter f w distribution, Table 3 and Ap-pendix B). This is contrary to the synthetic data resultsshown in Figure 7, where some of the ML estimators cor-respond to a disk-like morphology. On the other hand,the mean images for the 2009–2012 RG model fits showa significant flux density interior to the ring, indicat-ing that these data sets are consistent with a symmetricGaussian source model, having no central flux depres-sion. This is a consequence of the resolution being lim-ited by the lack of long baselines prior to 2013. Short onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT f w > , becoming a Gaussianin the limit of σ (cid:29) d .The slash parameter β can be measured to be . ± . for the 2013 data set, which is consistent with the fits tothe 2017 data sets that give β ∼ . (RT) or β ∼ . (RG). Fits to the 2012 data set indicate preference to-ward more symmetric brightness distribution, 2009 and2011 do not provide meaningful constraints on β .5.3. Quality of the visibility amplitude fits
The quality of fits and their behavior in terms of χ n are similar to the synthetic data sets (see the commentsin Section 4.3 and Table 4). In Figures 11-12 we explic-itly give the number of independent visibility amplitudeobservations for each data set N ob and the number ofindependent visibilities on nonzero (intersite) baselines N nz . Note that the latter is larger than the number ofdetections on nonzero, nonredundant baselines given inTable 1, as some detections are independent but redun-dant. We also provide the number of explicitly modeledamplitude gains N g for each data set (see Section 3.2and Section 4.3). Given the pathologies in the χ n met-ric described in Section 4.3, we characterize the qualityof the ML estimator fits to data using the two followingmetrics ¯ e = 1 N nz N nz (cid:88) i =1 (cid:12)(cid:12)(cid:12)(cid:12) | V i | − | ¯ V i | σ i (cid:12)(cid:12)(cid:12)(cid:12) , (12) ˆ e = 1 N nz N nz (cid:88) i =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) | V i | − | ˆ V i | σ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (13)In Equations 12-13 we follow the notation of Equation 6,that is, ¯ V i represents the visibilities of the geometricmodel while ˆ V i corresponds to the model modified byapplying the estimated gains, representing the final fitto observations V i . Uncertainties σ i correspond to thethermal error budget. We only account for nonzero base-lines, which describe the compact source properties. Inthe bottom rows of Figures 11-12 we indicate two errorbars. Black error bars correspond to the thermal un-certainties σ i , while the red ones correspond to inflateduncertainties s i = (cid:113) σ i + 2 · (0 . V i ) , (14)approximately capturing the uncertainty related to theamplitude gains. For all 2009–2013 data sets, the flexi-bility of the full model is sufficient to fit the sparse datato within the thermal uncertainty level with a best-fitML estimator. 5.4. Consistency with the prior analysis
In order to assess model-related biases and verify towhat extent our simplified models recover geometric pa-rameters consistent with the ones reported by EHT, i.e.,image domain results given in EHTC IV Tables 5 and 7,and geometric modeling results given in EHTC VI Ta-bles 2 and 3, we gather these results in Table 4. For thedetails of the methods and algorithms, see explanationsand references in EHTC IV and EHTC VI. We noticethat (1) differences between methods may be as large as30 ◦ for the same data set, (2) models considered in thiswork measure diameter and position angle consistentlywith more complex crescent models (GC) and with theimage domain methods within the expected intermodelvariation, (3) RT and RG models are too simplistic tofully capture the properties of the 2017 data sets, result-ing in underfitting, indicated by higher values of χ n , and(4) posteriors of the RT and RG models are narrower forthe 2017 data sets than the GC posteriors as an effectof the underfitting. Table 4.
Comparison between parameters extraction resultsreported in This Paper, EHTC IV, and EHTC VI.Source Method d ( µ as) φ B or η (deg) a χ n . +0 . − . . ± . . RG . ± .
07 154 . +1 . − . . EHTC VI
Themis . ± .
14 153 . +2 . − . . (GC) dynesty . +0 . − . . +1 . − . . EHTC IV
DIFMAP . ± . . ± . eht-imaging . ± . . ± . SMILI . ± . . ± . . . ± .
05 184 . +0 . − . . RG . +0 . − . . ± . . EHTC VI
Themis . +0 . − . . +2 . − . . (GC) dynesty . +0 . − . . +2 . − . . EHTC IV
DIFMAP . ± . . ± . . (image eht-imaging . ± . . ± . . domain) SMILI . ± . . ± . . a results from this work and EHTC VI correspond to visibil-ity domain-based estimator φ B , while results from EHTC IVcorrespond to the image domain estimator η , similar to theone given by Equation 116. DISCUSSION AND CONCLUSIONSProcessing five independent observations of M87 ∗ ina unified framework offers important insights into the0 − A m p li t ud e ( J y ) N ob = 38 , N nz = 38 N g = 46 N ob = 232 , N nz = 120 N g = 238 N ob = 67 , N nz = 53 N g = 52 − . . . R e s i du a l ( J y ) ¯ e = 0 . e = 0 .
20 0 1 2 3 4Baseline (G λ )¯ e = 0 . e = 0 .
28 0 1 2 3 4¯ e = 1 . e = 0 . Figure 11.
Top row: visibility amplitudes measured in 2009–2012 are shown as black diamonds with error bars correspondingto thermal uncertainties. Blue shaded regions correspond to the range of visibility amplitudes of the asymmetric ring RTmodel ML estimators, shown in the first row of Figure 9. The total number of observed visibility amplitudes N ob is given,along with the number of nonzero baseline visibility amplitudes N nz , and the number of modeled gains N g . Bottom row: differences between measured amplitudes | V i | and the geometric model amplitudes | ¯ V i | . Black error bars correspond to thermaluncertainties, while red ones correspond to error budget inflated by adding systematics approximately capturing the gainsuncertainties (Equation 14). Two fit quality metrics, defined with Equations 12-13, are provided for each data set. − − A m p li t ud e ( J y ) N ob = 187 , N nz = 146 N g = 176 N ob = 274 , N nz = 242 N g = 126 N ob = 216 , N nz = 189 N g = 104 − . . . R e s i du a l ( J y ) ¯ e = 2 . e = 0 .
26 0 2 4 6 8Baseline (G λ )¯ e = 15 . e = 2 .
36 0 2 4 6 8¯ e = 9 . e = 1 . Figure 12.
Same as Figure 11, but for the 2013–2017 data sets. source morphology and stability. While the constrain-ing power of the 2009–2013 data sets is rather weak incomparison to the 2017 observations, we find evidencefor a persistent ring structure that shows modest struc- tural variability. Both the persistence and variabilityoffer important constraints for models of M87 ∗ . We willnow discuss our evidence and theoretical implications forthe presence of the shadow feature in 2009–2017 (Sec- onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT ◦ ◦ ◦ ◦ ◦ ◦ Figure 13.
Top row: a single GRMHD simulation snapshot, ray-traced with total compact flux density adjusted to valuesbetween 0.1 Jy and 1.1 Jy. The brightness distribution in each panel has been normalized by the brightness maximum. EHTobservations of M87 ∗ found a total compact flux density of 0.8-0.9 Jy in 2009–2012 and 0.5-0.6 Jy in 2013 and 2017, see thebottom panel of Figure 3. Bottom row: same snapshots blurred to the resolution of 15 µ as. The 42 µ as diameter ring is indicatedwith a dashed circle. The position angle convention φ B , exp is shown with a dashed blue bar. The gray bar indicates the imagedomain position angle η calculated using Equation 11, and values of η obtained are given in the bottom-right corner of eachpanel. tion 6.1), the persistence of the source geometry as anargument for its theoretical interpretation (Section 6.2),and the variability of the source geometry (Section 6.3).We also summarize the limitations and caveats of thetheoretical interpretation within the GRMHD frame-work (Section 6.4).6.1. Presence of the shadow feature
In Section 5 we have discussed the fits of the asym-metric ring models to the M87 ∗ observations, indicatingthat all data sets are consistent with such a geometry.Within the framework of a ring model, a key question forthe archival observations is whether we unambiguouslydetect the inner flux depression seen in the 2017 results,the expected feature of a black hole shadow. While max-imum likelihood estimators clearly indicate this featurein all cases (Figure 8), a detailed inspection of the pos-terior distributions shown in Appendix B allows us toconclude that only the 2013 archival data set providesa robust detection of the central flux depression, con-straining the relative thickness of the ring f w to less than0.5 (see Table 3). For the 2009–2012 data sets, the RTmodel posteriors indicate some preference toward thepresence of a flux depression; however, a disk-like filled-in geometry cannot be excluded with a high degree ofconfidence (Appendix B). This is not entirely surpris-ing, as the 2009–2012 observations lack baselines withprojected lengths > . G λ , probing spatial frequencieshigher than the visibility null b (Section 2.4), whichare more sensitive to differences between an empty ring,a filled-in disk, and a Gaussian. However, the presence of the shadow feature is sensi-tive to changes in the optical depth. The total compactflux density in 2009–2012 was measured to be . − . Jy,significantly higher than the . − . Jy observed in2013 and 2017 (Figure 3, bottom panel). Mildly ele-vated levels of X-ray emission from the nucleus of M87before 2016 were also reported (Sun et al. 2018). Thesemeasurements suggest a higher mass accretion rate in2009–2012 and hence a larger density scale, in turn in-creasing the opacities and the optical depth, possiblychanging the appearance of the black hole shadow (Moś-cibrodzka et al. 2012; Chan et al. 2015). Is it possiblethat the shadow feature had been obscured by a moreoptically thick medium in 2009–2012 than in the caseof the more recent 2013 and 2017 observations? Whilethe fully general answer to that question would requireextensive testing of a variety of GRMHD models, weaddress this concern by analyzing a representative ex-ample from the library of simulated M87 ∗ images. Weconsider a random snapshot from a SANE simulationwith spin a ∗ = 0 . , and electron temperature param-eter R high = 20 . The R low parameter is equal to 1,as it is throughout this paper. The snapshot is thenrepeatedly ray-traced with its density scale (and henceopacities and emissivities) adjusted to give a total com-pact flux density equal to 0.1, 0.3, 0.5, 0.7, 0.9, and1.1 Jy. The resulting images, normalized by the bright-ness maximum, are shown in Figure 13. These findingsindicate that variation in the total compact flux densitybetween 0.1 and 1.1 Jy does not eliminate the centralbrightness depression. Judging from the similarity be-tween blurred images, the flux density scaling is also2not expected to influence the estimated diameter or theposition angle appreciably. However, it is likely to influ-ence measures of asymmetry, such as the slash parame-ter β in the RT/RG models, discussed in Section 3.1. In-terestingly, we see a preference toward more symmetricsource geometry (larger β ) in the 2012 posteriors (Ap-pendix B). We conclude that the central flux depressionwas most likely present throughout the 2009–2017 obser-vations and that a lack of a high-confidence detection ofthis feature in 2009–2012 is presumably a consequenceof the very limited ( u, v ) -coverage.None of the EHT observations so far took place dur-ing an unambiguous flaring activity period of the M87nucleus, such as the events discussed in Abramowskiet al. (2012). We note that such an event could po-tentially influence the image morphology more stronglythan the moderate increase of total brightness consid-ered in Figure 13. Future simultaneous multiwavelengthobservational campaigns will shed light on the structuralchanges in the M87 ∗ compact radio image in relation tothe enhanced activity in different parts of the spectrum,allowing the site of particle acceleration to be localized.6.2. Black hole shadow or a transient feature
As discussed in Section 3, the 2009–2013 data setscan be successfully modeled with both an asymmetricGaussian and a ring model, with each giving a simi-lar fit quality. However, the maximum likelihood Gaus-sian models are very inconsistent in size, shape, andorientation across different years, see Figure 14. In con-trast, the best-fitting ring models, as seen in Figure 8,are similar in size over all epochs under the priors de-scribed in Sections 3.2 and 3.3. Even when consideredseparately from the 2017 results, this consistency sup-ports our choice to interpret the 2009–2013 morphologyas a ring and to draw conclusions from the results offitting asymmetric ring models to the data. The 95%confidence intervals of all 2009–2017 diameter posteri-ors lie within 40 ± µ as, while ML estimators of the di-ameter from all years lie within 43 ± µ as, see Table 3.These values are also consistent with the M87 ∗ diam-eter measurement of 42 ± µ as reported by EHTC VIand with the expected size of an observed black holeof mass/distance corresponding to the stellar dynam-ics measurement (Blakeslee et al. 2009; Gebhardt et al.2011). All 2009–2017 diameter measurements are incon-sistent with the gas dynamics mass measurement (Walshet al. 2013), which predicts a diameter roughly half aslarge. The consistency in the diameter across multipleobservational epochs supports the interpretation of thering-like feature as emission from the immediate sur-roundings of the supermassive black hole.The four EHT observations of M87 ∗ in 2017 spannedabout a week, corresponding to a timescale of ∼ M .With such a short time span, we cannot exclude a tran-sient origin for the source morphology using 2017 dataalone. However, such a feature would need to attain µ as Figure 14.
Year to year consistency of the best-fitting (ML)asymmetric Gaussian model to the M87 ∗ data sets. the geometry and apparent size expected of a shadowof a black hole (of independently measured mass-to-distance ratio) through an unusual coincidence. More-over, transient features are unlikely to persist with a sim-ilar geometry throughout multiple years of observations,corresponding to ∼ − M timescales (the totalspan of our analyzed observations is M ). For ex-ample, a lensed background source would need a lowtransverse velocity of v (cid:46)
40 km s − to travel (cid:46) M be-tween 2009 and 2017. This is much smaller than thegas velocities seen in the nucleus of M87 (e.g., Mac-chetto et al. 1997). A bright feature moving with thejet ( v ∼ c ) should travel (cid:38) pc on that timescale, a fac-tor of roughly larger than the physical size of thering itself. A bright knot or other stationary jet fea-ture would need to persist with a similar location, fluxdensity, and ring morphology to remain consistent withthese results. The 8 yr span of the 2009–2017 monitoringis also much longer than the typical variability timescaleof the M87 nucleus observed at 230 GHz, which is ∼ days (Bower et al. 2015). While features remaining sta-tionary for many years in otherwise rapidly flowing jetshave been reported and interpreted as standing recolli-mation shocks (Lister et al. 2009), such a configurationwould constitute one more unusual coincidence. Thus,we conclude that with multiple years of observations re-maining consistent with a ∼ µ as ring model, it ishighly unlikely that the origin of the observed geometrycould be a transient feature.6.3. Time variability of the source geometry
In addition to conclusions from the persistence of thering structure, we can also draw inferences from the vari-ability observed in the ring structure across the 2009–2017 data sets. In particular, the spread of the diam-eter and brightness position angle estimates (Table 3)are significantly larger than the spread for correspond-ing static synthetic data sets (Table 2). As a specificexample, the circular standard deviation of the ML po-sition angle estimators given in Table 2 is equal to 19 ◦ ,19 ◦ , 11 ◦ , 4 ◦ for GRMHD1, GRMHD2, MODEL1, andMODEL2, respectively. For the observational data (Ta-ble 3) the circular standard deviation is equal to 48 ◦ .This larger spread suggests that we are detecting intrin-sic structural variability despite the large uncertaintiesin the parameters estimated with pre-2017 observations.Moreover, unambiguous signatures of intrinsic variabil- onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT a ∗ , the electron temper-ature parameter R high (see Section 4.1) and the accre-tion state – strongly magnetized magnetically arresteddisk (MAD, Narayan et al. 2003), or standard and nor-mal evolution (SANE) flow, such as those considered inSections 4.1 and 6.1. Other parameters (e.g., total com-pact flux density of 0.5 Jy, inclination, jet position angle)are adjusted to match the observed properties of M87 ∗ .For our exploratory study, we utilize the following foursimulations: (S1) MAD, a ∗ = 0 . , R high = 10 ; (S2)SANE, a ∗ = 0 . , R high = 10 ; (S3) MAD, a ∗ = − . , R high = 10 ; (S4) MAD, a ∗ = 0 . , R high = 160 . For eachsimulation, we take 500 ray-traced snapshots with M separation in time. For each snapshot we calculate theposition angle of the bright component η using Equa-tion 11. We then construct histograms of η for each ofthe four simulations, shown in Figure 15.Each of the simulation parameters influences the dis-tribution of η , both in terms of its mean and spread.Some of these differences can be readily understood,for instance, in the case of prograde accretion ontoa spinning black hole, the radiation is boosted bothwith Doppler and with frame-dragging effects. The po-sition angle of the bright component is thus expectedto be relatively more influenced by the geometry (as-sumed to be fixed) and not by the stochastic componentthan the retrograde case, in which Doppler effect andframe-dragging counteract and the geometry becomesrelatively less important. Irrespective of the mecha-nism, some variants of simulations have great difficultiesexplaining either source orientation on 2017 April 6th(bright component too far clockwise), or in 2013 (toofar counterclockwise), as indicated in Figure 15. Thus,continued EHT observations, with tight constraints on η spaced over multiple years, will constrain these typesof models on the basis of variability in η .6.4. Limitations of the current approach
There are caveats to the simple analysis outlinedabove that should be addressed in more focused futurestudies. The simulations that we consider do not cap-ture the full extent of the physics relevant for M87 ∗ . Inparticular, the electron temperature could be calculatedin a more self-consistent fashion than via a temperatureratio prescription of Equation 10, by separately evolvingthe energy of ions and electrons (Sądowski et al. 2017).Then, one could evolve a population of nonthermal elec-trons (Chael et al. 2017), determine the electrons heat- . . . . a ∗ = 0 . R high = 10 MADSANE φ B , e x p φ j e t . . . . P r o b a b ili t y d e n s i t y MAD R high = 10 a ∗ = 0 . a ∗ = − . φ B , e x p φ j e t
150 200 250 300 350 η (deg)0 . . . . . MAD a ∗ = 0 . R high = 10 R high = 160 φ B , e x p φ j e t Figure 15.
Histograms of the brightness position angle η ,measured in the image domain for 500 ray-traced snapshotsof GRMHD simulations from the EHT M87 ∗ Simulation Li-brary (EHTC V). In each panel, the blue histogram repre-sents the same fixed model S1: MAD, a ∗ = 0 . , R high = 10 ,and the red histogram represents a model in which a singleparameter has been altered with respect to the fixed model(S2-S4). Orientations measured in the 2009–2017 observa-tional data sets with an ML estimator are indicated, with68% confidence intervals (Table 3, two results from 2017 areshown without their very narrow error bars). The gray areaaround φ jet corresponds to the observed jet position anglevariation (Walker et al. 2018). ing with a well-motivated subgrid prescription (Chaelet al. 2018a; Davelaar et al. 2019), or even employ non-ideal MHD to model nonthermal emission caused by theparticles accelerated through magnetic reconnection ina more self-consistent manner (Ripperda et al. 2019). Inthe current analysis we also make an assumption of notilt between the plane of accretion and the black holespin. If the tilt is present, an additional degree of free-dom (“camera longitude”) corresponding to a positionangle of the black hole spin (misaligned with the jet po-sition angle φ jet ) in the image plane will influence theobserved crescent orientation (Chatterjee et al. 2020).4In that case, analysis of the position angle distribu-tions could place joint constraints on the tilt magnitude,the longitude, and other parameters of the simulation.Large-scale parameter surveys with these extensions tothe GRMHD setup are currently precluded by the im-mense computational costs.A separate concern is whether the orientation of cres-cent models fitted to the VLBI data is consistent withthe image domain η (Equation 11). In the case of thesynthetic data sets considered in Section 4, the twoGRMHD data sets exhibited quite large biases whilethe two MODEL data sets showed a high level of consis-tency, so this issue requires further study. Characteriz-ing GRMHD simulations in terms of VLBI observablesis the subject of continued work. SUMMARYWe have performed geometric modeling of the 2009–2017 EHT observations of M87 ∗ . Motivated by EHTimaging and modeling results using the 2017 observa-tions and the stability of fits across the archival obser-vations, we have used a simple asymmetric ring model.We found that the fitted ring diameter is stable through-out these observations, which strongly argues in favor ofits association with the shadow of a supermassive blackhole. We observe indications of modest intrinsic vari-ability in the total flux density of the ring and in itsposition angle.Specifically, we find the brightness asymmetry alongthe east-west direction in the previously unpublished2013 observations, while all other data sets are consis-tent with the north-south asymmetry direction seen inthe 2017 EHT images. This degree of position anglevariation is seen in some GRMHD simulations of M87 ∗ ,while others do not show position angle variations asbroad as those observed between 2013 and 2017. Thus,the source variability over these observations providesnew constraints on the simulation parameters, includingthe black hole spin, accretion flow magnetization, andelectron heating model. As an example, the GRMHDMAD model with spin a ∗ = 0 . and R high = 160 (lastpanel of Figure 15), which was determined to be viableby EHTC V, is inconsistent with the presented positionangle measurements. We expect that unmodeled phys-ical effects such as black hole and accretion flow spinmisalignment may also be important in interpreting thisvariability.Our results extend the temporal span of EHT con-straints on the ring morphology by nearly three ordersof magnitude, from ∼ M over the 2017 observationsto ∼ M between the 2009 and 2017 campaigns. Be-cause the correlation timescale for M87 ∗ is expected tobe at least a few tens of M , the longer span is critical fordecoupling stable image features such as the black holeshadow from transient features associated with the tur-bulent accretion flow. As continued EHT observationsbecome available, the variation of the estimated posi- tion angle should allow us to discriminate between viableGRMHD models, providing constraints on the physicalparameters of M87 ∗ and opening an exciting new av-enue for quantitative time-domain studies of structuralvariability in M87 ∗ .ACKNOWLEDGEMENTSThe authors of the present paper thank the fol-lowing organizations and programs: the Academy ofFinland (projects 274477, 284495, 312496); the Ad-vanced European Network of E-infrastructures for As-tronomy with the SKA (AENEAS) project, supportedby the European Commission Framework ProgrammeHorizon 2020 Research and Innovation action undergrant agreement 731016; the Alexander von Hum-boldt Stiftung; the Black Hole Initiative at HarvardUniversity, through a grant (60477) from the JohnTempleton Foundation; the China Scholarship Coun-cil; Comisión Nacional de Investigación Científica yTecnológica (CONICYT, Chile, via PIA ACT172033,Fondecyt projects 1171506 and 3190878, BASAL AFB-170002, ALMA-conicyt 31140007); Consejo Nacionalde Ciencia y Tecnología (CONACYT, Mexico, projects104497, 275201, 279006, 281692); the Delaney Fam-ily via the Delaney Family John A. Wheeler Chair atPerimeter Institute; Dirección General de Asuntos delPersonal Académico-—Universidad Nacional Autónomade México (DGAPA-—UNAM, project IN112417); theEuropean Research Council Synergy Grant "BlackHole-Cam: Imaging the Event Horizon of Black Holes"(grant 610058); the Generalitat Valenciana postdoc-toral grant APOSTD/2018/177 and GenT Program(project CIDEGENT/2018/021); the Gordon and BettyMoore Foundation (grants GBMF-3561, GBMF-5278);the Istituto Nazionale di Fisica Nucleare (INFN)sezione di Napoli, iniziative specifiche TEONGRAV;the International Max Planck Research School forAstronomy and Astrophysics at the Universities ofBonn and Cologne; the Jansky Fellowship program ofthe National Radio Astronomy Observatory (NRAO);the Japanese Government (Monbukagakusho: MEXT)Scholarship; the Japan Society for the Promotion of Sci-ence (JSPS) Grant-in-Aid for JSPS Research Fellowship(JP17J08829); the Key Research Program of FrontierSciences, Chinese Academy of Sciences (CAS, grantsQYZDJ-SSW-SLH057, QYZDJ-SSW-SYS008, ZDBS-LY-SLH011); the Leverhulme Trust Early Career Re-search Fellowship; the Max-Planck-Gesellschaft (MPG);the Max Planck Partner Group of the MPG and theCAS; the MEXT/JSPS KAKENHI (grants 18KK0090,JP18K13594, JP18K03656, JP18H03721, 18K03709,18H01245, 25120007); the MIT International Scienceand Technology Initiatives (MISTI) Funds; the Min-istry of Science and Technology (MOST) of Taiwan(105-2112-M-001-025-MY3, 106-2112-M-001-011, 106-2119-M-001-027, 107-2119-M-001-017, 107-2119-M-001-020, and 107-2119-M-110-005); the National Aero- onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT Abramowski, A., Acero, F., Aharonian, F., et al. 2012, ApJ,746, 151, doi: 10.1088/0004-637X/746/2/151Akiyama, K., Lu, R.-S., Fish, V. L., et al. 2015, ApJ, 807,150, doi: 10.1088/0004-637X/807/2/150Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010,arXiv e-prints, arXiv:1012.3754.https://arxiv.org/abs/1012.3754Arras, P., Frank, P., Haim, P., et al. 2020, arXiv e-prints,arXiv:2002.05218. https://arxiv.org/abs/2002.05218Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214,doi: 10.1086/170270Bardeen, J. M. 1973, in Black Holes (Les Astres Occlus),ed. C. Dewitt & B. S. Dewitt, 215–239Benkevitch, L., Akiyama, K., Lu, R., Doeleman, S., & Fish,V. 2016, arXiv e-prints, arXiv:1609.00055.https://arxiv.org/abs/1609.00055Blackburn, L., Pesce, D. W., Johnson, M. D., et al. 2020,ApJ, 894, 31, doi: 10.3847/1538-4357/ab8469Blackburn, L., Chan, C.-k., Crew, G. B., et al. 2019, ApJ,882, 23, doi: 10.3847/1538-4357/ab328dBlakeslee, J. P., Jordán, A., Mei, S., et al. 2009, ApJ, 694,556, doi: 10.1088/0004-637X/694/1/556Bower, G. C., Dexter, J., Markoff, S., et al. 2015, ApJL,811, L6, doi: 10.1088/2041-8205/811/1/L6Broderick, A. E., & Loeb, A. 2009, ApJ, 697, 1164,doi: 10.1088/0004-637X/697/2/1164Broderick, A. E., Fish, V. L., Johnson, M. D., et al. 2016,ApJ, 820, 137, doi: 10.3847/0004-637X/820/2/137Broderick, A. E., Gold, R., Karami, M., et al. 2020, ApJ,897, 139, doi: 10.3847/1538-4357/ab91a4Chael, A., Rowan, M., Narayan, R., Johnson, M., & Sironi,L. 2018a, MNRAS, 478, 5209,doi: 10.1093/mnras/sty1261 Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018b,ApJ, 857, 23, doi: 10.3847/1538-4357/aab6a8Chael, A. A., Johnson, M. D., Narayan, R., et al. 2016,ApJ, 829, 11, doi: 10.3847/0004-637X/829/1/11Chael, A. A., Narayan, R., & Sadowski, A. 2017, MNRAS,470, 2367, doi: 10.1093/mnras/stx1345Chan, C.-k., Psaltis, D., Özel, F., et al. 2015, ApJ, 812,103, doi: 10.1088/0004-637X/812/2/103Chatterjee, K., Younsi, Z., Liska, M., et al. 2020, arXive-prints, arXiv:2002.08386.https://arxiv.org/abs/2002.08386Davelaar, J., Olivares, H., Porth, O., et al. 2019, A&A, 632,A2, doi: 10.1051/0004-6361/201936150Dexter, J., Agol, E., Fragile, P. C., & McKinney, J. C. 2010,ApJ, 717, 1092, doi: 10.1088/0004-637X/717/2/1092Doeleman, S., Agol, E., Backer, D., et al. 2009, inastro2010: The Astronomy and Astrophysics DecadalSurvey, Vol. 2010, 68Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012,Science, 338, 355, doi: 10.1126/science.1224768Event Horizon Telescope Collaboration, Akiyama, K.,Alberdi, A., et al. 2019a, ApJL, 875, L1,doi: 10.3847/2041-8213/ab0ec7—. 2019b, ApJL, 875, L2, doi: 10.3847/2041-8213/ab0c96—. 2019c, ApJL, 875, L3, doi: 10.3847/2041-8213/ab0c57—. 2019d, ApJL, 875, L4, doi: 10.3847/2041-8213/ab0e85—. 2019e, ApJL, 875, L5, doi: 10.3847/2041-8213/ab0f43—. 2019f, ApJL, 875, L6, doi: 10.3847/2041-8213/ab1141Falcke, H., Melia, F., & Agol, E. 2000, ApJL, 528, L13,doi: 10.1086/312423Fish, V. L., Doeleman, S. S., Beaudoin, C., et al. 2011,ApJL, 727, L36, doi: 10.1088/2041-8205/727/2/L36 onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT Fish, V. L., Johnson, M. D., Doeleman, S. S., et al. 2016,ApJ, 820, 90, doi: 10.3847/0004-637X/820/2/90Fromm, C. M., Younsi, Z., Baczko, A., et al. 2019, A&A,629, A4, doi: 10.1051/0004-6361/201834724Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ,589, 444, doi: 10.1086/374594Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ,729, 119, doi: 10.1088/0004-637X/729/2/119Goddi, C., Martí-Vidal, I., Messias, H., et al. 2019, PASP,131, 075003, doi: 10.1088/1538-3873/ab136aHonma, M., Akiyama, K., Uemura, M., & Ikeda, S. 2014,PASJ, 66, 95, doi: 10.1093/pasj/psu070Issaoun, S., Folkers, T. W., Marrone, D. P., et al. 2017,EHT Memo Series, 2017-CE-03Janssen, M., Goddi, C., van Bemmel, I. M., et al. 2019,A&A, 626, A75, doi: 10.1051/0004-6361/201935181Johnson, M. D., Fish, V. L., Doeleman, S. S., et al. 2015,Science, 350, 1242, doi: 10.1126/science.aac7087Kamruddin, A. B., & Dexter, J. 2013, MNRAS, 434, 765,doi: 10.1093/mnras/stt1068Kim, J., Marrone, D. P., Chan, C.-K., et al. 2016, ApJ, 832,156, doi: 10.3847/0004-637X/832/2/156Kim, J. Y., Krichbaum, T. P., Lu, R. S., et al. 2018, A&A,616, A188, doi: 10.1051/0004-6361/201832921Liddle, A. R. 2007, MNRAS, 377, L74,doi: 10.1111/j.1745-3933.2007.00306.xLister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ,138, 1874, doi: 10.1088/0004-6256/138/6/1874Lu, R.-S., Fish, V. L., Akiyama, K., et al. 2013, ApJ, 772,13, doi: 10.1088/0004-637X/772/1/13Lu, R.-S., Krichbaum, T. P., Roy, A. L., et al. 2018, ApJ,859, 60, doi: 10.3847/1538-4357/aabe2eLuminet, J.-P. 1979, A&A, 75, 228Macchetto, F., Marconi, A., Axon, D. J., et al. 1997, ApJ,489, 579, doi: 10.1086/304823Matthews, L. D., Crew, G. B., Doeleman, S. S., et al. 2018,PASP, 130, 015002, doi: 10.1088/1538-3873/aa9c3dMościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A,586, A38, doi: 10.1051/0004-6361/201526630Mościbrodzka, M., & Gammie, C. F. 2018, MNRAS, 475,43, doi: 10.1093/mnras/stx3162 Mościbrodzka, M., Shiokawa, H., Gammie, C. F., &Dolence, J. C. 2012, ApJL, 752, L1,doi: 10.1088/2041-8205/752/1/L1Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A.2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K.2012, MNRAS, 426, 3241,doi: 10.1111/j.1365-2966.2012.22002.xPorth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS,243, 26, doi: 10.3847/1538-4365/ab29fdRipperda, B., Bacchini, F., Porth, O., et al. 2019, ApJS,244, 10, doi: 10.3847/1538-4365/ab3922Rogers, A. E. E., Doeleman, S. S., & Moran, J. M. 1995,AJ, 109, 1391, doi: 10.1086/117371Sądowski, A., Narayan, R., Penna, R., & Zhu, Y. 2013,MNRAS, 436, 3856, doi: 10.1093/mnras/stt1881Sądowski, A., Wielgus, M., Narayan, R., et al. 2017,MNRAS, 466, 705, doi: 10.1093/mnras/stw3116Sun, X.-N., Yang, R.-Z., Rieger, F. M., Liu, R.-Y., &Aharonian, F. 2018, A&A, 612, A106,doi: 10.1051/0004-6361/201731716Thompson, A. R., Moran, J. M., & Swenson, George W., J.2017, Interferometry and Synthesis in Radio Astronomy,3rd Edition (Springer), doi: 10.1007/978-3-319-44431-4van Cittert, P. H. 1934, Physica, 1, 201,doi: /10.1016/S0031-8914(34)90026-4Vincent, F. H., Wielgus, M., Abramowicz, M. A., et al.2020, arXiv e-prints, arXiv:2002.09226.https://arxiv.org/abs/2002.09226Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., &Junor, W. 2018, ApJ, 855, 128,doi: 10.3847/1538-4357/aaafccWalsh, J. L., Barth, A. J., Ho, L. C., & Sarzi, M. 2013,ApJ, 770, 86, doi: 10.1088/0004-637X/770/2/86White, C. J., Dexter, J., Blaes, O., & Quataert, E. 2020,ApJ, 894, 14, doi: 10.3847/1538-4357/ab8463Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529,doi: 10.1146/annurev-astro-082812-141003Zernike, F. 1938, Physica, 5, 785,doi: 10.1016/S0031-8914(38)80203-2 A. ML ESTIMATORS FOR THE RG MODELFigure 16 presents RG models best-fitting the synthetic data sets, following the procedures described in Section 4.
GRMHD1 I N P U T µ as GRMHD2 MODEL1 MODEL2
Figure 16.
Maximum likelihood estimators corresponding to the fits to four synthetic images, shown in the first row (noblurring). Estimators were obtained through synthetic VLBI observations with the ( u, v ) -coverage and uncertainties identicalto those of the real observations performed in 2009–2017. The blurred ring model (RG) was used, and the presented images ofML estimators were blurred to a µ as resolution. Blue dashed lines indicate the convention for the expected position angle ofthe bright component φ B , exp . The gray bar represents the ML estimate of φ B . The dashed circles correspond to a diameter of42 µ as. onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT B. CORNER PLOTS FOR THE RT AND THE RG MODELSWe present the posterior probability distributions corresponding to fitting the 2009–2013 data sets with the RTmodel, Figure 17, and with the RG model, Figure 18. Similarly, for the 2017 data sets, we show the posteriordistributions obtained for the RT model, Figure 19, and for the RG model, Figure 20. φ B ( d e g ) . . . . . f w . . . . I ( J y )
32 40 48 56 d ( µ as) . . . . . β
80 160 240 320 φ B (deg) . . . . . f w . . . . I (Jy) . . . . . β φ B ( d e g ) . . . . . f w . . . . . I ( J y )
30 36 42 48 d ( µ as) . . . . . β
80 160 240 320 φ B (deg) . . . . . f w .
70 0 .
75 0 .
80 0 .
85 0 . I (Jy) . . . . . β − φ B ( d e g ) . . . . . f w . . . . . I ( J y )
30 33 36 39 42 d ( µ as) . . . . . β −
80 0 80 160 240 φ B (deg) . . . . . f w . . . . . I (Jy) . . . . . β φ B ( d e g ) . . . . f w . . . . I ( J y ) . . . . . d ( µ as) . . . . . β
80 160 240 320 φ B (deg) .
08 0 .
16 0 .
24 0 . f w .
51 0 .
54 0 .
57 0 . I (Jy) . . . . . β Figure 17.
Posterior distributions of RT model parameters resulting from fitting to M87 ∗ Themis . The maximum likelihood estimate is indicated with red lines. Contours indicate 0.2, 0.5, 0.8 and 0.95 of the posteriorvolume. φ B ( d e g ) f w . . . . . I ( J y )
15 30 45 60 75 d ( µ as) . . . . . β
80 160 240 320 φ B (deg) f w . . . . . I (Jy) . . . . . β φ B ( d e g ) . . . . f w . . . . I ( J y )
20 30 40 50 d ( µ as) . . . . . β
80 160 240 320 φ B (deg) . . . . f w .
76 0 .
80 0 .
84 0 . I (Jy) . . . . . β − φ B ( d e g ) f w . . . . . I ( J y )
16 24 32 40 d ( µ as) . . . . . β −
80 0 80 160 240 φ B (deg) f w . . . . . I (Jy) . . . . . β φ B ( d e g ) . . . . f w . . . . I ( J y ) . . . . . d ( µ as) . . . . . β
80 160 240 320 φ B (deg) .
05 0 .
10 0 .
15 0 . f w .
51 0 .
54 0 .
57 0 . I (Jy) . . . . . β Figure 18.
Same as Figure 17, but for the RG model. onitoring the Morphology of M87 ∗ in 2009–2017 with the EHT φ B ( d e g ) . . . . f w . . . . . I ( J y ) . . . . . d ( µ as) . . . . β
152 154 156 158 φ B (deg) .
300 0 .
325 0 .
350 0 . f w .
33 0 .
36 0 .
39 0 .
42 0 . I (Jy) .
225 0 .
240 0 .
255 0 . β . . . . φ B ( d e g ) . . . f w . . . . . I ( J y ) . . . . . d ( µ as) . . . . . β . . . . φ B (deg) .
405 0 .
420 0 . f w . . . . . I (Jy) .
17 0 .
18 0 .
19 0 .
20 0 . β Figure 19.
Same as Figure 17, but for the M87 ∗ observations on 2017 April 6th and April 11th fitted with the RT model. φ B ( d e g ) . . . . . f w . . . . I ( J y ) .
40 41 .
55 41 .
70 41 . d ( µ as) . . . . β
150 152 154 156 158 φ B (deg) . . . . . f w .
36 0 .
39 0 .
42 0 . I (Jy) .
360 0 .
375 0 .
390 0 . β . . . φ B ( d e g ) . . . f w . . . . I ( J y ) . . . . d ( µ as) . . . . β . . . φ B (deg) .
285 0 .
300 0 . f w .
42 0 .
45 0 .
48 0 . I (Jy) .
33 0 .
34 0 .
35 0 . β Figure 20.
Same as Figure 17, but for the M87 ∗∗