The Event Horizon General Relativistic Magnetohydrodynamic Code Comparison Project

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

The equations describe particle number conservation:

$\begin{eqnarray}&&{\partial }_{t}(\sqrt{-g}\rho {u}^{t})=-{\partial }_{i}(\sqrt{-g}\rho {u}^{i}),\end{eqnarray} \tag{ 1 }$

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

$\begin{eqnarray}&&{\partial }_{t}\left(\sqrt{-g}\,{{T}^{t}}_{\nu }\right)=-{\partial }_{i}\left(\sqrt{-g}\,{{T}^{i}}_{\nu }\right)+\sqrt{-g}\,{{T}^{\kappa }}_{\lambda }\,{{{\rm{\Gamma }}}^{\lambda }}_{\nu \kappa },\end{eqnarray} \tag{ 2 }$

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

$\begin{eqnarray}&&{T}_{\mathrm{MHD}}^{\mu \nu }=(\rho +u+p+{b}^{2}){u}^{\mu }{u}^{\nu }+\left(p+\displaystyle \frac{1}{2}{b}^{2}\right){g}^{\mu \nu }-{b}^{\mu }{b}^{\nu },\end{eqnarray} \tag{ 3 }$

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

$\begin{eqnarray}&&{b}^{t}={B}^{i}{u}^{\mu }{g}_{i\mu },\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&{b}^{i}=({B}^{i}+{b}^{t}{u}^{i})/{u}^{t};\end{eqnarray} \tag{ 5 }$

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

$\begin{eqnarray}&&{\partial }_{t}(\sqrt{-g}{B}^{i})=-{\partial }_{j}(\sqrt{-g}\,({b}^{j}{u}^{i}-{b}^{i}{u}^{j})).\end{eqnarray} \tag{ 6 }$

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

$\begin{eqnarray}&&\displaystyle \frac{1}{\sqrt{-g}}{\partial }_{i}(\sqrt{-g}{B}^{i})=0,\end{eqnarray} \tag{ 7 }$

which follows from the time component of the homogeneous Maxwell equations. For closure, we adopt the equation of state of an ideal gas. This takes the form $p=(\hat{\gamma }-1)u$, where $\hat{\gamma }$ is the adiabatic index. More in-depth discussions of the ideal GRMHD system of equations can be found in the various publications of the authors, e.g., Gammie et al. (2003), Anninos et al. (2005), Del Zanna et al. (2007), White et al. (2016), and Porth et al. (2017).

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

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.

Zoom In Zoom Out Reset image size

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

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

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

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

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

