The Event Horizon General Relativistic Magnetohydrodynamic Code Comparison Project

TL;DR

This paper compares nine general relativistic magnetohydrodynamics (GRMHD) codes on a common accretion flow problem, demonstrating their maturity, consistency, and improved agreement with higher resolution.

Contribution

It provides a comprehensive comparison of nine GRMHD codes, establishing their reliability and consistency for modeling astrophysical phenomena near compact objects.

Findings

Code agreement improves with resolution

Community of GRMHD codes is mature and capable

Codes produce consistent results on test problems

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 (EHT) 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 between 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.

Figures

Tables

Table 1. Table 1: Algorithmic details of the code comparison runs for the low-resolution realizations. 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 ​ π ] 0 2 𝜋 [0,2\pi] in all cases.

Code	Timestepper	Riemann s.	Rec.	Mag. field	$(\Delta r_p,\Delta {\theta }_p,\Delta {\varphi }_p)$	Domain: r	Domain: $\theta$
Athena++	2nd Order	HLL	PPM	CT (Gardiner & Stone, 2005)	$(0.506,0.393,0.788)$	$[1.152,50]$	$[0,\pi ]$
BHAC	2nd Order	LLF	PPM	UCT (Del Zanna et al., 2007)	$(0.402,0.295,0.788)$	$[1.185,3333]$	$[0,\pi ]$
BHAC Cart.	2nd Order	LLF	PPM	UCT (Del Zanna et al., 2007)	Cartesian AMR $(0.176,0.163,0.163)$	See Sect. 4.1.2.	–
Cosmos++	SSPRK, 3rd Order	HLL	PPM	(Fragile et al., 2012)	$(0.508,0.375,0.788)$	$[1.25,354]$	$[0.070,3.072]$
ECHO	3rd Order	HLL	PPM	UCT (Del Zanna et al., 2007)	$(0.460,0.382,0.752)$	$[1.048,2500]$	$[0.060,3.082]$
H-AMR	2nd Order	HLL	PPM	UCT (Gardiner & Stone, 2005)	$(0.523,0.379,0.785)$	$[1.169,500]$	$[0,\pi ]$
HARM-Noble	2nd Order	LLF	PPM	PPM+Flux-CT (Tóth, 2000), (Noble et al., 2009)	$(0.578,0.182,0.784)$	$[1.090,80]$	$[0.053,3.088]$
iharm3D	2nd Order	LLF	PLM	Flux-CT (Tóth, 2000)	$(0.519,0.118,0.788)$	$[1.073,50]$	$[0,\pi ]$
IllinoisGRMHD	RK4	HLL	PPM	Vector potential-based PPM+Flux-CT (Tóth, 2000), (Etienne et al., 2010), (Etienne et al., 2012b)	Cartesian AMR $(0.246,0.228,0.229)$	See Sect. 4.1.8.	–
KORAL	2nd order	LLF	PPM	Flux-CT (Tóth, 2000)	$(0.478,0.266,0.785)$	$[1.15,50]$	$[0,\pi ]$

Table 2. Table 2: Quantifications – Time series data. Quantities in angular brackets ⟨ ⋅ ⟩ delimited-⟨⟩ ⋅ \langle\cdot\rangle denote time averages in the interval t ∈ [ 5000 , 10000 ] ​ M 𝑡 5000 10000 M t\in[5000,10000]\rm M with error given by the standard deviation. For further convergence testing, another two BHAC runs ( 256 ∗ superscript 256 256^{*} using 256 × 256 × 128 256 256 128 256\times 256\times 128 cells and 384 × 384 × 384 384 384 384 384\times 384\times 384 cells) and a H-AMR run with 1608 × 1056 × 1024 1608 1056 1024 1608\times 1056\times 1024 cells are listed here, where the resolutions correspond to effective N r × N θ × N ϕ subscript 𝑁 𝑟 subscript 𝑁 𝜃 subscript 𝑁 italic-ϕ N_{r}\times N_{\theta}\times N_{\phi} in the disk region.

$N_{\theta }$	Code	${\dot{M}}_{\mathrm{Peak}}$	$⟨\dot{M}⟩$	${\frac{{\Phi }_{\mathrm{BH}}}{\sqrt{\dot{M}}}|}_{\mathrm{Peak}}$	$⟨{\Phi }_{\mathrm{BH}}/\sqrt{\dot{M}}⟩$	${\dot{L}/\dot{M}|}_{\mathrm{Peak}}$	$⟨\dot{L}/\dot{M}⟩$	${\frac{\dot{E}-\dot{M}}{\dot{M}}|}_{\mathrm{Peak}}$	$⟨\frac{\dot{E}-\dot{M}}{\dot{M}}⟩$
96	Athena++	$0.461$	$0.041\pm 0.019$	$5.951$	$4.164\pm 0.869$	$1.95$	$1.276\pm 0.326$	$0.332$	$0.18\pm 0.044$
	BHAC	$0.507$	$0.074\pm 0.021$	$3.099$	$1.426\pm 0.16$	$1.964$	$1.865\pm 0.033$	$0.33$	$0.112\pm 0.017$
	Cosmos++	$0.602$	$0.142\pm 0.036$	$0.481$	$0.239\pm 0.065$	$2.268$	$2.207\pm 0.025$	$0.232$	$0.025\pm 0.009$
	ECHO	$0.331$	$0.04\pm 0.018$	$2.949$	$2.021\pm 0.377$	$2.005$	$1.211\pm 0.289$	$0.33$	$0.124\pm 0.022$
	H-AMR	$0.536$	$0.091\pm 0.048$	$2.908$	$1.249\pm 0.131$	$2.066$	$1.985\pm 0.028$	$0.336$	$0.046\pm 0.018$
	HARM-Noble	$0.797$	$0.125\pm 0.05$	$2.236$	$0.786\pm 0.119$	$2.206$	$1.964\pm 0.051$	$0.327$	$0.067\pm 0.017$
	iharm3D	$0.645$	$0.136\pm 0.06$	$3.705$	$1.239\pm 0.346$	$2.083$	$1.958\pm 0.087$	$0.331$	$0.067\pm 0.011$
	KORAL	$0.738$	$0.157\pm 0.051$	$3.13$	$0.458\pm 0.072$	$2.1$	$2.027\pm 0.026$	$0.328$	$0.05\pm 0.007$
	max/min	2.406	3.903	12.374	17.408	1.163	1.823	1.448	7.096
