# The Event Horizon General Relativistic Magnetohydrodynamic Code Comparison Project

Oliver Porth, Koushik Chatterjee, Ramesh Narayan, Charles F. Gammie, Yosuke Mizuno, Peter Anninos, John G. Baker, Matteo Bugli, Chi-Kwan Chan, Jordy Davelaar, et al.

► To cite this version:

Oliver Porth, Koushik Chatterjee, Ramesh Narayan, Charles F. Gammie, Yosuke Mizuno, et al. The Event Horizon General Relativistic Magnetohydrodynamic Code Comparison Project. The Astrophysical Journal Supplement, 2019, 243 (2), pp.26. 10.3847/1538-4365/ab29fd hal-02129669

HAL Id: hal-02129669

https://hal.science/hal-02129669v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Event Horizon General Relativistic Magnetohydrodynamic Code Comparison Project

Oliver Porth<fnref n="1" />,<fnref n="2" /> ID, Koushik Chatterjee<fnref n="1" />, Ramesh Narayan<fnref n="3" />,<fnref n="4" /> ID, Charles F. Gammie<fnref n="5" /> ID, Yosuke Mizuno<fnref n="2" /> ID, Peter Anninos<fnref n="6" />, John G. Baker<fnref n="7" />,<fnref n="8" />, Matteo Bugli<fnref n="9" /> ID, Chi-kwan Chan<fnref n="10" />,<fnref n="11" /> ID, Jordy Davelaar¹² ID, Luca Del Zanna¹³,¹⁴, Zachariah B. Etienne¹⁵,¹⁶, P. Chris Fragile¹⁷ ID, Bernard J. Kelly<fnref n="7" />,¹⁸,¹⁹ ID, Matthew Liska<fnref n="1" />, Sera Markoff²⁰ ID, Jonathan C. McKinney²¹, Bhupendra Mishra²², Scott C. Noble<fnref n="7" />,²³ ID, Héctor Olivares<fnref n="2" /> ID, Ben Prather²⁴ ID, Luciano Rezzolla<fnref n="2" /> ID, Benjamin R. Ryan²⁵,²⁶ ID, James M. Stone²⁷, Niccolò Tomei¹³,¹⁴, Christopher J. White²⁸, and Ziri Younsi<fnref n="2" />,²⁹ ID

Kazunori Akiyama<fnref n="4" />,³⁰,³¹,³² ID, Antxon Alberdi³³ ID, Walter Alef³⁴, Keiichi Asada³⁵, Rebecca Azulay³⁴,³⁶,³⁷ ID, Anne-Kathrin Baczko³⁴ ID, David Ball<fnref n="10" />, Mislav Balokovic<fnref n="3" />,<fnref n="4" /> ID, John Barrett³¹ ID, Dan Bintley³⁸, Lindy Blackburn<fnref n="3" />,<fnref n="4" /> ID, Wilfred Boland³⁹, Katherine L. Bouman<fnref n="3" />,<fnref n="4" />,⁴⁰ ID, Geoffrey C. Bower⁴¹ ID, Michael Bremer⁴², Christiaan D. Brinkerink¹² ID, Roger Brissenden<fnref n="3" />,<fnref n="4" /> ID, Silke Britzen³⁴ ID, Avery E. Broderick⁴³,⁴⁴,⁴⁵ ID, Dominique Broguiere⁴², Thomas Bronzwaer¹², Do-Young Byun⁴⁶,⁴⁷ ID, John E. Carlstrom⁴⁸,⁴⁹,⁵⁰,⁵¹, Andrew Chael<fnref n="3" />,<fnref n="4" /> ID, Shami Chatterjee⁵² ID, Ming-Tang Chen⁴¹, Yongjun Chen (陈永军)⁵³,⁵⁴ ID, Ilje Cho⁴⁶,⁴⁷ ID, Pierre Christian<fnref n="3" />,<fnref n="10" /> ID, John E. Conway⁵⁵ ID, James M. Cordes⁵² ID, Geoffrey, B. Crew³¹ ID, Yuzhu Cui⁵⁶,⁵⁷ ID, Mariafelicia De Laurentis<fnref n="2" />,⁵⁸,⁵⁹ ID, Roger Deane⁶⁰,⁶¹ ID, Jessica Dempsey³⁸ ID, Gregory Desvignes³⁴ ID, Sheperd S. Doeleman<fnref n="3" />,<fnref n="4" /> ID, Ralph P. Eatough³⁴ ID, Heino Falcke¹² ID, Vincent L. Fish³¹ ID, Ed Fomalont³⁰ ID, Raquel Fraga-Encinas¹² ID, Bill Freeman⁶²,⁶³, Per Friberg³⁸, Christian M. Fromm<fnref n="2" />, José L. Gómez³³ ID, Peter Galison<fnref n="4" />,⁶⁴,⁶⁵ ID, Roberto García⁴², Olivier Gentaz⁴², Boris Georgiev⁴⁴ ID, Ciriaco Goddi¹²,⁶⁶, Roman Gold<fnref n="2" /> ID, Minfeng Gu (顾敏峰)⁵³,⁶⁷ ID, Mark Gurwell<fnref n="3" /> ID, Kazuhiro Hada⁵⁶,⁵⁷ ID, Michael H. Hecht³¹, Ronald Hesper⁶⁸ ID, Luis C. Ho (何子山)⁶⁹,⁷⁰ ID, Paul Ho³⁵, Mareki Honma⁵⁶,⁵⁷ ID, Chih-Wei L. Huang³⁵ ID, Lei Huang (黄磊)⁵³,⁶⁷ ID, David H. Hughes⁷¹, Shiro Ikeda³²,⁷²,⁷³,⁷⁴ ID, Makoto Inoue³⁵, Sara Issaoun¹² ID, David J. James<fnref n="3" />,<fnref n="4" /> ID, Buell T. Jannuzi<fnref n="10" />, Michael Janssen¹² ID, Britton Jeter⁴⁴ ID, Wu Jiang (江悟)⁵³ ID, Michael D. Johnson<fnref n="3" />,<fnref n="4" /> ID, Svetlana Jorstad⁷⁵,⁷⁶ ID, Taehyun Jung⁴⁶,⁴⁷ ID, Mansour Karami⁴³,⁷⁷ ID, Ramesh Karuppusamy³⁴ ID, Tomohisa Kawashima³² ID, Garrett K. Keating<fnref n="3" /> ID, Mark Kettenis⁷⁸ ID, Jae-Young Kim³⁴ ID, Junhan Kim<fnref n="10" /> ID, Jongsoo Kim⁴⁶, Motoki Kino³²,⁷⁹ ID, Jun Yi Koay³⁵ ID, Patrick, M. Koch³⁵ ID, Shoko Koyama³⁵ ID, Michael Kramer³⁴ ID, Carsten Kramer⁴² ID, Thomas P. Krichbaum³⁴ ID, Cheng-Yu Kuo⁸⁰, Tod R. Lauer⁸¹ ID, Sang-Sung Lee⁴⁶ ID, Yan-Rong Li (李彦荣)⁸³,⁸⁰ ID, Zhiyuan Li (李志远)⁸³,⁸⁴ ID, Michael Lindqvist⁵⁵, Kuo Liu³⁴ ID, Elisabetta Liuzzo⁸⁵ ID, Wen-Ping Lo³⁵,⁸⁶, Andrei P. Lobanov³⁴, Laurent Loinard⁸⁷,⁸⁸ ID, Colin Lonsdale³¹, Ru-Sen Lu (路如淼)³⁴,⁵³ ID, Nicholas R. MacDonald³⁴ ID, Jirong Mao (毛基荣)⁸⁹,⁹⁰,⁹¹ ID, Daniel P. Marrone<fnref n="10" /> ID, Alan P. Marscher⁷⁵ ID, Iván Martí-Vidal⁵⁵,⁹² ID, Satoki Matsushita³⁵, Lynn D. Matthews³¹ ID, Lia Medeiros<fnref n="10" />,⁹³ ID, Karl M. Menten³⁴ ID, Izumi Mizuno³⁸ ID, James M. Moran<fnref n="3" />,<fnref n="4" /> ID, Kotaro Moriyama⁵⁶ ID, Monika Moscibrodzka¹² ID, Cornelia Müller¹²,³⁴ ID, Hiroshi Nagai⁵⁷,⁹⁴ ID, Neil M. Nagar⁹⁵ ID, Masanori Nakamura³⁵ ID, Gopal Narayanan⁹⁶, Iniyan Natarajan⁶¹ ID, Roberto Neri⁴² ID, Chunchong Ni⁴⁴ ID, Aristeidis Noutsos³⁴ ID, Hiroki Okino⁵⁶,⁹⁷, Tomoaki Oyama⁵⁶, Feryal Özel<fnref n="10" />, Daniel C. M. Palumbo<fnref n="3" />,<fnref n="4" /> ID, Nimesh Patel<fnref n="3" /> ID, Ue-Li Pen⁴³,⁹⁸,⁹⁹,¹⁰⁰ ID, Dominic W. Pesce<fnref n="3" />,<fnref n="4" /> ID, Vincent Piétu⁴², Richard Plambeck¹⁰¹, Aleksandar PopStefanija⁹⁶, Jorge A. Preciado-López⁴³ ID, Dimitrios Psaltis<fnref n="10" /> ID, Hung-Yi Pu⁴³ ID, Venkatessh Ramakrishnan⁹⁵ ID, Ramprasad Rao⁴¹ ID, Mark G. Rawlings³⁸, Alexander W. Raymond<fnref n="3" />,<fnref n="4" /> ID, Bart Ripperda<fnref n="2" /> ID, Freek Roelofs¹² ID, Alan Rogers³¹, Eduardo Ros³⁴ ID, Mel Rose<fnref n="10" /> ID, Arash Roshanineshat<fnref n="10" />, Helge Rottmann³⁴, Alan L. Roy³⁴ ID, Chet Ruszczyk³¹ ID, Kazi L. J. Rygl⁸⁵ ID, Salvador Sánchez¹⁰², David Sánchez-Arguelles⁷¹,¹⁰³ ID, Mahito Sasada⁵⁶,¹⁰⁴ ID, Tuomas Savolainen<fnref n="3" />,<fnref n="4" />,¹⁰⁵,¹⁰⁶ ID, F. Peter Schloerb⁹⁶, Karl-Friedrich Schuster⁴², Lijing Shao³⁴,⁷⁰ ID, Zhiqiang Shen (沈志强)⁵³,⁵⁴ ID, Des Small⁷⁸ ID, Bong Won Sohn⁴⁶,⁴⁷,¹⁰⁷ ID, Jason SooHoo³¹ ID, Fumie Tazaki⁵⁶ ID, Paul Tiede⁴³,⁴⁴ ID, Remo P. J. Tilanus¹²,⁶⁶,¹⁰⁸ ID, Michael Titus³¹ ID, Kenji Toma¹⁰⁹,¹¹⁰ ID, Pablo Torne³⁴,¹⁰² ID, Tyler Trent<fnref n="10" />, Sascha Trippe¹¹¹ ID, Shuichiro Tsuda¹¹², Ilse van Bemmel⁷⁸ ID, Huib Jan van Langevelde⁷⁸,¹¹³ ID, Daniel R. van Rossum¹² ID, Jan Wagner³⁴, John Wardle¹¹⁴ ID, Jonathan Weintroub<fnref n="3" />,<fnref n="4" /> ID, Norbert Wex³⁴, Robert Wharton³⁴ ID, Maciek Wielgus<fnref n="3" />,<fnref n="4" /> ID, George N. Wong²⁴ ID, Qingwen Wu (吴庆文)¹¹⁵ ID, Ken Young<fnref n="3" /> ID, André Young¹² ID, Feng Yuan (袁峰)⁵³,⁶⁷,¹¹⁶ ID, Ye-Fei Yuan (袁业飞)¹¹⁷ ID, J. Anton Zensus³⁴ ID, Guangyao Zhao⁴⁶ ID, Shan-Shan Zhao¹²,⁸³ ID, and Ziyan Zhu⁶⁵

(The Event Horizon Telescope Collaboration)

<footnote n="1">¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH, Amsterdam, The Netherlands; o.porth@uva.nl</footnote>

<footnote n="2">² Institut für Theoretische Physik, Goethe-Universität Frankfurt, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany</footnote>

<footnote n="3">³ Center for Astrophysics|Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="4">⁴ Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="5">⁵ Department of Astronomy and Department of Physics, University of Illinois, Urbana, IL 61801, USA</footnote>

<footnote n="6">⁶ Lawrence Livermore National Laboratory, Livermore, CA 94550, USA</footnote>

<footnote n="7">⁷ Gravitational Astrophysics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="8">⁸ Joint Space-Science Institute, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="9">⁹ IRFU/Département d'Astrophysique, Laboratoire AIM, CEA/DRF-CNRS-Université Paris Diderot, CEA-Saclay F-91191, France</footnote>

<footnote n="10">¹⁰ Steward Observatory and Department of Astronomy, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA</footnote>

<footnote n="11">¹¹ Data Science Institute, University of Arizona, 1230 N. Cherry Avenue, Tucson, AZ 85721, USA</footnote>
¹² Department of Astrophysics, Institute for Mathematics, Astrophysics and Particle Physics (IMAPP), Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

¹³ Dipartimento di Fisica e Astronomia, Università di Firenze e INFN—Sez. di Firenze, via G. Sansone 1, I-50019 Sesto F.no, Italy

¹⁴ INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

¹⁵ Department of Mathematics, West Virginia University, Morgantown, WV 26506, USA

¹⁶ Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA

¹⁷ Department of Physics & Astronomy, College of Charleston, 66 George Street, Charleston, SC 29424, USA

¹⁸ Department of Physics, University of Maryland, Baltimore County, Baltimore, MD 21250, USA

¹⁹ CRESST, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

²⁰ Anton Pannekoek Institute for Astronomy & GRAPPA, University of Amsterdam, Postbus 94249, 1090GE Amsterdam, The Netherlands
²¹ H2O.ai, 2307 Leghorn Street, Mountain View, CA 94043, USA

²² JILA, University of Colorado and National Institute of Standards and Technology, 440 UCB, Boulder, CO 80309-0440, USA

²³ Department of Physics and Engineering Physics, The University of Tulsa, Tulsa, OK 74104, USA

²⁴ Department of Physics, University of Illinois, 1110 West Green Street, Urbana, IL 61801, USA

²⁵ CCS-2, Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, USA

²⁶ Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

²⁷ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

²⁸ Kayli Institute for Theoretical Physics, University of California Santa Barbara, Kohn Hall, Santa Barbara, CA 93107, USA

²⁹ Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey, RH5 6NT, UK

³⁰ National Radio Astronomy Observatory, 520 Edgemont Road, 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

³³ Instituto de Astrofísica de Andalucía—CSIC, Glorieta de la Astronomía s/n, E-18008 Granada, Spain

³⁴ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

³⁵ 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

³⁸ East Asian Observatory, 660 N. A'ohoku Place, Hilo, HI 96720, USA

³⁹ Nederlandse Onderzoeksgool voor Astronomie (NOVA), P.O. Box 9513, 2300 RA Leiden, the Netherlands, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

⁴⁰ California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

⁴¹ Institute of Astronomy and Astrophysics, Academia Sinica, 645 N. A'ohoku Place, Hilo, HI 96720, USA

⁴² Institut de Radioastronomie Millimétrique, 300 rue de la Piscine, F-38406 Saint Martin d'Hères, France

⁴³ 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, ON N2L 3G1, Canada

⁴⁶ Korea Astronomy and Space Science Institute, Daedok-daero 776, Yuseong-gu, Daejeon 34055, Republic of Korea

⁴⁷ University of Science and Technology, Gajeong-ro 217, Yuseong-gu, Daejeon 34113, Republic of Korea

⁴⁸ Kayli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA

⁴⁹ Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA

⁵⁰ Department of Physics, University of Chicago, Chicago, IL 60637, USA

⁵¹ Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA

⁵² Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, NY 14853, USA

⁵³ Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

⁵⁴ 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-439 92 Onsala, Sweden

⁵⁶ 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, Tokyo 181-8588, Japan

⁵⁸ Dipartimento di Fisica “E. Pancini,” Università di Napoli “Federico II,” Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy

⁵⁹ INFN Sez. di Napoli, Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy

⁶⁰ 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, Grahamstown 6140, South Africa

⁶² Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, 32-D476, 77 Massachusetts Avenue, Cambridge, MA 02142, USA

⁶³ Google Research, 355 Main Street, Cambridge, MA 02142, USA

⁶⁴ 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

⁶⁷ 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

⁷⁰ Kayli 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

⁷² The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa, Tokyo, 190-8562, Japan

⁷³ Department of Statistical Science, The Graduate University for Advanced Studies (SOKENDAI), 10-3 Midori-cho, Tachikawa, Tokyo 190-8562, Japan

⁷⁴ Kayli Institute for the Physics and Mathematics of the Universe, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, 277-8583, Japan

⁷⁵ 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

⁷⁷ University of Waterloo, 200 University Avenue West, Waterloo, Ontario N2L 3G1, Canada

⁷⁸ Joint Institute for VLBI ERIC (JIVE), Oude Hoogeveensdijk 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 Road, Kaosiung City 80424, Taiwan, R.O.C.

⁸¹ National Optical Astronomy Observatory, 950 North Cherry Avenue, 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 Road, Taipei 10617, Taiwan, R.O.C

⁸⁷ 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

⁸⁹ 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's Republic of China

⁹¹ Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese Academy of Sciences, 650011 Kunming, People's Republic of China

⁹² Centro Astronómico de Yebes (IGN), Apartado 148, E-19180 Yebes, Spain

⁹³ Department of Physics, Broida Hall, University of California Santa Barbara, Santa Barbara, CA 93106, USA

⁹⁴ National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan

⁹⁵ 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

¹⁰¹ Radio Astronomy Laboratory, University of California, Berkeley, CA 94720, USA

¹⁰² Instituto de Radioastronomía Milimétrica, IRAM, Avenida Divina Pastora 7, Local 20, E-18012, Granada, Spain

¹⁰³ Consejo Nacional de Ciencia y Tecnología, Av. Insurgentes Sur 1582, 03940, Ciudad de México, México

¹⁰⁴ 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

¹¹² Mizusawa VLBI Observatory, National Astronomical Observatory of Japan, Hoshigaoka 2-12, Mizusawa-ku, Oshu-shi, Iwate 023-0861, Japan

¹¹³ Leiden Observatory, Leiden University, Postbus 2300, 9513 RA Leiden, The Netherlands

¹¹⁴ Physics Department, Brandeis University, 415 South Street, Waltham, MA 02453, USA

¹¹⁵ School of Physics, Huazhong University of Science and Technology, Wuhan, Hubei, 430074, People's Republic of China

¹¹⁶ 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

Received 2019 April 5; revised 2019 May 27; accepted 2019 May 28; published 2019 August 1

## Abstract

Recent developments in compact object astrophysics, especially the discovery of merging neutron stars by LIGO, the imaging of the black hole in M87 by the Event Horizon Telescope, and high-precision astrometry of the Galactic Center at close to the event horizon scale by the GRAVITY experiment motivate the development of numerical source models that solve the equations of general relativistic magnetohydrodynamics (GRMHD). Here we compare GRMHD solutions for the evolution of a magnetized accretion flow where turbulence is promoted by the magnetorotational instability from a set of nine GRMHD codes: Athena++, BHAC, Cosmos++, ECHO, H-AMR, iharm3D, HARM-Noble, IllinoisGRMHD, and KORAL. Agreement among the codes improves as resolution increases, as measured by a consistently applied, specially developed set of code performance metrics. We conclude that the community of GRMHD codes is mature, capable, and consistent on these test problems.

**Key words:** black hole physics – magnetic fields – magnetohydrodynamics (MHD) – methods: numerical – relativistic processes

## 1. Introduction

Fully general relativistic models of astrophysical sources are in high demand, not only since the discovery of gravitational waves emitted from merging stellar mass black holes (Abbott et al. 2016). The need for an accurate description of the interplay between strong gravity, matter, and electromagnetic fields is further highlighted by the recent detection of electromagnetic counterpart radiation to the coalescence of a neutron star binary (Abbott et al. 2017). Our own effort is motivated by the Event Horizon Telescope (EHT) project, which allows direct imaging of hot, luminous plasma near a black hole (hereafter BH) event horizon (Doeleman et al. 2008; Goddi et al. 2017; Event Horizon Telescope Collaboration et al. 2019a). The main targets of the EHT are the BH at the center of the Milky Way (also known by the name of the associated radio source, Sgr A*; e.g., Lu et al. 2018) and the BH at the center of the galaxy M87 with the associated central radio source M87* (Doeleman et al. 2012;

Akiyama et al. 2015). In order to extract information on the dynamics of the plasma that lies within a few $GM/c^2$ of the event horizon ($M \equiv \text{BH mass}$) as well as information about the black hole's gravitational field, it is necessary to develop models of the accretion flow, associated winds and relativistic jets, and the emission properties in each of the components.

Earlier semianalytic works (Narayan & Yi 1995; Narayan et al. 1998; Yuan et al. 2002) have provided the general parameter regime of the Galactic Center by exploiting spectral information. For example, Mahadevan & Quataert (1997) demonstrated that electrons and ions are only weakly collisionally coupled and unlikely to be in thermal equilibrium. Also, key parameters like the accretion rate are typically estimated based on simple one-dimensional models (Marrone et al. 2007). They have solidified the notion that the accretion rate in Sgr A* is far below the Eddington limit $\dot{M}_{\text{Edd}} = L_{\text{Edd}}/(0.1c^2) \simeq 2 M/(10^8 M_{\odot}) M_{\odot} \text{yr}^{-1}$, where $L_{\text{Edd}} = 4\pi GMc/\sigma_T$ is the Eddington luminosity (with $\sigma_T$ being the Thomson electron cross section). New observational
capabilities like millimeter and IR interferometry, provided by the EHT and GRAVITY (Gravity Collaboration et al. 2018) collaborations, now allow us to go much closer to the source, which requires a description of general relativistic and dynamical (hence time-dependent) effects.

The most common approach to dynamical relativistic source modeling uses the ideal general relativistic magnetohydrodynamic (GRMHD) approximation. It is worth reviewing the nature and quality of the two approximations inherent in the GRMHD model. First, the plasma is treated as a fluid rather than as a collisionless plasma. Second, the exchange of energy between the plasma and the radiation field is ignored.

The primary EHT sources Sgr A\* and M87\* fall in the class of low-luminosity active galactic nuclei (AGN) and accrete with $\dot{M}/\dot{M}_{\text{Edd}} \lesssim 10^{-6}$ (Marrone et al. 2007) and $\dot{M}/\dot{M}_{\text{Edd}} \lesssim 10^{-5}$ (Kuo et al. 2014) far below the Eddington limit. In both cases, the accretion flow is believed to form an optically thin disk that is geometrically thick and therefore has temperature comparable to the virial temperature (see Yuan & Narayan 2014 for a review). The plasma is at sufficiently high temperature and low density that it is collisionless: ions and electrons can travel $\gg GM/c^2$ along magnetic field lines before being significantly deflected by Coulomb scattering, while the effective mean free path perpendicular to field lines is the gyroradius, which is typically $\ll GM/c^2$. A rigorous description of the accreting plasma would thus naively require integrating the Boltzmann equation at far greater expense than integrating the fluid equations. Full Boltzmann treatments of accretion flows have so far been limited to the study of localized regions within the source (e.g., Hoshino 2015; Kunz et al. 2016). Global models that incorporate nonideal effects using PIC-inspired closure models suggest, however, that the effects of thermal conduction and pressure anisotropy (viscosity) are small (Chandra et al. 2015, 2017; Foucart et al. 2017), and thus that one would not do too badly with an ideal fluid prescription.

For Sgr A\*, radiative cooling is negligible (Dibi et al. 2012). For M87, radiative cooling is likely important (e.g., Mościbrodzka et al. 2011; Ryan et al. 2018; Chael et al. 2019). Cooling through the synchrotron process and via inverse Compton scattering primarily affects the electrons, which are weakly coupled to the ions and therefore need not be in thermal equilibrium with them. To properly treat the radiation field for the nonlocal process of Compton scattering requires solving the Boltzmann equation for photons (the radiative transport equation) in full (e.g., Ryan et al. 2015) or in truncated form with “closure.” A commonly employed closure is to assume the existence of a frame in which the radiation field can be considered isotropic, yielding the “M1” closure (Levermore 1984) for which a general relativistic derivation is shown, for example, in Sądowski et al. (2013). As expected, the computational demands imposed by the additional “radiation fluid” are considerable. It may, however, be possible to approximate the effects of cooling by using a suitable model to assign an energy density (or temperature) to the electrons (Mościbrodzka et al. 2016b). Again, an ideal fluid description, which automatically satisfies energy, momentum, and particle number conservation laws, is not a bad place to start.