In RIAF accretion, a special role is played by the horizon-penetrating magnetic flux ${{\rm{\Phi }}}_{\mathrm{BH}}$: normalized by the accretion rate $\phi := {{\rm{\Phi }}}_{\mathrm{BH}}/\sqrt{\dot{M}}$, it was shown that a maximum for the magnetic flux ${\phi }_{\max }\approx 15$ (in our system of units) exists, which depends only mildly on BH spin, but somewhat on the disk scale height (with taller disks being able to hold more magnetic flux; Tchekhovskoy et al. 2012). Once the magnetic flux reaches ${\phi }_{\max }$, accretion is brought to a near-stop by the accumulation of magnetic field near the BH (Tchekhovskoy et al. 2011; McKinney et al. 2012), leading to a fundamentally different dynamic of the accretion flow and maximal energy extraction via the Blandford & Znajek (1977) process. This state is commonly referred to as Magnetically Arrested Disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1976; Narayan et al. 2003) to contrast with the Standard and Normal Evolution (SANE), where accretion is largely unaffected by the BH magnetosphere (here $\phi \sim \mathrm{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.

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

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

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. ${\mathtt{Cosmos}}++$ utilizes unstructured meshes with adaptive (h-) refinement (Anninos et al. 2005), moving-mesh (r-refinement; Anninos et al. 2012), and adaptive order (p-refinement; Anninos et al. 2017) capabilities, enabling it to evolve fluid systems over a wide range of spatial scales with targeted precision. It includes numerous hydrodynamics solvers (conservative and nonconservative), magnetic fields (ideal and nonideal), radiative cooling and transport, geodesic transport, generic tracer fields, and full Navier–Stokes viscosity (Fragile et al. 2018). For this work, we utilize the High Resolution Shock-Capturing scheme with staggered magnetic fields and Constrained Transport as described in Fragile et al. (2012). Code verification is described in Anninos et al. (2005).

The origin of the ${\mathtt{Eulerian}}\,{\mathtt{Conservative}}\,{\mathtt{High}}\,-{\mathtt{Order}}$ (ECHO) code dates back to the year 2000 (Londrillo & Del Zanna 2000, 2004), when it was first proposed as a shock-capturing scheme for classical MHD based on high-order, finite-difference reconstruction routines, one-wave or two-wave Riemann solvers, and a rigorous enforcement of the solenoidal constraint for staggered electromagnetic field components (the Upwind Constraint Transport, UCT). The GRMHD version of ECHO used in the present paper is described in Del Zanna et al. (2007) and preserves the same basic characteristics. Important extensions of the code were later presented for dynamical spacetimes (Bucciantini & Del Zanna 2011) and nonideal Ohm equations (Bucciantini & Del Zanna 2013; Del Zanna et al. 2016; Del Zanna & Bucciantini 2018). Specific recipes for the simulation of accretion tori around Kerr black holes can be found in Bugli et al. (2014, 2018). Further references and applications may be found at www.astro.unifi.it/echo. Code verification is described in Del Zanna et al. (2007).

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

The HARM-Noble code (Gammie et al. 2003; Noble et al. 2006, 2009) is a flux-conservative, high-resolution shock-capturing GRMHD code that originated from the 2D GRMHD code called HARM (Gammie et al. 2003; Noble et al. 2006) and the 3D version (C. F. Gammie 2006, private communication) of the 2D Newtonian MHD code called HAM (Guan & Gammie 2008; Gammie & Guan 2012). Because of its shared history, HARM-Noble is very similar to the iharm3D 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 that a physical set of primitive variables is always recovered. HARM-Noble was also written to be agnostic to coordinate and spacetime choices, making it in a sense generally covariant. This feature was most extensively demonstrated when dynamically warped coordinate systems were implemented (Zilhão & Noble 2014) and time-dependent spacetimes were incorporated (e.g., Noble et al. 2012; Bowen et al. 2018).

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

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

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

Although KORAL is suitable for pure GRMHD simulations such as the ones discussed in this paper, the code was developed with the goal of simulating general relativistic flows with radiation (Sa̧dowski et al. 2013, 2014) and multispecies fluid. Radiation is handled as a separate fluid component via a moment formalism using M1 closure (Levermore 1984). A radiative viscosity term is included (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.1.1. Athena++

4.1.2. BHAC

4.1.3. Cosmos++

4.1.4. ECHO

4.1.5. H-AMR

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

$\begin{eqnarray}&&r({x}^{{1}^{{\prime} }})=\exp \left({x}^{{1}^{{\prime} }}\right),\quad \theta ({x}^{{2}^{{\prime} }})={\theta }_{c}+h\sin \left(2\pi {x}^{{2}^{{\prime} }}\right),\end{eqnarray} \tag{ 11 }$

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 }\to -{B}^{\theta }$, ${v}^{\theta }\to -{v}^{\theta }$).

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

Recovery of the primitive variables from the conserved variables was performed with the "2D" and "1${{\rm{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 equation of motion (EOM) replaces the energy EOM in the primitive variable method whenever it fails to converge; the resultant primitive variables are unphysical, or $u\lt {10}^{-2}\,{B}^{2}$.¹²⁰

4.1.7. iharm3D

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

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

$\begin{eqnarray}&&\theta =N\,y\left(1+\displaystyle \frac{{(y/{x}_{t})}^{\alpha }}{\alpha +1}\right)+\displaystyle \frac{\pi }{2},\end{eqnarray} \tag{ 12 }$

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

iharm3d imposes floors on rest-mass density ρ and internal energy u to enforce $\rho \gt {10}^{-5}{r}^{-2}$ and $u\gt {10}^{-7}{r}^{-2\hat{\gamma }}$, where $\hat{\gamma }$ is the adiabatic index. It also requires that $\rho \gt {B}^{2}/50$, $u\,\gt {B}^{2}/2500$, and subsequently that $\rho \gt 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 $\lt 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¹ and X³ fluxes of the X² component of the magnetic field are antiparallel across the boundary. The X³ boundaries are periodic.

4.1.8. IllinoisGRMHD

4.1.9. KORAL

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

  $\begin{eqnarray}&&\dot{M}:= {\int }_{0}^{2\pi }{\int }_{0}^{\pi }\rho {u}^{r}\sqrt{-g}\,d\theta \,d\phi ,\end{eqnarray} \tag{ 13 }$

  $\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{BH}}:= \displaystyle \frac{1}{2}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{| }^{* }{F}^{{rt}}| \sqrt{-g}\,d\theta \,d\phi ,\end{eqnarray} \tag{ 14 }$

  $\begin{eqnarray}&&\dot{L}:= {\int }_{0}^{2\pi }{\int }_{0}^{\pi }{T}_{\phi }^{r}\sqrt{-g}\,d\theta \,d\phi ,\end{eqnarray} \tag{ 15 }$

  $\begin{eqnarray}&&\dot{E}:= {\int }_{0}^{2\pi }{\int }_{0}^{\pi }(-{T}_{t}^{r})\sqrt{-g}\,d\theta \,d\phi ,\end{eqnarray} \tag{ 16 }$

  where all quantities are evaluated at the outer horizon ${r}_{{\rm{h}}}$. A cadence of $1M$ or less is chosen. In practice, these quantities are nondimensionalized with the accretion rate, yielding, for example, the normalized magnetic flux $\phi ={{\rm{\Phi }}}_{\mathrm{BH}}/\sqrt{\dot{M}}$, also known as the "MAD parameter," which, for spin a = 0.9375 and torus scale height $H/R\approx 0.3$, has the critical value ${\phi }_{\max }\approx 15$ (within the units adopted here; Tchekhovskoy et al. 2012).