128	Athena++	$0.847$	$0.023\pm 0.012$	$5.385$	$2.995\pm 0.694$	$2.082$	$1.743\pm 0.208$	$0.332$	$0.108\pm 0.029$
	BHAC	$0.66$	$0.078\pm 0.046$	$2.492$	$1.74\pm 0.407$	$1.974$	$1.829\pm 0.074$	$0.331$	$0.123\pm 0.015$
	BHAC Cart.	$0.57$	$0.163\pm 0.091$	$2.137$	$1.123\pm 0.277$	$1.947$	$1.864\pm 0.061$	$0.331$	$0.097\pm 0.009$
	Cosmos++	$0.642$	$0.177\pm 0.045$	$0.484$	$0.245\pm 0.032$	$2.133$	$2.073\pm 0.019$	$0.232$	$0.078\pm 0.01$
	ECHO	$0.507$	$0.059\pm 0.021$	$1.575$	$1.056\pm 0.182$	$2.028$	$1.876\pm 0.083$	$0.331$	$0.073\pm 0.012$
	H-AMR	$0.588$	$0.117\pm 0.039$	$2.378$	$1.729\pm 0.204$	$2.04$	$1.961\pm 0.031$	$0.336$	$0.061\pm 0.008$
	HARM-Noble	$1.28$	$0.128\pm 0.051$	$3.23$	$0.61\pm 0.179$	$2.259$	$2.093\pm 0.054$	$0.328$	$0.052\pm 0.017$
	iharm3D	$0.841$	$0.132\pm 0.035$	$3.464$	$1.019\pm 0.135$	$2.089$	$1.964\pm 0.037$	$0.325$	$0.076\pm 0.01$
	IllinoisGRMHD	$0.6$	$0.177\pm 0.041$	$1.477$	$0.779\pm 0.08$	$2.141$	$2.053\pm 0.044$	$0.361$	$0.026\pm 0.009$
	KORAL	$1.221$	$0.199\pm 0.062$	$2.351$	$0.517\pm 0.097$	$2.084$	$2.017\pm 0.025$	$0.332$	$0.056\pm 0.008$
	max/min	2.525	8.604	11.126	12.24	1.16	1.201	1.552	4.631
192	Athena++	$0.875$	$0.134\pm 0.071$	$2.739$	$1.405\pm 0.39$	$2.109$	$2.025\pm 0.047$	$0.332$	$0.065\pm 0.008$
	BHAC	$0.825$	$0.238\pm 0.042$	$2.754$	$0.773\pm 0.107$	$2.121$	$2.072\pm 0.021$	$0.331$	$0.059\pm 0.009$
	ECHO	$0.843$	$0.167\pm 0.047$	$1.698$	$0.506\pm 0.057$	$2.139$	$2.047\pm 0.033$	$0.332$	$0.05\pm 0.008$
	H-AMR	$0.671$	$0.237\pm 0.104$	$2.798$	$0.691\pm 0.049$	$2.086$	$2.019\pm 0.026$	$0.265$	$0.048\pm 0.009$
	HARM-Noble	$0.994$	$0.201\pm 0.116$	$3.593$	$0.972\pm 0.216$	$2.224$	$2.097\pm 0.044$	$0.33$	$0.057\pm 0.013$
	iharm3D	$1.106$	$0.18\pm 0.091$	$3.161$	$1.292\pm 0.19$	$2.066$	$1.871\pm 0.046$	$0.32$	$0.057\pm 0.01$
	KORAL	$1.01$	$0.183\pm 0.081$	$2.256$	$1.217\pm 0.338$	$2.082$	$1.987\pm 0.041$	$0.331$	$0.046\pm 0.007$
	max/min	1.649	1.772	2.116	2.777	1.077	1.121	1.25	1.405
${256}^{\ast }$	BHAC	$0.697$	$0.203\pm 0.065$	$2.552$	$1.021\pm 0.091$	$2.097$	$2.034\pm 0.019$	$0.332$	$0.064\pm 0.007$
384	BHAC	$0.91$	$0.221\pm 0.067$	$3.424$	$0.946\pm 0.128$	$2.12$	$2.043\pm 0.021$	$0.331$	$0.049\pm 0.006$
1056	H-AMR	$1.143$	$0.152\pm 0.056$	$2.324$	$1.341\pm 0.068$	$1.975$	$1.901\pm 0.028$	$0.336$	$0.05\pm 0.006$

Table 3. Table 3: Quantifications – Time series data, linear reconstruction. Quantities in angular brackets ⟨ ⋅ ⟩ delimited-⟨⟩ ⋅ \langle\cdot\rangle denote time averages in the interval t ∈ [ 5000 , 10000 ] ​ M 𝑡 5000 10000 M t\in[5000,10000]\rm M with error given by the standard deviation.

$N_{\theta }$	Code	${\dot{M}}_{\mathrm{Peak}}$	$⟨\dot{M}⟩$	${\frac{{\Phi }_{\mathrm{BH}}}{\sqrt{\dot{M}}}|}_{\mathrm{Peak}}$	$⟨{\Phi }_{\mathrm{BH}}/\sqrt{\dot{M}}⟩$	${\dot{L}/\dot{M}|}_{\mathrm{Peak}}$	$⟨\dot{L}/\dot{M}⟩$	${\frac{\dot{E}-\dot{M}}{\dot{M}}|}_{\mathrm{Peak}}$	$⟨\frac{\dot{E}-\dot{M}}{\dot{M}}⟩$
96	Athena++	$0.209$	$0.049\pm 0.026$	$4.762$	$3.308\pm 0.732$	$1.929$	$1.512\pm 0.225$	$0.332$	$0.145\pm 0.024$
	H-AMR	$0.201$	$0.044\pm 0.015$	$1.341$	$0.914\pm 0.15$	$1.931$	$1.473\pm 0.187$	$0.446$	$0.241\pm 0.049$
	KORAL	$0.863$	$0.09\pm 0.077$	$5.247$	$3.226\pm 1.022$	$2.146$	$1.734\pm 0.286$	$0.328$	$0.033\pm 0.015$
	max/min	4.296	2.062	3.913	3.62	1.112	1.177	1.358	7.258
128	Athena++	$0.266$	$0.035\pm 0.011$	$6.553$	$3.845\pm 0.664$	$1.918$	$1.387\pm 0.272$	$0.371$	$0.167\pm 0.041$
	H-AMR	$0.244$	$0.089\pm 0.023$	$1.211$	$0.707\pm 0.105$	$1.907$	$1.738\pm 0.099$	$0.336$	$0.098\pm 0.012$
	KORAL	$0.65$	$0.179\pm 0.066$	$4.671$	$0.881\pm 0.143$	$2.245$	$2.118\pm 0.046$	$0.332$	$0.041\pm 0.013$
	max/min	2.665	5.165	5.41	5.439	1.177	1.527	1.117	4.049
192	Athena++	$0.42$	$0.069\pm 0.038$	$3.108$	$1.73\pm 0.311$	$2.064$	$1.947\pm 0.047$	$0.332$	$0.094\pm 0.007$
	H-AMR	$0.612$	$0.04\pm 0.018$	$1.349$	$0.916\pm 0.192$	$2.022$	$1.487\pm 0.22$	$0.336$	$0.194\pm 0.045$
	KORAL	$0.789$	$0.202\pm 0.036$	$3.357$	$0.776\pm 0.059$	$2.084$	$2.025\pm 0.02$	$0.337$	$0.065\pm 0.008$
	max/min	1.88	5.027	2.489	2.23	1.031	1.362	1.016	2.967