It is possible to write the GRMHD equations in conservation form. This enables one to evolve the GRMHD equations using techniques developed to evolve other conservation laws such as those describing nonrelativistic fluids and magnetized fluids. Over the last decades, a number of GRMHD codes have been developed, most using conservation form, and applied to a

large variety of astrophysical scenarios (Hawley et al. 1984; Koide et al. 1999; De Villiers & Hawley 2003; Gammie et al. 2003; Anninos et al. 2005; Baiotti et al. 2005; Duez et al. 2005; Antón et al. 2006; Mizuno et al. 2006; Del Zanna et al. 2007; Giacomazzo & Rezzolla 2007; Tchekhovskoy et al. 2011; Radice & Rezzolla 2013; McKinney et al. 2014; Radice et al. 2014; Sądowski et al. 2014; Etienne et al. 2015; Zanotti & Dumbser 2015; Meliani et al. 2016; White et al. 2016; Liska et al. 2018a).

Despite the conceptual simplicity of the MHD equations, the nonlinear properties, which allow for shocks and turbulence, render their treatment difficult. This is particularly true for the case study considered here: in state-of-the-art simulations of BH accretion, angular momentum transport is provided by Maxwell and Reynolds stresses of the orbiting plasma. MHD turbulence is seeded by the magnetorotational instability (MRI) in the differentially rotating disk (Balbus & Hawley 1991, 1998) and gives rise to chaotic behavior, which hinders strict convergence of the solutions. Nonetheless, it can be demonstrated that certain global properties of the solutions exhibit signs of convergence.

Another challenge is posed by the “funnel” region near the polar axis where low angular momentum material will be swallowed up by the BH (e.g., McKinney 2006). The strong magnetic fields that permeate the BH (held in place by the equatorial accretion flow) are able to extract energy in the form of Poynting flux from a rotating BH, giving rise to a relativistic “jet” (Blandford & Znajek 1977; Takahashi et al. 1990). The ensuing near-vacuum and magnetic dominance are traditionally difficult to handle for fluid-type simulations, but analytic calculations (e.g., Pu et al. 2015) and novel kinetic approaches (Parfrey et al. 2019) can be used to validate the flow in this region.

Due to their robustness when dealing with, e.g., supersonic jet outflows, current production codes typically employ high-resolution shock-capturing schemes in a finite-volume or finite-difference discretization (Font 2008; Rezzolla & Zanotti 2013; Martí & Müller 2015). However, new schemes, for example based on discontinuous Galerkin finite-element methods, are continuously being developed (Anninos et al. 2017; Fambri et al. 2018).

Given the widespread and increasing application of GRMHD simulations, it is critical for the community to evaluate the underlying systematics due to different numerical treatments and demonstrate the general robustness of the results. Furthermore, at the time of writing, all results on the turbulent BH accretion problem under consideration are obtained without explicitly resolving dissipative scales in the system (called the implicit large eddy simulation (ILES) technique). Hence, differences are in fact expected to prevail even for the highest resolutions achievable. Quantifying how large the systematic differences are is one of the main objectives in this first comprehensive code comparison of an accreting BH scenario relevant to the EHT science case. This work has directly informed the generation of the simulation library utilized in the modeling of the EHT 2017 observations (Event Horizon Telescope Collaboration et al. 2019b). We use independent production codes that differ widely in their algorithms and implementation. In particular, the codes that are being compared are Athena++, BHAC, Cosmos++, ECHO, H-AMR, iharm3D, HARM-Noble, IllinoisGRMHD, and KORAL. These codes are described further in Section 3 below.
The structure of this paper is as follows: in Section 2.1, we introduce the GRMHD system of equations, define the notation used throughout this paper, and briefly discuss the astrophysical problem. Code descriptions with references are given in Section 3. The problem setup is described in Section 4, where code-specific choices are also discussed. Results are presented in Section 5, and we close with a discussion in Section 6.

## 2. Astrophysical Problem

Let us first give a brief overview of the problem under investigation and the main characteristics of the accretion flow.

### 2.1. GRMHD System of Equations

For clarity and consistency of notation, we give a brief summary here of the ideal GRMHD equations in a coordinate basis $(t, x^i)$ with four-metric $(g_{\mu\nu})$ and metric determinant $g$. As is customary, Greek indices run through [0, 1, 2, 3] while Roman indices span [1, 2, 3]. The equations are solved in (geometric) code units (i.e., setting the gravitational constant and speed of light to unity $G = c = 1$), where, compared to the standard Gauss cgs system, the factor $1/\sqrt{4\pi}$ is further absorbed in the definition of the magnetic field. Hence, in the following, we will report times and radii simply in units of the mass of the central object $M$.

The equations describe particle number conservation:

$$ \partial_t(\sqrt{-g} \rho u^i) = -\partial_i(\sqrt{-g} \rho u^i), \qquad (1) $$

where $\rho$ is the rest-mass density and $u^i$ is the four-velocity; conservation of energy-momentum:

$$ \partial_t(\sqrt{-g} T^i_\nu) = -\partial_i(\sqrt{-g} T^i_\nu) + \sqrt{-g} T^{\kappa\lambda} \Gamma^\lambda_{\nu\kappa}, \qquad (2) $$

where $\Gamma^\lambda_{\nu\kappa}$ is the metric connection; the definition of the stress-energy tensor for ideal MHD:

$$ T_{\text{MHD}}^{\mu\nu} = (\rho + u + p + b^2)u^\mu u^\nu + \left( p + \frac{1}{2}b^2 \right) g^{\mu\nu} - b^\mu b^\nu, \qquad (3) $$

where $u$ is the fluid internal energy, $p$ the fluid pressure, and $b^\mu$ the magnetic field four-vector; the definition of $b^\mu$ in terms of magnetic field variables $B^i$, which are commonly used as “primitive” or dependent variables:

$$ b^t = B^i u^\mu g_{i\mu}, \qquad (4) $$

$$ b^i = (B^i + b^t u^i) / u^t; \qquad (5) $$

and the evolution equation for $B^i$, which follows from the source-free Maxwell equations:

$$ \partial_t(\sqrt{-g} B^i) = -\partial_j(\sqrt{-g} (b^j u^i - b^i u^j)). \qquad (6) $$

The system is subject to the additional no-monopoles constraint,

$$ \frac{1}{\sqrt{-g}} \partial_t(\sqrt{-g} B^i) = 0, \qquad (7) $$

which follows from the time component of the homogeneous Maxwell equations. For closure, we adopt the equation of state of an ideal gas. This takes the form $p = (\hat{\gamma} - 1)u$, where $\hat{\gamma}$ is the adiabatic index. More in-depth discussions of the ideal

GRMHD system of equations can be found in the various publications of the authors, e.g., Gammie et al. (2003), Anninos et al. (2005), Del Zanna et al. (2007), White et al. (2016), and Porth et al. (2017).

To establish a common notation for use in this paper, we note the following definitions: the magnetic field strength as measured in the fluid frame is given by $B := \sqrt{b^\alpha b_\alpha}$. This leads to the definition of the magnetization $\sigma := B^2/\rho$ and the plasma $\beta \beta := 2p/B^2$. In addition, we denote with $\Gamma$ the Lorentz factor with respect to the normal observer frame.

### 2.2. The Magnetized Black Hole Accretion Problem

We will now discuss the most important features of the problem at hand and introduce the jargon that has developed over the years. A schematic overview with key aspects of the accretion flow is given in Figure 1.

At a very low Eddington rate $\dot{M} \sim 10^{-6} \dot{M}_{\text{Edd}}$, the radiative cooling timescale becomes longer than the accretion timescale. In such radiatively inefficient accretion flows (RIAF), dynamics and radiation emission effectively decouple. For the primary EHT targets, Sgr A* and M87*, this is a reasonable first approximation, and hence, purely nonradiative GRMHD simulations ignoring cooling can be used to model the data. For an RIAF, the protons assume temperatures close to the virial one, which leads to an extremely “puffed-up” appearance of the tenuous accretion disk.

In the polar regions of the BH, plasma is either sucked in or expelled in an outflow, leaving behind a highly magnetized region called the funnel. The magnetic field of the funnel is held in place by the dynamic and static pressure of the disk. Because in ideal MHD, plasma cannot move across magnetic field lines (due to the frozen-in condition), there is no way to resupply the funnel with material from the accretion disk, and hence, the funnel would be completely devoid of matter if no pairs were created locally. In state-of-the-art GRMHD calculations, this is the region where numerical floor models inject a finite amount of matter to keep the hydrodynamic evolution viable.

The general morphology is separated into the following components: (i) the disk, which contains the bound matter, (ii) the evacuated funnel extending from the polar caps of the BH, and the (iii) jet sheath, which is the remaining outflowing matter. In Figure 1, the regions are indicated by commonly used discriminators in a representative simulation snapshot: the blue contour shows the bound/unbound transition defined via the geometric Bernoulli parameter $u_t = -1$,<fnref n="118" /> the red contour demarcates the funnel boundary $\sigma = 1$, and the green contour the equipartition $\beta = 1$, which is close to the bound/unbound line along the disk boundary (consistent with McKinney & Gammie 2004). In McKinney & Gammie (2004) a disk corona was also introduced for the material with $\beta \in [1, 3]$; however, as this choice is arbitrary, there is no compelling reason to label the corona as a separate entity in the RIAF scenario.

Because plasma is evacuated within the funnel, it has been suggested that unscreened electric fields in the charge-starved region can lead to particle acceleration which might fill the

<footnote n="118">¹¹⁸ There are various ways to define the unbound material in the literature. For example, McKinney & Gammie (2004) used the geometric Bernoulli parameter $u_t$. The hydrodynamic Bernoulli parameter used, for example, by Mościbrodzka et al. (2016a) is given by $-hu_t$, where $h = c^2 + p + u$ denotes the specific enthalpy. Chael et al. (2019) and Narayan et al. (2012) also included magnetic and radiative contributions. The geometric and hydrodynamic prescriptions certainly underestimate the amount of outflowing material.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Views of the radiatively inefficient turbulent black hole accretion problem at $t_{KS} = 10,000 M$ against the Kerr–Schild coordinates (subscript KS). Left: logarithmic rest-frame density (hue) and rendering of the magnetic field structure using line-integral convolution (luminance), showing the ordered field in the funnel region and turbulence in the disk. Center: the logarithm of the magnetization with colored contours indicating characteristics of the flow. The magnetized funnel is demarcated by $\sigma = 1$, (red lines), the disk is indicated by $\beta = 1$ (green lines), and the geometric Bernoulli criterion ($u_t = -1$) is given by the blue solid line in the region outside of the funnel. Right: schematic of the main components. In these plots, the BH horizon is the black disk and the ergosphere is shown as a black contour. The snapshot was obtained from a simulation with BHAC.</figcaption></figure>

magnetosphere via a pair cascade (e.g., Blandford & Znajek 1977; Beskin et al. 1992; Hirotani & Okamoto 1998; Levinson & Rieger 2011; Broderick & Tchekhovskoy 2015). The most promising alternative mechanism to fill the funnel region is by pair creation via $\gamma\gamma$ collisions of seed photons from the accretion flow itself (e.g., Stepney & Guilbert 1983; Phinney 1995). Neither of these processes is included in current state-of-the-art GRMHD simulations; however, the efficiency of pair formation via $\gamma\gamma$ collisions can be evaluated in postprocessing as demonstrated by Mościbrodzka et al. (2011).

Turning back to the morphology of the RIAF accretion, Figure 1, one can see that between the evacuated funnel demarcated by the funnel wall (red) and the bound disk material (blue), there is a strip of outflowing material often also referred to as the jet sheath (Dexter et al. 2012; Mościbrodzka & Falcke 2013; Mościbrodzka et al. 2016a; Davelaar et al. 2018). As argued by Hawley & Krolik (2006, who coined the alternative term funnel-wall jet for this region), this flow emerges as plasma from the disk is driven against the centrifugal barrier by magnetic and thermal pressure. In current GRMHD-based radiation models utilized, e.g., by the Event Horizon Telescope Collaboration et al. (2019b), as the density in the funnel region is dominated by the artificial floor model, the funnel is typically excised from the radiation transport. The denser region outside the funnel wall remains, which naturally leads to a limb-brightened structure of the observed M87 “jet” at radio frequencies (e.g., Mościbrodzka et al. 2016a; Chael et al. 2019; Davelaar et al. 2019). In the millimeter band (Event Horizon Telescope Collaboration et al. 2019a), the horizon-scale emission originates either from the body of the disk or from the region close to the funnel wall, depending on the assumptions on the electron temperatures (Event Horizon Telescope Collaboration et al. 2019b).

In RIAF accretion, a special role is played by the horizon-penetrating magnetic flux $\Phi_{BH}$: normalized by the accretion rate $\phi := \Phi_{BH}/\sqrt{M}$, it was shown that a maximum for the magnetic

flux $\phi_{max} \approx 15$ (in our system of units) exists, which depends only mildly on BH spin, but somewhat on the disk scale height (with taller disks being able to hold more magnetic flux; Tchekhovskoy et al. 2012). Once the magnetic flux reaches $\phi_{max}$, accretion is brought to a near-stop by the accumulation of magnetic field near the BH (Tchekhovskoy et al. 2011; McKinney et al. 2012), leading to a fundamentally different dynamic of the accretion flow and maximal energy extraction via the Blandford & Znajek (1977) process. This state is commonly referred to as Magnetically Arrested Disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1976; Narayan et al. 2003) to contrast with the Standard and Normal Evolution (SANE), where accretion is largely unaffected by the BH magnetosphere (here $\phi \sim$ few). While the MAD case is certainly of great scientific interest, in this initial code comparison, we focus on the SANE case for two reasons: (i) the SANE case is already extensively discussed in the literature and hence provides the natural starting point, (ii) the MAD dynamics poses additional numerical challenges (and remedies), which render it ill suited to establish a baseline agreement of GMRHD accretion simulations.

### 3. Code Descriptions

In this section, we give a brief, alphabetically ordered overview of the codes participating in this study, with notes on development history and target applications. Links to public release versions are provided, if applicable.

#### 3.1. Athena++

Athena++ is a general-purpose, finite-volume astrophysical fluid dynamics framework, based on a complete rewrite of Athena (Stone et al. 2008). It allows for adaptive mesh refinement in numerous coordinate systems, with additional physics added in a modular fashion. It evolves magnetic fields via the staggered mesh constrained transport algorithm of
Gardiner & Stone (2005), based on the ideas of Evans & Hawley (1988), exactly maintaining the divergence-free constraint. The code can use a number of different time integration and spatial reconstruction algorithms. Athena++ can run GRMHD simulations in arbitrary stationary spacetimes using a number of different Riemann solvers (White et al. 2016). Code verification is described in White et al. (2016), and a public release version can be obtained from https://github.com/PrincetonUniversity/athena-public-version.

### 3.2. BHAC

The BlackHoleAccretionCode (BHAC) first presented by Porth et al. (2017) is a multidimensional GRMHD module for the MPI-AMRVAC framework (Keppens et al. 2012; Porth et al. 2014; Xia et al. 2018). BHAC has been designed to solve the equations of GRMHD in arbitrary spacetimes/coordinates and exploits adaptive mesh refinement techniques with an oct-tree block-based approach. The algorithm is centered on second-order finite-volume methods, and various schemes for the treatment of the magnetic field update have been implemented on ordinary and staggered grids. More details on the various $\nabla \cdot \mathbf{B}$ preserving schemes and their implementation in BHAC can be found in Olivares et al. (2018, 2019). Originally designed to study BH accretion in ideal MHD, BHAC has been extended to incorporate nuclear equations of state, neutrino leakage, charged and purely geodetic test particles (Bacchini et al. 2018, 2019), and nonblack-hole fully numerical metrics. In addition, a nonideal resistive GRMHD module has been developed (e.g., Ripperda et al. 2019a, 2019b). Code verification is described in Porth et al. (2017) and a public release version can be obtained from https://bhac.science.

### 3.3. Cosmos++

Cosmos++ (Anninos et al. 2005; Fragile et al. 2012, 2014) is a parallel, multidimensional, fully covariant, modern object-oriented (C++) radiation hydrodynamics and MHD code for both Newtonian and general relativistic astrophysical and cosmological applications. Cosmos++ utilizes unstructured meshes with adaptive (h-) refinement (Anninos et al. 2005), moving-mesh (r-refinement; Anninos et al. 2012), and adaptive order (p-refinement; Anninos et al. 2017) capabilities, enabling it to evolve fluid systems over a wide range of spatial scales with targeted precision. It includes numerous hydrodynamics solvers (conservative and nonconservative), magnetic fields (ideal and nonideal), radiative cooling and transport, geodesic transport, generic tracer fields, and full Navier-Stokes viscosity (Fragile et al. 2018). For this work, we utilize the High Resolution Shock-Capturing scheme with staggered magnetic fields and Constrained Transport as described in Fragile et al. (2012). Code verification is described in Anninos et al. (2005).

### 3.4. ECHO

The origin of the Eulerian Conservative High-Order (ECHO) code dates back to the year 2000 (Londrillo & Del Zanna 2000, 2004), when it was first proposed as a shock-capturing scheme for classical MHD based on high-order, finite-difference reconstruction routines, one-wave or two-wave Riemann solvers, and a rigorous enforcement of the solenoidal constraint for staggered electromagnetic field components (the Upwind Constraint Transport, UCT). The GRMHD version of ECHO used in the present paper is described in Del Zanna et al. (2007) and preserves the same basic characteristics. Important

extensions of the code were later presented for dynamical spacetimes (Bucciantini & Del Zanna 2011) and nonideal Ohm equations (Bucciantini & Del Zanna 2013; Del Zanna et al. 2016; Del Zanna & Bucciantini 2018). Specific recipes for the simulation of accretion tori around Kerr black holes can be found in Bugli et al. (2014, 2018). Further references and applications may be found at www.astro.unifi.it/echo. Code verification is described in Del Zanna et al. (2007).

### 3.5. H-AMR