• 2.  
  Disk-averaged quantities. We compare averages of the basic flow variables $q\in \{\rho ,p,{u}^{\phi },\sqrt{{b}_{\alpha }{b}^{\alpha }},{\beta }^{-1}\}$. These are defined similarly to Shiokawa et al. (2012), White et al. (2016), and Beckwith et al. (2008). Hence, for a quantity $q(r,\theta ,\phi ,t)$, the shell average is defined as

  $\begin{eqnarray}&&\langle q\rangle (r,t):= \displaystyle \frac{{\int }_{0}^{2\pi }{\int }_{\theta \min }^{\theta \max }q(r,\theta ,\phi ,t)\sqrt{-g}\,d\phi \,d\theta }{{\int }_{0}^{2\pi }{\int }_{\theta \min }^{\theta \max }\sqrt{-g}\,d\phi \,d\theta },\end{eqnarray} \tag{ 17 }$

  which is then further averaged over the time interval ${t}_{\mathrm{KS}}\in [5000,{\rm{10,000}}]M$ to yield $\langle q\rangle (r)$ (note that we omit the weighting with the density as done by Shiokawa et al. 2012; White et al. 2016). The limits ${\theta }_{\min }=\pi /3$, ${\theta }_{\max }=2\pi /3$ ensure that only material from the disk is taken into account in the averaging.
• 3.  
  Emission proxy. To get a feeling for the code-to-code variations in synthetic, optically thin light curves, we also integrate the pseudo-emissivity for thermal synchrotron radiation following an appropriate nonlinear combination of flow variables:

  $\begin{eqnarray}&&{ \mathcal L }(t):= {\int }_{0}^{2\pi }{\int }_{\theta \min }^{\theta \max }{\int }_{{r}_{{\rm{h}}}}^{{r}_{\max }}j(B,p,\rho )\sqrt{-g}\,d\phi \,d\theta \,{dr},\end{eqnarray} \tag{ 18 }$

  where again ${\theta }_{\min }=\pi /3$, ${\theta }_{\max }=2\pi /3$, and ${r}_{\max }=50M$ are used. The emissivity j is here defined as follows: $j={\rho }^{3}{p}^{-2}\exp {(-C({\rho }^{2}/({{Bp}}^{2}))}^{1/3})$, which captures the high-frequency limit of the thermal synchrotron emissivity $\nu \gg {\nu }_{c}{{\rm{\Theta }}}_{e}^{2}\sin (\theta )$ (compare with Equation (57) of Leung et al. 2011). The constant C is chosen such that the radiation is cut off after a few gravitational radii, resembling the expected millimeter emission from the Galactic Center, that is, C = 0.2 in geometrized units. This emission proxy is optimized for the science case of the EHT where near optically thin synchrotron emission is expected.¹²¹ In order to compare the variability, again a cadence of $1M$ or less is chosen in most cases.
• 4.  
  t, ϕ averages. Finally, we compare temporally and azimuthally averaged data for a more global impression of the disk and jet system:

  $\begin{eqnarray}&&\langle q\rangle (r,\theta ):= \displaystyle \frac{{\int }_{{t}_{\mathrm{beg}}}^{{t}_{\mathrm{end}}}{\int }_{0}^{2\pi }q(r,\theta ,\phi ,t)\,d\phi \,{dt}}{2\pi ({t}_{\mathrm{end}}-{t}_{\mathrm{beg}})}\end{eqnarray} \tag{ 19 }$

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