Table 4. Table 4: Quantifications – Time series data, run to run variation. Quantities in angular brackets ⟨ ⋅ ⟩ delimited-⟨⟩ ⋅ \langle\cdot\rangle denote time averages in the interval t ∈ [ 5000 , 10000 ] ​ M 𝑡 5000 10000 M t\in[5000,10000]\rm M with error given by the standard deviation.

$N_{\theta }$	Code	${\dot{M}}_{\mathrm{Peak}}$	$⟨\dot{M}⟩$	${\frac{{\Phi }_{\mathrm{BH}}}{\sqrt{\dot{M}}}|}_{\mathrm{Peak}}$	$⟨{\Phi }_{\mathrm{BH}}/\sqrt{\dot{M}}⟩$	${\dot{L}/\dot{M}|}_{\mathrm{Peak}}$	$⟨\dot{L}/\dot{M}⟩$	${\frac{\dot{E}-\dot{M}}{\dot{M}}|}_{\mathrm{Peak}}$	$⟨\frac{\dot{E}-\dot{M}}{\dot{M}}⟩$
96	KORAL Odyssey	$0.821$	$0.408\pm 0.129$	$3.029$	$0.791\pm 0.362$	$2.126$	$2.031\pm 0.038$	$0.33$	$0.061\pm 0.008$
	KORAL Stampede2	$0.742$	$0.418\pm 0.112$	$2.818$	$1.068\pm 0.16$	$2.105$	$2.014\pm 0.043$	$0.33$	$0.061\pm 0.01$
	max/min	1.106	1.024	1.075	1.349	1.01	1.008	1.0	1.009
128	KORAL Odyssey	$0.859$	$0.318\pm 0.102$	$2.899$	$1.279\pm 0.306$	$2.073$	$1.978\pm 0.049$	$0.331$	$0.066\pm 0.007$
	KORAL Stampede2	$1.217$	$0.361\pm 0.099$	$2.988$	$1.171\pm 0.17$	$2.051$	$1.963\pm 0.032$	$0.331$	$0.062\pm 0.008$
	max/min	1.417	1.134	1.031	1.093	1.011	1.008	1.0	1.053
192	KORAL Odyssey	$1.067$	$0.29\pm 0.132$	$2.459$	$1.254\pm 0.231$	$2.036$	$1.953\pm 0.023$	$0.331$	$0.056\pm 0.007$
	KORAL Stampede2	$0.933$	$0.327\pm 0.116$	$2.58$	$1.013\pm 0.1$	$2.047$	$1.981\pm 0.025$	$0.331$	$0.056\pm 0.008$
	max/min	1.144	1.128	1.049	1.237	1.006	1.014	1.0	1.014

Equations

$${\partial }_t(\sqrt{-g}\rho u^t)=-{\partial }_i(\sqrt{-g}\rho u^i)$$

$${\partial }_t\left(\sqrt{-g}{T^t}_{\nu }\right)=-{\partial }_i\left(\sqrt{-g}{T^i}_{\nu }\right)+\sqrt{-g}{T^{\kappa }}_{\lambda }{{\Gamma }^{\lambda }}_{\nu \kappa },$$

$$T_{\mathrm{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 }$$

$$b^t=B^i{u}^{\mu }{g}_{i\mu },$$

$$b^i=(B^i+b^t{u}^i)/u^t;$$

$${\partial }_t(\sqrt{-g}{B}^i)=-{\partial }_j(\sqrt{-g}(b^j{u}^i-b^i{u}^j)).$$

$$\frac{1}{\sqrt{-g}}{\partial }_i(\sqrt{-g}{B}^i)=0,$$

$$A_{\varphi }\propto \max (\rho /{\rho }_{\max }-0.2,0),$$

$$p^{\ast }=p(1+X_p)$$

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

$$r(x^{1^{\prime }})=\exp \left(x^{1^{\prime }}\right),\theta (x^{2^{\prime }})={\theta }_c+h\sin \left(2\pi x^{2^{\prime }}\right)$$

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

$$\dot{M}$$

$${\Phi }_{\mathrm{BH}}$$

$$\dot{L}$$

$$\dot{E}$$

$$⟨q⟩(r,t):=\frac{{\int }_0^{2\pi }{\int }_{\theta \min }^{\theta \max }q(r,\theta ,\varphi ,t)\sqrt{-g}d\varphi d\theta }{{\int }_0^{2\pi }{\int }_{\theta \min }^{\theta \max }\sqrt{-g}d\varphi d\theta },$$

$$\mathcal{L}(t):={\int }_0^{2\pi }{\int }_{\theta \min }^{\theta \max }{\int }_{r_{\mathrm{h}}}^{r_{\max }}j(B,p,\rho )\sqrt{-g}d\varphi d\theta dr$$

$$⟨q⟩(r,\theta ):=\frac{{\int }_{t_{\mathrm{beg}}}^{t_{\mathrm{end}}}{\int }_0^{2\pi }q(r,\theta ,\varphi ,t)d\varphi dt}{2\pi (t_{\mathrm{end}}-t_{\mathrm{beg}})}$$

$$\Psi (r,\theta )={\left(\frac{r}{r_{\mathrm{h}}}\right)}^{\kappa }(1-\cos (\theta ))$$

$$H/R(r):=\frac{{\int }_{t_{\mathrm{beg}}}^{t_{\mathrm{end}}}{\int }_0^{2\pi }{\int }_0^{\pi }|\pi /2-{\theta }_{\mathrm{KS}}|\rho (r,\theta ,\varphi ,t)\sqrt{-g}dtd\varphi d\theta }{{\int }_{t_{\mathrm{beg}}}^{t_{\mathrm{end}}}{\int }_0^{2\pi }{\int }_0^{\pi }\rho (r,\theta ,\varphi ,t)\sqrt{-g}dtd\varphi d\theta },$$

$$⟨w^{r\varphi }⟩(r,t):=⟨-b^r\sqrt{{\gamma }_{rr}}{b}^{\varphi }\sqrt{{\gamma }_{\varphi \varphi }}⟩$$

$$\alpha (r,t):=\frac{⟨w^{r\varphi }⟩(r,t)}{⟨P_{\mathrm{tot}}⟩(r,t)}$$

$$⟨r⟩(t):=\frac{{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_{r_{\mathrm{h}}}^{r_{\max }}{r}_{\mathrm{KS}}\rho (r,\theta ,\varphi ,t)\sqrt{-g}d\varphi d\theta dr}{{\int }_0^{2\pi }{\int }_0^{\pi }{\int }_{r_{\mathrm{h}}}^{r_{\max }}\rho (r,\theta ,\varphi ,t)\sqrt{-g}d\varphi d\theta dr}$$

$${\lambda }^{(\theta )}:=\frac{2\pi }{\sqrt{(\rho h+b^2})\Omega }{b}^{\mu }{e}_{\mu }^{(\theta )}$$

$$\Delta x^{(\theta )}:=\Delta {\theta }^{\mu }{e}_{\mu }^{(\theta )}$$

$$⟨w^{r\varphi }⟩(r,t)$$

$$⟨w^{(\bar{r})(\bar{\varphi })}⟩(r,t)$$

Full text

black hole physics - magnetic fields - magnetohydrodynamics (MHD) - methods: numerical - relativistic processes

††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, (Sa̧dowski et al., 2013, 2014)

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 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 black hole 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 black hole at the center of the galaxy M87 with 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$ black hole 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 semi-analytic works (Narayan & Yi, 1995; Narayan et al., 1998; Yuan et al., 2002) have provided with the general parameter regime of the galactic center by exploiting spectral information. For example, Mahadevan & Quataert (1997) demonstrated that the electrons and ions are only weakly collisionally coupled and unlikely 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}_{\rm Edd}=L_{\rm Edd}/(0.1c^{2})\simeq 2\,M/(10^{8}M_{\odot})\,M_{\odot}\rm yr^{-1}$ where $L_{\rm Edd}=4\pi GMc/\sigma_{\rm T}$ is the Eddington luminosity (with $\sigma_{\rm T}$ being the Thomson electron cross section). New observational capabilities like mm- and IR- interferometry, as provided by the EHT and GRAVITY (Gravity Collaboration et al., 2018) collaborations now allow 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 a collisionless plasma. Second, the exchange of energy between the plasma and the radiation field is neglected.

The primary EHT sources Sgr A* and M87* fall in the class of low-luminosity active galactive nuclei (AGN) and accrete with $\dot{M}/\dot{M}_{\rm Edd}\lesssim 10^{-6}$ (Marrone et al., 2007) and $\dot{M}/\dot{M}_{\rm 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 are so far 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., 2018). 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 non-local 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 Sa̧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; Baiotti et al., 2005; Duez et al., 2005; Anninos 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; Radice et al., 2014; Sa̧dowski et al., 2014; McKinney et al., 2014; Etienne et al., 2015; White et al., 2016; Zanotti & Dumbser, 2015; Meliani et al., 2016; Liska et al., 2018a).