H-AMR is a 3D GRMHD code that builds upon HARM (Gammie et al. 2003; Noble et al. 2006) and the public code HARM-PI (https://github.com/atchekho/harmpi), and has been extensively rewritten to increase the code's speed and add new features (Liska et al. 2018a; Chatterjee et al. 2019). H-AMR makes use of GPU acceleration in a natively developed hybrid CUDA-OpenMP-MPI framework with adaptive mesh refinement (AMR) and locally adaptive time-stepping (LAT) capability. LAT is superior to the “standard” hierarchical time-stepping approach implemented in other AMR codes because the spatial and temporal refinement levels are decoupled, giving much greater speedups on logarithmically spaced spherical grids. These advancements bring GRMHD simulations with hereto unachieved grid resolutions for durations exceeding $10^5 M$ within the realm of possibility.

### 3.6. HARM-Noble

The HARM-Noble code (Gammie et al. 2003; Noble et al. 2006, 2009) is a flux-conservative, high-resolution shock-capturing GRMHD code that originated from the 2D GRMHD code called HARM (Gammie et al. 2003; Noble et al. 2006) and the 3D version (C. F. Gammie 2006, private communication) of the 2D Newtonian MHD code called HARM (Guan & Gammie 2008; Gammie & Guan 2012). Because of its shared history, HARM-Noble is very similar to the iharm3D code.<fnref n="119" /> Numerous features and changes were made from these original sources, though. Some additions include piecewise parabolic interpolation for the reconstruction of primitive variables at cell faces and the electric field at the cell edges for the constrained transport scheme, and new schemes for ensuring that a physical set of primitive variables is always recovered. HARM-Noble was also written to be agnostic to coordinate and spacetime choices, making it in a sense generally covariant. This feature was most extensively demonstrated when dynamically warped coordinate systems were implemented (Zilhão & Noble 2014) and time-dependent spacetimes were incorporated (e.g., Noble et al. 2012; Bowen et al. 2018).

### 3.7. iharm3D

The iharm3D code (Gammie et al. 2003; Noble et al. 2006, 2009; also J. Dolence 2019, private communication) is a conservative, 3D GRMHD code. The equations are solved on a logically Cartesian grid in arbitrary coordinates. Variables are zone centered, and the divergence constraint is enforced using the FluxCT constrained transport algorithm of Tóth (2000). Time integration uses a second-order predictor-corrector step. Several spatial reconstruction options are available, although

<footnote n="119">¹¹⁹ We use the name HARM-Noble to differentiate this code from the other HARM-related codes referenced herein. However, HARM-Noble is more commonly referred to as HARM3d in papers using the code.</footnote>
linear and WENO5 algorithms are preferred. Fluxes are evaluated using the local Lax-Friedrichs (LLF) method (Rusanov 1961). Parallelization is achieved with a hybrid MPI/OpenMP domain decomposition scheme. iharm3D has demonstrated convergence at second order on a suite of problems in Minkowski and Kerr spacetimes (Gammie et al. 2003). Code verification is described in Gammie et al. (2003) and a public release 2D version of the code (which differs from that used here) can be obtained from http://horizon.astro.illinois.edu/codes.

### 3.8. IllinoisGRMHD

The IllinoisGRMHD code (Etienne et al. 2015) is an open-source, vector potential-based Cartesian AMR code in the Einstein Toolkit (Löffler et al. 2012), used primarily for dynamical spacetime GRMHD simulations. For the simulation presented here, spacetime and grid dynamics are disabled. IllinoisGRMHD exists as a complete rewrite of (yet agrees to roundoff precision with) the long-standing GRMHD code (Duez et al. 2005; Etienne et al. 2010, 2012b) developed by the Illinois Numerical Relativity group to model a large variety of dynamical spacetime GRMHD phenomena (see, e.g., Etienne et al. 2006; Paschalidis et al. 2011, 2012, 2015; Gold et al. 2014 for a representative sampling). Code verification is described in Etienne et al. (2015) and a public release can be obtained from https://math.wvu.edu/~zetienne/ILGRMHD/.

### 3.9. KORAL

KORAL (Sądowski et al. 2013, 2014) is a multidimensional GRMHD code that closely follows, in its treatment of the MHD conservation equations, the methods used in the iharm3D code (Gammie et al. 2003; Noble et al. 2006; see description above). KORAL can be run with various first-order reconstruction schemes (Minmod, Monotonized Central, Superbee) or with the higher order piecewise parabolic method (PPM) scheme. Fluxes can be computed using either the Local Lax Friedrich (LLF) or the Harten Lax van Leer (HLL) method. There is an option to include an artificial magnetic dynamo term in the induction equation (Sądowski et al. 2015), which is helpful for running 2D axisymmetric simulations for long durations (not possible without this term because the MRI dies away in 2D).

Although KORAL is suitable for pure GRMHD simulations such as the ones discussed in this paper, the code was developed with the goal of simulating general relativistic flows with radiation (Sądowski et al. 2013, 2014) and multispecies fluid. Radiation is handled as a separate fluid component via a moment formalism using M1 closure (Levermore 1984). A radiative viscosity term is included (Sądowski et al. 2015) to mitigate “radiation shocks” which can sometimes occur with M1 in optically thin regions, especially close to the symmetry axis. Radiative transfer includes continuum opacity from synchrotron free-free and atomic bound-free processes, as well as Comptonization (Sądowski & Narayan 2015; Sądowski et al. 2017). In addition to radiation density and flux (which are the radiation moments considered in the M1 scheme), the code also separately evolves the photon number density, thereby retaining some information on the effective temperature of the radiation field. Apart from radiation, KORAL can handle two-temperature plasmas, with separate evolution equations (thermodynamics, heating, cooling, energy transfer) for the two particle species (Sądowski et al. 2017), and can also

evolve an isotropic population of nonthermal electrons (Chael et al. 2017). Code verification is described in Sądowski et al. (2014).

## 4. Setup Description

As the initial condition for our 3D GRMHD simulations, we consider a torus in hydrodynamic equilibrium threaded by a single weak ($\beta \gg 1$) poloidal magnetic field loop. The particular equilibrium torus solution with constant angular momentum $l := u_\phi u^\prime$ was first presented by Fishbone & Moncrief (1976) and Kozlowski et al. (1978) and is now a standard test for GRMHD simulations (see, e.g., Font & Daigne 2002; Gammie et al. 2003; Zanotti et al. 2003; Antón et al. 2006; Rezzolla & Zanotti 2013; White et al. 2016; Porth et al. 2017). Note that there exist two possible choices for the constant angular momentum, the alternative being $-u_\phi u^\prime$, which was used, e.g., by Kozlowski et al. (1978) throughout most of their work. For ease of use, the coordinates $(r, \theta, \phi)$ noted in the following are standard spherical Kerr-Schild coordinates; however, each code might employ different coordinates internally. To facilitate cross-comparison, we set the initial conditions in the torus close to those adopted by the more recent works of Shiokawa et al. (2012) and White et al. (2016). Hence, the spacetime is a Kerr BH with dimensionless spin parameter $a = 0.9375$. The inner radius of the torus is set to $r_{\text{in}} = 6 M$ and the density maximum is located at $r_{\text{max}} = 12 M$. We adopt an ideal gas equation of state with an adiabatic index of $\hat{\gamma} = 4/3$. A weak single magnetic field loop defined by the vector potential

$$A_{\phi} \propto \max(\rho/\rho_{\text{max}} - 0.2, 0) \quad (8)$$

is added to the stationary solution. The field strength is set such that $2p_{\text{max}}/(B^2)_{\text{max}} = 100$, where the global maxima of pressure $p_{\text{max}}$ and magnetic field strength $B_{\text{max}}^2$ do not necessarily coincide. With this choice of initial magnetic field geometry and strength, the simulations are anticipated to progress according to the SANE regime, although this can only be verified a posteriori.

In order to excite the MRI inside the torus, the thermal pressure is perturbed by white noise of amplitude 4%. More precisely,

$$p^* = p(1 + X_p), \quad (9)$$

and $X_p$ is a uniformly distributed random variable between -0.02 and 0.02.

To avoid density and pressures dropping to zero in the funnel region, floor models are customarily employed in fluid codes. Because the strategies differ significantly between implementations, only a rough guideline to the floor model was given. The following floor model was suggested: $\rho_{\text{fl}} = 10^{-5} r^{-3/2}$ and $p_{\text{fl}} = 1/3 \times 10^{-7} r^{-5/2}$, which corresponds to the one used by McKinney & Gammie (2004) in the power-law indices. Thus, for all cells that satisfy $\rho \le \rho_{\text{fl}}$, set $\rho = \rho_{\text{fl}}$; in addition, if $p \le p_{\text{fl}}$, set $p = p_{\text{fl}}$. It is well known that occasionally, unphysical cells are encountered with, e.g., negative pressures and high Lorentz factor in the funnel. For example, it can be beneficial to enforce that the Lorentz factor stay within reasonable bounds. This delimits the momentum and energy of faulty cells and thus aids in keeping the error localized. The various fail-safes and checks of each code are described in more detail in Section 4.1. The implications of the different choices will be discussed in Sections 5.4 and 6.
In terms of coordinates and gridding, we deliberately gave only loose guidelines. The reasoning is twofold: first, this way, the results can inform us about the typical setup used with a particular code, thus allowing us to get a feeling for how existing results compare. The second reason is purely utilitarian, as settling for a common grid setup would incur extra work and likely introduce unnecessary bugs. For spherical setups, which are the majority of the participants, a form of logarithmically stretched Kerr–Schild coordinates with optional warping/cylindrification in the polar region was hence suggested.

Similarly, the positioning and nature of the boundary conditions has been left free for each group with only the guideline to capture the domain of interest $r \in [r_h, 50]$, $\theta \in [0, \pi]$, $\phi \in [0, 2\pi]$. The implications of the different choices will be discussed in Sections 5.3 and 6. Three rounds of resolutions are suggested in order to judge the convergence between codes. These are low-res: 96³, mid-res: 128³, and high-res: 192³, where the resolution corresponds to the domain of interest mentioned above.

To make sure the initial data is set up correctly in all codes, a stand-alone Fortran 90 program was supplied and all participants have provided radial cuts in the equatorial region. This has proven to be a very effective way to validate the initial configuration.

An overview of the algorithms employed for the various codes can be found in Table 1. Here, the resolutions ($\Delta r_p$, $\Delta\theta_p$, $\Delta\phi_p$) refer to the proper distance between grid cells at the density maximum in the equatorial plane for the low-resolution realization (typically 96 × 96 × 96). Specifically, we define

$$ \Delta r_p := \sqrt{g_{rr}(12M, \pi/2)} \Delta r, \quad (10) $$

and analog for the other directions. For the two Cartesian runs, we report the proper grid spacings at the same position (in the xz-plane) for the x-, y-, and z-directions, respectively, and treat $\Delta z_p$ as representative of the out-of-plane resolution $\Delta\theta_p$ in the following sections.

## 4.1. Code-specific Choices

### 4.1.1. Athena++

For these simulations, Athena++ uses the second-order van Leer integrator (van Leer 2006) with third-order PPM reconstruction (Colella & Woodward 1984). Magnetic fields are evolved on a staggered mesh as described in Gardiner & Stone (2005) and generalized to general relativity in White et al. (2016). The two-wave HLL approximate Riemann solver is used to calculate fluxes (Harten et al. 1983). The coordinate singularity at the poles is treated by communicating data across the pole in order to apply reconstruction and by using the magnetic fluxes at the pole to help construct edge-centered electric fields in order to properly update magnetic fields near the poles. Mass and/or internal energy is added in order to ensure $\rho > 10^{-5}r^{-3/2}$, $p > 1/3 \times 10^{-7}r^{-5/2}$, $\sigma < 100$, and $\beta > 10^{-3}$. Additionally, the normal-frame Lorentz factor is kept under 50 by reducing the velocity if necessary.

All Athena++ simulations are done in Kerr–Schild coordinates. The grids are logarithmically spaced in $r$ and uniformly spaced in $\theta$ and $\phi$. They use the fiducial resolution but then employ static mesh refinement to derefine in all coordinates toward the poles, stepping by a factor of 2 each time. The 96³

grid achieves the fiducial resolution for $\pi/4 < \theta < 3\pi/4$ at all radii and for $\pi/12 < \theta < 11\pi/12$ when $r > 19.48$, derefining twice; the 128³ grid achieves this resolution for $3\pi/16 < \theta < 13\pi/16$ at all radii and for $\pi/16 < \theta < 15\pi/16$ when $r > 19.68$, derefining twice; and the 192³ grid achieves this resolution for $7\pi/24 < \theta < 17\pi/24$ when $r > 1.846$ and for $\pi/8 < \theta < 7\pi/8$ when $r > 31.21$, derefining three times. The outer boundary is always at $r = 50$, where the material is kept at the background initial conditions. The inner boundaries are at radii of 1.152, 1.2, and 1.152, respectively, ensuring that exactly one full cell at the lowest refinement level is inside the horizon. Here, the material is allowed to freely flow into the BH, with the velocity zeroed if it becomes positive.

### 4.1.2. BHAC

In BHAC, we employ the LLF fluxes in combination with PPM reconstruction (Colella & Woodward 1984) and a two-step predictor-corrector scheme. The setup analyzed here was run with the staggered upwind constrained transport (UCT) scheme of Del Zanna et al. (2007). The simulations are performed in modified Kerr–Schild coordinates (McKinney & Gammie 2004) with $\theta$-coordinate stretching parameter $h = 0.75$. In the staggered case, two to three grid levels are utilized (three for the high-resolution run) with static mesh refinement chosen such that the polar axis at small radii is fully derefined and the torus is fully resolved on the highest grid level. This allows us to significantly increase the time step, which is otherwise dominated by the cells close to the singular axis. Hence, compared to the brute-force uniform grid setup, the time step in the three-level run is increased by a factor of 16. We adopt a floor model in the Eulerian frame as suggested; however, because staggered fields need to be interpolated to the cell centers, which introduces an additional error, we have to increase the floors to $\rho_{fl} = 10^{-4}r^{-3/2}$ and $p_{fl} = 1/3 \times 10^{-6} r^{-5/2}$. Floors that are lower by an order of magnitude were no problem for centered field FluxCT runs (Tóth 2000).

On the singular axis, we explicitly zero out the fluxes as well as the $\phi$- and $r$-components of the electric fields. Furthermore, we employ $\pi$-periodic boundary conditions and hence fill the axial ghost cells with their diagonal counterpart. No further pole fix was required for the stable evolution of the setup.

To increase the variety of the setups considered in the comparison and to introduce a case for which the difficulties related to the polar axis are not present, we also carry out a simulation in Cartesian Kerr–Schild (CKS) coordinates. For this run (referred to as BHAC Cart. in the following), we use a combination of adaptive and fixed mesh refinement in an attempt to simultaneously resolve the MRI turbulence within the disk and follow the propagation of the jet. Adaptive mesh refinement is triggered by variations in the plasma magnetization $\sigma$ and the density, quantified by means of the Löhner scheme (Löhner 1987). We structure the domain as a set of nested boxes such that the highest AMR level achievable for each box increases inwards, and the highest level in the simulation can be achieved only by the innermost box, containing the event horizon. The simulation employs eight such levels and a base resolution of $N_x \times N_y \times N_z = 96 \times 96 \times 192$, and extends over $x, y \in [-500 M, 500 M]$ and $z \in [-1000 M, 1000 M]$. In order to prevent an unphysical outflow from the BH interior, we apply cutoff values for the density and pressure in the region $r < 0.5(r_{H-} + r_{H+})$, where $r_{H-}$ and $r_{H+}$ are the locations of the inner and outer event horizons. Specifically, we set $\rho_{cut} = 10^{-2}$ or $p_{cut} = 10^{-4}$. Other
<table><caption>**Table 1**

Algorithmic Details of the Code Comparison Runs for the Low-resolution Realizations</caption><thead><tr><th>Code</th><th>Time Stepper</th><th>Riemann s.</th><th>Rec.</th><th>Mag. Field</th><th>(Δr<sub>p</sub>, Δθ<sub>p</sub>, Δφ<sub>p</sub>)</th><th>Domain: r</th><th>Domain: θ</th></tr></thead><tbody><tr><td>Athena++</td><td>Second Order</td><td>HLL</td><td>PPM</td><td>CT (Gardiner & Stone 2005)</td><td>(0.506, 0.393, 0.788)</td><td>[1.152, 50]</td><td>[0, π]</td></tr><tr><td>BHAC</td><td>Second Order</td><td>LLF</td><td>PPM</td><td>UCT (Del Zanna et al. 2007)</td><td>(0.402, 0.295, 0.788)</td><td>[1.185, 3333]</td><td>[0, π]</td></tr><tr><td>BHAC Cart.</td><td>Second Order</td><td>LLF</td><td>PPM</td><td>UCT (Del Zanna et al. 2007)</td><td>Cartesian AMR (0.176, 0.163, 0.163)</td><td>See Section 4.1.2.</td><td>...</td></tr><tr><td>Cosmos++</td><td>SSPRK, third Order</td><td>HLL</td><td>PPM</td><td>(Fragile et al. 2012)</td><td>(0.508, 0.375, 0.788)</td><td>[1.25, 354]</td><td>[0.070, 3.072]</td></tr><tr><td>ECHO</td><td>Third Order</td><td>HLL</td><td>PPM</td><td>UCT (Del Zanna et al. 2007)</td><td>(0.460, 0.382, 0.752)</td><td>[1.048, 2500]</td><td>[0.060, 3.082]</td></tr><tr><td>H-AMR</td><td>Second Order</td><td>HLL</td><td>PPM</td><td>UCT (Gardiner & Stone 2005)</td><td>(0.523, 0.379, 0.785)</td><td>[1.169, 500]</td><td>[0, π]</td></tr><tr><td>HARM-Noble</td><td>Second Order</td><td>LLF</td><td>PPM</td><td>PPM+FluxCT (Tóth 2000), (Noble et al. 2009)</td><td>(0.578, 0.182, 0.784)</td><td>[1.090, 80]</td><td>[0.053, 3.088]</td></tr><tr><td>iharm3D</td><td>Second Order</td><td>LLF</td><td>PLM</td><td>FluxCT (Tóth 2000)</td><td>(0.519, 0.118, 0.788)</td><td>[1.073, 50]</td><td>[0, π]</td></tr><tr><td>IllinoisGRMHD</td><td>RK4</td><td>HLL</td><td>PPM</td><td>Vector potential-based PPM+FluxCT (Tóth 2000), (Etienne et al. 2010), (Etienne et al. 2012b)</td><td>Cartesian AMR (0.246, 0.228, 0.229)</td><td>See Section 4.1.8.</td><td>...</td></tr><tr><td>KORAL</td><td>Second order</td><td>LLF</td><td>PPM</td><td>FluxCT (Tóth 2000)</td><td>(0.478, 0.266, 0.785)</td><td>[1.15, 50]</td><td>[0, π]</td></tr></tbody></table>

**Note.** The columns are code name, order of the time integration, approximate Riemann solver, spatial reconstruction scheme, scheme for magnetic field evolution, proper distance between grid cells at the midplane, and radial and meridional extents of the computational domain. The azimuthal direction spans [0, 2π] in all cases.
than that, the evolution in the interior of the event horizon is followed with the same algorithm as in the exterior; in particular, the magnetic field update is obtained by the UCT algorithm providing zero divergence of the magnetic field throughout.

For all simulations, to more accurately treat the highly magnetized polar region, we employ the entropy strategy discussed in Porth et al. (2017), that is, each time the plasma β drops below the threshold value of 10⁻², the advected entropy is used for primitive variable recovery instead of the conserved energy.

### 4.1.3. Cosmos++

In Cosmos++, we use the HLL Riemann solvers with PPM reconstruction (Colella & Woodward 1984) and a five-stage, strong-stability-preserving Runge-Kutta time-integration scheme (Spiteri & Ruuth 2002), set to third order. The magnetic fields were evolved using the staggered constrained transport scheme described in Fragile et al. (2012). A uniform θ-coordinate was used with small cutouts near the poles (θ < 4°) to keep the time step reasonable. The outer radial boundary was placed at (50)¹.⁵M to reduce boundary effects. We then increased the number of zones in the radial direction by N_r¹.⁵ to maintain the desired resolution in the region of interest. Outflow boundary conditions (copying scalar fields to ghost zones, while ensuring the velocity component normal to the boundary points outward) were used on the radial and polar boundaries. For the primitive inversion step, we primarily used the “2D” option from Noble et al. (2006), with a 5D numerical inversion scheme (similar to the 9D inversion described in Fragile et al. 2014) as a backup. In cases where both solvers failed, we ignored the conserved energy and instead used the entropy to recover the primitive variables. Otherwise, the default suggestions as laid out in Section 4 were used.

### 4.1.4. ECHO

The time integration performed in ECHO uses the third-order accurate IMplicit-EXplicit (IMEX) strong-stability-preserving scheme (Pareschi & Russo 2005), which in the case of ideal GRMHD reduces effectively to a third-order Runge-Kutta. The upwind fluxes are computed with the HLL Riemann solver with PPM reconstruction (Colella & Woodward 1984), using the UCT scheme of Del Zanna et al. (2007) for the treatment of the magnetic field.

The numerical grid is logarithmically stretched in radius and uniform in both θ and φ, excluding from the polar angle domain the regions close to the rotation axis to avoid excessively small time steps. At the radial and polar boundaries, we impose outflow boundary conditions, i.e., we copy the value of the primitive variables and set the velocity normal to the boundary to zero whenever it has a positive (negative) value at the inner (outer) boundary.

As suggested, we adopt the floor model for the rest-mass density and pressure used by McKinney & Gammie (2004). For the primitive variable recovery, we first use the three-dimensional scheme described in Bucciantini & Del Zanna (2013), and in the case of failure, we apply the 1D scheme from Porth et al. (2017). Should none of these routines be successful, we then attempt to retrieve the primitives using the advected entropy instead of the total energy. In case of persisting failures, we finally reset the value of density and pressure to the

atmospheric floor values and set the Eulerian velocity to zero, without modifying the magnetic field.

### 4.1.5. H-AMR

Like HARM (Gammie et al. 2003), H-AMR evolves the GRMHD equations in arbitrary (fixed) spacetimes. H-AMR is third-order accurate in space by using the PPM (Colella & Woodward 1984) reconstruction of primitive variables at cell faces, and second-order accurate in time. The fluxes at cell faces are calculated using an approximate HLL Riemann solver, while the magnetic field is evolved on a staggered grid, where the electric fields have been velocity upwinded to add dissipation (Gardiner & Stone 2005). Because the funnel is devoid of matter, H-AMR artificially inserts mass in the drift frame (Ressler et al. 2017). This does not lead to a runaway in velocity, which occurs when mass is inserted in the fluid frame (Gammie et al. 2003), or to artificial drag on the field lines, which occurs when mass is inserted in the zero angular momentum observer (ZAMO) frame (McKinney et al. 2012). We enforce floor values on the rest-mass density ρ<sub>fl</sub> = MAX[B²/20, 10⁻⁵(r/M)⁻³/²] and internal energy u<sub>fl</sub> = MAX[B²/750, 3.33 × 10⁻⁸(r/M)⁻⁵/²]. To provide a backup inversion method for primitive variable recovery if all other primary inversion method(s) fail (Noble et al. 2006; Hamlin & Newman 2013), H-AMR also advects the conserved entropy (Noble et al. 2009). We typically use three to four levels of local adaptive time stepping to speed up the code by an additional factor of 3–5 to 3–5 × 10⁸ zone cycles/s/GPU (Chatterjee et al. 2019).

We use a (close to) uniformly spaced logarithmic spherical grid combined with outflow boundary conditions in the radial direction, transmissive boundary conditions in the θ-direction, and periodic boundary conditions in the φ-direction. Because cells get squeezed near the pole, the time step in all spherical grids is reduced by an additional factor proportional to the resolution in the φ-direction. To remedy this issue, we stretch out cells immediately adjacent to the pole in the θ-direction (Tchekhovskoy et al. 2011) and use multiple levels of static mesh derefinement in the φ-direction to keep the cell’s aspect ratio close to uniform at high latitudes. As an example, the very high-resolution (effectively, 1608 × 1056 × 1024) run in this work uses a φ-resolution of 64 cells for 0° < θ < 3.75, 128 cells for 3.75 < θ < 7.5, 256 cells for 7.5 < θ < 15°, 512 cells for 15° < θ < 30°, and the full 1024 cells for 30° < θ < 90°. This method prevents the squeezing of cells near the pole from reducing the global time step, while maintaining high accuracy in all three dimensions.

### 4.1.6. HARM-Noble

The results using HARM-Noble given here used the LF approximate Riemann solver as defined in Gammie et al. (2003), RK2 time integration, the FluxCT method of Tóth (2000), and PPM reconstruction (Colella & Woodward 1984) for the primitive variables at cell faces and the electric fields at cell edges (for the sake of the FluxCT method). We used a “modified Kerr-Schild” coordinate system specified by

$$r(x^l) = \exp(x^l), \quad \theta(x^2) = \theta_c + h \sin(2\pi x^2), \quad (11)$$

with θ_c = 0.017π and h = 0.25. Zeroth-order extrapolation was used to set values in the ghost zones at the radial and poloidal boundaries; outflow conditions were enforced at the inner and outer radial boundaries. The poloidal vector
components were reflected across the poloidal boundary (e.g., $B^\theta \to -B^\theta$, $v^\theta \to -v^\theta$).

Instead of using 96 cells (for instance for the low-resolution run) within $r = 50M$ as many of the others used, the HARM-Noble runs used this number over the entire radial extent it used, out to $r = 80M$. This means that a HARM-Noble run has a lower radial resolution than another code's run with the same cell count.

Recovery of the primitive variables from the conserved variables was performed with the “2D” and “1D<sub>w</sub>” methods of Noble et al. (2006), where the latter method is used if the former one fails to converge to a sufficiently accurate physical solution.

The HARM-Noble runs also used the so-called “entropy fix” of Balsara & Spicer (1999) as described in Noble et al. (2009), wherein the entropy equation of motion (EOM) replaces the energy EOM in the primitive variable method whenever it fails to converge; the resultant primitive variables are unphysical, or $u < 10^{-2} B^2$.<fnref n="120" />

### 4.1.7. iharm3D

iharm3D is an unsplit method-of-lines scheme that is spatially and temporally second order. A predictor-corrector scheme is used for time stepping. For models presented here, spatial reconstruction was carried out with a piecewise linear method using the monotonized central slope limiter. The divergence-free condition on the magnetic field was maintained to machine precision with the FluxCT scheme of Tóth (2000).

The simulations used “funky” modified Kerr-Schild (FMKS) coordinates $t, X^1, X^2, X^3$ that are similar to the MKS coordinates of McKinney & Gammie (2004), with $R_0 = 0$ and $h = 0.3$, but with an additional modification to MKS to enlarge grid zones near the polar axis at small radii and increase the time step. For FMKS, then, $X^2$ is chosen to smoothly interpolate from KS $\theta$ to an expression that deresolves the poles:

$$ \theta = N y \left( 1 + \frac{(y/x_t)^\alpha}{\alpha + 1} \right) + \frac{\pi}{2}, \quad (12) $$

where $N$ is a normalization factor, $y \equiv 2X^2 - 1$, and we choose $x_t = 0.82$ and $\alpha = 14$.

iharm3d imposes floors on rest-mass density $\rho$ and internal energy $u$ to enforce $\rho > 10^{-5}r^{-2}$ and $u > 10^{-7}r^{-2}\hat{\gamma}$, where $\hat{\gamma}$ is the adiabatic index. It also requires that $\rho > B^2/50$, $u > B^2/2500$, and subsequently that $\rho > u/50$ (Ressler et al. 2017). At high magnetizations, mass and energy created by the floors are injected in the drift frame (Ressler et al. 2017); otherwise, they are injected in the normal observer frame. The fluid velocity primitive variables are rescaled until the Lorentz factor in the normal observer frame is <50. If inversion from conserved to primitive variables fails, the primitive variables are set to an average over a stencil of neighboring cells in which the inversion has not failed.

The radial boundaries use outflow-only conditions. The polar axis boundaries are purely reflecting, and the $X^1$ and $X^3$ fluxes

of the $X^2$ component of the magnetic field are antiparallel across the boundary. The $X^3$ boundaries are periodic.

### 4.1.8. IllinoisGRMHD

IllinoisGRMHD simulations presented here adopt a Cartesian FMR grid (using the Cactus/Carpet infrastructure), in which four overlapping cubes with half-side length 27.34375M are centered on the $x-y$ plane at positions $(x/M, y/M) = (25, 25)$; $(25, -25)$; $(-25, 25)$; $(-25, -25)$. These cubes, all at resolution $\Delta x = \Delta y = \Delta z \approx M/4.388571$, constitute the highest refinement level, and a total of six factor-of-2 derefinements (with approximate resolutions $M/2.19$, $M/1.10$, $1.82M$, $3.65M$, $7.29M$, and $14.583M$ [exact]) are placed outside this high-resolution level so that the outer boundary is a cube with half-side length $1750M$, centered at $(x, y, z) \approx (M/8.78, M/8.78, M/8.78)$. This ensures full cell-centering on the finest refinement level, which maximally avoids the $z$-axis and $r=0$ coordinate singularities when mapping initial data to the Cartesian grids.

IllinoisGRMHD adopts the same formalism as iharm3D (Gammie et al. 2003) for the GRMHD field equations, with the exception of the magnetic field evolution; IllinoisGRMHD evolves the vector potential directly (Etienne et al. 2010, 2012b). Evolving the vector potential enables any interpolation scheme to be used at AMR refinement boundaries, and the formulation IllinoisGRMHD adopts reduces to the standard staggered FluxCT scheme on uniform-resolution grids. As for GRMHD field reconstruction, IllinoisGRMHD makes use of the PPM algorithm and HLL for its approximate Riemann solver. The conservative-to-primitives solver in IllinoisGRMHD is based on the open-source Noble et al. code (Noble et al. 2006), but has been extended to adjust conservative variables that violate physicality inequalities prior to the solve (Etienne et al. 2012a).

### 4.1.9. KORAL

The simulations presented here used PPM reconstruction and Lax-Friedrichs fluxes. Modified Kerr-Schild coordinates were employed, using the technique developed in Tchekhovskoy et al. (2011), whereby the $\theta$-grid was concentrated modestly toward the equator and was moved away from the poles at small radii to avoid very small time steps. The following floors and ceilings were applied for numerical stability (they mostly affect the polar low-density regions, which are not the focus of the present paper): $10^{-8} \le u/\rho c^2 \le 10^2$, $B^2/u \le 10^5$, $B^2/\rho c^2 \le 50$, $\Gamma \le 20$. Outflow boundary conditions were used at the outer radius and reflecting boundary conditions at the poles.

## 5. Results

A consistent set of diagnostics focusing both on horizon-scale quantities and on the global evolution of the accretion flow is performed with all codes. For ease of implementation, all diagnostics are performed in the standard Kerr-Schild coordinates. We now first describe the quantifications and then move on to the intercode comparison.

### 5.1. Diagnostics

1. Horizon-penetrating fluxes. The relevant quantities: mass, magnetic, angular momentum, and (reduced) energy fluxes

<footnote n="120">¹²⁰ During the production stage of publication, Noble discovered the HARM-Noble simulations errantly used the same seed across all MPI processes for the random number generator used to perturb the initial conditions. Because the azimuthal dimension was decomposed and because the other dimensions are decomposed uniformly, using the same seed meant the initial conditions were azimuthally periodic with period equal to $\Delta\phi = 2\pi/N$, where $N$ is the number of azimuthal decompositions. Due to the symmetry-preserving implementation of HARM-Noble, these modes do not appreciably grow out of the accumulated roundoff error during the course of the simulation. In other words, the HARM-Noble simulations are missing the $m = 1, \dots, (N-1)$ azimuthal modes of dynamics, where $N = (3, 6, 6)$ for runs $96^3, 128^3, 192^3$, respectively.</footnote>
are defined as follows:

$$ \dot{M} := \int_{0}^{2\pi} \int_{0}^{\pi} \rho u^{r} \sqrt{-g} \, d\theta \, d\phi, \qquad (13) $$

$$ \Phi_{BH} := \frac{1}{2} \int_{0}^{2\pi} \int_{0}^{\pi} |*F^{rt}| \sqrt{-g} \, d\theta \, d\phi, \qquad (14) $$

$$ \dot{L} := \int_{0}^{2\pi} \int_{0}^{\pi} T_{\phi}^{r} \sqrt{-g} \, d\theta \, d\phi, \qquad (15) $$

$$ \dot{E} := \int_{0}^{2\pi} \int_{0}^{\pi} (-T_{t}^{r}) \sqrt{-g} \, d\theta \, d\phi, \qquad (16) $$

where all quantities are evaluated at the outer horizon $r_h$. A cadence of 1M or less is chosen. In practice, these quantities are nondimensionalized with the accretion rate, yielding, for example, the normalized magnetic flux $\phi = \Phi_{BH}/\sqrt{\dot{M}}$, also known as the “MAD parameter,” which, for spin $a = 0.9375$ and torus scale height $H/R \approx 0.3$, has the critical value $\phi_{max} \approx 15$ (within the units adopted here; Tchekhovskoy et al. 2012).

2. Disk-averaged quantities. We compare averages of the basic flow variables $q \in \{\rho, p, u^{\phi}, \sqrt{b_{\alpha}b^{\alpha}}, \beta^{-1}\}$. These are defined similarly to Shiokawa et al. (2012), White et al. (2016), and Beckwith et al. (2008). Hence, for a quantity $q(r, \theta, \phi, t)$, the shell average is defined as

$$ \langle q \rangle (r, t) := \frac{\int_{0}^{2\pi} \int_{\theta_{\min}}^{\theta_{\max}} q(r, \theta, \phi, t) \sqrt{-g} \, d\phi \, d\theta}{\int_{0}^{2\pi} \int_{\theta_{\min}}^{\theta_{\max}} \sqrt{-g} \, d\phi \, d\theta}, \qquad (17) $$

which is then further averaged over the time interval $t_{KS} \in [5000, 10,000]M$ to yield $\langle q \rangle (r)$ (note that we omit the weighting with the density as done by Shiokawa et al. 2012; White et al. 2016). The limits $\theta_{min} = \pi/3$, $\theta_{max} = 2\pi/3$ ensure that only material from the disk is taken into account in the averaging.

3. Emission proxy. To get a feeling for the code-to-code variations in synthetic, optically thin light curves, we also integrate the pseudo-emissivity for thermal synchrotron radiation following an appropriate nonlinear combination of flow variables:

$$ \mathcal{L}(t) := \int_{0}^{2\pi} \int_{\theta_{\min}}^{\theta_{\max}} \int_{r_{h}}^{r_{\max}} j(B, p, \rho) \sqrt{-g} \, d\phi \, d\theta \, dr, \qquad (18) $$

where again $\theta_{min} = \pi/3$, $\theta_{max} = 2\pi/3$, and $r_{max} = 50M$ are used. The emissivity $j$ is here defined as follows: $j = \rho^3 p^{-2} \exp(-C(\rho^2/(Bp^2))^{1/3})$, which captures the high-frequency limit of the thermal synchrotron emissivity $\nu \gg \nu_c \Theta_e^2 \sin(\theta)$ (compare with Equation (57) of Leung et al. 2011). The constant $C$ is chosen such that the radiation is cut off after a few gravitational radii, resembling the expected millimeter emission from the Galactic Center, that is, $C = 0.2$ in geometrized units. This emission proxy is optimized for the science case of the EHT where near optically thin synchrotron emission is expected.<fnref n="121" /> In order to compare the variability, again a cadence of 1M or less is chosen in most cases.

4. t, φ averages. Finally, we compare temporally and azimuthally averaged data for a more global impression

of the disk and jet system:

$$ \langle q \rangle (r, \theta) := \frac{\int_{t_{\text{beg}}}^{t_{\text{end}}} \int_{0}^{2\pi} q(r, \theta, \phi, t) \, d\phi \, dt}{2\pi(t_{\text{end}} - t_{\text{beg}})} \qquad (19) $$

for the quantities $q \in \{\rho, \beta^{-1}, \sigma\}$ with the averaging interval ranging from 5000 M to 10,000 M.

## 5.2. Time Series

Time series data of horizon-penetrating fluxes is presented in Figures 2–4 for low, medium, and high resolutions respectively. Because the mass accretion governs the behavior of these fluxes, the data have been appropriately normalized by the accretion rate. All codes capture accurately the linear phase of the MRI leading to an onset of accretion to the BH at $t \simeq 300M$. While there is still a good correspondence of $\dot{M}$ for early times $<1000M$, the chaotic nature of the problem fully asserts itself after 2000M. At low resolution, there exist order-of-magnitude variations in the data, most notably for the normalized horizon-penetrating magnetic fluxes $\phi$ and the energy fluxes. The low value of $\phi$ for the Cosmos++ data is caused by the choice of boundary conditions near the polar region as will be discussed in Section 5.4 in more detail. As indicated by $\phi \sim 0.5-6 < \phi_{max}$, all simulations are in the SANE regime of radiatively inefficient BH accretion. As a measure for the variance between codes, in Table 2 we quantify the peak fluxes as well as the average values far in the nonlinear regime.

It is worthwhile to point out the good agreement in angular momentum and energy fluxes at high resolution for the codes employing spherical grids, where the highest average values differ from the lowest ones by 12% and 40%.

In Figure 3, we additionally show the Cartesian realizations next to the medium-resolution cases. It should be kept in mind though that the resolution in the Cartesian cases is typically much worse near the horizon but better at larger radii, which makes it hard to directly compare the simulations. Qualitatively, there is very good agreement of the Cartesian runs with the spherical data in all quantities except for the energy fluxes where BHAC Cart. is systematically higher and IllinoisGRMHD is slightly lower than all spherical codes.

Because meshes with various amounts of compression in the midplane were employed by the different groups, a comparison purely based on the number of grid cells is, however, hardly fair even in the spherical cases. Hence, in Figure 5 we show the data of Table 2 against the proper grid spacing at the location of the initial density maximum $\Delta\theta_p$. Ordered by $\Delta\theta_p$, the trends in the data become clearer: for example, the midplane resolution of the iharm3D run at $N_\theta = 96$ is higher than the resolutions employed, e.g., in BHAC and KORAL at $N_\theta = 192$. Hence, even the lowest resolution iharm3D case is able to properly resolve the MRI at late times (see Section 5.6), yielding consistent results across all runs. Taken for itself, KORAL and HARM-Noble also show very good agreement among the three resolution instances (with the exception of horizon-penetrating magnetic flux for KORAL); however, the horizon-penetrating flux in Cosmos++ is consistently off by a factor of ~4 compared to the distribution average. The trend can be explained by noting that Cosmos++ used a polar cutout with open outflow boundary conditions. This allows magnetic flux

<footnote n="121">¹²¹ In this, the proxy is distinct from the prescription often used for optically thick, geometrically thin disks which is based on an estimation of the turbulent dissipation, (e.g., Hubeny & Hubeny 1998; Armitage & Reynolds 2003).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time series at the resolution of 96³ in the domain of interest for all codes. From top to bottom, the panels show the mass accretion rate, horizon-penetrating magnetic flux, accretion of angular momentum and energy, as well as the pseudo-luminosity according to Equation (18). The data in panels 2–5 have been nondimensionalized with the accretion rate $\dot{M}$. For reference, the H-AMR $N_{\theta} = 1056$ solution was added (thick dashed red lines).</figcaption></figure>

to effectively leak out of the domain, which further reduces the magnetization of the funnel region (see also Section 5.4).

### 5.3. Disk Profiles

The disk-averaged profiles of the relevant quantities are presented in Figures 6–8. At the late times under consideration, viscous spreading and accretion have significantly transformed the initial distributions; for example, the peak densities are now an order of magnitude below the initial state. There is a fair spread in some quantities between codes, most notably in the profiles of density and (inverse) plasma beta, while all codes capture very well the rotation law of the disk, which can be approximated by a power law with slope $r^{-1.75}$, somewhat steeper than the Keplerian case. Interestingly, the profiles of the magnetic field are also captured quite accurately and tend to agree very well with a very high-resolution reference case computed with the H-AMR code (dashed red curves). The approximate scaling of the disk magnetic field as $r^{-1}$ indicates a dominance of the toroidal field component (e.g., Hirose et al. 2004).

At low and medium resolution, the density profile in BHAC (blue curves) demonstrates an inner peak or “mini torus,” which goes away in high resolution. As the additional analysis

of the MRI quality reveals (see Section 5.6), the low- and mid-resolution BHAC runs do not properly resolve the fastest growing wavelengths $\lambda_r$, $\lambda_\theta$, and hence, we consider the stable mini torus a numerical artifact, due to the lack of angular momentum transport in the underresolved simulations.

Taking the inverse plasma $\beta$ (bottom-right panel) as a proxy for the Maxwell stresses, the decrease in torus density is consistent with an increase in the turbulent transport of angular momentum, due to the increase in magnetization. Inspection of the torus density profiles shows another trend: we obtain systematically smaller tori in setups where the outer boundary is near the torus edge (Athena++, iharm3D, KORAL). This is expected, as the outflow boundary prohibits matter from falling back from large radii as would otherwise occur, and hence, mass is effectively removed from the simulation.

Boundary effects are also visible in the run of magnetic field and pressures: again, the three small domain models have the steepest pressure profiles, with the Athena++ showing clear kinks in several quantities at the boundary and the iharm3D measurement for $\langle\beta^{-1}\rangle$ rising sharply at the boundary (due to a few high $\beta^{-1}$ cells in the equatorial region close to the boundary; see, for example, also Section 5.4 and Figure 11).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** As Figure 2 for the 128³ data and the Cartesian runs "IllinoisGRMHD" and "BHAC Cart."</figcaption></figure>

Turning to the high-resolution data, the spread in all quantities becomes dramatically reduced. Starting at a few gravitational radii, the magnetization of the disk is nearly constant, with values in the range 0.04–0.1. Here, the Cartesian runs are the exception; their magnetization decreases within $r \sim 10M$ and reaches a minimum of $\sim 0.01$. We suppose that this stems from increased numerical dissipation, due to the comparatively low resolution in the inner regions of the Cartesian grids. In fact, we checked how many resolution elements capture the fastest growing MRI mode (quality factor) for BHAC Cart. and confirm that these fall below a critical value of six for $r \sim 10M$.

Taking advantage of the scale freedom in the ideal MHD approximation, we have also performed a comparison with density, pressure, and magnetic field variables rescaled to the H-AMR $N_\theta = 1056$ solution. This is exemplified in Appendix C, which shows the expected better match for the density and pressure profiles. However, as can already be anticipated from the difference in plasma $\beta$, the variance in the magnetic field is bound to increase.

### 5.4. Axisymmetrized Data

To gain an overall impression of the solutions in the quasi-stationary state, we compare averages over the $t$- and $\phi$-coordinates

in Figures 9–12. The different panels depict the rest-frame density, inverse plasma $\beta$, and the magnetization obtained for each code at the highest resolution available.

Though qualitatively very similar between codes, the density maps portray correlations between torus size and position of the outflow boundary, as noted already in Section 5.3. Visually, the boundary is most pronounced in the Athena++ run where the low-density region is spread out over a large polar angle.

As a large variety of floor models is employed, maps of plasma $\beta$ and magnetization exhibit substantial spread in the funnel region. Aside from the difference in absolute values reached in the funnel (which are merely related to the floor level), there are also qualitative differences concerning the cells near the axis: depending on the pole treatment, the inverse plasma $\beta$ decreases toward the axis for some codes (Cosmos+, KORAL, HARM-Noble, Athena++, H-AMR) while others see an increase of $\beta^{-1}$ (BHAC, iharm3D) or a near-constant behavior (BHAC Cart.). We should note that the Cartesian runs do not require any pole treatment at all and as such give perhaps the most reliable answer in this region. However, Cartesian grids are also far from perfect: in the IllinoisGRMHD run, we observe some artifacts in the form of horizontal features in $\langle \beta^{-1} \rangle$ which are associated with resolution jumps in the computational grid.
<figure><img src="image_6.png" /><figcaption>Figure 4. As Figure 2 for the 192³ data.</figcaption></figure>

Significant leakage of magnetic flux leads to a complete loss of the magnetically dominated region in the medium-resolution *Cosmos++* run. A similar behavior is observed with the medium-resolution *HARM-Noble* case (not shown). In both cases, the polar axis was excised, which necessitates the use of boundary conditions at the polar cutout. If no special precautions are taken, the magnetic field simply leaks through the boundary, leaving less flux on the BH (see Section 5.2). The open outflow boundaries adopted in the *Cosmos++* runs are particularly prone to this effect.

Although not the focus of this work (reflected, e.g., in the choice of gridding by the various groups), it is interesting to examine the jet–disk boundary defined as $σ = 1$ (e.g., McKinney & Gammie 2004; Nakamura et al. 2018) in more detail. Figure 13 illustrates that as the resolution is increased, the contour is more faithfully recovered across codes: whereas some low-resolution runs do not capture the highly magnetized funnel, at high resolution, all runs show a clearly defined jet–disk boundary. Despite the large variances in floor treatment, at $N_θ = 192$ cells, the difference is reduced to five degrees and the polar angle of the disk–jet boundary ranges between 10° and 15° at $r = 50M$. For illustrative purposes, we also overplot flux–functions of the approximate force-free solutions

discussed in Tchekhovskoy et al. (2008) and Nakamura et al. (2018):

$$ \Psi(r, \theta) = \left(\frac{r}{r_h}\right)^{\kappa} (1 - \cos(\theta)). \qquad (20) $$

In the recent 2D simulations of Nakamura et al. (2018), the jet boundary given by $σ = 1$ was accurately described by the choice $κ = 0.75$ for a wide range of initial conditions, BH spins, and horizon-penetrating flux: $Φ_{BH}/\sqrt{M} \in [5, 10]$. This results in a field line shape $z \propto R^{1.6}$, which matches well the shape derived from VLBI observations on the scale of 10–10⁵*M* (e.g., Asada & Nakamura 2012; Hada et al. 2013). In our 3D SANE models, the line $κ = 0.75$ is recovered at low resolution by ECHO and iharm3D, and the more collimated genuinely paraboloid shape $κ = 1$ resembling the solution of Blandford & Znajek (1977) seems to be a better match for the funnel-wall shape at higher resolutions. Incidentally, the match with the $N_θ = 1056$ H-AMR run with the paraboloid shape is particularly good. We suspect that the narrower jet profile compared to Nakamura et al. (2018) is due to the lower $φ \simeq 2–3$ of our benchmark configuration (at high resolution).
<table><caption>**Table 2**
Quantifications—Time Series Data</caption><thead><tr><th>N<sub>θ</sub></th><th>Code</th><th>M&#775;<sub>Peak</sub></th><th>⟨M&#775;⟩</th><th>Φ<sub>BH</sub>/√M&#775;<sub>Peak</sub></th><th>⟨Φ<sub>BH</sub>/√M&#775;⟩</th><th>L&#775;/M&#775;<sub>Peak</sub></th><th>⟨L&#775;/M&#775;⟩</th><th>E&#775;/M&#775;<sub>Peak</sub></th><th>⟨E&#775;/M&#775;⟩</th></tr></thead><tbody><tr><td rowspan="8">96</td><td>Athena++</td><td>0.461</td><td>0.041 ± 0.019</td><td>5.951</td><td>4.164 ± 0.869</td><td>1.95</td><td>1.276 ± 0.326</td><td>0.332</td><td>0.18 ± 0.044</td></tr><tr><td>BHAC</td><td>0.507</td><td>0.074 ± 0.021</td><td>3.099</td><td>1.426 ± 0.16</td><td>1.964</td><td>1.865 ± 0.033</td><td>0.33</td><td>0.112 ± 0.017</td></tr><tr><td>Cosmos++</td><td>0.602</td><td>0.142 ± 0.036</td><td>0.481</td><td>0.239 ± 0.065</td><td>2.268</td><td>2.207 ± 0.025</td><td>0.232</td><td>0.025 ± 0.009</td></tr><tr><td>ECHO</td><td>0.331</td><td>0.04 ± 0.018</td><td>2.949</td><td>2.021 ± 0.377</td><td>2.005</td><td>1.211 ± 0.289</td><td>0.33</td><td>0.124 ± 0.022</td></tr><tr><td>H-AMR</td><td>0.536</td><td>0.091 ± 0.048</td><td>2.908</td><td>1.249 ± 0.131</td><td>2.066</td><td>1.985 ± 0.028</td><td>0.336</td><td>0.046 ± 0.018</td></tr><tr><td>HARM-Noble</td><td>0.797</td><td>0.125 ± 0.05</td><td>2.236</td><td>0.786 ± 0.119</td><td>2.206</td><td>1.964 ± 0.051</td><td>0.327</td><td>0.067 ± 0.017</td></tr><tr><td>iharm3D</td><td>0.645</td><td>0.136 ± 0.06</td><td>3.705</td><td>1.239 ± 0.346</td><td>2.083</td><td>1.958 ± 0.087</td><td>0.331</td><td>0.067 ± 0.011</td></tr><tr><td>KORAL</td><td>0.738</td><td>0.157 ± 0.051</td><td>3.13</td><td>0.458 ± 0.072</td><td>2.1</td><td>2.027 ± 0.026</td><td>0.328</td><td>0.05 ± 0.007</td></tr><tr><td>max/min</td><td></td><td>2.406</td><td>3.903</td><td>12.374</td><td>17.408</td><td>1.163</td><td>1.823</td><td>1.448</td><td>7.096</td></tr><tr><td rowspan="10">128</td><td>Athena++</td><td>0.847</td><td>0.023 ± 0.012</td><td>5.385</td><td>2.995 ± 0.694</td><td>2.082</td><td>1.743 ± 0.208</td><td>0.332</td><td>0.108 ± 0.029</td></tr><tr><td>BHAC</td><td>0.66</td><td>0.078 ± 0.046</td><td>2.492</td><td>1.74 ± 0.407</td><td>1.974</td><td>1.829 ± 0.074</td><td>0.331</td><td>0.123 ± 0.015</td></tr><tr><td>BHAC Cart.</td><td>0.57</td><td>0.163 ± 0.091</td><td>2.137</td><td>1.123 ± 0.277</td><td>1.947</td><td>1.864 ± 0.061</td><td>0.331</td><td>0.097 ± 0.009</td></tr><tr><td>Cosmos++</td><td>0.642</td><td>0.177 ± 0.045</td><td>0.484</td><td>0.245 ± 0.032</td><td>2.133</td><td>2.073 ± 0.019</td><td>0.232</td><td>0.078 ± 0.01</td></tr><tr><td>ECHO</td><td>0.507</td><td>0.059 ± 0.021</td><td>1.575</td><td>1.056 ± 0.182</td><td>2.028</td><td>1.876 ± 0.083</td><td>0.331</td><td>0.073 ± 0.012</td></tr><tr><td>H-AMR</td><td>0.588</td><td>0.117 ± 0.039</td><td>2.378</td><td>1.729 ± 0.204</td><td>2.04</td><td>1.961 ± 0.031</td><td>0.336</td><td>0.061 ± 0.008</td></tr><tr><td>HARM-Noble</td><td>1.28</td><td>0.128 ± 0.051</td><td>3.23</td><td>0.61 ± 0.179</td><td>2.259</td><td>2.093 ± 0.054</td><td>0.328</td><td>0.052 ± 0.017</td></tr><tr><td>iharm3D</td><td>0.841</td><td>0.132 ± 0.035</td><td>3.464</td><td>1.019 ± 0.135</td><td>2.089</td><td>1.964 ± 0.037</td><td>0.325</td><td>0.076 ± 0.01</td></tr><tr><td>IllinoisGRMHD</td><td>0.6</td><td>0.177 ± 0.041</td><td>1.477</td><td>0.779 ± 0.08</td><td>2.141</td><td>2.053 ± 0.044</td><td>0.361</td><td>0.026 ± 0.009</td></tr><tr><td>KORAL</td><td>1.221</td><td>0.199 ± 0.062</td><td>2.351</td><td>0.517 ± 0.097</td><td>2.084</td><td>2.017 ± 0.025</td><td>0.332</td><td>0.056 ± 0.008</td></tr><tr><td>max/min</td><td></td><td>2.525</td><td>8.604</td><td>11.126</td><td>12.24</td><td>1.16</td><td>1.201</td><td>1.552</td><td>4.631</td></tr><tr><td rowspan="7">192</td><td>Athena++</td><td>0.875</td><td>0.134 ± 0.071</td><td>2.739</td><td>1.405 ± 0.39</td><td>2.109</td><td>2.025 ± 0.047</td><td>0.332</td><td>0.065 ± 0.008</td></tr><tr><td>BHAC</td><td>0.825</td><td>0.238 ± 0.042</td><td>2.754</td><td>0.773 ± 0.107</td><td>2.121</td><td>2.072 ± 0.021</td><td>0.331</td><td>0.059 ± 0.009</td></tr><tr><td>ECHO</td><td>0.843</td><td>0.167 ± 0.047</td><td>1.698</td><td>0.506 ± 0.057</td><td>2.139</td><td>2.047 ± 0.033</td><td>0.332</td><td>0.05 ± 0.008</td></tr><tr><td>H-AMR</td><td>0.671</td><td>0.237 ± 0.104</td><td>2.798</td><td>0.691 ± 0.049</td><td>2.086</td><td>2.019 ± 0.026</td><td>0.265</td><td>0.048 ± 0.009</td></tr><tr><td>HARM-Noble</td><td>0.994</td><td>0.201 ± 0.116</td><td>3.593</td><td>0.972 ± 0.216</td><td>2.224</td><td>2.097 ± 0.044</td><td>0.33</td><td>0.057 ± 0.013</td></tr><tr><td>iharm3D</td><td>1.106</td><td>0.18 ± 0.091</td><td>3.161</td><td>1.292 ± 0.19</td><td>2.066</td><td>1.871 ± 0.046</td><td>0.32</td><td>0.057 ± 0.01</td></tr><tr><td>KORAL</td><td>1.01</td><td>0.183 ± 0.081</td><td>2.256</td><td>1.217 ± 0.338</td><td>2.082</td><td>1.987 ± 0.041</td><td>0.331</td><td>0.046 ± 0.007</td></tr><tr><td>max/min</td><td></td><td>1.649</td><td>1.772</td><td>2.116</td><td>2.777</td><td>1.077</td><td>1.121</td><td>1.25</td><td>1.405</td></tr><tr><td>256*</td><td>BHAC</td><td>0.697</td><td>0.203 ± 0.065</td><td>2.552</td><td>1.021 ± 0.091</td><td>2.097</td><td>2.034 ± 0.019</td><td>0.332</td><td>0.064 ± 0.007</td></tr><tr><td>384</td><td>BHAC</td><td>0.91</td><td>0.221 ± 0.067</td><td>3.424</td><td>0.946 ± 0.128</td><td>2.12</td><td>2.043 ± 0.021</td><td>0.331</td><td>0.049 ± 0.006</td></tr><tr><td>1056</td><td>H-AMR</td><td>1.143</td><td>0.152 ± 0.056</td><td>2.324</td><td>1.341 ± 0.068</td><td>1.975</td><td>1.901 ± 0.028</td><td>0.336</td><td>0.05 ± 0.006</td></tr></tbody></table>

**Note.** Quantities in angular brackets ⟨⟩ denote time averages in the interval $t \in [5000, 10,000]M$ with error given by the standard deviation. For further convergence testing, another two BHAC runs (256* using 256 × 256 × 128 cells and 384 × 384 × 384 cells) and an H-AMR run with 1608 × 1056 × 1024 cells are listed here, where the resolutions correspond to effective $N_r \times N_\theta \times N_\phi$ in the disk region.

Comparing the torus density in more detail, in Figure 14 we illustrate contours of $\langle\rho\rangle/\langle\rho\rangle_{\text{max}}$ for the runs with $N_\theta = 192$ cells. Normalized in this way, the agreement in torus extent generally improves as it compensates for the spread in peak densities (see Figure 8). Nonetheless, there remains a large variance in the given contours.

For a more quantitative view, we compute the density scale height:

$$H/R(r) = \frac{\int_{t_{\text{beg}}}^{t_{\text{end}}} \int_{0}^{2\pi} \int_{0}^{\pi} |\pi/2 - \theta_{\text{KS}}| \rho(r, \theta, \phi, t) \sqrt{-g} \, dt \, d\phi \, d\theta}{\int_{t_{\text{beg}}}^{t_{\text{end}}} \int_{0}^{2\pi} \int_{0}^{\pi} \rho(r, \theta, \phi, t) \sqrt{-g} \, dt \, d\phi \, d\theta}, \quad (21)$$

where the averaging time is again taken between $t_{\text{beg}} = 5000M$ and $t_{\text{end}} = 10,000M$. As shown in Figure 15, a scale height of $H/R = 0.25-0.3$ between $n_{\text{KS}} \in [10, 50]$ is consistently recovered with all codes. The largest departures from the general trend are seen in the inner regions for the Cartesian runs that show a slight “flaring up,” which can also be observed in

Figures 9 and 11. Furthermore, the presence of the outflow boundary in Athena++ leads to a decrease of density in the outer equatorial region (see Figure 9), which shows up as kink in Figure 15.

### 5.5. Variability Analysis

The analysis of light curves from accreting compact objects is a powerful diagnostic of the dynamics of the inner accretion flow. To get a feeling for the predictive power of variability from GRMHD simulations in a turbulent regime, we here compare the salient features between codes.

In Sgr A*, episodic flares are detected in X-ray and near-infrared (NIR) on a roughly daily basis (Baganoff et al. 2001; Nowak et al. 2012; Mossoux & Grosso 2017), and the recent detection of orbital motion by the Gravity Collaboration et al. (2018) places the origin of an IR flare within 10 gravitational radii of the BH. In addition to large flares, a low-level continuous variability is present in NIR and radio (e.g., Dodds-Eden et al. 2011). While the nature of the strong coincident X-ray and IR flares is most likely associated with discrete events of particle heating or acceleration (Markoff et al. 2001;
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Resolution dependence of the averaged quantities for all runs of Table 2. The figure has been normalized to the mean of the distribution, allowing the overall spread to be judged, and errors mark the standard deviation over the averaging interval. Upon increased resolution, the spread in the quantities decreases and one can make out a critical resolution $\Delta\theta_p = 0.0125-0.25$, below which accretion rates and horizon-penetrating fluxes tend to converge to the same answer. The accreted angular momentum flux $\dot{L}/\dot{M}$ is very well captured at all resolutions with minimal spread. Note that the data of the magnetic fluxes (top-right panel) for *Cosmos++* lie below the plotting window.</figcaption></figure>

Dodds-Eden et al. 2009), weaker IR flares could also be caused by lensed thermal plasma in the turbulent accretion flow (Chan et al. 2015), and flares in the submillimeter band can be attributed solely to turbulent fluctuations in the accretion flow (Dexter et al. 2009, 2010). On short timescales (below the characteristic timescales $\tau \sim 4$ hr in IR and $\tau \sim 8$ hr in submillimeter), the continuous variability is well described by a red-noise spectrum with power-law slope of $\approx 2$ with no indication for a further break down to $\sim 2-10$ minutes, delimited by the typical sampling (Do et al. 2009; Dexter et al. 2014; Witzel et al. 2018). As pointed out by Witzel et al. (2018), the IR and submillimeter characteristic timescales are in fact statistically compatible, indicating a direct relation of the emitting regions. For low-luminosity AGN, Bower et al. (2015) found that the submillimeter characteristic timescale is consistent with a linear dependence on the BH mass as would be expected for optically thin emission in radiatively inefficient general relativistic models. Hence, in M87, the variability with its characteristic timescale of $45_{-24}^{+61}$ days (Bower et al. 2015) can be obtained by scaling the Sgr A* result.

Due to the potential impact of the spacetime on the variability of the emission,<fnref n="122" /> it is clear that the characterization of the spectrum can hold great merit. From theoretical grounds, a break in the power spectrum at the orbital frequency of the innermost emitting annulus (often identified with the innermost stable circular orbit, ISCO) is expected (see also the discussion in Armitage & Reynolds 2003). This high-frequency break is, however, close to the current sampling frequency of a few minutes for Sgr A*.

To get a feeling for the millimeter variability without actually subjecting each code's output to a rigorous general relativistic radiative transfer calculation (GRRT), an approximate optically thin light curve was calculated as a volume integral over the domain of interest according to Equation (18). Before we discuss the variability from this pseudo-emissivity, let us briefly compare its properties to an actual GRRT calculation of an exemplary simulation output. For this purpose, we employ the BHOSS code (Z. Younsi et al. 2019, in preparation) to compute the 1.3 mm emission corresponding to the EHT observing frequency from the high-resolution BHAC data scaled to Sgr A*.

The light curves are computed assuming a resolution of 1024 × 1024 pixels and a field of view of 100 M in the horizontal and vertical directions, i.e., [-50 M, 50 M]. We here ignore the effect of the finite light travel time through the domain (also known as the “fast light” approximation; see, e.g., Dexter et al. 2010 for a discussion). The chosen synchrotron emissivity is the photon pitch angle-dependent approximate formula from Leung et al. (2011; Equation (72) therein), and the corresponding self-absorption is included. The ion-to-electron temperature ratio used to obtain the local electron temperature as a function of the local ion temperature is fixed to 3 as in the intermediate model of Mościbrodzka et al. (2009). The domain within which the GRRT calculations are performed is set identically to that specified in the emission proxy, i.e., Equation (18). The mass and distance of Sgr A* adopt the standard fiducial values (Genzel et al. 2010), and the 1.3 mm flux is normalized to 3.4 Jy (Marrone et al. 2007) for an observer inclination of 60° using the fiducial snapshot at 10,000 M.

A comparison of the light curves for $i = 0°$, $i = 90°$, the (scaled) accretion rate, and the emission proxy is given in Figure 16. It demonstrates that both the accretion rate and the emission proxy follow reasonably well the trend of the GRRT

<footnote n="122">¹²² Numerical simulations by Dolence et al. (2012) have observed a steepening of the IR and X-ray spectral index to ~3 accompanied by quasi-periodic oscillations (QPOs) at the frequency of the ISCO (corresponding to $f_0 = 0.041/M$ in our case and thus periods of 8 minutes if scaled to Sgr A*).</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Disk-averaged quantities at the resolution of 96³ in the domain of interest for all codes. Data have been averaged over the $\theta$ and $\phi$ angles as well as over the time interval $t \in [5000M, 10,000M]$. Standard deviation over time is marked by the shaded region (where data were available) to judge the secular evolution within the averaging window. For reference, the H-AMR $N_\theta = 1056$ solution was added. See text for details.</figcaption></figure>

light curve on the scale of a few hours. However, the accretion rate shows more pronounced variability on small scales than the proxy or the GRRT calculations, which can be understood as a consequence of the extended size of the emission region in the latter two prescriptions. Doppler boosting adds additional variability to the light curve, which leads to larger fluctuations in the edge-on light curve compared to the face-on ($i = 0^\circ$) curve.

Having validated the emission proxy as a reasonable descriptor for the flux in the millimeter and, we now compute the power-spectrum densities (PSDs) of the emission proxy $\mathcal{L}$. The PSD is obtained by dividing the data into three non-overlapping segments between 5000 M and 10,000 M, and fast Fourier transformation of each Hamming-windowed segment, followed by Fibonacci binning in frequency space. The three binned PSDs are then averaged, and we also take note of the standard deviation between them. To validate this procedure, synthetic pure red-noise data with a PSD of $f^{-\beta}$ with $\beta \in \{2, 3\}$ was also created and is faithfully recovered by the algorithm.<fnref n="123" /> It is clear that as each window only contains 10 hr, the observed break at the characteristic frequency $\tau \sim 8$hr (for Sgr A*) is not measurable by this procedure.

The data is presented in Figure 17 for all available resolutions. At low and mid resolution, all light curves are

compatible with a red-noise spectrum with $\beta = 2$, and we obtain a steepening near $f = 0.1M^{-1}$. For high resolution, this break is less apparent; however, there is a slight trend for an overall steepening to $\beta = 3$, e.g., in HARM-Noble, iharm3D, and IllinoisGRMHD. It is reassuring that the light curve obtained by GRRT at $i = 0^\circ$ (light gray) shown in the lower-left panel displays a very similar behavior to the corresponding data obtained via the proxy (BHAC code, blue).

Because the accretion rates have been computed with all codes and generally are an easily accessible and relevant diagnostic in GRMHD simulations, which can give some guidance on the overall variability features (see, e.g., Reynolds & Miller 2009; Dexter et al. 2010), we perform the same analysis with the time series of $\dot{M}$. The right panels of Figure 17 indicates that the low-frequency power law seen in $\mathcal{L}$ is not recovered in $\dot{M}$. For frequencies $f \lesssim 0.01M^{-1}$, the power in all codes is definitely shallower than $f^{-2}$, approaching a flicker-type noise, $f^{-1}$, and we find indications for a spectral break in the vicinity of the ISCO frequency in several codes. The $\beta = 1$ slope for $f < 10^{-2}M$ is consistent with that of, e.g., Hogg & Reynolds (2016), who extracted the accretion rate at the ISCO itself. These authors also note that the PSD of an “emission proxy” is steeper than the one obtained from $\dot{M}$, due to the fact that additional low-frequency power is added from larger radii in the extended emission proxy.

To also compare the variability magnitude, we have computed the root mean square (rms) of the accretion rate after zeroing out

<footnote n="123">¹²³ Note that due to the spectral leakage of the finite time light curves, using a boxcar window resulted in the artificial flattening of the $\beta = 3$ spectrum to $\beta \sim 2.5$</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 7. As Figure 6 for the 128³ data.</figcaption></figure>

all Fourier amplitudes with frequency below 1/1000M. This serves as an effective detrending of the low-frequency secular evolution. In order to avoid edge effects, we compute the rms in the region $t \in [2500, 7500]M$. Exemplary data of the remaining high-frequency variability is shown for the high-resolution runs in the left panel of Figure 18, where we have normalized by the average accretion rate in the region of interest. One can visually make out differences in variability amplitude with the BHAC data on the low end and HARM-Noble and iharm3D on the high end. The (normalized) rms values shown against midplane resolution (right panel) quantify this further and indicate quite a universal behavior with values in the range 1.2–1.6 across all codes and resolutions.<fnref n="124" /> With increased resolution, all codes tend to be attracted to the point $k = \text{rms}(\dot{M})/\langle\dot{M}\rangle \simeq 1.3$. Repeating the same analysis with the emission proxy yielded essentially the same outcome. This quite striking result (after all, there is a scatter in mean accretion rates itself by a factor of ~8 in the sample) is a restatement of the rms–flux relationship, which is a ubiquitous feature of BH accretion across all mass ranges (Uttley & McHardy 2001; Heil et al. 2012). The quotient $k \sim 1.3$ is consistent with the recent simulations of Hogg & Reynolds (2016) who find an rms–flux relationship of the accretion rate with $k = 1.4 \pm 0.4$.

Summarizing the comparison of variabilities, although it remains a challenge to extract accurate power spectra from the

finite time series of the data, all codes agree quite well on the salient features:

1. A $\beta \simeq 2$ power-law slope is recovered for the emission proxy $\mathcal{L}$ with all codes and resolutions.

2. The PSD for $\mathcal{L}$ and $\dot{M}$ steepens at $f \simeq 0.1M^{-1}$.

3. The PSD of the accretion rate is shallower than the one for the proxy and is more accurately described by an index of $\beta \simeq 1$.

4. The “rms–flux” relationship is quite universally recovered for $\dot{M}$ and $\mathcal{L}$ for all codes and resolutions.

## 5.6. Further Analysis

To get a deeper understanding of the dynamical effect that increasing the grid resolution has, some further analysis was carried out with the results of BHAC and H-AMR.

### 5.6.1. Maxwell Stress and $\alpha$

We first compute the disk’s $\alpha$-parameter due to Maxwell stresses; hence, we are interested only in the term $-b^r b^\phi$ of the energy–momentum tensor (3). As noted, however, by Krolik et al. (2005), Beckwith et al. (2008), and Penna et al. (2013), the stress depends on the coordinate system and cannot unambiguously be defined in general relativity. The best one can do is to measure the stress in a locally flat frame comoving with the fluid. Even so, further ambiguities arise when defining the orientation of the comoving tetrad and hence the $b^{(r)}$ and $b^{(\phi)}$ directions. With the aim of merely comparing codes and

<footnote n="124">¹²⁴ Using the nonnormalized (raw) values of the rms yielded far less consistent results with differing resolution-dependent trends between codes. Furthermore, exceptional low-number statistics “flare” events dominate over raw light curves, as seen in the low-resolution iharm3D data in Figure 2, for example.</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 8. As Figure 6 for the 192³ data.</figcaption></figure>

convergence properties in global disk simulations, we here settle for the simpler diagnostic in the Kerr–Schild coordinate frame. As shown in Appendix D, for our case, this is a fair approximation to the fluid-frame stress in the most commonly used basis. Hence, using the same operation as for the disk profiles (see Equation (17)), we compute the shell average of

$$ \langle w^{r\phi} \rangle (r, t) := \langle -b^r \sqrt{\bar{\gamma}_{rr}} b^{\phi} \sqrt{\bar{\gamma}_{\phi\phi}} \rangle \quad (22) $$

and define the $\alpha(r, t)$ parameter due to Maxwell stresses as

$$ \alpha(r, t) := \frac{\langle w^{r\phi} \rangle (r, t)}{\langle P_{\text{tot}} \rangle (r, t)}. \quad (23) $$

This $\alpha(r, t)$ is further averaged in time to yield $\alpha(r)$ or volume-averaged to yield $\bar{\alpha}(t)$.

Figure 19 illustrates the resolution dependence of the $\alpha$-parameter for the H-AMR and BHAC runs. Overall, there is good agreement between BHAC and H-AMR for $N_\theta > 128$. Not surprisingly, the stress profiles resemble closely the run of $\langle \beta^{-1} \rangle$ already shown in Figures 6–8, but $\alpha(r)$ is lower by a factor of $\sim 3$. In particular, there is no “stress edge” at or within the ISCO (located at $r_{\text{KS}} = 2.04M$), where $\alpha$ drops to zero. This is consistent with the results of Krolik et al. (2005) for the highly spinning case and marginally consistent with Penna et al. (2013), who find that $\alpha(r)$ peaks at $\sim 2–3M$. The latter difference can be explained in part through the coordinate choice and in part through the higher spin adopted here (see also Appendix D). With increased resolution, the $\alpha$ in our global simulations increases, as

opposed to what is observed in local shearing boxes. For example, Bodo et al. (2014) and Ryan et al. (2017) reported $\bar{\alpha} \propto N^{-1/3}$ for stratified local simulations and Guan et al. (2009) found $\bar{\alpha} \propto N^{-1}$ for the unstratified case ($N$ is the number of grid cells per scale height). Over time, $\bar{\alpha}$ decreases, but less so as the grid resolution is increased. This secular trend and the large gap in resolution between the $N_\theta = 384$ BHAC run and the $N_\theta = 1056$ H-AMR run make a strict quantification of the convergence behavior difficult. At a comparatively high resolution of $N_\theta = 384$, though, it seems that $\bar{\alpha}$ has not converged yet.

### 5.6.2. Viscous Spreading

The viscosity induced by the MRI leads to the accretion of material and viscous spreading of the disk as a whole. Hence, with increased resolution and thus higher disk stresses, a global effect on the viscous spreading is expected. We quantify this by calculating the rest-frame density-weighted radius, $\langle r \rangle$, according to

$$ \langle r \rangle (t) := \frac{\int_{0}^{2\pi} \int_{0}^{\pi} \int_{r_{\text{h}}}^{r_{\text{max}}} r_{\text{KS}} \rho(r, \theta, \phi, t) \sqrt{-g} \, d\phi \, d\theta \, dr}{\int_{0}^{2\pi} \int_{0}^{\pi} \int_{r_{\text{h}}}^{r_{\text{max}}} \rho(r, \theta, \phi, t) \sqrt{-g} \, d\phi \, d\theta \, dr}, \quad (24) $$

where the outer radius of integration has again been set to the domain of interest $r_{\text{max}} = 50M$. As more mass is expelled beyond $50M$; however, a larger integration domain would yield
<figure><img src="image_11.png" /><figcaption>**Figure 9.** *t*- and *φ*-averaged data. We show averages of the rest-frame density, inverse plasma *β*, and the magnetization *σ*. From top to bottom: Athena++ (N<sub>θ</sub> = 192), BHAC (N<sub>θ</sub> = 192), and BHAC Cart. (see Table 1).</figcaption></figure>

somewhat different results. The corresponding temporal
evolution is displayed in Figure 20 for BHAC and H-AMR,
where linear and piecewise parabolic reconstruction has been

employed. It is directly apparent that the disk size increases
with resolution. Furthermore, the saturation of the low- and
medium-resolution PPM runs occurring at *t* ~ 6000 to values
<figure><img src="image_12.png" /><figcaption>**Figure 10.** As Figure 9 for the codes: Cosmos++ ($N_{\theta} = 128$), ECHO ($N_{\theta} = 192$), and H-AMR ($N_{\theta} = 1056$).</figcaption></figure>

of $\langle r \rangle \in [24, 27]$ are indications that the turbulence starts to decay at this time. It is noteworthy that the linear reconstruction experiments with H-AMR indicate a shrinking of the disk even at a high resolution of $N_{\theta} = 192$.

The behavior of $\langle r \rangle$ is consistent with the fact that higher resolution in the disk leads to larger Maxwell stresses, driving stronger disk winds and thus larger disk expansion. This lowers the overall disk density (as seen in Figures 6–8:
<figure><img src="image_13.png" /><figcaption>**Figure 11.** As Figure 9 for the codes: HARM-Noble, iharm3D (both $N_{\theta} = 192$), and IllinoisGRMHD (see Table 1).</figcaption></figure>

⟨ρ(r)⟩, but the disk’s magnetic field is affected to a lesser extent (see Figures 6–8: ⟨|B|(r)⟩), hence the disk’s α increases further in a nonlinear fashion. As shown in Section 5.7,

however, the disk radius starts to converge for $N_{\theta} > 192$, indicating that the global effect of MRI stresses is accurately captured.
<figure><img src="image_14.png" /><figcaption>Figure 12. As Figure 9 for KORAL with $N_\theta = 192$.</figcaption></figure>

### 5.6.3. MRI Quality Factors

These resolution-dependent effects can be understood by an inspection of the “MRI quality factors” (Q-factors for short), defined as the number of cells available to resolve the fastest growing MRI mode in a given direction. While the fastest growing mode might be well resolved in the initial condition, whether this remains to be the case throughout the simulation can only be established a posteriori. If the simulation fails to capture at least the fastest growing mode, turbulent field amplification cannot proceed, and the field will decay due to numerical dissipation, leading also to a cessation of mass accretion. Although strictly a criterion for sufficient resolution of the linear MRI, the Q-factors also inform on the adequacy in the nonlinear regime (see, e.g., Sano et al. 2004; Noble et al. 2010; Hawley et al. 2013). The consensus values for adequate resolution are given by $Q^{(z)} \ge 10$, $Q^{(\phi)} \ge 20-25$ of Hawley et al. (2011, 2013); however, lower toroidal resolution requirements have also been suggested, for example, $Q^{(z)} \ge 10-15$ for $Q^{(\phi)} \approx 10$ in Sorathia et al. (2012). Furthermore, these authors note that the toroidal and poloidal resolutions are coupled, and low $Q^{(z)}$ can be compensated by $Q^{(\phi)} \ge 25$. Hence, the product of the quality factors $Q^{(z)}Q^{(\phi)} \ge 200-250$ has also been suggested (Narayan et al. 2012; Dhang & Sharma 2019) as a quality metric.

As a basis for our MRI quality factors, we compare the wavelength of the fastest growing MRI mode evaluated in a fluid-frame tetrad basis $e_\mu^{(r)}$ with grid spacing. In particular, we orient the tetrads along the locally nonrotating reference frame (LNRF; see, e.g., Takahashi 2008 for the transformations). Take, for example, the $\theta$-direction:

$$ \lambda^{(\theta)} := \frac{2\pi}{\sqrt{(\rho h + b^2)\Omega}} b^{\mu} e_{\mu}^{(\theta)}, \qquad (25) $$

and adopt the corresponding grid resolutions as seen in the orthonormal fluid frame:

$$ \Delta x^{(\theta)} := \Delta \theta^{\mu} e_{\mu}^{(\theta)}, \qquad (26) $$

where $\Delta\theta^\mu = (0, 0, \Delta\theta, 0)$ is the distance between two adjacent grid lines, and $\Omega = u^\phi/u^\mu$ is the angular velocity of the fluid. The resulting $Q^{(\theta)} := \lambda^{(\theta)}/\Delta x^{(\theta)}$ are shown for the standard resolutions at the final simulation times in the top panel of Figure 21. It is clear that the low- and medium-resolution simulations are not particularly well resolved with regard to the poloidal MRI wavelengths. The time and spatial average values in the domain of interest are $Q_{96}^{(r,\theta,\phi)} \simeq 2.9, 3, 14$, $Q_{128}^{(r,\theta,\phi)} \simeq 3, 3.8, 17$, and $Q_{192}^{(r,\theta,\phi)} \simeq 7.3, 7.9, 19.7$. This is still somewhat below the nominal $Q^{(z,\phi)} \sim 10, 20$ of Hawley et al. (2011) but well above the six cells suggested by Sano et al. (2004) to adequately capture the saturation level of the MRI. However, it seems that only the highest resolution case keeps MRI-driven turbulence alive up to the end of the simulation. At $96^3, 128^3$, or $192^3$ resolution in the domain of interest, the bottleneck is clearly in the poloidal direction with $Q^{(r)} \simeq Q^{(\theta)}$ but the azimuthal factor $Q^{(\phi)}$ is approximately twice as large. Hence, one might surmise that the resolution in the $\phi$-direction will have little influence on the overall evolution. This has been confirmed by additional experiments with 96 instead of 192 azimuthal cells using BHAC and by runs where the $r$ direction was also reduced using the KORAL code. In both cases, the agreement was very good, with horizon-penetrating fluxes consistent with the full resolution case. As noted by Hawley et al. (2011), an increase in azimuthal resolution can also lead to an improved behavior of the poloidal quality factors. We find the same; however, the effect is not dramatic. The dotted curve in Figure 21 shows that with half the azimuthal resolution, $Q^{(r)}$ and $Q^{(\theta)}$ indeed slightly decrease (to $Q^{(r,\theta)} \simeq 5.3, 6.8$), and $Q^{(\phi)} \simeq 9.6$ is essentially halved. Upon increasing the resolution to $N_\theta = 384$, we obtain $Q_{384}^{(r,\theta,\phi)} \simeq 21.3, 24.5, 40.9$, and hence, a superlinear increase in the Q-factor compared to the $N_\theta = 192$ case. This behavior
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Contours of the jet–disk boundary $\sigma = 1$ for increasing resolution (left) and a zoom into the region $\theta \in [0^\circ, 45^\circ]$ (right). Only in the high-resolution case is the jet recovered with all codes. To guide the eye, we also show paraboloid curves according to Equation (20), corresponding to $z \propto R^2$ ($\kappa = 1$, thin dotted black) and $z \propto R^{1.6}$ ($\kappa = 0.75$, thin solid black). The overall spread of the jet angle (given by its $\theta$-coordinate at $r = 50M$) is on the $5^\circ$ level.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparison of the density contours at fixed ratios of the peak value for the high-resolution models. Contours are placed at $\langle\rho\rangle/\langle\rho\rangle_{\text{max}} \in [0.0078125, 0.125, 0.5]$.</figcaption></figure>

will be encountered again for the very high-resolution run and is discussed further in the next section.

These results give a framework to explain the resolution-dependent trends in the data. As an example, the formation of the peculiar “mini torus” in the low- and mid-resolution BHAC

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Scale height profiles of the high-resolution models ($N_\theta = 192$) according to Equation (21). As expected, all codes are consistent in this quantity, with the largest departures seen in the inner regions for the Cartesian runs.</figcaption></figure>

runs is likely a result of insufficiently resolved MRI in these runs. Furthermore, the convergence of accretion rates and MAD parameter $\phi$ shown in Figure 5 coincides with the point when Q-factors reach a sufficient level.

## 5.7. A Very High-resolution Run

With H-AMR, we are able to produce a global simulation with a resolution of $N_r \times N_\theta \times N_\phi = 1608 \times 1056 \times 1024$ with five levels of static mesh refinement in the $\phi$-direction (see Section 4.1.5), bringing us to local shearing box regimes (e.g., Davis et al. 2010; Ryan et al. 2017). In terms of the density scale height $H/R \simeq 0.25–0.3$, we obtain 85–100 cells per scale height and disk-averaged quality factors $Q_{1056}^{(r,\theta,\phi)} \ge 120, 100, 190$. This is more than sufficient for capturing MRI-driven turbulence as per the criteria ($Q^{(z,\phi)} \ge 10, 25$ in cylindrical coordinates). Indeed, Hawley et al. (2011) suggests a resolution of $600 \times 450 \times 200$ for a global simulation with $H/R = 0.1$, which we have achieved. Sufficiently resolved global simulations are essential for reproducing microphysical phenomena such as MRI in a macrophysical environment. Hawley et al. (2013) notes that even though MRI is a local instability, turbulent structures may form on larger scales (seen in large shearing boxes by Simon et al. 2012). Additionally, shearing box simulations inadvertently affect poloidal flux generation by implicitly conserving the total vertical
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Comparison of the light curves scaled to Sgr A* from the BHAC high-resolution run showing the typical short-term variability and long-term stability of low-luminosity AGN: emission proxy $\mathcal{L}_{\text{proxy}}$ (blue), accretion rate $\dot{M}$ (red), and two synthetic light curves obtained with the GRRT BHOSS code: inclination $i = 0^{\circ}$ ($\mathcal{L}_{\text{proxy},i=0^{\circ}}$; black, thick) and $i = 90^{\circ}$ ($\mathcal{L}_{\text{proxy},i=90^{\circ}}$; black, thin). The accretion rate and emission proxy have been scaled to match the mean flux of the $i = 0^{\circ}$ calculation, and time has been reset at $t_0 = 5000 M$ of the simulation. Both follow the overall trend of the GRRT light curve, though the accretion rate exhibits substantially more small-scale variability. One hour corresponds to $\sim 180 M$ simulation time.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Periodograms, compensated with $f^2$, of the emission proxy $\mathcal{L}$ (left) and for the accretion rate $\dot{M}$ (right) for increasing resolutions (top to bottom). Also, the PSD of a synthetic light curve at 230 GHz and inclination $i = 0^{\circ}$ is shown in the lower panel (light gray curve). To better visually separate the curves, they have been progressively shifted by factors of 10. The blue dotted vertical line denotes the orbital frequency of the ISCO. In the data for $\dot{M}$, we obtain a flatter low-frequency spectrum and a peak in the vicinity of the ISCO frequency of $f_0 = 0.04 M^{-1}$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Variability of the detrended accretion rates. The long time secular evolution has been removed by zeroing out all Fourier amplitudes below 1/1000M. An example of the high-resolution data is shown in the left panel. In the right, we compare the resolution dependence in the rms accretion rate normalized by its mean value. We obtain a clear “rms–flux” relation in all codes with quotient converging to $k \simeq 1.3$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Maxwell stress analysis. Left: profiles of the disk $\alpha(r)$ parameter. Shaded regions indicate the standard deviation over time. Right: Time evolution of the spatially averaged $\bar{\alpha}(t)$. The $\alpha(r)$ parameter is constant throughout most of the disk and behaves similarly to $\langle\beta^{-1}\rangle$ as shown in Figures 6–8.</figcaption></figure>

magnetic flux, and hence, may not be as reliable as global simulations in dealing with large-scale poloidal flux generation (Liska et al. 2018b), which is required for jet launching and propagation.

It is quite striking that the increase of the MRI quality factors for the very high-resolution run is beyond the factor of ~5.5, which would be gained by going from $N_\theta = 192$ to $N_\theta = 1056$. Instead, we find an extra factor of 2 increase in Q-factors (see Appendix A), consistent with the simultaneous rise in disk magnetization by the same amount. Overall trends of parameters over time shown in Figure 4 are quite similar to the iharm3D and KORAL high-resolution data, with some other 192³ PPM data within range.

This gives a first impression of the behavior of viscous angular momentum transport and the turbulent dissipation and generation of magnetic energy. As evidenced by Figure 22, there are signs of convergence in the disk radius at $N_\theta = 1056$,

though the disk magnetization continues to increase slightly with increased resolution. Due to the lack of intermediate data between $N_\theta = 192$ and $N_\theta = 1056$, it is, however, not possible to judge whether a saturation point has already been reached in between. The steady increase of the disk magnetization with resolution has also been noted by Shiokawa et al. (2012), who performed iharm3D simulations up to $N_\theta = 384$. Quantitatively, the values of $\beta \sim 15–20$ and the resolution-dependent trend reported in their Figure 3 seem consistent with our findings.

One notes that the magnetization reverses its trend in most codes as $\langle\beta^{-1}\rangle$ at first even decreases but then picks up at $N_\theta = 128$ and continues to increase up to the highest resolutions tested. This is suggestive of a “resolution threshold” above which turbulent amplification of the magnetic field starts to operate, leading to the higher average value of the magnetization. At the same time, the near-constant profile of
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Barycentric radii of the disk according to Equation (24). The data have been computed with *BHAC* and *H-AMR*, both using linear (PLM; dashed curves) and piecewise parabolic reconstruction (PPM; solid curves). The disk stalls its spreading early for low and medium resolutions, and even shrinks when linear reconstruction is employed.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** MRI quality factors for increasing resolution (left to right) obtained at $t = 10,000 M$ with the *BHAC* code (top row) and corresponding temporal evolution of the average quality factor in the above highlighted region (bottom row). In the high-resolution time series, dotted lines indicate the *Q*-factors for a run with half the azimuthal resolution (hence 96 cells). At high resolution, the late time decay of $Q^{(r)}$, $Q^{(\theta)}$ is halted and a mean value of ~8 is reached, indicating just about sufficient resolution to keep the MRI going.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Left: convergence in the disk spreading as quantified by the barycentric radius at $t = 10K M$. The data have been computed with *BHAC* and *H-AMR* both using linear (PLM) and piecewise parabolic reconstruction (PPM). We note indications for convergence against the very high-resolution run. Right: average magnetization of the disk for $r_{KS} \in [6, 40]M$ against the physical midplane resolution. Similar to the disk $\alpha$, a weak dependence on resolution is obtained in all codes, with $\langle \beta^{-1} \rangle$ increasing by $\sim 3$ from lowest to highest resolution in the set.</figcaption></figure>
$\langle\beta^{-1}(r)\rangle$ is established throughout the disk, reflected also in a decrease of the scatter in Figure 22.

# 6. Discussion and Conclusions

Using a standardized setup of radiatively inefficient BH accretion and a set of relevant diagnostics, we have compared the MRI-driven turbulent quasi-stationary state in nine GRMHD codes which are widely used in the community. Many of the codes have been developed independently; others share the same heritage but are completely rewritten. There are both differences and commonalities between the codes. Listing the most important common elements, we note that all apply overall second-order conservative schemes and preserve the divergence of the magnetic field to machine precision. All except iharm3D applied the PPM for spatial reconstruction. Regarding the algorithmic differences, the order of time integration, use of approximate Riemann solver (HLL and LLF), and specific integration of the induction equation are the most striking ones. Furthermore, two configurations have been computed in Cartesian coordinates, while mostly spherical coordinates are used. The results shed light on the systematics between codes and elucidate the resolution requirements to reach a certain level of agreement. Generally, we find that the level of agreement between codes improves for all diagnostics when the resolution is increased toward convergence. Once sufficient resolution is used, key parameters like the accretion rate of the saturated turbulent state agree within their temporal variations among all codes, and the spread in the mean accretion rate is given by a factor of ~1.7 between the lowest and the highest.

## 6.1. Resolution and Convergence

We have performed simulations at resolutions commonly used for BH accretion and forward modeling of EHT observations. In fact, early GRMHD models employed simulations that use a comparable number of or fewer cells than our low- and mid-resolution setups of 96³–128³ cells. For example, the fiducial runs of McKinney & Blandford (2009) use 256 × 128 × 32 cells (though with fourth-order reconstruction), and Fragile et al. (2007) used a standard effective resolution of 128³ for their tilted disks. Economizing on the spatial resolution is particularly important for large parameter surveys as in Mościbrodzka et al. (2016a), who applied HARM3D with 96 × 96 × 64 cells and for long time evolutions as in Penna et al. (2010), who ran a large suite of simulations with resolutions of effective 256 × 64 × 64 cells (accounting for the restricted Δφ = π). Although comparing cell numbers should be done with a grain of salt, taking into account gridding and initial conditions, these numbers give an indication of the recent state of affairs.

Balancing the desire to explore a large parameter space (as in the generation of a simulation library for EHT observations) against numerical satisfaction, one will end up performing simulations that are “just about” sufficiently resolved. Hence, it is interesting to note that with the choice of our resolutions, we have sampled a resolution edge in most codes. This is best seen in the runs of accretion rates and magnetic fluxes, which substantially reduce their scatter once midplane resolutions of Δθp ≃ 0.25–0.0125M are adopted. Below this threshold, the algorithmic differences still dominate, leading to large scatter with sometimes opposite trends between codes. The resolution

threshold was also confirmed by inspection of the MRI quality factors in BHAC and H-AMR. Only the high-resolution PPM runs would qualify as resolving the saturated MRI per the criteria of Hawley et al. (2011).

Using our sample of runs at different resolutions, we have also investigated the convergence of global quantities of interest. While it is clear that convergence in the strict sense cannot be achieved for ideal MHD simulations of turbulence as viscous and resistive length scales depend on the numerical resolution and numerical scheme (e.g., Kritsuk et al. 2011), convergence in certain physical quantities such as the stress parameter $\alpha$ or the “magnetic tilt angle” between the average poloidal field directions could still be achieved (see also the discussions in Sorathia et al. 2012; Hawley et al. 2013). Hence, we find that while the viscous spreading of the disk appears converged at accessible resolutions (~192³ cells), the disk $\alpha$ parameter and magnetization continue their weak resolution-dependent rise up to ~1024³ cells. The latter confirms and extends the findings of Shiokawa et al. (2012) by a resolution factor of ×2.75.

## 6.2. Systematics

Turning to the systematics between codes, we observe a veritable diversity in the averaged density profiles. Although the spread diminishes with increasing resolution, there is still a factor of ~3 difference in peak density at nominal 192³ resolution. As the magnetic field profiles in the disk are quite robust against codes and resolutions, the lower density and pressure values can lead to a nonlinear runaway effect as Maxwell stresses in the simulations with lower densities will increase, leading to additional viscous spreading of the disk and hence again lower peak densities.

When modeling radiative properties of accretion systems as, e.g., in the EHT workflow, the variance in absolute disk densities is somewhat mitigated by the scale freedom of the ideal MHD simulations, which allow the value of the density together with the other MHD variables to be rescaled. Not only do the absolute dissipative length scales depend on resolution, but their ratio, given through the effective magnetic Prandtl number $Pr_m$, also varies with the Riemann solver, reconstruction method, magnetic field solver, grid, and resolution (e.g., Kritsuk et al. 2011). Moreover, it is well known that transport properties of the nonlinear MRI depend critically on $Pr_m$ (Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009; Longaretti & Lesur 2010), and hence, some systematic differences in the saturated state are in fact expected to prevail even for the highest resolutions.

Another systematic is introduced by the axial boundary condition: the jet-disk boundary $\sigma = 1$ is faithfully recovered only at high resolution, and the jet can be particularly skimpy when a polar cutout is used. In the low- and medium-resolution HARM-Noble and Cosmos++ simulations, we found that the axial excision significantly decreases the magnetization of the funnel region. The reason for the loss of the magnetized funnel is well known: standard “soft” reflective boundaries at the polar cones allow a finite (truncation error) numerical flux of $B^r$ to leave the domain. This effect diminishes with resolution and can be counteracted by zero-flux “hard” boundaries at the polar cutout as noted by (Shiokawa et al. 2012). In Cosmos++, the use of outflow boundaries in the excised region leads to a particularly striking difference, resembling the situation described by Igumenshchev et al. (2003).
Inspecting the trends in the jet collimation profiles, we note that there is a close relationship between the jet width and the horizon-penetrating magnetic flux. Indeed, Table 2 and Figure 13 for the $N_\theta = 128$ data reveal that a one-to-one relation holds for the four models with the smallest $\langle\phi\rangle$ and the one with the absolute largest value of $\langle\phi\rangle$ where the differences in opening angle are most pronounced. Because the jet collimation plays an important role in the acceleration of the bulk flow beyond the light cylinder (e.g., Camenzind 1986; Komissarov 2007), this systematic will translate into an effect on the overall flow acceleration of the funnel jets.

### 6.3. The Floor Region

We have also calculated the 2D averaged maps of density, (inverse-) plasma $\beta$, and magnetization. These allow the variance of the jet-opening angles between codes to be quantified, and we find that the spread around a “mean simulation” is reduced to $\sim 2.5$ at $r = 50M$ at high resolutions. Interestingly, the jet boundaries up to $r_{KS} = 50M$ follow closely the paraboloid shape of the solution by Blandford & Znajek (1977). It will be insightful to see how this result changes quantitatively when the magnetic flux on the BH is increased toward the MAD case, which is known to exhibit a wider jet base. The maps of plasma $\beta$ and magnetization truthfully illustrate the variance in the jet region, which is introduced by the different floor treatments. Demonstrating the considerable diversity serves two purposes: (i) to underline that the plasma parameters of the Poynting-dominated jet region should not be taken at face value in current MHD treatments, and (ii) that despite this variance in the jet region itself, the effect on the dynamics (e.g., the jet-opening angle and profile) is minor.

### 6.4. Time Variability

Using an emission proxy tailored to optically thin synchrotron emission from electrons distributed according to a relativistic Maxwellian, we have computed light curves with all codes and resolutions and have compared their power spectra. The same analysis was performed with the accretion rates extracted at the BH horizon. In broad terms, all PSDs of the light curves are compatible with a red-noise spectrum up to $f \simeq 1/10M$, where a steepening is observed. Inspection of the PSD for the accretion rate yields flatter PSDs $\sim f^{-1}$ for $f < 1/100M$ and a similar steepening at $f \simeq 1/10M$. This is quite consistent across all codes and resolutions, and agrees with earlier results of Armitage & Reynolds (2003) and Hogg & Reynolds (2016). The steeper low-frequency PSDs of the proxy can be explained by noting that the integration over the disk adds additional low-frequency fluctuations from larger radii. Whether the high-frequency break occurs at the ISCO frequency of $f_0 = 0.041/M$ or at somewhat higher frequencies $f > 0.1M$ is not that clear cut in the data; however, the presence of a high-frequency break indicates an inner annulus of the emission at or within the ISCO. With the detrended time series of the accretion rates and emission proxy, we have computed the rms variability and found an ubiquitous relationship between the rms and the absolute value (rms–flux relationship) with slope of $k \simeq 1.2–1.6$ across all codes and resolutions. Upon increased resolution, $k$ tends to converge toward $\sim 1.3$ for our benchmark problem.

### 6.5. Implications

This first large GRMHD code comparison effort shows that simulation outcomes are quite robust against the numerical algorithm, implementation, and choice of grid geometry in current state-of-the-art codes. Once certain resolution standards are fulfilled (which might be reached at differing computational expense for different codes), we can find no preference favoring one solution against the other. In modeling the EHT observations, we find it beneficial to use several of the codes tested here interchangeably. In fact, a large simulation library comprising 43 well-resolved GRMHD simulations has been created for comparison to the observations, using *iharm3D*, *KORAL*, *H-AMR*, and *BHAC* (Event Horizon Telescope Collaboration et al. 2019b). Several parameter combinations have been run with multiple codes for added redundancy, and the diagnostics which were calibrated here are used for cross-checks.

To serve as benchmark for future developments, the results obtained here are freely available on the platform *cyverse*. This also includes the unprecedented high-resolution H-AMR run for which a thorough analysis might provide us with new general insights into rotating turbulence. See Section 7 for access instructions. Furthermore, animations of the simulation output can be found in the supporting material at 10.7910/DVN/UCFCLK.

### 6.6. Outlook

The benchmark problem of this work, a spin $a = 0.9375$, turbulent MHD torus with scale height $H/R \approx 0.3$, and normalized horizon-penetrating magnetic flux of $\phi \approx 2$ falls into the class of radiatively inefficient (Narayan et al. 1995) SANE disks and is perhaps the simplest case one might consider (its widespread use likely goes back to the initial conditions provided with the public HARM code; Gammie et al. 2003). To increase the challenge, one might consider MHD simulations of MAD (e.g., Narayan et al. 2003). The numerical problems and new dynamics introduced by the increase of magnetic flux are considerable. New violent interchange-type accretion modes and stronger magnetized disks (potentially with suppressed MRI) come into play (e.g., Tchekhovskoy et al. 2011; McKinney et al. 2012), presenting a physical scenario well suited to bringing GRMHD codes to their limits. A resolution study of the MAD scenario was recently presented by White et al. (2019), showing the difficulty of converging various quantities of interest, e.g., the synchrotron variability. How different GRMHD codes fare with the added challenge shall be studied in a future effort, where also the jet dynamics will be more in the focus. This shall include the dynamics in the funnel, e.g., the properties of the stagnation surface (Nakamura et al. 2018) and the acceleration/collimation profiles.

Another area where code comparison will prove useful is in the domain of nonideal GRMHD modeling, e.g., radiative MHD (e.g., Fragile et al. 2012; McKinney et al. 2014; Sądowski et al. 2014; Ryan et al. 2015; Takahashi et al. 2016) and resistive MHD (e.g., Bugli et al. 2014; Qian et al. 2017; Ripperda et al. 2019b). This would be particularly informative as such codes are not yet widely used in the community, e.g., no public versions have been released to date.

One of the prevailing systematics in GRMHD modeling lies in the ILES approximation, which is typically applied.
Development of explicit general relativistic treatments of magnetic diffusivity and viscosity (e.g., Fragile et al. 2018; Fujibayashi et al. 2018) will soon allow direct numerical simulations of turbulent BH accretion covering both the low and high $Pr_m$ regimes, which are expected to be present in such systems (Balbus & Henri 2008), to be performed.

## 7. Supplementary Information

Animations of the quantities given in Figures 9–12 were prepared with BHAC, ECHO, and H-AMR for meridional slices and can be found at doi:10.7910/DVN/UCFCLK.

Processed data used to create all figures in this manuscript as well as raw snapshot data of several high-resolution runs (including the $N_\theta = 1056$ H-AMR run) is available through the *cyverse Data Commons* (http://datacommons.cyverse.org/). The data are freely accessible at http://datacommons.cyverse.org/browse/iplant/home/shared/eht/2019/GRMHDCodeComparisonProject. Please consult the README.txt and LICENSE.txt for usage instructions.

R.N. thanks the National Science Foundation (NSF; grants OISE-1743747, AST-1816420) and acknowledges computational support from the NSF via XSEDE resources (grant TG-AST080026N). L.D.Z. acknowledges support from the PRIN-MIUR project Multi-scale Simulations of High-Energy Astrophysical Plasmas (Prot. 2015L5EE2Y) and from the INFN—TEONGRAV initiative. C.J.W. made use of the Extreme Science and Engineering Discovery Environment (XSEDE) Comet at the San Diego Supercomputer Center through allocation AST170012. The H-AMR high-resolution simulation was made possible by NSF PRAC award Nos. 1615281 and OAC-1811605 at the Blue Waters sustained-petascale computing project and supported in part under grant No. NSF PHY-1125915 (PI A. Tchekhovskoy). K.C. and S.M. are supported by the Netherlands Organization for Scientific Research (NWO) VICI grant (No. 639.043.513); M.L. is supported by the NWO Spinoza Prize (PI M.B.M. van der Klis). The HARM-Noble simulations were made possible by NSF PRAC award No. NSF OAC-1515969, OAC-1811228 at the Blue Waters sustained-petascale computing project, and supported in part under grant No. NSF PHY-1125915. The BHAC CKS-GRMHD simulations were performed on the Dutch National Supercomputing cluster Cartesius and are funded by the NWO computing grant 16431. S.C.N. was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center administered by USRA through a contract with NASA. Y.M., H.O., O.P., and L.R. acknowledge support from the ERC synergy grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (grant No. 610058). M.B. acknowledges support from the European Research Council (grant No. 715368—MagBURST) and from the Gauss Centre for Supercomputing e.V. (www.Gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (www.lrz.de). P.C.F. was supported by NSF grant AST-1616185 and used resources from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant No. ACI-1548562. Work by P.A. was performed in part under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344.

The authors of the present paper further thank the following organizations and programs: the Academy of Finland (projects 274477, 284495, 312496); the Advanced European Network of E-infrastructures for Astronomy with the SKA (AENEAS) project, supported by the European Commission Framework Programme Horizon 2020 Research and Innovation action under grant agreement 731016; the Alexander von Humboldt Stiftung; the Black Hole Initiative at Harvard University, through a grant (60477) from the John Templeton Foundation; the China Scholarship Council; Comisión Nacional de Investigación Científica y Tecnológica (CONICYT, Chile, via PIA ACT172033, Fondecyt 1171506, BASAL AFB-170002, ALMA-conicyt 31140007); Consejo Nacional de Ciencia y Tecnología (CONACYT, Mexico, projects 104497, 275201, 279006, 281692); the Delaney Family via the Delaney Family John A. Wheeler Chair at Perimeter Institute; Dirección General de Asuntos del Personal Académico-Universidad Nacional Autónoma de México (DGAPA-UNAM, project IN112417); the European Research Council Synergy Grant “BlackHoleCam: Imaging the Event Horizon of Black Holes” (grant 610058); the Generalitat Valenciana postdoctoral grant APOSTD/2018/177; the Gordon and Betty Moore 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 for Astronomy and Astrophysics at the Universities of Bonn and Cologne; the Jansky Fellowship program of the National Radio Astronomy Observatory (NRAO); the Japanese Government (Monbukagakusho: MEXT) Scholarship; the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for JSPS Research Fellowship (JP17J08829); the Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS, grants QYZDJ-SSW-SLH057, QYZDJ-SSW-SYS008); the Leverhulme Trust Early Career Research Fellowship; the Max-Planck-Gesellschaft (MPG); the Max Planck Partner Group of the MPG and the CAS; the MEXT/JSPS KAKENHI (grants 18KK0090, JP18K13594, JP18K03656, JP18H03721, 18K03709, 18H01245, 25120007); the MIT International Science and Technology Initiatives (MISTI) Funds; the Ministry 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 Aeronautics and Space Administration (NASA, Fermi Guest Investigator grant 80NSSC17K0649); the National Institute of Natural Sciences (NINS) of Japan; the National Key Research and Development Program of China (grant 2016YFA0400704, 2016YFA0400702); the National Science Foundation (NSF, grants AST-0096454, AST-0352953, AST-0521233, AST-0705062, AST-0905844, AST-0922984, AST-1126433, AST-1140030, DGE-1144085, AST-1207704, AST-1207730, AST-1207752, MRI-1228509, OPP-1248097, AST-1310896, AST-1312651, AST-1337663, AST-1440254, AST-1555365, AST-1715061, AST-1615796, AST-1716327, OISE-1743747, AST-1816420); the Natural Science Foundation of China (grants 11573051, 11633006, 11650110427, 10625314, 11721303, 11725312); the Natural Sciences and Engineering Research Council of Canada (NSERC, including a Discovery Grant and the NSERC Alexander Graham Bell Canada Graduate Scholarships Doctoral Program); the National Youth Thousand Talents Program of China; the National Research Foundation of Korea (the Global PhD Fellowship Grant: grants NRF-2015H1A2A1033752, 2015-R1D1A1A01056807, the Korea
Research Fellowship Program: NRF-2015H1D3A1066561); the Netherlands Organization for Scientific Research (NWO) VICI award (grant 639.043.513) and Spinoza Prize SPI 78-409; the New Scientific Frontiers with Precision Radio Interferometry Fellowship awarded by the South African Radio Astronomy Observatory (SARAO), which is a facility of the National Research Foundation (NRF), an agency of the Department of Science and Technology (DST) of South Africa; the Onsala Space Observatory (OSO) national infrastructure, for the provisioning of its facilities/observational support (OSO receives funding through the Swedish Research Council under grant 2017-00648); the Perimeter Institute for Theoretical Physics (research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development, and by the Province of Ontario through the Ministry of Research, Innovation and Science); the Russian Science Foundation (grant 17-12-01029); the Spanish Ministerio de Economía y Competitividad (grants AYA2015-63939-C2-1-P, AYA2016-80889-P); the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709); the Toray Science Foundation; the US Department of Energy (USDOE) through the Los Alamos National Laboratory (operated by Triad National Security, LLC, for the National Nuclear Security Administration of the USDOE (Contract 89233218CNA000001)); the Italian Ministero dell'Istruzione Università e Ricerca through the grant Progetti Premiali 2012-iALMA (CUP C52113000140001); the European Union's Horizon 2020 research and innovation programme under grant agreement No. 730562 RadioNet; ALMA North America Development Fund; the Academia Sinica; and Chandra TM6-17006X. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. The simulations were performed in part on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, and on the HazelHen cluster at the HLRS in Stuttgart. This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca), and Compute Canada (http://www.computecanada.ca).

Software: Athena++ (White et al. 2016); BHAC (Porth et al. 2017); Cosmos++, (Anninos et al. 2005; Fragile et al. 2012, 2014); ECHO (Del Zanna et al. 2007); H-AMR (Liska et al. 2018a; Chatterjee et al. 2019); HARM-Noble (Gammie et al. 2003; Noble et al. 2006, 2009); iharm3D (Gammie et al. 2003;

Noble et al. 2006, 2009); IllinoisGRMHD (Etienne et al. 2015); KORAL, (Sądowski et al. 2013, 2014).

## Appendix A
### Comparison with Linear Reconstruction

Using H-AMR, KORAL, and Athena++, additional simulations with piecewise linear reconstruction (PLM) were performed. To be more specific, in order of increasing numerical diffusivity, Athena++ used the modified “van Leer” limiter of Mignone (2014), H-AMR the MC<sub>β</sub> limiter (van Leer 1977) with β = 1.5, and KORAL MC<sub>β</sub> with β = 1.9, where MC<sub>β</sub> reduces to the most diffusive “minmod-limiter” in the case of β = 1 and to “MC” for β = 2. It has been observed before that the higher order reconstruction of PPM compensates for the large dissipation in the LLF and HLL approximate Riemann solvers which are applied throughout this work. For example, Flock et al. (2010) noted that with PPM reconstruction, the HLL and LLF methods properly resolve the growth of the linear MRI with 10 cells per mode as opposed to 16 cells for the PLM reconstruction.

For completeness, we list the values of horizon-penetrating fluxes from these runs in Table 3. It is striking that even with a resolution as high as 192³, the PLM runs can give a very different answer: the time-averaged accretion rate in the H-AMR and Athena++ PLM runs falls short of their PPM counterpart by a factor of 2 to almost 6. For KORAL on the other hand, we obtain a more consistent behavior between the reconstruction techniques.

Apart from a slightly higher midplane resolution in the KORAL case, there are significant algorithmic differences between the codes. Whereas the implementation of KORAL follows closely the one by Gammie et al. (2003), where a cell-centered representation of the magnetic field along with the arithmetic averaging of the electric fields (ACT) is chosen, H-AMR and Athena++ both implement a staggered UCT scheme following Gardiner & Stone (2005). In contrast to ACT methods, the UCT scheme reproduces the correct solution for plane-parallel grid-aligned fields; however, the amount of dissipation of the scheme is effectively doubled.

Inspection of the MRI quality factors for H-AMR as shown in Figure 23 reveals that the turbulence is not sufficiently resolved in the 96–192 PLM runs. Not only does $\langle Q_0 \rangle$ stay below 10 during the entire run, it is also decreasing over time, indicating net dissipation of magnetic energy.

These results show that even when the overall order of the scheme is fixed to second order, the PPM method can significantly reduce the dissipation of the numerical integration and improve the results for a given resolution.
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Domain-averaged and density-weighted quality factors with the H-AMR code using linear (PLM) and parabolic reconstruction (PPM). The values are consistent with the ones reported in Section 5.6 with the BHAC code. All PLM runs remain below ~8 and decrease further in the late time evolution.</figcaption></figure>

<table><caption>**Table 3**

Quantifications—Time Series Data, Linear Reconstruction</caption><thead><tr><th>N<sub>θ</sub></th><th>Code</th><th>M&#775;<sub>Peak</sub></th><th>⟨M&#775;⟩</th><th>Φ<sub>BH</sub>/√M&#775;<sub>Peak</sub></th><th>⟨Φ<sub>BH</sub>/√M&#775;⟩</th><th>L&#775;/M&#775;<sub>Peak</sub></th><th>⟨L&#775;/M&#775;⟩</th><th>E&#775;/M&#775;<sub>Peak</sub></th><th>⟨E&#775;/M&#775;⟩</th></tr></thead><tbody><tr><td rowspan="4">96</td><td>Athena++</td><td>0.209</td><td>0.049 ± 0.026</td><td>4.762</td><td>3.308 ± 0.732</td><td>1.929</td><td>1.512 ± 0.225</td><td>0.332</td><td>0.145 ± 0.024</td></tr><tr><td>H-AMR</td><td>0.201</td><td>0.044 ± 0.015</td><td>1.341</td><td>0.914 ± 0.15</td><td>1.931</td><td>1.473 ± 0.187</td><td>0.446</td><td>0.241 ± 0.049</td></tr><tr><td>KORAL</td><td>0.863</td><td>0.09 ± 0.077</td><td>5.247</td><td>3.226 ± 1.022</td><td>2.146</td><td>1.734 ± 0.286</td><td>0.328</td><td>0.033 ± 0.015</td></tr><tr><td>max/min</td><td>4.296</td><td>2.062</td><td>3.913</td><td>3.62</td><td>1.112</td><td>1.177</td><td>1.358</td><td>7.258</td></tr><tr><td rowspan="4">128</td><td>Athena++</td><td>0.266</td><td>0.035 ± 0.011</td><td>6.553</td><td>3.845 ± 0.664</td><td>1.918</td><td>1.387 ± 0.272</td><td>0.371</td><td>0.167 ± 0.041</td></tr><tr><td>H-AMR</td><td>0.244</td><td>0.089 ± 0.023</td><td>1.211</td><td>0.707 ± 0.105</td><td>1.907</td><td>1.738 ± 0.099</td><td>0.336</td><td>0.098 ± 0.012</td></tr><tr><td>KORAL</td><td>0.65</td><td>0.179 ± 0.066</td><td>4.671</td><td>0.881 ± 0.143</td><td>2.245</td><td>2.118 ± 0.046</td><td>0.332</td><td>0.041 ± 0.013</td></tr><tr><td>max/min</td><td>2.665</td><td>5.165</td><td>5.41</td><td>5.439</td><td>1.177</td><td>1.527</td><td>1.117</td><td>4.049</td></tr><tr><td rowspan="4">192</td><td>Athena++</td><td>0.42</td><td>0.069 ± 0.038</td><td>3.108</td><td>1.73 ± 0.311</td><td>2.064</td><td>1.947 ± 0.047</td><td>0.332</td><td>0.094 ± 0.007</td></tr><tr><td>H-AMR</td><td>0.612</td><td>0.04 ± 0.018</td><td>1.349</td><td>0.916 ± 0.192</td><td>2.022</td><td>1.487 ± 0.22</td><td>0.336</td><td>0.194 ± 0.045</td></tr><tr><td>KORAL</td><td>0.789</td><td>0.202 ± 0.036</td><td>3.357</td><td>0.776 ± 0.059</td><td>2.084</td><td>2.025 ± 0.02</td><td>0.337</td><td>0.065 ± 0.008</td></tr><tr><td>max/min</td><td>1.88</td><td>5.027</td><td>2.489</td><td>2.23</td><td>1.031</td><td>1.362</td><td>1.016</td><td>2.967</td></tr></tbody></table>

Note. Quantities in angular brackets ⟨⟩ denote time averages in the interval $t \in [5000, 10,000]M$ with error given by the standard deviation.

## Appendix B

### Run-to-run Variation

Due to the chaotic nature of the turbulence at late times, different initial random perturbations can accumulate to large differences in the realization of the dynamics. Likewise, because the compiler version and optimization (e.g., order of execution) influences the roundoff error, a similar effect can be observed if the same physical scenario is run on two machines with differing compiler and runtime configurations. In order to judge the impact this makes compared to the intercode

differences, the KORAL runs were repeated on two clusters: Harvard's Odyssey machine and Stampede2 of the Texas Advanced Computing Center. Here, both the initial random perturbations and the machine architecture differed.

Again, we list the standard quantifications of the time-series data in Table 4. It shows that with a few exceptions of the peak values (and normalized magnetic flux), the differences are typically in the low percent regime. Hence, the differences promoted by the various adopted algorithms are larger than the run to run variation of the KORAL code.
<table><caption>**Table 4**
Quantifications—Time Series Data, Run to Run Variation</caption><thead><tr><th>N<sub>θ</sub></th><th>Code</th><th>M&#775;<sub>Peak</sub></th><th>⟨M&#775;⟩</th><th>Φ<sub>BH</sub>/√M&#775;<sub>Peak</sub></th><th>⟨Φ<sub>BH</sub>/√M&#775;⟩</th><th>L&#775;/M&#775;<sub>Peak</sub></th><th>⟨L&#775;/M&#775;⟩</th><th>E&#775;/M&#775;<sub>Peak</sub></th><th>⟨E&#775;/M&#775;⟩</th></tr></thead><tbody><tr><td rowspan="3">96</td><td>KORAL Odyssey</td><td>0.821</td><td>0.408 ± 0.129</td><td>3.029</td><td>0.791 ± 0.362</td><td>2.126</td><td>2.031 ± 0.038</td><td>0.33</td><td>0.061 ± 0.008</td></tr><tr><td>KORAL Stampede2</td><td>0.742</td><td>0.418 ± 0.112</td><td>2.818</td><td>1.068 ± 0.16</td><td>2.105</td><td>2.014 ± 0.043</td><td>0.33</td><td>0.061 ± 0.01</td></tr><tr><td>max/min</td><td>1.106</td><td>1.024</td><td>1.075</td><td>1.349</td><td>1.01</td><td>1.008</td><td>1.0</td><td>1.009</td></tr><tr><td rowspan="3">128</td><td>KORAL Odyssey</td><td>0.859</td><td>0.318 ± 0.102</td><td>2.899</td><td>1.279 ± 0.306</td><td>2.073</td><td>1.978 ± 0.049</td><td>0.331</td><td>0.066 ± 0.007</td></tr><tr><td>KORAL Stampede2</td><td>1.217</td><td>0.361 ± 0.099</td><td>2.988</td><td>1.171 ± 0.17</td><td>2.051</td><td>1.963 ± 0.032</td><td>0.331</td><td>0.062 ± 0.008</td></tr><tr><td>max/min</td><td>1.417</td><td>1.134</td><td>1.031</td><td>1.093</td><td>1.011</td><td>1.008</td><td>1.0</td><td>1.053</td></tr><tr><td rowspan="3">192</td><td>KORAL Odyssey</td><td>1.067</td><td>0.29 ± 0.132</td><td>2.459</td><td>1.254 ± 0.231</td><td>2.036</td><td>1.953 ± 0.023</td><td>0.331</td><td>0.056 ± 0.007</td></tr><tr><td>KORAL Stampede2</td><td>0.933</td><td>0.327 ± 0.116</td><td>2.58</td><td>1.013 ± 0.1</td><td>2.047</td><td>1.981 ± 0.025</td><td>0.331</td><td>0.056 ± 0.008</td></tr><tr><td>max/min</td><td>1.144</td><td>1.128</td><td>1.049</td><td>1.237</td><td>1.006</td><td>1.014</td><td>1.0</td><td>1.014</td></tr></tbody></table>

Note. Quantities in angular brackets ⟨⟩ denote time averages in the interval $t \in [5000, 10,000]M$ with error given by the standard deviation.

## Appendix C
Rescaled Disk Profiles

With the scale freedom allowed by the test-fluid assumption, the density, pressure, and magnetic fields of a given simulation can in principle be rescaled by a constant factor (respectively, its square root for the magnetic fields), for example, to perform spectral fitting of observations. We here exploit this freedom and rescale the simulations to a reference case for which we use the $N_{\theta} = 1056$ H-AMR simulation. In particular, we match the

density at the initial density maximum of the disk, $r = 12M$. The result of this procedure is exemplified in Figure 24 for the high-resolution data. It is obvious that this improves the overlap in the inner regions $r < 12M$ in densities and pressure, but the spread in magnetic fields is worsened. This is merely a consequence of the nonconvergence of magnetization, for if the density is made to match, the spread in magnetic field increases with the reference run having roughly a magnetization twice as large as most 192³ runs.
<figure><img src="image_26.png" /><figcaption>Figure 24. Disk profiles rescaled to match a “reference solution” at the point $r = 12M$. We take the H-AMR very high-resolution run for reference.</figcaption></figure>

## Appendix D
### Measuring the Maxwell Stress

Maxwell stresses play an important role in disk accretion; however, their definition is ambiguous in general relativity, due to the difficulty of identifying spatial directions in one frame with the ones in another. In particular, due to the local nature of the magnetorotational instability, the physical interpretation is best guided by stresses measured in a frame corotating with the disk. Hence, to relate to the classical disk model of Novikov & Thorne (1973), Krolik et al. (2005) devised a set of comoving tetrads most closely preserving the Boyer–Lindquist azimuthal and radial directions. The expressions for the basis are written explicitly in Appendix A of Beckwith et al. (2008) and need not be reproduced here. Note that in order to use them, the four-vectors $b^{\mu}$ (here assumed to be in Kerr–Schild coordinates) first need to be converted to $b'^{\mu}$ in Boyer–Lindquist coordinates.

To validate the simplified approach taken in Section 5.6.1, we here compare Maxwell stresses obtained with the following two prescriptions:

$$ \langle w^{r\phi} \rangle (r, t) := - \frac{\int_{0}^{2\pi} d\phi \int_{\theta_{\min}}^{\theta_{\max}} d\theta \sqrt{-g} \sqrt{g_{rr}} b^r \sqrt{g_{\phi\phi}} b^{\phi}}{\int_{0}^{2\pi} d\phi \int_{\theta_{\min}}^{\theta_{\max}} d\theta \sqrt{-g}}, \quad (27) $$

$$ \langle w^{(\bar{r})}(\bar{\phi}) \rangle (r, t) := - \frac{\int_{\bar{\phi}_{\min}}^{\bar{\phi}_{\max}} dx^{(\bar{\phi})} \int_{\bar{\theta}_{\min}}^{\bar{\theta}_{\max}} dx^{(\bar{\theta})} b^{(\bar{r})} b^{(\bar{\phi})}}{\int_{\bar{\phi}_{\min}}^{\bar{\phi}_{\max}} dx^{(\bar{\phi})} \int_{\bar{\theta}_{\min}}^{\bar{\theta}} dx^{(\bar{\theta})}}, \quad (28) $$

where $\langle w^{r\phi} \rangle (r, t)$ is the Kerr–Schild frame measurement used in Section 5.6.1 and $\langle w^{(\bar{r})}(\bar{\phi}) \rangle (r, t)$ is the stress in a frame comoving with the local fluid velocity as in Krolik et al. (2005). Integration is carried out over the equatorial wedge $\theta_{KS} \in [\pi/3, 2\pi/3]$ and the full azimuthal range. Note that we have taken the fluid-frame integration over the comoving coordinate increments $dx^{(\bar{r})}$ which result from the transformation of $dx^{\mu} := (0, \Delta r, \Delta\theta, \Delta\phi)$ as appropriate (see the discussion in Noble et al. 2010).

The resulting stress profiles and volume-integrated time series (nondimensionalized by the time- and volume- averaged total pressure $\langle P_{tot} \rangle$) are shown in Figure 25 for the five runs with the BHAC code. Significant departures of the two measurements only occur within $r_{KS} \simeq 2M$, where the fluid-frame stress flattens out and shows indications for a maximum for the two highest resolution runs. For reference, the ISCO is located at $2.04M$, which roughly coincides with the change of slope. The overall stress differs only by a few percent between diagnostics (27) and (28) as demonstrated in the time series. This is much smaller than differences between resolutions. We hence conclude that the simple coordinate-frame measurement is appropriate for our purpose of studying the convergence and robustness of $\alpha$ as carried out in Section 5.6.1. It is important, however, to stress that the good correspondence between fluid- and coordinate-frame stresses does not necessarily hold for all values of the BH spin. In the Schwarzschild case, for example, while the fluid-frame stress drops to zero at the horizon (Krolik et al. 2005; Noble et al. 2010), the coordinate-frame stress remains monotonic.
<figure><img src="image_27.png" /><figcaption>**Figure 25.** Comparison of the Maxwell stress in different frames. Solid curves show the fluid-frame measurement and dotted curves the Kerr–Schild coordinate-frame stress. There is very good correspondence between the profiles up to ~2M. Consequently, the time evolution of the volume-averaged stresses agrees within a few percent.</figcaption></figure>

## ORCID iDs

Oliver Porth https://orcid.org/0000-0002-4584-2557

Ramesh Narayan https://orcid.org/0000-0002-1919-2730

Charles F. Gammie https://orcid.org/0000-0001-7451-8935

Yosuke Mizuno https://orcid.org/0000-0002-8131-6730

Matteo Bugli https://orcid.org/0000-0002-7834-0422

Chi-kwan Chan https://orcid.org/0000-0001-6337-6126

Jordy Davelaar https://orcid.org/0000-0002-2685-2434

P. Chris Fragile https://orcid.org/0000-0002-5786-186X

Bernard J. Kelly https://orcid.org/0000-0002-3326-4454

Sera Markoff https://orcid.org/0000-0001-9564-0876

Scott C. Noble https://orcid.org/0000-0003-3547-8306

Héctor Olivares https://orcid.org/0000-0001-6833-7580

Ben Prather https://orcid.org/0000-0002-0393-7734

Luciano Rezzolla https://orcid.org/0000-0002-1330-7103

Benjamin R. Ryan https://orcid.org/0000-0001-8939-4461

Ziri Younsi https://orcid.org/0000-0001-9283-1191

Kazunori Akiyama https://orcid.org/0000-0002-9475-4254

Antxon Alberdi https://orcid.org/0000-0002-9371-1033

Rebecca Azulay https://orcid.org/0000-0002-2200-5393

Anne-Kathrin Baczko https://orcid.org/0000-0003-3090-3975

Mislav Baloković https://orcid.org/0000-0003-0476-6647

John Barrett https://orcid.org/0000-0002-9290-0764

Lindy Blackburn https://orcid.org/0000-0002-9030-642X

Katherine L. Bouman https://orcid.org/0000-0003-0077-4367

Geoffrey C. Bower https://orcid.org/0000-0003-4056-9982

Christiaan D. Brinkerink https://orcid.org/0000-0002-2322-0749

Roger Brissenden https://orcid.org/0000-0002-2556-0894

Silke Britzen https://orcid.org/0000-0001-9240-6734

Avery E. Broderick https://orcid.org/0000-0002-3351-760X

Do-Young Byun https://orcid.org/0000-0003-1157-4109

Andrew Chael https://orcid.org/0000-0003-2966-6220

Shami Chatterjee https://orcid.org/0000-0002-2878-1502

Ilje Cho https://orcid.org/0000-0001-6083-7521

Pierre Christian https://orcid.org/0000-0001-6820-9941

John E. Conway https://orcid.org/0000-0003-2448-9181

James M. Cordes https://orcid.org/0000-0002-4049-1882

Geoffrey, B. Crew https://orcid.org/0000-0002-2079-3189

Yuzhu Cui https://orcid.org/0000-0001-6311-4345

Mariafelicia De Laurentis https://orcid.org/0000-0002-9945-682X

Roger Deane https://orcid.org/0000-0003-1027-5043

Jessica Dempsey https://orcid.org/0000-0003-1269-9667

Gregory Desvignes https://orcid.org/0000-0003-3922-4055

Sheperd S. Doeleman https://orcid.org/0000-0002-9031-0904

Ralph P. Eatough https://orcid.org/0000-0001-6196-4135

Heino Falcke https://orcid.org/0000-0002-2526-6724

Vincent L. Fish https://orcid.org/0000-0002-7128-9345

Ed Fomalont https://orcid.org/0000-0002-9036-2747

Raquel Fraga-Encinas https://orcid.org/0000-0002-5222-1361

José L. Gómez https://orcid.org/0000-0003-4190-7613

Peter Galison https://orcid.org/0000-0002-6429-3872

Boris Georgiev https://orcid.org/0000-0002-3586-6424

Roman Gold https://orcid.org/0000-0003-2492-1966

Minfeng Gu (顾敏峰) https://orcid.org/0000-0002-4455-6946

Mark Gurwell https://orcid.org/0000-0003-0685-3621

Kazuhiro Hada https://orcid.org/0000-0001-6906-772X

Ronald Hesper https://orcid.org/0000-0003-1918-6098

Luis C. Ho (何子山) https://orcid.org/0000-0001-6947-5846

Mareki Honma https://orcid.org/0000-0003-4058-9000

Chih-Wei L. Huang https://orcid.org/0000-0001-5641-3953

Shiro Ikeda https://orcid.org/0000-0002-2462-1448

Sara Issaoun https://orcid.org/0000-0002-5297-921X

David J. James https://orcid.org/0000-0001-5160-4486

Michael Janssen https://orcid.org/0000-0001-8685-6544

Britton Jeter https://orcid.org/0000-0003-2847-1712

Wu Jiang (江悟) https://orcid.org/0000-0001-7369-3539

Michael D. Johnson https://orcid.org/0000-0002-4120-3029

Svetlana Jorstad https://orcid.org/0000-0001-6158-1708

Taehyun Jung https://orcid.org/0000-0001-7003-8643

Mansour Karami https://orcid.org/0000-0001-7387-9333

Ramesh Karuppusamy https://orcid.org/0000-0002-5307-2919

Tomohisa Kawashima https://orcid.org/0000-0001-8527-0496

Garrett K. Keating https://orcid.org/0000-0002-3490-146X

Mark Kettenis https://orcid.org/0000-0002-6156-5617

Jae-Young Kim https://orcid.org/0000-0001-8229-7183

Junhan Kim https://orcid.org/0000-0002-4274-9373

Motoki Kino https://orcid.org/0000-0002-2709-7338

Jun Yi Koay https://orcid.org/0000-0002-7029-6658
Patrick, M. Koch https://orcid.org/0000-0003-2777-5861

Shoko Koyama https://orcid.org/0000-0002-3723-3372

Michael Kramer https://orcid.org/0000-0002-4175-2271

Carsten Kramer https://orcid.org/0000-0002-4908-4925

Thomas P. Krichbaum https://orcid.org/0000-0002-4892-9586

Tod R. Lauer https://orcid.org/0000-0003-3234-7247

Sang-Sung Lee https://orcid.org/0000-0002-6269-594X

Yan-Rong Li (李彦荣) https://orcid.org/0000-0001-5841-9179

Zhiyuan Li (李志远) https://orcid.org/0000-0003-0355-6437

Kuo Liu https://orcid.org/0000-0002-2953-7376

Elisabetta Liuzzo https://orcid.org/0000-0003-0995-5201

Laurent Loinard https://orcid.org/0000-0002-5635-3345

Ru-Sen Lu (路如森) https://orcid.org/0000-0002-7692-7967

Nicholas R. MacDonald https://orcid.org/0000-0002-6684-8691

Jirong Mao (毛基荣) https://orcid.org/0000-0002-7077-7195

Daniel P. Marrone https://orcid.org/0000-0002-2367-1080

Alan P. Marscher https://orcid.org/0000-0001-7396-3332

Iván Martí-Vidal https://orcid.org/0000-0003-3708-9611

Lynn D. Matthews https://orcid.org/0000-0002-3728-8082

Lia Medeiros https://orcid.org/0000-0003-2342-6728

Karl M. Menten https://orcid.org/0000-0001-6459-0669

Izumi Mizuno https://orcid.org/0000-0002-7210-6264

James M. Moran https://orcid.org/0000-0002-3882-4414

Kotaro Moriyama https://orcid.org/0000-0003-1364-3761

Monika Moscibrodzka https://orcid.org/0000-0002-4661-6332

Cornelia Müller https://orcid.org/0000-0002-2739-2994

Hiroshi Nagai https://orcid.org/0000-0003-0292-3645

Neil M. Nagar https://orcid.org/0000-0001-6920-662X

Masanori Nakamura https://orcid.org/0000-0001-6081-2420

Iniyan Natarajan https://orcid.org/0000-0001-8242-4373

Roberto Neri https://orcid.org/0000-0002-7176-4046

Chunchong Ni https://orcid.org/0000-0003-1361-5699

Aristeidis Noutsos https://orcid.org/0000-0002-4151-3860

Daniel C. M. Palumbo https://orcid.org/0000-0002-7179-3816

Ue-Li Pen https://orcid.org/0000-0003-2155-9578

Dominic W. Pesce https://orcid.org/0000-0002-5278-9221

Jorge A. Preciado-López https://orcid.org/0000-0002-4146-0113

Dimitrios Psaltis https://orcid.org/0000-0003-4058-2837

Hung-Yi Pu https://orcid.org/0000-0001-9270-8812

Venkatessh Ramakrishnan https://orcid.org/0000-0002-9248-086X

Ramprasad Rao https://orcid.org/0000-0002-1407-7944

Alexander W. Raymond https://orcid.org/0000-0002-5779-4767

Bart Ripperda https://orcid.org/0000-0002-7301-3908

Freek Roelofs https://orcid.org/0000-0001-5461-3687

Eduardo Ros https://orcid.org/0000-0001-9503-4892

Mel Rose https://orcid.org/0000-0002-2016-8746

Alan L. Roy https://orcid.org/0000-0002-1931-0135

Chet Ruszczyk https://orcid.org/0000-0001-7278-9707

Kazi L. J. Rygl https://orcid.org/0000-0003-4146-9043

David Sánchez-Arguelles https://orcid.org/0000-0002-7344-9920

Mahito Sasada https://orcid.org/0000-0001-5946-9960

Tuomas Savolainen https://orcid.org/0000-0001-6214-1085

Lijing Shao https://orcid.org/0000-0002-1334-8853

Zhiqiang Shen (沈志强) https://orcid.org/0000-0003-3540-8746

Des Small https://orcid.org/0000-0003-3723-5404

Bong Won Sohn https://orcid.org/0000-0002-4148-8378

Jason SooHoo https://orcid.org/0000-0003-1938-0720

Fumie Tazaki https://orcid.org/0000-0003-0236-0600

Paul Tiede https://orcid.org/0000-0003-3826-5648

Remo P. J. Tilanus https://orcid.org/0000-0002-6514-553X

Michael Titus https://orcid.org/0000-0001-9001-3275

Kenji Toma https://orcid.org/0000-0002-7114-6010

Pablo Torne https://orcid.org/0000-0001-8700-6058

Sascha Trippe https://orcid.org/0000-0003-0465-1559

Ilse van Bemmel https://orcid.org/0000-0001-5473-2950

Huib Jan van Langevelde https://orcid.org/0000-0002-0230-5946

Daniel R. van Rossum https://orcid.org/0000-0001-7772-6131

John Wardle https://orcid.org/0000-0002-8960-2942

Jonathan Weintroub https://orcid.org/0000-0002-4603-5204

Robert Wharton https://orcid.org/0000-0002-7416-5209

Maciek Wielgus https://orcid.org/0000-0002-8635-4242

George N. Wong https://orcid.org/0000-0001-6952-2147

Qingwen Wu (吴庆文) https://orcid.org/0000-0003-4773-4987

Ken Young https://orcid.org/0000-0002-3666-4920

André Young https://orcid.org/0000-0003-0000-2682

Feng Yuan (袁峰) https://orcid.org/0000-0003-3564-6437

Ye-Fei Yuan (袁业飞) https://orcid.org/0000-0002-7330-4756

J. Anton Zensus https://orcid.org/0000-0001-7470-3321

Guangyao Zhao https://orcid.org/0000-0002-4417-1659

Shan-Shan Zhao https://orcid.org/0000-0002-9774-3606

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, PhRvL, 116, 241103

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, PhRvL, 119, 161101

Akiyama, K., Lu, R.-S., Fish, V. L., et al. 2015, ApJ, 807, 150

Anninos, P., Bryant, C., Fragile, P. C., et al. 2017, ApJS, 231, 17

Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723

Anninos, P., Fragile, P. C., Wilson, J., & Murray, S. D. 2012, ApJ, 759, 132

Antón, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296

Armitage, P. J., & Reynolds, C. S. 2003, MNRAS, 341, 1041

Asada, K., & Nakamura, M. 2012, ApJL, 745, L28

Bacchini, F., Ripperda, B., Chen, A. Y., & Sironi, L. 2018, ApJS, 237, 6

Bacchini, F., Ripperda, B., Porth, O., & Sironi, L. 2019, ApJS, 240, 40

Baganoff, F. K., Bautz, M. W., Brandt, W. N., et al. 2001, Natur, 413, 45

Baiotti, L., Hawke, I., Montero, P. J., et al. 2005, PhRvD, 71, 024035

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A., & Hawley, J. F. 1998, RvMP, 70, 1

Balbus, S. A., & Henri, P. 2008, ApJ, 674, 408

Balsara, D. S., & Spicer, D. 1999, JCoPh, 148, 133

Beckwith, K., Hawley, J. F., & Krolik, J. H. 2008, MNRAS, 390, 21

Beskin, V. S., Istomin, Y. N., & Parev, V. I. 1992, SvA, 36, 642

Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1976, Ap&SS, 42, 401

Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433

Bodo, G., Cattaneo, F., Mignone, A., & Rossi, P. 2014, ApJL, 787, L13

Bowen, D. B., Mewes, V., Campanelli, M., et al. 2018, ApJL, 853, L17

Bower, G. C., Markoff, S., Dexter, J., et al. 2015, ApJ, 802, 69

Broderick, A. E., & Tchekhovskoy, A. 2015, ApJ, 809, 97

Bucciantini, N., & Del Zanna, L. 2011, A&A, 528, A101

Bucciantini, N., & Del Zanna, L. 2013, MNRAS, 428, 71

Bugli, M., Del Zanna, L., & Bucciantini, N. 2014, MNRAS, 440, L41

Bugli, M., Guilet, J., Müller, E., et al. 2018, MNRAS, 475, 108

Carnemzind, M. 1986, A&A, 162, 32

Chael, A., Narayan, R., & Johnson, M. D. 2019, MNRAS, 486, 2873

Chael, A. A., Narayan, R., & Sałowski, A. 2017, MNRAS, 470, 2367

Chan, C.-k., Psaltis, D., Özel, F., et al. 2015, ApJ, 812, 103

Chandra, M., Foucart, F., & Gammie, C. F. 2017, ApJ, 837, 92
Chandra, M., Gammie, C. F., Foucart, F., & Quataert, E. 2015, ApJ, 810, 162

Chatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff, S. B. 2019, arXiv:1904.03243

Colella, P., & Woodward, P. R. 1984, JCoPh, 54, 174

Davelaar, J., Mościbrodzka, M., Bronzwaer, T., & Falcke, H. 2018, A&A, 612, A34

Davelaar, J., Olivares, H., Porth, O., et al. 2019, A&A, submitted (arXiv:1906.10065)

Davis, S. W., Stone, J. M., & Pessah, M. E. 2010, ApJ, 713, 52

De Villiers, J.-P., & Hawley, J. F. 2003, ApJ, 589, 458

Del Zanna, L., & Bucciantini, N. 2018, MNRAS, 479, 657

Del Zanna, L., Papini, E., Landi, S., Bugli, M., & Bucciantini, N. 2016, MNRAS, 460, 3753

Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11

Dexter, J., Agol, E., & Fragile, P. C. 2009, ApJL, 703, L142

Dexter, J., Agol, E., Fragile, P. C., & McKinney, J. C. 2010, ApJ, 717, 1092

Dexter, J., Kelly, B., Bower, G. C., et al. 2014, MNRAS, 442, 2797

Dexter, J., McKinney, J. C., & Agol, E. 2012, MNRAS, 421, 1517

Dhang, P., & Sharma, P. 2019, MNRAS, 482, 848

Dibi, S., Drappeau, S., Fragile, P. C., Markoff, S., & Dexter, J. 2012, MNRAS, 426, 1928

Do, T., Ghez, A. M., Morris, M. R., et al. 2009, ApJ, 691, 1021

Dodds-Eden, K., Gillessen, S., Fritz, T. K., et al. 2011, ApJ, 728, 37

Dodds-Eden, K., Porquet, D., Trap, G., et al. 2009, ApJ, 698, 676

Doeleman, S. S., Fish, V. L., Schenck, D. E., et al. 2012, Sci, 338, 355

Doeleman, S. S., Weintroub, J., Rogers, A. E. E., et al. 2008, Natur, 455, 78

Dolence, J. C., Gammie, C. F., Shiokawa, H., & Noble, S. C. 2012, ApJL, 746, L10

Duez, M. D., Liu, Y. T., Shapiro, S. L., & Stephens, B. C. 2005, PhRvD, 72, 024028

Etienne, Z. B., Liu, Y. T., Paschalidis, V., & Shapiro, S. L. 2012a, PhRvD, 85, 064029

Etienne, Z. B., Liu, Y. T., & Shapiro, S. L. 2006, PhRvD, 74, 044030

Etienne, Z. B., Liu, Y. T., & Shapiro, S. L. 2010, PhRvD, 82, 084031

Etienne, Z. B., Paschalidis, V., Haas, R., Moesta, P., & Shapiro, S. L. 2015, CQGra, 32, 175009

Etienne, Z. B., Paschalidis, V., Liu, Y. T., & Shapiro, S. L. 2012b, PhRvD, 85, 024013

Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659

Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJL, 875, L1

Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019b, ApJL, 875, L5

Fambri, F., Dumbser, M., Köppel, S., Rezzolla, L., & Zanotti, O. 2018, MNRAS, 477, 4543

Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962

Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26

Font, J. A. 2008, LRR, 11, 7

Font, J. A., & Daigne, F. 2002, MNRAS, 334, 383

Foucart, F., Chandra, M., Gammie, C. F., Quataert, E., & Tchekhovskoy, A. 2017, MNRAS, 470, 2240

Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417

Fragile, P. C., Etheridge, S. M., Anninos, P., Mishra, B., & Kluźniak, W. 2018, ApJ, 857, 1

Fragile, P. C., Gillespie, A., Monahan, T., Rodriguez, M., & Anninos, P. 2012, ApJS, 201, 9

Fragile, P. C., Olejar, A., & Anninos, P. 2014, ApJ, 796, 22

Fromang, S., Papalozou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123

Fujibayashi, S., Kiuchi, K., Nishimura, N., Sekiguchi, Y., & Shibata, M. 2018, ApJ, 860, 64

Gammie, C. F., & Guan, X. 2012, HAM2D: 2D Shearing Box Model, Astrophysics Source Code Library, ascl:1210.022

Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444

Gardiner, T. A., & Stone, J. M. 2005, JCoPh, 205, 509

Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, RvMP, 82, 3121

Giacomazzo, B., & Rezzolla, L. 2007, CQGra, 24, 235

Goddi, C., Falcke, H., Kramer, M., et al. 2017, IJMPD, 26, 1730001

Gold, R., Paschalidis, V., Ruiz, M., et al. 2014, PhRvD, 90, 104030

Gravity Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 618, L10

Guan, X., & Gammie, C. F. 2008, ApJS, 174, 145

Guan, X., Gammie, C. F., Simon, J. B., & Johnson, B. M. 2009, ApJ, 694, 1010

Hada, K., Kino, M., Doi, A., et al. 2013, ApJ, 775, 70

Hamlin, N. D., & Newman, W. I. 2013, PhRvE, 87, 043101

Harten, A., Lax, P. D., & Leer, B. v. 1983, SIAMR, 25, 35

Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84

Hawley, J. F., & Krolik, J. H. 2006, ApJ, 641, 103

Hawley, J. F., Richers, S. A., Guan, X., & Krolik, J. H. 2013, ApJ, 772, 102

Hawley, J. F., Smarr, L. L., & Wilson, J. R. 1984, ApJ, 277, 296

Heil, L. M., Vaughan, S., & Uttley, P. 2012, MNRAS, 422, 2620

Hirose, S., Krolik, J. H., De Villiers, J.-P., & Hawley, J. F. 2004, ApJ, 606, 1083

Hirotani, K., & Okamoto, I. 1998, ApJ, 497, 563

Hogg, J. D., & Reynolds, C. S. 2016, ApJ, 826, 40

Hoshino, M. 2015, PhRvL, 114, 061101

Hubeny, I., & Hubeny, V. 1998, ApJ, 505, 558

Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042

Keppens, R., Meliani, Z., van Marle, A., et al. 2012, JCoPh, 231, 718

Koide, S., Shibata, K., & Kudoh, T. 1999, ApJ, 522, 727

Komissarov, S. S. 2007, MNRAS, 382, 995

Kozłowski, M., Jaroszynski, M., & Abramowicz, M. A. 1978, A&A, 63, 209

Kritsuk, A. G., Nordlund, Å., Collins, D., et al. 2011, ApJ, 737, 13

Krolik, J. H., Hawley, J. F., & Hirose, S. 2005, ApJ, 622, 1008

Kunz, M. W., Stone, J. M., & Quataert, E. 2016, PhRvL, 117, 235101

Kuo, C. Y., Asada, K., Rao, R., et al. 2014, ApJL, 783, L33

Lesur, G., & Longaretti, P.-Y. 2007, MNRAS, 378, 1471

Leung, P. K., Gammie, C. F., & Noble, S. C. 2011, ApJ, 737, 21

Levermore, C. D. 1984, JQSRT, 31, 149

Levinson, A., & Rieger, F. 2011, ApJ, 730, 123

Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018a, MNRAS, 474, L81

Liska, M. T. P., Tchekhovskoy, A., & Quataert, E. 2018b, arXiv:1809.04608

Löffler, F., Faber, J., Bentivegna, E., et al. 2012, CQGra, 29, 115001

Löhner, R. 1987, CMAME, 61, 323

Londrillo, P., & Del Zanna, L. 2000, ApJ, 530, 508

Londrillo, P., & Del Zanna, L. 2004, JCoPh, 195, 17

Longaretti, P.-Y., & Lesur, G. 2010, A&A, 516, A51

Lu, R.-S., Krichbaum, T. P., Roy, A. L., et al. 2018, ApJ, 859, 60

Mahadevan, R., & Quataert, E. 1997, ApJ, 490, 605

Markoff, S., Falcke, H., Yuan, F., & Biermann, P. L. 2001, A&A, 379, L13

Marrone, D. P., Moran, J. M., Zhao, J., & Rao, R. 2007, ApJL, 654, L57

Martí, J. M., & Müller, E. 2015, LRCA, 1, 3

McKinney, J. C. 2006, MNRAS, 367, 1797

McKinney, J. C., & Blandford, R. D. 2009, MNRAS, 394, L126

McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977

McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083

McKinney, J. C., Tchekhovskoy, A., Sadowski, A., & Narayan, R. 2014, MNRAS, 441, 3177

Meliani, Z., Grandclément, P., Casse, F., et al. 2016, CQGra, 33, 155010

Mignone, A. 2014, JCoPh, 270, 784

Mizuno, Y., Nishikawa, K.-I., Koide, S., Hardee, P., & Fishman, G. J. 2006, arXiv:astro-ph/0609004

Mościbrodzka, M., & Falcke, H. 2013, A&A, 559, L3

Mościbrodzka, M., Falcke, H., & Noble, S. 2016a, A&A, 596, A13

Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016b, A&A, 586, A38

Mościbrodzka, M., Gammie, C. F., Dolence, J. C., & Shiokawa, H. 2011, ApJ, 735, 9

Mościbrodzka, M., Gammie, C. F., Dolence, J. C., Shiokawa, H., & Leung, P. K. 2009, ApJ, 706, 497

Mossoux, E., & Grosso, N. 2017, A&A, 604, A85

Nakamura, M., Asada, K., Hada, K., et al. 2018, ApJ, 868, 146

Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69

Narayan, R., Mahadevan, R., Grindlay, J. E., Popham, R. G., & Gammie, C. 1998, ApJ, 492, 554

Narayan, R., Sadowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241

Narayan, R., & Yi, I. 1995, ApJ, 452, 710

Narayan, R., Yi, I., & Mahadevan, R. 1995, Natur, 374, 623

Noble, S. C., Gammie, C. F., McKinney, J. C., & Del Zanna, L. 2006, ApJ, 641, 626

Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411

Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959

Noble, S. C., Mundim, B. C., Nakano, H., et al. 2012, ApJ, 755, 51

Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlusus), ed. C. DeWitt & B. DeWitt (New York: Gordon and Breach), 343

Nowak, M. A., Neilsen, J., Markoff, S. B., et al. 2012, ApJ, 759, 95

Olivares, H., Porth, O., Davelaar, J., et al. 2019, A&A, in press (arXiv:1906.10795)

Olivares, H., Porth, O., & Mizuno, Y. 2018, arXiv:1802.00860
Pareschi, L., & Russo, G. 2005, JSCom, 25, 129

Parfrey, K., Philippov, A., & Cerutti, B. 2019, PhRvL, 122, 035101

Paschalidis, V., Etienne, Z. B., & Shapiro, S. L. 2012, PhRvD, 86, 064032

Paschalidis, V., Liu, Y. T., Etienne, Z. B., & Shapiro, S. L. 2011, PhRvD, 84, 104032

Paschalidis, V., Ruiz, M., & Shapiro, S. L. 2015, ApJL, 806, L14

Penna, R. F., McKinney, J. C., Narayan, R., et al. 2010, MNRAS, 408, 752

Penna, R. F., Sądowski, A., Kulkarni, A. K., & Narayan, R. 2013, MNRAS, 428, 2255

Phinney, E. S. 1995, BAAS, 27, 1450

Porth, O., Olivares, H., Mizuno, Y., et al. 2017, ComAC, 4, 1

Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4

Pu, H.-Y., Nakamura, M., Hirotani, K., et al. 2015, ApJ, 801, 56

Qian, Q., Fendt, C., Noble, S., & Bugli, M. 2017, ApJ, 834, 29

Radice, D., & Rezzolla, L. 2013, in ASP Conf. Ser. 474, Numerical Modeling of Space Plasma Flows (ASTRONUM2012), ed. N. V. Pogorelov, E. Audit, & G. P. Zank (San Francisco, CA: ASP), 25

Radice, D., Rezzolla, L., & Galeazzi, F. 2014, MNRAS L., 437, L46

Ressler, S. M., Tchekhovskoy, A., Quataert, E., & Gammie, C. F. 2017, MNRAS, 467, 3604

Reynolds, C. S., & Miller, M. C. 2009, ApJ, 692, 869

Rezzolla, L., & Zanotti, O. 2013, Relativistic Hydrodynamics (Oxford: Oxford Univ. Press)

Ripperda, B., Bacchini, F., Porth, O., et al. 2019a, ApJS, submitted (arXiv:1907.07197)

Ripperda, B., Porth, O., Sironi, L., & Keppens, R. 2019b, MNRAS, 485, 299

Rusanov, V. V. 1961, ZVMMF, 1, 267

Ryan, B. R., Dolence, J. C., & Gammie, C. F. 2015, ApJ, 807, 31

Ryan, B. R., Gammie, C. F., Fromang, S., & Kestener, P. 2017, ApJ, 840, 6

Ryan, B. R., Ressler, S. M., Dolence, J. C., Gammie, C., & Quataert, E. 2018, ApJ, 864, 126

Sądowski, A., & Narayan, R. 2015, MNRAS, 454, 2372

Sądowski, A., Narayan, R., McKinney, J. C., & Tchekhovskoy, A. 2014, MNRAS, 439, 503

Sądowski, A., Narayan, R., Tchekhovskoy, A., et al. 2015, MNRAS, 447, 49

Sądowski, A., Narayan, R., Tchekhovskoy, A., & Zhu, Y. 2013, MNRAS, 429, 3533

Sądowski, A., Wielgus, M., Narayan, R., et al. 2017, MNRAS, 466, 705

Sano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M. 2004, ApJ, 605, 321

Shiokawa, H., Dolence, J. C., Gammie, C. F., & Noble, S. C. 2012, ApJ, 744, 187

Simon, J. B., Beckwith, K., & Armitage, P. J. 2012, MNRAS, 422, 2685

Simon, J. B., & Hawley, J. F. 2009, ApJ, 707, 833

Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 2012, ApJ, 749, 189

Spiteri, R. J., & Ruuth, S. J. 2002, SJNA, 40, 469

Stepney, S., & Guilbert, P. W. 1983, MNRAS, 204, 1269

Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137

Takahashi, H. R., Ohsuga, K., Kawashima, T., & Sekiguchi, Y. 2016, ApJ, 826, 23

Takahashi, M., Nitta, S., Tatematsu, Y., & Tomimatsu, A. 1990, ApJ, 363, 206

Takahashi, R. 2008, MNRAS, 383, 1155

Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2008, MNRAS, 388, 551

Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, J. Phys. Conf. Ser, 372, 012040

Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79

Tóth, G. 2000, JCoPh, 161, 605

Uttley, P., & McHardy, I. M. 2001, MNRAS, 323, L26

van Leer, B. 1977, JCoPh, 23, 276

van Leer, B. J. 2006, CCoPh, 1, 192

White, C. J., Stone, J. M., & Gammie, C. F. 2016, ApJS, 225, 22

White, C. J., Stone, J. M., & Quataert, E. 2019, ApJ, 874, 168

Witzel, G., Martinez, G., Hora, J., et al. 2018, ApJ, 863, 15

Xia, C., Teunissen, J., El Mellah, I., Chané, E., & Keppens, R. 2018, ApJS, 234, 30

Yuan, F., Markoff, S., & Falcke, H. 2002, A&A, 383, 854

Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529

Zanotti, O., & Dumbser, M. 2015, CoPhC, 188, 110

Zanotti, O., Rezzolla, L., & Font, J. A. 2003, MNRAS, 341, 832

Zilhão, M., & Noble, S. C. 2014, CQGra, 31, 065013