Time series data of horizon-penetrating fluxes is presented in Figures 2–4 for low, medium, and high resolutions respectively. Because the mass accretion governs the behavior of these fluxes, the data have been appropriately normalized by the accretion rate. All codes capture accurately the linear phase of the MRI leading to an onset of accretion to the BH at $t\simeq 300M$. While there is still a good correspondence of $\dot{M}$ for early times $\lt 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 ϕ and the energy fluxes. The low value of ϕ 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\mbox{--}6\lt {\phi }_{\max }$, all simulations are in the SANE regime of radiatively inefficient BH accretion. As a measure for the variance between codes, in Table 2 we quantify the peak fluxes as well as the average values far in the nonlinear regime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Quantifications—Time Series Data

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

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

Download table as:  ASCIITypeset image

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

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

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At low and medium resolution, the density profile in BHAC (blue curves) demonstrates an inner peak or "mini torus," which goes away in high resolution. As the additional analysis of the MRI quality reveals (see Section 5.6), the low- and mid-resolution BHAC runs do not properly resolve the fastest growing wavelengths ${\lambda }_{r},{\lambda }_{\theta }$, and hence, we consider the stable mini torus a numerical artifact, due to the lack of angular momentum transport in the underresolved simulations.

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

Boundary effects are also visible in the run of magnetic field and pressures: again, the three small domain models have the steepest pressure profiles, with the Athena++ showing clear kinks in several quantities at the boundary and the iharm3D measurement for $\langle {\beta }^{-1}\rangle$ rising sharply at the boundary (due to a few high ${\beta }^{-1}$ cells in the equatorial region close to the boundary; see, for example, also Section 5.4 and Figure 11).

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

Taking advantage of the scale freedom in the ideal MHD approximation, we have also performed a comparison with density, pressure, and magnetic field variables rescaled to the H-AMR N_θ = 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 β, the variance in the magnetic field is bound to increase.

To gain an overall impression of the solutions in the quasi-stationary state, we compare averages over the t- and $\phi$-coordinates in Figures 9–12. The different panels depict the rest-frame density, inverse plasma β, and the magnetization obtained for each code at the highest resolution available.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

As a large variety of floor models is employed, maps of plasma β and magnetization exhibit substantial spread in the funnel region. Aside from the difference in absolute values reached in the funnel (which are merely related to the floor level), there are also qualitative differences concerning the cells near the axis: depending on the pole treatment, the inverse plasma β decreases toward 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, Cartesian grids are also far from perfect: in the IllinoisGRMHD run, we observe some artifacts in the form of horizontal features in $\langle {\beta }^{-1}\rangle$ which are associated with resolution jumps in the computational grid.

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

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

$\begin{eqnarray}&&{\rm{\Psi }}(r,\theta )={\left(\displaystyle \frac{r}{{r}_{{\rm{h}}}}\right)}^{\kappa }(1-\cos (\theta )).\end{eqnarray} \tag{ 20 }$