Despite the conceptual simplicity of the MHD equations, the non-linear 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 black hole 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 black hole (e.g. McKinney, 2006). The strong magnetic fields that permeate the black hole (held in place by the equatorial accretion flow) are able to extract energy in the form of Poynting flux from a rotating black hole, 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., 2018) can be used to validate the flow in this region.

Due to their robustness when dealing e.g. with super-sonic 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 black hole 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 black hole 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 here give a brief summary of the ideal GRMHD equations in a coordinate basis $(t,x^{i})$ with four-metric $(g_{\mu\nu})$ and metric determinant $g$. As 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 $\rm M$.

The equations describe particle number conservation;

[TABLE]

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

[TABLE]

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

[TABLE]

where $u$ is the fluid internal energy, $p$ the fluid pressure, and $b^{\mu}$ is 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:

[TABLE]

[TABLE]

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

[TABLE]

The system is subject to the additional no-monopoles constraint

[TABLE]

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); 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 magnetised 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 very low Eddington rate $\dot{M}\sim 10^{-6}\dot{M}_{\rm 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 non-radiative GRMHD simulations neglecting cooling can be used to model the data. For a 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 black hole, 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. Since in ideal MHD, plasma cannot move across magnetic field lines (due to the frozen-in condition), there is no way to re-supply 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 components of i) the disk which contains the bound matter ii) the evacuated funnel extending from the polar caps of the black hole 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$ 111There 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. (2018); Narayan et al. (2012) included also magnetic and radiative contributions. Certainly the geometric and hydrodynamic prescriptions underestimate the amount of outflowing material. , 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) also a disk-corona was 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 separate entity in the RIAF scenario.

Since 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 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 post-processing as demonstrated by Mościbrodzka et al. (2011).