In the recent 2D simulations of Nakamura et al. (2018), the jet boundary given by σ = 1 was accurately described by the choice κ = 0.75 for a wide range of initial conditions, BH spins, and horizon-penetrating flux: ${{\rm{\Phi }}}_{\mathrm{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\mbox{--}{10}^{5}M$ (e.g., Asada & Nakamura 2012; Hada et al. 2013). In our 3D SANE models, the line κ = 0.75 is recovered at low resolution by ECHO and iharm3D, and the more collimated genuinely paraboloid shape κ = 1 resembling the solution of Blandford & Znajek (1977) seems to be a better match for the funnel-wall shape at higher resolutions. Incidentally, the match with the N_θ = 1056 H-AMR run with the paraboloid shape is particularly good. We suspect that the narrower jet profile compared to Nakamura et al. (2018) is due to the lower $\phi \simeq 2\mbox{--}3$ of our benchmark configuration (at high resolution).

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

$\begin{eqnarray}&&\begin{array}{l}H/R\,(r)\\ := \,\displaystyle \frac{{\int }_{{t}_{\mathrm{beg}}}^{{t}_{\mathrm{end}}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\ | \pi /2-{\theta }_{\mathrm{KS}}| \ \rho (r,\theta ,\phi ,t)\sqrt{-g}\,{dt}\,d\phi \,d\theta }{{\int }_{{t}_{\mathrm{beg}}}^{{t}_{\mathrm{end}}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\rho (r,\theta ,\phi ,t)\sqrt{-g}\,{dt}\,d\phi \,d\theta },\end{array}\end{eqnarray} \tag{ 21 }$

where the averaging time is again taken between ${t}_{\mathrm{beg}}=5000M$ and ${t}_{\mathrm{end}}={\rm{10,000}}M$. As shown in Figure 15, a scale height of $H/R=0.25\mbox{--}0.3$ between ${r}_{\mathrm{KS}}\in [10,50]$ is consistently recovered with all codes. The largest departures from the general trend are seen in the inner regions for the Cartesian runs that show a slight "flaring up," which can also be observed in Figures 9 and 11. Furthermore, the presence of the outflow boundary in Athena++ leads to a decrease of density in the outer equatorial region (see Figure 9), which shows up as kink in Figure 15.

Zoom In Zoom Out Reset image size

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

In Sgr A*, episodic flares are detected in X-ray and near-infrared (NIR) on a roughly daily basis (Baganoff et al. 2001; Nowak et al. 2012; Mossoux & Grosso 2017), and the recent detection of orbital motion by the Gravity Collaboration et al. (2018) places the origin of an IR flare within 10 gravitational radii of the BH. In addition to large flares, a low-level continuous variability is present in NIR and radio (e.g., Dodds-Eden et al. 2011). While the nature of the strong coincident X-ray and IR flares is most likely associated with discrete events of particle heating or acceleration (Markoff et al. 2001; Dodds-Eden et al. 2009), weaker IR flares could also be caused by lensed thermal plasma in the turbulent accretion flow (Chan et al. 2015), and flares in the submillimeter band can be attributed solely to turbulent fluctuations in the accretion flow (Dexter et al. 2009, 2010). On short timescales (below the characteristic timescales $\tau \sim 4\,\mathrm{hr}$ in IR and $\tau \sim 8\,\mathrm{hr}$ in submillimeter), the continuous variability is well described by a red-noise spectrum with power-law slope of $\approx 2$ with no indication for a further break down to ∼2–10 minutes, delimited by the typical sampling (Do et al. 2009; Dexter et al. 2014; Witzel et al. 2018). As pointed out by Witzel et al. (2018), the IR and submillimeter characteristic timescales are in fact statistically compatible, indicating a direct relation of the emitting regions. For low-luminosity AGN, Bower et al. (2015) found that the submillimeter characteristic timescale is consistent with a linear dependence on the BH mass as would be expected for optically thin emission in radiatively inefficient general relativistic models. Hence, in M87, the variability with its characteristic timescale of ${45}_{-24}^{+61}$ days (Bower et al. 2015) can be obtained by scaling the Sgr A* result.

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

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

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

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

Zoom In Zoom Out Reset image size

Having validated the emission proxy as a reasonable descriptor for the flux in the millimeter and, we now compute the power-spectrum densities (PSDs) of the emission proxy ${ \mathcal L }$. The PSD is obtained by dividing the data into three non-overlapping segments between 5000 M and 10,000 M, and fast Fourier transformation of each Hamming-windowed segment, followed by 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, synthetic pure red-noise data with a PSD of ${f}^{-\beta }$ with $\beta \in \{2,3\}$ was also created and is faithfully recovered by the algorithm.¹²³ It is clear that as each window only contains $10\,\mathrm{hr}$, the observed break at the characteristic frequency $\tau \sim 8\mathrm{hr}$ (for Sgr A*) is not measurable by this procedure.

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

Zoom In Zoom Out Reset image size

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

To also compare the variability magnitude, we have computed the root mean square (rms) of the accretion rate after zeroing out all Fourier amplitudes with frequency below $1/1000M$. This serves as an effective detrending of the low-frequency secular evolution. In order to avoid edge effects, we compute the rms in the region $t\in [2500,7500]M$. Exemplary data of the remaining high-frequency variability is shown for the high-resolution runs in the left panel of Figure 18, where we have normalized by the average accretion rate in the region of interest. One can visually make out differences in variability amplitude with the BHAC data on the low end and HARM-Noble and iharm3D on the high end. The (normalized) rms values shown against midplane resolution (right panel) quantify this further and indicate quite a universal behavior with values in the range 1.2–1.6 across all codes and resolutions.¹²⁴ 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 ∼8 in the sample) is a restatement of the rms–flux relationship, which is a ubiquitous feature of BH accretion across all mass ranges (Uttley & McHardy 2001; Heil et al. 2012). The quotient $k\sim 1.3$ is consistent with the recent simulations of Hogg & Reynolds (2016) who find an rms–flux relationship of the accretion rate with $k=1.4\pm 0.4$.

Zoom In Zoom Out Reset image size

Summarizing the comparison of variabilities, although it remains a challenge to extract accurate power spectra from the finite time series of the data, all codes agree quite well on the salient features:

• 1.  
  A $\beta \simeq 2$ power-law slope is recovered for the emission proxy ${ \mathcal L }$ with all codes and resolutions.
• 2.  
  The PSD for ${ \mathcal L }$ and $\dot{M}$ steepens at $f\simeq 0.1{M}^{-1}$.
• 3.  
  The PSD of the accretion rate is shallower than the one for the proxy and is more accurately described by an index of $\beta \simeq 1$.
• 4.  
  The "rms–flux" relationship is quite universally recovered for $\dot{M}$ and ${ \mathcal L }$ for all codes and resolutions.

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

5.6.1. Maxwell Stress and α

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

$\begin{eqnarray}&&\langle {w}^{r\phi }\rangle (r,t):= \langle -{b}^{r}\sqrt{{\gamma }_{{rr}}}{b}^{\phi }\sqrt{{\gamma }_{\phi \phi }}\rangle \end{eqnarray} \tag{ 22 }$

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

$\begin{eqnarray}&&\alpha (r,t):= \displaystyle \frac{\langle {w}^{r\phi }\rangle (r,t)}{\langle {P}_{\mathrm{tot}}\rangle (r,t)}.\end{eqnarray} \tag{ 23 }$

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

Figure 19 illustrates the resolution dependence of the α-parameter for the H-AMR and BHAC runs. Overall, there is good agreement between BHAC and H-AMR for ${N}_{\theta }\gt 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 ∼3. In particular, there is no "stress edge" at or within the ISCO (located at ${r}_{\mathrm{KS}}=2.04M$), where α drops to zero. This is consistent with the results of Krolik et al. (2005) for the highly spinning case and marginally consistent with Penna et al. (2013), who find that $\alpha (r)$ peaks at ∼2–3M. The latter difference can be explained in part through the coordinate choice and in part through the higher spin adopted here (see also Appendix D). With increased resolution, the α in our global simulations increases, as opposed to what is observed in local shearing boxes. For example, Bodo et al. (2014) and Ryan et al. (2017) reported $\bar{\alpha }\propto {N}^{-1/3}$ for stratified local simulations and Guan et al. (2009) found $\bar{\alpha }\propto {N}^{-1}$ for the unstratified case (N is the number of grid cells per scale height). Over time, $\bar{\alpha }$ decreases, but less so as the grid resolution is increased. This secular trend and the large gap in resolution between the N_θ = 384 BHAC run and the N_θ = 1056 H-AMR run make a strict quantification of the convergence behavior difficult. At a comparatively high resolution of N_θ = 384, though, it seems that $\bar{\alpha }$ has not converged yet.

Zoom In Zoom Out Reset image size

5.6.2. Viscous Spreading

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

$\begin{eqnarray}&&\langle r\rangle (t):= \displaystyle \frac{{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{{r}_{{\rm{h}}}}^{{r}_{\max }}\ {r}_{\mathrm{KS}}\ \rho (r,\theta ,\phi ,t)\sqrt{-g}\,d\phi \,d\theta \,{dr}}{{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{{r}_{{\rm{h}}}}^{{r}_{\max }}\rho (r,\theta ,\phi ,t)\sqrt{-g}\,d\phi \,d\theta \,{dr}}\,,\end{eqnarray} \tag{ 24 }$

where the outer radius of integration has again been set to the domain of interest ${r}_{\max }=50M$. As more mass is expelled beyond $50M$; however, 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 has been employed. It is directly apparent that the disk size increases with resolution. Furthermore, the saturation of the low- and medium-resolution PPM runs occurring at $t\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_θ = 192.

Zoom In Zoom Out Reset image size

The behavior of $\langle r\rangle$ is consistent with the fact that higher resolution in the disk leads to larger Maxwell stresses, driving stronger disk winds and thus larger disk expansion. This lowers the overall disk density (as seen in Figures 6–8: $\langle \rho (r)\rangle$), but the disk's magnetic field is affected to a lesser extent (see Figures 6–8: $\langle | B| (r)\rangle$), hence the disk's α increases further in a nonlinear fashion. As shown in Section 5.7, however, the disk radius starts to converge for ${N}_{\theta }\gt 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" (Q-factors for short), defined as the number of cells available to resolve the fastest growing MRI mode in a given direction. While the fastest growing mode might be well resolved in the initial condition, whether this remains to be the case throughout the simulation can only be established a posteriori. If the simulation fails to capture at least the fastest growing mode, turbulent field amplification cannot proceed, and the field will decay due to numerical dissipation, leading also to a cessation of mass accretion. Although strictly a criterion for sufficient resolution of the linear MRI, the Q-factors also inform on the adequacy in the nonlinear regime (see, e.g., Sano et al. 2004; Noble et al. 2010; Hawley et al. 2013). The consensus values for adequate resolution are given by ${Q}^{(z)}\geqslant 10$, ${Q}^{(\phi )}\geqslant 20\mbox{--}25$ of Hawley et al. (2011, 2013); however, lower toroidal resolution requirements have also been suggested, for example, ${Q}^{(z)}\geqslant 10\mbox{--}15$ for ${Q}^{(\phi )}\approx 10$ in Sorathia et al. (2012). Furthermore, these authors note that the toroidal and poloidal resolutions are coupled, and low ${Q}^{(z)}$ can be compensated by ${Q}^{(\phi )}\geqslant 25$. Hence, the product of the quality factors ${Q}^{(z)}{Q}^{(\phi )}\geqslant 200\mbox{--}250$ has also been suggested (Narayan et al. 2012; Dhang & Sharma 2019) as a quality metric.

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

$\begin{eqnarray}&&{\lambda }^{(\theta )}:= \displaystyle \frac{2\pi }{\sqrt{(\rho h+{b}^{2}}){\rm{\Omega }}}{b}^{\mu }{e}_{\mu }^{(\theta )},\end{eqnarray} \tag{ 25 }$

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

$\begin{eqnarray}&&{\rm{\Delta }}{x}^{(\theta )}:= {\rm{\Delta }}{\theta }^{\mu }{e}_{\mu }^{(\theta )}\,,\end{eqnarray} \tag{ 26 }$

where ${\rm{\Delta }}{\theta }^{\mu }=(0,0,{\rm{\Delta }}\theta ,0)$ is the distance between two adjacent grid lines, and ${\rm{\Omega }}={u}^{\phi }/{u}^{t}$ is the angular velocity of the fluid. The resulting ${Q}^{(\theta )}:= {\lambda }^{(\theta )}/{\rm{\Delta }}{x}^{(\theta )}$ are shown for the standard resolutions at the final simulation times in the top panel of Figure 21. It is clear that the low- and medium-resolution simulations are not particularly well resolved with regard to the poloidal MRI wavelengths. The time and spatial average values in the domain of interest are ${Q}_{96}^{(r,\theta ,\phi )}\,\simeq 2.9,3,14$, ${Q}_{128}^{(r,\theta ,\phi )}\simeq 3,3.8,17$, and ${Q}_{192}^{(r,\theta ,\phi )}\simeq 7.3,7.9,19.7$. This is still somewhat below the nominal ${Q}^{(z,\phi )}\sim 10,20$ of Hawley et al. (2011) but well above the six cells suggested by Sano et al. (2004) to adequately capture the saturation level of the MRI. However, it seems that only the highest resolution case keeps MRI-driven turbulence alive up to the end of the simulation. At 96³, 128³, or 192³ resolution in the domain of interest, the bottleneck is clearly in the poloidal direction with ${Q}^{(r)}\simeq {Q}^{(\theta )}$ but the azimuthal factor ${Q}^{(\phi )}$ is approximately twice as large. Hence, one might surmise that the resolution in the $\phi$-direction will have little influence on the overall evolution. This has been confirmed by additional experiments with 96 instead of 192 azimuthal cells using BHAC and by runs where the r direction was also reduced using the KORAL code. In both cases, the agreement was very good, with horizon-penetrating fluxes consistent with the full resolution case. As noted by Hawley et al. (2011), an increase in azimuthal resolution can also lead to an improved behavior of the poloidal quality factors. We find the same; however, the effect is not dramatic. The dotted curve in Figure 21 shows that with half the azimuthal resolution, ${Q}^{(r)}$ and ${Q}^{(\theta )}$ indeed slightly decrease (to ${Q}^{(r,\theta )}\simeq 5.3,6.8$), and ${Q}^{(\phi )}\simeq 9.6$ is essentially halved. Upon increasing the resolution to N_θ = 384, we obtain ${Q}_{384}^{(r,\theta ,\phi )}\simeq 21.3,24.5,40.9$, and hence, a superlinear increase in the $Q-$factor compared to the N_θ = 192 case. This behavior will be encountered again for the very high-resolution run and is discussed further in the next section.

Zoom In Zoom Out Reset image size

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 ϕ shown in Figure 5 coincides with the point when Q-factors reach a sufficient level.

With H-AMR, we are able to produce a global simulation with a resolution of ${N}_{r}\times {N}_{\theta }\times {N}_{\phi }$ = 1608 × 1056 × 1024 with five levels of static mesh refinement in the ϕ-direction (see Section 4.1.5), bringing us to local shearing box regimes (e.g., Davis et al. 2010; Ryan et al. 2017). In terms of the density scale height $H/R\simeq 0.25\mbox{--}0.3$, we obtain 85–100 cells per scale height and disk-averaged quality factors ${Q}_{1056}^{(r,\theta ,\phi )}\geqslant 120,100,190$. This is more than sufficient for capturing MRI-driven turbulence as per the criteria (${Q}^{(z,\phi )}\geqslant 10,25$ in cylindrical coordinates). Indeed, Hawley et al. (2011) suggests a resolution of 600 × 450 × 200 for a global simulation with $H/R=0.1$, which we have achieved. Sufficiently resolved global simulations are essential for reproducing microphysical phenomena such as MRI in a macrophysical environment. Hawley et al. (2013) notes that even though MRI is a local instability, turbulent structures may form on larger scales (seen in large shearing boxes by Simon et al. 2012). Additionally, shearing box simulations inadvertently affect poloidal flux generation by implicitly conserving the total vertical magnetic flux, and hence, may not be as reliable as global simulations in dealing with large-scale poloidal flux generation (Liska et al. 2018b), which is required for jet launching and propagation.

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

This gives a first impression of the behavior of viscous angular momentum transport and the turbulent dissipation and generation of magnetic energy. As evidenced by Figure 22, there are signs of convergence in the disk radius at N_θ = 1056, though the disk magnetization continues to increase slightly with increased resolution. Due to the lack of intermediate data between N_θ = 192 and N_θ = 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_θ = 384. Quantitatively, the values of β ∼ 15–20 and the resolution-dependent trend reported in their Figure 3 seem consistent with our findings.

Zoom In Zoom Out Reset image size

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.

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

Balancing the desire to explore a large parameter space (as in the generation of a simulation library for EHT observations) against numerical satisfaction, one will end up performing simulations that are "just about" sufficiently resolved. Hence, it is interesting to note that with the choice of our resolutions, we have sampled a resolution edge in most codes. This is best seen in the runs of accretion rates and magnetic fluxes, which substantially reduce their scatter once midplane resolutions of ${\rm{\Delta }}{\theta }_{{\rm{p}}}\simeq 0.25\mbox{--}0.0125M$ are adopted. Below this threshold, the algorithmic differences still dominate, leading to large scatter with sometimes opposite trends between codes. The resolution threshold was also confirmed by inspection of the MRI quality factors in BHAC and H-AMR. Only the high-resolution PPM runs would qualify as resolving the saturated MRI per the criteria of Hawley et al. (2011).

Using our sample of runs at different resolutions, we have also investigated the convergence of global quantities of interest. While it is clear that convergence in the strict sense cannot be achieved for ideal MHD simulations of turbulence as viscous and resistive length scales depend on the numerical resolution and numerical scheme (e.g., Kritsuk et al. 2011), convergence in certain physical quantities such as the stress parameter α 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 α parameter and magnetization continue their weak resolution-dependent rise up to ∼1024³ cells. The latter confirms and extends the findings of Shiokawa et al. (2012) by a resolution factor of ×2.75.

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

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

Another systematic is introduced by the axial boundary condition: the jet–disk boundary σ = 1 is faithfully recovered only at high resolution, and the jet can be particularly skimpy when a polar cutout is used. In the low- and medium- resolution HARM-Noble and Cosmos++ simulations, we found that the axial excision significantly decreases the magnetization of the funnel region. The reason for the loss of the magnetized funnel is well known: standard "soft" reflective boundaries at the polar cones allow a finite (truncation error) numerical flux of B^r to leave the domain. This effect diminishes with resolution and can be counteracted by zero-flux "hard" boundaries at the polar cutout as noted by (Shiokawa et al. 2012). In Cosmos++, the use of outflow boundaries in the excised region leads to a particularly striking difference, resembling the situation described by Igumenshchev et al. (2003).

Inspecting the trends in the jet collimation profiles, we note that there is a close relationship between the jet width and the horizon-penetrating magnetic flux. Indeed, Table 2 and Figure 13 for the ${N}_{\theta }=128$ data reveal that a one-to-one relation holds for the four models with the smallest $\langle \phi \rangle$ and the one with the absolute largest value of $\langle \phi \rangle$ where the differences in opening angle are most pronounced. Because the jet collimation plays an important role in the acceleration of the bulk flow beyond the light cylinder (e.g., Camenzind 1986; Komissarov 2007), this systematic will translate into an effect on the overall flow acceleration of the funnel jets.

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

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

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

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

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

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

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