Turning back to the morphology of the RIAF accretion, Figure 1, one can see that between evacuated funnel demarcated by the funnel wall (red) and 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), this flow emerges as plasma from the disk is driven against the centrifugal barrier by magnetic and thermal pressure (which coined the alternative term funnel wall jet for this region). In current GRMHD based radiation models as utilized e.g. in 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., 2018; Davelaar et al., 2019). In the mm-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_{\rm BH}$: normalized by the accretion rate $\phi:=\Phi_{\rm BH}/\sqrt{\dot{M}}$, it was shown that a maximum for the magnetic flux $\phi_{\rm max}\approx 15$ (in our system of units) exists which depends only mildly on black hole 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_{\rm max}$, accretion is brought to a near-stop by the accumulation of magnetic field near the black hole (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 black hole magnetosphere (here $\phi\sim\rm 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 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 general-relativistic magnetohydrodynamics in arbitrary spacetimes/coordinates and exploits adaptive mesh refinement techniques with an oct-tree block-based approach. The algorithm is centred 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 $\mathbf{\nabla\cdot B}$ preserving schemes and their implementation in BHAC can be found in Olivares et al. (2018, 2019). Originally designed to study black hole 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 non-black hole fully numerical metrics. In addition, a non-ideal resistive GRMHD module is under development (e.g. Ripperda et al., 2019a, b). 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 non-conservative), magnetic fields (ideal and non-ideal), 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; Londrillo & Del Zanna, 2004), when it was first proposed a shock-capturing scheme for classical MHD based on high-order finite-differences reconstruction routines, one-wave or two-waves 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 non-ideal 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 which 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 timestepping (LAT) capability. LAT is superior to the ’standard’ hierarchical timestepping approach implemented in other AMR codes since the spatial and temporal refinement levels are decoupled, giving much greater speedups on logarithmic 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 (Gammie 2006, private communication) of the 2D Newtonian MHD code called HAM (Guan & Gammie, 2008; Gammie & Guan, 2012). Because of its shared history, HARM-Noble is very similar to the iharm3D code. 222We 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. 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 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, 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 Flux-CT constrained transport algorithm of Tóth (2000). Time integration uses a second-order predictor-corrector step. Several spatial reconstruction options are available, although 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 (Sa̧dowski et al., 2013, 2014) is a multidimensional GRMHD code which 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 PPM scheme. Fluxes can be computed using either the LLF or the HLL method. There is an option to include an artifical magnetic dynamo term in the induction equation (Sa̧dowski et al., 2015), which is helpful for running 2D axisymmetric simulations for long durations (not possible without this term since 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 (Sa̧dowski et al., 2013, 2014) and multi-species fluid. Radiation is handled as a separate fluid component via a moment formalism using M1 closure (Levermore, 1984). A radiative viscosity term is included (Sa̧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 (Sa̧dowski & Narayan, 2015; Sa̧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 (Sa̧dowski et al., 2017), and can also evolve an isotropic population of nonthermal electrons (Chael et al., 2017). Code verification is described in Sa̧dowski et al. (2014).

4 Setup description

As initial condition for our 3D GRMHD simulations, we consider a hydrodynamic equilibrium torus 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^{t}$ 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; Zanotti et al., 2003; Gammie 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_{t}$ 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); White et al. (2016). Hence the spacetime is a Kerr black hole (hereafter BH) with dimensionless spin parameter $a=0.9375$. The inner radius of the torus is set to $r_{\rm in}=6\,M$ and the density maximum is located at $r_{\rm 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

[TABLE]

is added to the stationary solution. The field strength is set such that $2p_{\rm max}/(B^{2})_{\rm max}=100$ , where global maxima of pressure $p_{\rm max}$ and magnetic field strength $B^{2}_{\rm max}$ do not necessarily coincide. With this choice for 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:

[TABLE]

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. Since the strategies differ significantly between implementations, only a rough guideline on the floor model was given. The following floor model was suggested: $\rho_{\mathrm{fl}}=10^{-5}r^{-3/2}$ and $p_{\mathrm{fl}}=1/3\times 10^{-7}\ r^{-5/2}$ which corresponds to the one used by McKinney & Gammie (2004) in the powerlaw indices. Thus for all cells which satisfy $\rho\leq\rho_{\mathrm{fl}}$, set $\rho=\rho_{\mathrm{fl}}$, in addition if $p\leq p_{\mathrm{fl}}$, set $p=p_{\mathrm{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 to keep the error localized. The various failsafes 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 two-fold: first, this way the results can inform on the typical setup used with a particular code, thus allowing 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^{3}$, mid-res: $128^{3}$ and high-res: $192^{3}$ where the resolution corresponds to the domain of interest mentioned above.

To make sure the initial data is setup 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_{\rm p},\Delta\theta_{\rm p},\Delta\phi_{\rm p})$ refer to the proper distance between grid cells at the density maximum in the equatorial plane for the low resolution realization (typically $96\times 96\times 96$). Specifically, we define

[TABLE]

and analogue 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 the $\Delta z_{\rm p}$ as representative for the out-of-plane resolution $\Delta\theta_{\rm 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 GR 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 are 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^{3}$ 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^{3}$ 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^{3}$ 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 black hole, 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 de-refined and the torus is fully resolved on the highest grid level. This allows to significantly increase the timestep which is otherwise dominated by the cells close to the singular axis. Hence compared to the brute force uniform grid setup, the timestep in the 3-level run is increased by a factor of 16. We adopt a floor model in the Eulerian frame as suggested, however, since staggered fields need to be interpolated to the cell centers which introduces an additional error, we have to increase the floors to $\rho_{\mathrm{fl}}=10^{-4}r^{-3/2}$ and $p_{\mathrm{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 hence fill the axial ghost-cells with their diagonal counterpart. No further pole-fix was required for 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 8 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 black hole interior, we apply cut-off values for the density and pressure in the region $r<0.5(r_{\rm H-}+r_{\rm H+})$ where $r_{\rm H-}$ and $r_{\rm H+}$ are the location of the inner and outer event horizons. Specifically, we set $\rho_{\rm cut}=10^{-2}$ or $p_{\rm cut}=10^{-4}$. Other 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 magnetised polar region, we employ the entropy strategy discussed in Porth et al. (2017), that is, each time the plasma-$\beta$ drops below the threshold value of $10^{-2}$, 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 $\theta$-coordinate was used with small cut-outs near the poles ($\theta<4^{\circ}$) to keep the timestep reasonable. The outer radial boundary was placed at $(50)^{1.5}M$ to reduce boundary effects. We then increased the number of zones in the radial direction by $N_{r}^{1.5}$ 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 upwind constrained transport 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 $\theta$ and $\phi$, 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 rest-mass density and pressure used by McKinney & Gammie (2004). For the primitive variable recovery we first use three-dimensional scheme described in Bucciantini & Del Zanna (2013), and in 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 $3^{rd}$ order accurate in space, by using PPM (Colella & Woodward, 1984) reconstruction of primitive variables at cell faces, and $2^{nd}$ 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). Since 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 ZAMO frame (McKinney et al., 2012). We enforce floor values on the rest-mass density $\rho_{\rm fl}=$MAX$[B^{2}/20,10^{-5}(r/M)^{-3/2}]$ and internal energy $u_{\rm fl}=$MAX$[B^{2}/750,3.33\times 10^{-8}(r/M)^{-5/2}]$. 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 $3-4$ levels of local adaptive timestepping to speed up the code by an additional factor $3-5$ to $3-5\times 10^{8}$ 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 $\theta$-direction and periodic boundary conditions in the $\phi$-direction. Since cells get squeezed near the pole, the timestep in all spherical grids is reduced by an additional factor proportional to the resolution in the $\phi$-direction. To remedy this issue, we stretch out cells immediately adjacent to the pole in the $\theta$-direction (Tchekhovskoy et al., 2011) and use multiple levels of static mesh de-refinement in the $\phi$-direction to keep the cell’s aspect ratio close to uniform at high latitudes. As an example, the very high resolution (effectively, $1608\times 1056\times 1024$) run in this work uses a $\phi$-resolution of 64 cells for $0^{\circ}<\theta<3.75^{\circ}$, 128 cells for $3.75^{\circ}<\theta<7.5^{\circ}$, 256 cells for $7.5^{\circ}<\theta<15^{\circ}$, 512 cells for $15^{\circ}<\theta<30^{\circ}$ and the full 1024 cells for $30^{\circ}<\theta<90^{\circ}$. This method prevents the squeezing of cells near the pole from reducing the global timestep, while maintaining high accuracy in all 3 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

[TABLE]

with $\theta_{c}=0.017\pi$, 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}\rightarrow-B^{\theta}$, $v^{\theta}\rightarrow-v^{\theta}$).

Instead of using $96$ cells (for instance for the low resolution run) within $r=50\rm M$ as many of the others used, the HARM-Noble runs used this number over the entire radial extent it used, out to $r=80\rm M$. This means that a HARM-Noble run has 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 “1$\mathrm{D}_{W}$” 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 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}$. 333During 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. Since 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 accumulated round-off error during the course of the simulation. In other words, the HARM-Noble simulations are missing the $m=1,\ldots,\left(N-1\right)$ azimuthal modes of dynamics, where $N=\left(3,6,6\right)$ for runs $96^{3},128^{3},192^{3}$, respectively.

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 timestepping. 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 Flux-CT 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$, $h=0.3$, but with an additional modification to MKS to enlarge grid zones near the polar axis at small radii and increase the timestep. For FMKS, then, $X^{2}$ is chosen to smoothly interpolate from KS $\theta$ to an expression that deresolves the poles:

[TABLE]

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-sidelength 27.34375$M$ 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-two derefinements (with approximate resolutions $M/2.19$, $M/1.10$, $1.82M$, $3.65M$, $7.29M$, and $14.58\bar{3}M$ [exact]) are placed outside this high-resolution level so that the outer boundary is a cube with half-sidelength $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 Flux-CT 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 towards 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}\leq u/\rho c^{2}\leq 10^{2}$, $B^{2}/u\leq 10^{5}$, $B^{2}/\rho c^{2}\leq 50$, $\Gamma\leq 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 inter-code comparison.

5.1 Diagnostics

• •

Horizon penetrating fluxes. The relevant quantities: mass, magnetic, angular momentum and (reduced) energy fluxes are defined as follows:

[TABLE]

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

• •

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 similar to Shiokawa et al. (2012); White et al. (2016); Beckwith et al. (2008). Hence for a quantity $q(r,\theta,\phi,t)$, the shell average is defined as

[TABLE]

which is then further averaged over the time interval $t_{\rm KS}\in[5\,000,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_{\rm min}=\pi/3$, $\theta_{\rm max}=2\pi/3$ ensure that only material from the disk is taken into account in the averaging.

• •

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 non-linear combination of flow variables.

[TABLE]

where again $\theta_{\rm min}=\pi/3$ and $\theta_{\rm max}=2\pi/3$ are and $r_{\rm 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 cutoff after a few gravitational radii, resembling the expected $mm$-emission from the galactic center, that is $C=0.2$ in geometrised units. This emission proxy is optimized for the science case of the EHT where near optically thin synchrotron emission is expected.444In 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)). In order to compare the variability, again a cadence of $1M$ or less is chosen in most cases.

• •

$t,\phi$** - averages.** Finally, we compare temporally and azimuthally averaged data for a more global impression of the disk and jet system.

[TABLE]

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

5.2 Time series

Time series data of horizon penetrating fluxes is presented in Figures 2 to 4 for low, medium and high resolutions respectively. Since the mass accretion governs the behavior of these fluxes, the data has 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_{\rm max}$, all simulations are in the SANE regime of radiatively inefficient black hole 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 non-linear 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 slightly lower than all spherical codes.

Since meshes with various amounts of compression in the mid-plane 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_{\rm p}$. Ordered by $\Delta\theta_{\rm p}$, the trends in the data becomes clearer: for example, the mid-plane 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 also KORAL and HARM-Noble 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++ are consistently off by a factor of $\sim 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 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 relevant quantities are presented in Figures 6 to 8. At the late times under consideration, viscous spreading and accretion has 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 powerlaw with slope $r^{-1.75}$, somewhat steeper than the Keplerian case. Interestingly, also the profiles of magnetic field are 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 as $r^{-1}$ of the disk magnetic field 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 the 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 artefact due to lack of angular momentum transport in the under-resolved 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 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 to fall 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 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).

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 10\rm M$ 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 10\rm M$.

Taking advantage of the scale-freedom in the ideal MHD approximation, we have also performed a comparison with density, pressure and magnetic field variables re-scaled 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 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 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, 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. Besides 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 towards the axis for some codes (Cosmos++, KORAL, HARM-Noble, Athena++, H-AMR) while others see a 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, also cartesian grids are far from perfect: in the IllinoisGRMHD run we observe some artefacts in the form of horizontal features in $\langle\beta^{-1}\rangle$ which are associated with resolution jumps in the computational grid.

Significant leakage of magnetic flux leads to a complete loss of the magnetically dominated region in the medium resolution Cosmos++ run. A similar behaviour 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, magnetic field simply leaks through the boundary, leaving less flux on the black hole (cf. 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 $\sigma=1$ (e.g. McKinney & Gammie, 2004; Nakamura et al., 2018) in more detail. Figure 13 illustrates that as resolution is increased, the contour is more faithfully recovered across codes: whereas some low resolution runs do not capture the highly magnetised funnel, at high resolution all runs show a clearly defined jet-disk boundary. Despite the large variances in floor treatment, at $N_{\theta}=192$ cells, the difference is reduced to five degrees and the polar angle of the disk-jet boundary ranges between $10^{\circ}$ and $15^{\circ}$ at $r=50\rm M$. For illustrative purposes, we also overplot flux-functions of the approximate force-free solutions discussed in Tchekhovskoy et al. (2008); Nakamura et al. (2018):

[TABLE]

In the recent 2D simulations of Nakamura et al. (2018) the jet boundary given by $\sigma=1$ was accurately described by the choice $\kappa=0.75$ for a wide range of initial conditions, black hole spins and horizon penetrating flux $\Phi_{\rm BH}/\sqrt{\dot{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^{5}\rm M$ (e.g. Asada & Nakamura, 2012; Hada et al., 2013). In our 3D SANE models, the line $\kappa=0.75$ is recovered at low resolution by ECHO and iharm3D and the more collimated genuinely paraboloidal shape $\kappa=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_{\theta}=1056$ H-AMR run with the paraboloidal shape is particularly good. We suspect that the narrower jet profile compared to Nakamura et al. (2018) is due to the lower $\phi\simeq 2-3$ of our benchmark configuration (at high resolution).

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

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

[TABLE]

where the averaging-time is again taken between $t_{\rm beg}=5000\rm M$ and $t_{\rm end}=10\,000\rm M$. As shown in Figure 15, a scale height of $H/R=0.25-0.3$ between $r_{\rm 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 which 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 (cf. Figure 9) which shows up as kink in Figure 15.

5.5 Variability analysis

The analysis of lightcurves 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 Gravity Collaboration et al. (2018) places the origin of an IR flare within 10 gravitational radii of the black hole. 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; 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 sub-mm 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$ hrs in IR and $\tau\sim 8$ hrs in sub-mm) the continuous variability is well described by a red-noise spectrum with powerlaw slope of $\approx 2$ with no indication for a further break down to $\sim 2-10$ min, 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 sub-mm characteristic timescales are in fact statistically compatible indicating a direct relation of the emitting regions. For low-luminosity AGN, Bower et al. (2015) have found that the sub-mm characteristic timescale is consistent with a linear dependence on the black hole 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^{+61}_{-24}$ days (Bower et al., 2015) can be obtained by scaling of the Sgr A* result.

Due to the potential impact of the space-time on the variability of the emission555Numerical simulations by Dolence et al. (2012) have observed a steepening of the IR and X-ray spectral index to $\sim 3$ accompanied by quasi periodic oscillations (QPOs) at the frequency of the innermost stable circular orbit (ISCO) (corresponding to $f_{0}=0.041/M$ in our case and thus periods of $8$ minutes if scaled to Sgr A*)., it is clear that the characterisation 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 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 mm- variability without actually subjecting each code’s output to a rigorous general-relativistic radiative transfer calculation (GRRT), an approximate optically thin lightcurve was calculated as a volume integral over the domain of interest according to Eq. (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 (Younsi & et al., 2019 in prep.) to compute the 1.3 mm emission corresponding to the EHT observing frequency from the high resolution BHAC data scaled to Sgr A*.

The lightcurves are computed assuming a resolution of $1024\times 1024$ pixels and a field-of-view of 100 M in the horizontal and vertical directions, i.e., $[-50~{}{\rm M},50~{}{\rm M}]$. We here neglect 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) [eqn. (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 identical to that specified in the emission proxy, i.e., eqn. (18). The mass and distance of Sgr A* adopt the standard fiducial values (Genzel et al., 2010), and the 1.3 mm flux is normalised to $3.4$ Jy (Marrone et al., 2007) for an observer inclination of $60^{\circ}$ using the fiducial snapshot at 10 000 M.

A comparison of the lightcurves for $i=0^{\circ}$, $i=90^{\circ}$, 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 lightcurve 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 lightcurve which leads to larger fluctuations in the edge-on lightcurve compared to the face on ($i=0^{\circ}$) curve.

Having validated the emission proxy as a reasonable descriptor for the flux in the mm-band, 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 transform of each Hamming-windowed segment, followed by Fibonnacci 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, also synthetic pure red-noise data with a PSD of $f^{-\beta}$ with $\beta\in\{2,3\}$ was created and is faithfully recovered by the algorithm.666Note that due to the spectral leakage of the finite time lightcurves, using a boxcar window resulted in artificially flattening of the $\beta=3$ spectrum to $\beta\sim 2.5$ It is clear that as each window only contains $10\,\rm hrs$, the observed break at the characteristic frequency $\tau\sim 8\rm hrs$ (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 lightcurves 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 lightcurve obtained by GRRT at $i=0^{\circ}$ (light gray) shown in the lower left panel displays a very similar behavior as the corresponding data obtained via the proxy (BHAC code, blue).

Since 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 powerlaw seen in $\mathcal{L}$ is not recovered in $\dot{M}$. For frequencies $f\lesssim 0.01\rm M^{-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}\rm M$ is consistent with e.g. Hogg & Reynolds (2016) who extracted the accretion rate at the ISCO itself. The latter 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 all Fourier amplitudes with frequency below $1/1000\rm M$. 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]\rm 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, iharm3D on the high end. The (normalized) RMS values shown against mid-plane resolution (right panel) quantify this further and indicate quite a universal behaviour with values in the range $1.2-1.6$ across all codes and resolutions.777Using the non-normalized (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 lightcurves, as for example seen in the low-resolution iharm3D data in Figure 2. With increased resolution, all codes tend to be attracted to the point $k=$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 $\sim 8$ in the sample) is a re-statement of the RMS-flux relationship which is a ubiquitous feature of black hole 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$.

Summarising the comparison of variability, 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:

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

The PSD for $\mathcal{L}$ and $\dot{M}$ steepens at $f\simeq 0.1\rm M^{-1}$. 3. 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. 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 results of BHAC and H-AMR.

5.6.1 Maxwell stress and $\alpha$

We first compute the disks $\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); 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 co-moving with the fluid. Even so, further ambiguities arise when defining the orientation of the co-moving tetrad and hence the $b^{(r)}$ and $b^{(\phi)}$ directions. With the aim of merely comparing codes and 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

[TABLE]

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

[TABLE]

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_{\rm KS}=2.04\rm M$) where $\alpha$ drops to zero. This is consistent with 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-3\rm M$. 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, $\alpha$ in our global simulations increases, as opposed to what is observed in local shearing boxes. For example, Bodo et al. (2014); Ryan et al. (2017) report $\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 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 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

[TABLE]

where the outer radius of integration has again been set to the domain of interest $r_{\rm max}=50\rm M$. As more mass is expelled beyond $50\rm M$ though, a larger integration domain would yield somewhat different results. The corresponding temporal evolution is displayed in Figure 20 for BHAC and H-AMR where linear and piecewise parabolic reconstruction have 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\sim 6000$ to values 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 Fig. 6 - 8: $\langle\rho(r)\rangle$), but the disks magnetic field is affected to a lesser extent (cf. Fig. 6 - 8: $\langle|B|(r)\rangle$), hence the disks $\alpha$ increases further in a non-linear 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.

5.6.3 MRI quality factors

These resolution dependent effects can be understood by an inspection of the “MRI quality factors” (in short $Q$-factors), 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 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 non-linear 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)}\geq 10$, $Q^{(\phi)}\geq 20-25$ of Hawley et al. (2011, 2013), however, lower toroidal resolution requirements have also been suggested, for example $Q^{(z)}\geq 10-15$ for $Q^{(\phi)}\approx 10$ of Sorathia et al. (2012). Furthermore, the latter authors note that the toroidal and poloidal resolutions are coupled and low $Q^{(z)}$ can be compensated by a $Q^{(\phi)}\geq 25$. Hence the product of quality factors $Q^{(z)}Q^{(\phi)}\geq 200-250$ has also been suggested (Narayan et al., 2012; Dhang & Sharma, 2018) as quality metric.

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

[TABLE]

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

[TABLE]

where $\Delta\theta^{\mu}=(0,0,\Delta\theta,0)$ is the distance between two adjacent grid lines and $\Omega=u^{\phi}/u^{t}$ 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 regarding the poloidal MRI wavelengths. The time- and spatial average values in the domain of interest are $Q^{(r,\theta,\phi)}_{96}\simeq 2.9,3,14$, $Q^{(r,\theta,\phi)}_{128}\simeq 3,3.8,17$, $Q^{(r,\theta,\phi)}_{192}\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 bottle-neck is clearly in the poloidal direction with $Q^{(r)}\simeq Q^{(\theta)}$ but the azimuthal factor $Q^{(\phi)}$ is approximately twice larger. Hence one might surmise that the resolution in $\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 also the $r$ direction was 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 azimuthal resolution, indeed $Q^{(r)}$ and $Q^{(\theta)}$ slightly decrease (to $Q^{(r,\theta)}\simeq 5.3,6.8$) and $Q^{(\phi)}\simeq 9.6$ essentially halves. Upon increasing the resolution to $N_{\theta}=384$, we obtain $Q^{(r,\theta,\phi)}_{384}\simeq 21.3,24.5,40.9$ and hence a super-linear increase in $Q-$factor compared to the $N_{\theta}=192$ case. This behavior 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 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 5 levels of static mesh refinement in the $\phi$ direction (see Sec. 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^{(r,\theta,\phi)}_{1056}\geq 120,100,190$. This is more than sufficient for capturing MRI-driven turbulence as per the criteria ($Q^{(z,\phi)}\geq 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 micro-physical phenomena such as MRI in a macro-physical 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 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 increase of the MRI quality factors for the very high resolution run is beyond the factor of $\sim 5.5$ which would be gained by going from $N_{\theta}=192$ to $N_{\theta}=1056$. Instead we find an extra factor of two increase in $Q$-factors (cf. Appendix A), consistent with the simultaneous raise in disk-magnetization by the same amount. Overall trends of parameters over time shown in Fig. 4 are quite similar to the iharm3D and KORAL high-resolution data, with some other $192^{3}$ PPM data within range.

This gives a first impression on the behavior of viscous angular momentum transport and the turbulent dissipation and generation of magnetic energy. As evidenced by Fig. 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 $\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 standardised setup of radiatively inefficient black hole 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 re-written. There are both differences and commonalities between the codes. Listing the most important common elements, we note that all apply overall 2nd order conservative schemes and preserve the divergence of the magnetic field to machine precision. All except iharm3D applied the piecewise-parabolic-method (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 towards 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 $\sim 1.7$ between the lowest and the highest.

6.1 Resolution and convergence

We have performed simulations at resolutions commonly used for black hole accretion and forward modeling of EHT observations. In fact, early GRMHD models employed simulations that use comparable or less cells than our low- and mid-resolution setups of $96^{3}-128^{3}$ cells. For example, the fiducial runs of McKinney & Blandford (2009) use $256\times 128\times 32$ cells (though with fourth-order reconstruction) and Fragile et al. (2007) used a standard effective resolution of $128^{3}$ 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\times 96\times 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\times 64\times 64$ cells (accounting for the restricted $\Delta\phi=\pi$). Although comparing cell numbers should be done with a fair grain of salt, taking into account gridding and initial conditions, these numbers give an indication of the recent state of affairs.

Balancing the desire of exploring 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 mid-plane resolutions of $\Delta\theta_{\rm p}\simeq 0.25-0.0125\rm M$ 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 as 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 ($\sim 192^{3}$ cells), the disk $\alpha$ parameter and magnetization continue their weak resolution dependent rise up to $\sim 1024^{3}$ cells. The latter confirms and extends the findings of Shiokawa et al. (2012) by a resolution factor of $\times 2.75$.

6.2 Systematics

Turning to the systematics between codes, we observe a veritable diversity in averaged density profiles. Although the spread diminishes with increasing resolution, there is still a factor of $\sim 3$ difference in peak density at nominal $192^{3}$ 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 non-linear runnaway 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 to re-scale the value of density together with the other MHD variables. Not only do the absolute dissipative length-scales depend on resolution, also their ratio given through the effective magnetic Prandtl number $Pr_{\rm m}$ varies with 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 non-linear MRI depend critically on $Pr_{\rm m}$ (Lesur & Longaretti, 2007; Simon & Hawley, 2009; Fromang et al., 2007; 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 counter-acted by zero-flux “hard” boundaries at the polar cutout as noticed 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 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. Since 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 2D averaged maps of density, (inverse-) plasma-$\beta$ and the magnetization. These allow to quantify the variance of jet opening angles between codes and we find that the spread around a ’mean simulation’ is reduced to $\sim 2.5^{\circ}$ at $r=50\rm M$ at high resolutions. Interestingly, the jet boundaries up to $r_{\rm KS}=50\rm M$ follow closely the paraboloidal 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 black hole is increased towards the MAD case which is known to exhibit a wider jet base. The maps of plasma-$\beta$ and magnetization illustrate truthfully 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. jet opening angle and profile) is minor.

6.4 Time variability

Using an emission proxy taylored to optically thin synchrotron emission from electrons distributed according to a relativistic Maxwellian, we have computed lightcurves with all codes and resolutions and have compared their power spectra. The same analysis was performed with the accretion rates extracted at the black hole horizon. In broad terms, all PSDs of the lightcurves are compatible with a red-noise spectrum up to $f\simeq 1/10\rm M$ 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/10\rm M$. This is quite consistent across all codes and resolutions and agrees with earlier results of Armitage & Reynolds (2003); 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.1\rm M$ 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 de-trended timeseries of accretion rates and emission proxy, we have computed the root-mean-square (RMS) variability and find 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 towards $\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 with new general insights into rotating turbulence. See Section 7 for access instructions. Furthermore, animations of simulation output can be found in the supporting material at 10.7910/DVN/UCFCLK (catalog 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 magnetically arrested disks (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 magnetised 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 bring 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 non-ideal GRMHD modeling e.g. radiative MHD (e.g. Fragile et al., 2012; Sa̧dowski et al., 2014; McKinney et al., 2014; Ryan et al., 2015; Takahashi et al., 2016) and resistive MHD (e.g. Bugli et al., 2014; Qian et al., 2016; Ripperda et al., 2019a). 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 to perform direct numerical simulations of turbulent black hole accretion covering both the low and high $Pr_{\rm m}$ regimes which are expected to be present in such systems (Balbus & Henri, 2008).

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 10.7910/DVN/UCFCLK (catalog 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 cyverse Data Commons (http://datacommons.cyverse.org/). The data is freely accessible at http://datacommons.cyverse.org/browse/iplant/home/shared/eht/2019/GRMHDCodeComparisonProject. Please consult README.txt and LICENSE.txt for usage instructions.

RN thanks the National Science Foundation (NSF, grants OISE-1743747, AST-1816420) and acknowledges computational support from the NSF via XSEDE resources (grant TG-AST080026N). LDZ acknowledges support from the PRIN-MIUR project Multi-scale Simulations of High-Energy Astrophysical Plasmas (Prot. 2015L5EE2Y) and from the INFN - TEONGRAV initiative. CJW 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 awards no. 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). KC and SM are supported by the Netherlands Organization for Scientific Research (NWO) VICI grant (no. 639.043.513), ML 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). MB 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 number ACI-1548562. Work by P.A. was performed in part under the auspices of the U. S. 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 C52I13000140001); the European Union’s Horizon 2020 research and innovation programme under grant agreement No 730562 RadioNet; ALMA North America Development Fund; the Academia Sinica; 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).

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_{\beta}$ limiter (van Leer, 1977) with $\beta=1.5$ and KORAL $MC_{\beta}$ with $\beta=1.9$ where $MC_{\beta}$ reduces to the most diffusive “minmod-limiter” in case of $\beta=1$ and to “MC” for $\beta=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) note 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^{3}$, 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 two to almost six. For KORAL on the other hand, we obtain a more consistent behavior between the reconstruction techniques.

Apart from a slightly higher mid-plane 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 arithmetic averaging of the electric fields, (ACT) is chosen, H-AMR and Athena++ both implement a staggered upwind constrained transport scheme (UCT) 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_{\theta}\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.

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, since the compiler version and optimization (e.g. order of execution) influences the round-off error, a similar effect can be observed if the same physical scenario is run on two machines with differing compiler and run-time configurations. In order to judge the impact this makes compared to the inter-code 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.

Appendix C Rescaled disk profiles

With the scale-freedom allowed by the test-fluid assumption, density, pressure and the magnetic fields of a given simulation can in principle be re-scaled by a constant factor (respecitvely its square root for the magnetic fields), for example to perform spectral fitting of observations. We here exploit this freedom and re-scale 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=12\rm M$. 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<\rm 12M$ in densities and pressure, but the spread in magnetic fields is worsened. This is merely a consequence of the non-convergence of magnetization, for if the density is made to match, the spread in magnetic field increases with the reference run having roughly a twice larger magnetization than most $192^{3}$ runs.

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 magneto-rotational instability, the physical interpretation is best guided by stresses measured in a frame co-rotating with the disk. Hence to relate to the classical disk model of Novikov & Thorne (1973), Krolik et al. (2005) devised a set of co-moving 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^{\prime}_{\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:

[TABLE]

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 co-moving with the local fluid velocity as in Krolik et al. (2005). Integration is carried out over the equatorial wedge $\theta_{\rm KS}\in[\pi/3,2\pi/3]$ and the full azimuthal range. Note that we have taken the fluid frame integration over the co-moving coordinate increments $dx^{(\bar{\mu})}$ which result from the transformation of $dx^{\mu}:=(0,\Delta r,\Delta\theta,\Delta\phi)$ as appropriate (see discussion in Noble et al. 2010).

The resulting stress profiles and volume integrated time-series (non-dimensionalised by time- and volume- averaged total pressure $\langle P_{\rm 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_{\rm KS}\simeq 2\rm M$ 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.04\rm M$ which roughly coincides with the change of slope. The overall stress differs only by a few percent between diagnostics (D1) and (D2) 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 black hole 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 monotonous.