The Black Hole Accretion Code

We present the black hole accretion code (BHAC), a new multidimensional general-relativistic magnetohydrodynamics module for the MPI-AMRVAC framework. BHAC has been designed to solve the equations of ideal general-relativistic magnetohydrodynamics in arbitrary spacetimes and exploits adaptive mesh refinement techniques with an efficient block-based approach. Several spacetimes have already been implemented and tested. We demonstrate the validity of BHAC by means of various one-, two-, and three-dimensional test problems, as well as through a close comparison with the HARM3D code in the case of a torus accreting onto a black hole. The convergence of a turbulent accretion scenario is investigated with several diagnostics and we find accretion rates and horizon-penetrating fluxes to be convergent to within a few percent when the problem is run in three dimensions. Our analysis also involves the study of the corresponding thermal synchrotron emission, which is performed by means of a new general-relativistic radiative transfer code, BHOSS. The resulting synthetic intensity maps of accretion onto black holes are found to be convergent with increasing resolution and are anticipated to play a crucial role in the interpretation of horizon-scale images resulting from upcoming radio observations of the source at the Galactic Center.


Introduction
Accreting black holes (BHs) are amongst the most powerful astrophysical objects in the Universe. A substantial fraction of the gravitational binding energy of the accreting gas is released within tens of gravitational radii from the BH, and this energy supplies the power for a rich phenomenology of astrophysical systems including active galactic nuclei, X-ray binaries and gamma-ray bursts. Since the radiated energy originates from the vicinity of the BH, a fully general-relativistic treatment is essential for the modelling of these objects and the flows of plasma in their vicinity.
Depending on the mass accretion rate, a given system can be found in various spectral states, with different radiation mechanisms dominating and varying degrees of coupling between radiation and gas [1,2]. Some supermassive BHs, including the primary targets of observations by the Event-Horizon-Telescope Collaboration (EHTC [1] ), i.e., Sgr A* and M87, are accreting well below the Eddington accretion rate [3,4]. In this regime, the accretion flow advects most of the viscously released energy into the BH rather than radiating it to infinity. Such optically thin, radiatively inefficient and geometrically thick flows are termed advection-dominated accretion flows (ADAFs, see [5][6][7][8]) and can be modelled without radiation feedback. Next to the ADAF, two additional radiatively inefficient accretion flows (RIAFs) [1] http://www.eventhorizontelescope.org arXiv:1611.09720v2 [gr-qc] 18 Jun 2017 exist: The advection-dominated inflow-outflow solution (ADIOS) [9,10] and the convection-dominated accretion flow (CDAF) [11,12], which include respectively, the physical effects of outflows and convection. Analytical and semi-analytical approaches are reasonably successful in reproducing the main features in the spectra of ADAFs [see, e.g., 13]. However, numerical general-relativistic magnetohydrodynamic (GRMHD) simulations are essential to gain an understanding of the detailed physical processes at play in the Galactic Centre and other low-luminosity compact objects.
Modern BH accretion-disk theory suggests that angular momentum transport is due to MHD turbulence driven by the magnetorotational instability (MRI) within a differentially rotating disk [14,15]. Recent non-radiative GRMHD simulations of BH accretion systems in an ADAF regime have resolved these processes and reveal a flow structure that can be decomposed into a disk, a corona, a disk-wind and a highly magnetized polar funnel [see, e.g., [16][17][18][19]. The simulations show complex time-dependent behaviour in the disk, corona and wind. Depending on BH spin, the polar regions of the flow contain a nearly force-free, Poynting-flux-dominated jet [see, e.g., 17,18,20,21].
In addition to having to deal with highly nonlinear dynamics that spans a large range in plasma parameters, the numerical simulations also need to follow phenomena that occur across multiple physical scales. For example, in the MHD paradigm, jet acceleration is an intrinsically inefficient process that requires a few thousand gravitational radii to reach equipartition of the energy fluxes [22,23] (purely hydrodynamical mechanisms can however be far more efficient [24]). Jet-environment interactions like the prominent HST-1 feature of the radio-galaxy M87 [25][26][27] occur on scales of ∼ 5 × 10 5 gravitational radii. Hence, for a self-consistent picture of accretion and ejection, jet formation and recollimation due to interaction with the environment [see, e.g., 28], numerical simulations must capture horizon-scale processes, as well as parsec-scale interactions with an overall spatial dynamic range of ∼ 10 5 . The computational cost of such large-scale grid-based simulations quickly becomes prohibitive. Adaptive mesh refinement (AMR) techniques promise an effective solution for problems where it is necessary to resolve small and large scale dynamics simultaneously.
Another challenging scenario is presented by radiatively efficient geometrically thin accretion disks that mandate extreme resolution in the equatorial plane in order to resolve the growth of MRI instabilities. Typically this is dealt with by means of stretched grids that concentrate resolution where needed [29,30]. However, when the disk is additionally tilted with respect to the spin axis of the BH [31,32], lack of symmetry forbids such an approach. Here an adaptive mesh that follows the warping dynamics of the disk can be of great value. The list of scenarios where AMR can have transformative qualities due to the lack of symmetries goes on, the modelling of stardisk interactions [33], star-jet interactions [34], tidal disruption events [35], complex shock geometries [36,37], and intermittency in driven-turbulence phenomena [38,39], will benefit greatly from adaptive mesh refinement.
Qualitative aspects of BH accretion simulations are code-independent [see, e.g., 16,46,49], but quantitative variations raise questions regarding numerical convergence of the observables [58,67]. In preparation for the upcoming EHTC observations, a large international effort, whose European contribution is represented in part by the BlackHoleCam project [2] [68], is concerned with forward modelling of the future event horizon-scale interferometric observations of Sgr A* and M87 at submillimeter (EHTC; [69]) and near-infrared wavelengths (VLTI GRAVITY; [70]). To this end, GRMHD simulations have been coupled to general-relativistic radiative transfer (GRRT) calculations [see, e.g., [71][72][73][74][75][76]. In order to assess the credibility of these radiative models, it is necessary to assess the quantitative convergence of the underlying GRMHD simulations. In order to demonstrate the utility of BHAC for the EHTC science-case, we therefore validate the results obtained with BHAC against the HARM3D code [46,77] and investigate the convergence of the GRMHD simulations and resulting observables obtained with the GRRT post-processing code BHOSS [78].
The structure of the paper is as follows. In Sect. 2 we describe the governing equations and numerical methods. In Sect. 3 we show numerical tests in specialrelativistic and general-relativistic MHD. In Sect. 4 the results of 2D and 3D GRMHD simulations of magnetised accreting tori are presented. In Sect. 5 we briefly describe the GRRT post-processing calculation and the resulting image maps from the magnetised torus simulation shown in Sect. 4. In Sect. 6 we present our conclusions and outlook.
Throughout this paper, we adopt units where the speed of light, c = 1, the gravitational constant, G = 1, and the gas mass is normalised to the central compact object mass. Greek indices run over space and time, i.e., (0, 1, 2, 3), and Roman indices run over space only, i.e., (1,2,3). We assume a (−, +, +, +) signature for the spacetime metric. Self-gravity arising from the gas is neglected.

GRMHD formulation and numerical methods
In this section we briefly describe the covariant GRMHD equations, introduce the notation used throughout this paper, and present the numerical approach taken in our solution of the GRMHD system. The computational infrastructure underlying BHAC is the versatile open-source MPI-AMRVAC toolkit [79,80].
In-depth derivations of the covariant fluid-and magneto-fluid dynamical equations can be found in the textbooks by [40,81,82]. We follow closely the derivation of the GRMHD equations by [52]. This is very similar to the "Valencia formulation", cf. [40] and [50]. The general considerations of the "3+1" split of spacetime are discussed in greater detail in [83][84][85].
We start from the usual set of MHD equations in covariant notation ∇ µ (ρu µ ) = 0 , which respectively constitute mass conservation, conservation of the energymomentum tensor T µν , and the homogeneous Faraday's law. The Faraday tensor F µν may be constructed from the electric and magnetic fields E α , B α as measured in a generic frame U α as where η µνλδ is the fully-antisymmetric symbol (see, e.g., [40]) and g the determinant of the spacetime four-metric. The dual Faraday tensor * F µν := 1 2 (−g) −1/2 η µνλδ F λδ is then * We are interested only in the ideal MHD limit of vanishing electric fields in the fluid frame u µ , hence such that the inhomogeneous Faraday's law is not required and electric fields are dependent functions of velocities and magnetic fields. To eliminate the electric fields from the equations it is convenient to introduce vectors in the fluid frame and therefore we define the corresponding electric and magnetic field four-vectors as where e µ = 0 and we obtain the constraint u µ b µ = 0. The Faraday tensor is then and we can write the total energy-momentum tensor in terms of the vectors u µ and b µ alone [86] as Here the total pressure p tot = p + b 2 /2 was introduced, as well as the total specific enthalpy h tot = h + b 2 /ρ. In addition, we define the scalar b 2 := b ν b ν , denoting the square of the fluid frame magnetic field strength as b 2 = B 2 − E 2 .

3+1 split of spacetime
We proceed to split spacetime into 3+1 components by introducing a foliation into space-like hyper-surfaces Σ t defined as iso-surfaces of a scalar time function t. This leads to the timelike unit vector normal to the slices Σ t [40,85] n µ := −α∇ µ t , where α is the so-called lapse-function. The four-velocity n µ defines the frame of the Eulerian observer. If g µν is the metric associated with the four-dimensional manifold, we can define the metric associated with each timelike slice as This also allows us to introduce the spatial projection operator such that γ µ ν n µ = 0 and through which we can project any four-vector V µ (or tensor) into its temporal and spatial components. Introducing a coordinate system adapted to the foliation Σ t , the line element is given in 3+1 form [87] as where the spatial vector β µ is called the shift vector. Written in terms of coordinates, it describes the motion of coordinate lines as seen by an Eulerian observer More explicitly, we write the metric g µν and its inverse g µν as From (13) we find the following useful relation between the determinants of the 3-metric and 4-metric In a coordinate system specified by (11), the four-velocity of the Eulerian observer reads n µ = (−α, 0 i ), It is easy to verify that this normalised vector is indeed orthogonal to any spacelike vector on the foliation Σ t . Given a fluid element with four-velocity u µ , the Lorentz factor with respect to the Eulerian observer is [3] Γ := −u µ n µ = αu 0 and we introduce the three-vectors which denote the fluid three-velocity. In the following, purely spatial vectors (e.g., v 0 = 0) are denoted by Roman indices.
Further useful three-vectors are the electric and magnetic fields in the Eulerian frame which differ by a factor α from the definitions used in [46,88]. Writing the general Faraday tensor (2) in terms of quantities in the Eulerian frame, the ideal MHD condition (4) leads to the well known relation or put simply: E = B × v (here η ijk is the standard Levi-Civita antisymmetric symbol). Combining (6) with (17), one obtains the transformation between b µ and B µ as which enables the dual Faraday tensor (6) to be expressed in terms of the Eulerian fields * Equation (1)  or put more simply: ∇ · B = 0. Following (19) we obtain the scalar b 2 as where Using the spatial projection operator, the GRMHD equations (1) can be decomposed into spatial and temporal components. We skip ahead over the involved algebra [see e.g., 52] and directly state the final conservation laws with the conserved variables U and fluxes F i defined as where we define the transport velocity V i := αv i −β i . Hence we solve for conservation of quantities in the Eulerian frame: the density D := −ρu ν n ν , the covariant threemomentum S j , the rescaled energy density τ = U − D [4] (where U denotes the total energy density as seen by the Eulerian observer), and the Eulerian magnetic three-fields B j . The conserved energy density U is given by The purely spatial variant of the stress-energy tensor W ij was introduced for example in (23). It reads just as in special relativity Correspondingly, the covariant three-momentum density in the Eulerian frame is as usual. For the sources S we employ the convenient Valencia formulation without Christoffel symbols, yielding which is valid for stationary spacetimes that are considered for the remainder of this work (Cowling approximation). Following the definitions (23) and (30), all vectors and tensors are now specified through their purely spatial variants and thus apart from the occurrence of the lapse function α and the shift vector β i , the equations take on a form identical to the special-relativistic MHD (SRMHD) equations. This [4] Using τ = U −D instead of U improves accuracy in regions of low energy and enables one to consistently recover the Newtonian limit. fact allows for a straightforward transformation from the SRMHD physics module of MPI-AMRVAC into a full GRMHD code.
In addition to the set of conserved variables U , knowledge of the primitive variables P (U ) is required for the calculation of fluxes and source terms. They are given by While the transformation U (P ) is straightforward, the inversion P (U ) is a nontrivial matter which will be discussed further in Sect 2.10. Note that just like in MPI-AMRVAC , we do not store the primitive variables P but extend the conserved variables by the set of auxiliary variables where ξ := Γ 2 ρh. Knowledge of A allows for quick transformation of P (U ). The issue of inversion then becomes a matter of finding an A consistent with both P and U .

Finite volume formulation
Since BHAC solves the equations in a finite volume formulation, we take the integral of equation (22) over the spatial element of each cell dx 1 dx 2 dx 3 This can be written (cf. [89]) as with the volume averages defined as and We next define also the "surfaces" ∆S i and corresponding surface-averaged fluxes Considering that ∆V is assumed constant in time, this leads to the evolution equation We aim to achieve second-order accuracy and represent the interface-averaged flux, e.g.,F 1 ∂V (x 1 +∆x 1 /2) , with the value at the midpoint, change to an intuitive index notation F 1 i+1/2,j,k , and then arrive at a semi-discrete equation for the average state in the cell (i, j, k) as Here the source term S i,j,k is also evaluated at the cell barycenter to second-order accuracy [90]. Barycenter coordinatesx i are straightforwardly defined as This finite volume form is readily solved with the MPI-AMRVAC toolkit. For ease of implementation, we pre-compute all static integrals yielding cell volumes ∆V , Surfaces ∆S i and barycenter coordinates. The integrations are performed numerically at the phase of initialisation using a fourth-order Simpson's rule.
For the temporal update, we interpret the semi-discrete form (40) as an ordinary differential equation in time for each cell and employ a multi-step Runge-Kutta scheme to evolve the average state in the cellŪ i,j,k , a procedure also known as "method of lines". At each sub-step, the point-wise interface fluxes F i are obtained by performing a limited reconstruction operation of the cell-averaged stateŪ to the interfaces (see Sect. 2.8) and employing approximate Riemann solvers, e.g., HLL or TVDLF (Sect. 2.9).
Several temporal update schemes are available: simple predictor-corrector, thirdorder Runge-Kutta (RK) RK3 [91] and the strong-stability preserving s-step, pthorder RK schemes SSPRK(s,p) schemes: SSPRK(4,3), SSPRK(5,4) due to [92]. [5] 2.3 Metric data-structure The metric data-structure is built to be optimal in terms of storage while remaining convenient to use. Since the metric and its derivatives are often sparsely populated, the data is ultimately stored using index lists. For example, each element in the index list for the four-metric g µν holds the indices of the non-zero element together with a Fortran90 array of the corresponding metric coefficient for the grid block. A summation over indices, e.g., "lowering" can then be cast as a loop over entries in the index-list only. For convenience, all elements can also be accessed directly over intuitive identifiers which point to the storage in the index list, e.g., m%g(mu,nu)%elem yields the grid array of the g µν metric coefficients as expected. Similarly, the lowertriangular indices point to the transposed indices in the presence of symmetries. In addition, one block of zeros is allocated in the metric data-structure and all zero elements are set to point towards it. An overview of the available identifiers is given in Table 1.
As a consequence, only 14 grid functions are required for the Schwarzschild coordinates and 29 grid functions need to be allocated in the Kerr-Schild (KS) case. This is still less than half of the 68 grid functions which a brute-force approach would yield. The need for efficient storage management becomes apparent when we consider that the metric is required in the barycenter as well as on the interfaces, thus multiplying the required grid functions by a factor of four for three-dimensional simulations (yielding 116 grid functions in the KS case).
In order to eliminate the error-prone process of implementing complicated functions for metric derivatives, BHAC can obtain derivatives by means of an accurate complex-step numerical differentiation [93]. This elegant method takes advantage of the Cauchy-Riemann differential equations for complex derivatives and achieves full double-precision accuracy, thereby avoiding the stepsize dilemma of common finite-differencing formulae [94]. The small price to pay is that at the initialisation stage, metric elements are provided via functions of the complexified coordinates. However, the intrinsic complex arithmetic of Fortran90 allows for seamless implementation. [5] For implementation details, see [80].
To promote full flexibility in the spacetime, we always calculate the inverse metric γ ij using the standard LU decomposition technique [95]. As a result, GRMHD simulations on any metric can be performed after providing only the non-zero elements of the three-metric γ ij (x 1 , x 2 , x 3 ), the lapse function α(x 1 , x 2 , x 3 ) and the shift vector β i (x 1 , x 2 , x 3 ). As an additional convenience, BHAC can calculate the required elements and their derivatives entirely from the four-metric g µν (x 0 , x 1 , x 2 , x 3 ).

Equations of state
For closure of the system (1)-(4), an equation of state (EOS) connecting the specific enthalpy h with the remaining thermodynamic variables h(ρ, p) is required [40]. The currently implemented closures are , where the relativistic temperature is given by Θ = p/ρ and K n denotes the modified Bessel function of the second kind. In fact, we use an approximation to the previous expression that does not contain Bessel functions [see 79,96]. • Isentropic flow : Assumes an ideal gas with the additional constraint p = κργ, where the pseudo-entropy κ may be chosen arbitrarily. This allows one to omit the energy equation entirely and only the reduced set P = {ρ, v j , B j } is solved.
As long as h(ρ, p) is analytic, its implementation in BHAC is straightforward.

Divergence cleaning and augmented Faraday's law
To control the ∇ · B = 0 constraint on AMR grids, we have adopted a constraint dampening approach customarily used in Newtonian MHD [97]. In this approach, which is usually referred as Generalized Lagrangian Multiplier (GLM) of the Maxwell equations (but is also known as the "divergence-cleaning" approach), we extend the usual Faraday tensor by the scalar φ, such that the homogeneous Maxwell equation reads and the scalar φ follows from contraction φ = ( * F µν − φg µν )n µ n ν . Naturally, for φ → 0, the usual set of Maxwell equations is recovered. It is straightforward to show [see, e.g., 98] that (42) leads to a telegraph equation for the constraint violation parameter φ which becomes advected at the speed of light and decays on a timescale 1/κ. With the modification (42), the time-component of Maxwell's equation now becomes an evolution equation for φ. After some algebra (see Appendix A), we obtain Equivalently, the modified evolution equations for B i (see Appendix B) read Now equation (44) replaces the usual Faraday's law and (43) is evolved alongside the modified MHD system. Due to the term ∂ i φ on the right hand side of equation (44), the new equation is non-hyperbolic. Hence, numerical stability can be a more involved issue than for hyperbolic equations. We find that the numerical stability of the system is enhanced when using an upwinded discretisation for ∂ i φ. Note that Equations (43) and (44) are in agreement with [63] when accounting for ∂i √ γ √ γ = 1 2 γ lm ∂ i γ lm and taking the ideal MHD limit.

Flux-interpolated Constrained Transport
As an alternative to the GLM approach, the ∇ · B = 0 constraint can be enforced using a cell-centred version of Flux-interpolated Constrained Transport (FCT) consistent with the finite volume scheme used to evolve the hydrodynamic variables. Constrained Transport (CT) schemes aim to keep to zero at machine precision the sum of the magnetic fluxes through all surfaces bounding a cell, and therefore (in the continuous limit) the divergence of the magnetic field inside the cell. In the original version [99] this is achieved by evolving the magnetic flux through the cell faces and computing the circulation of the electric field along the edges bounding each face. Since each edge appears with opposite signs in the time update of two faces belonging to the same cell, the total magnetic flux leaving each cell is conserved during evolution. The magnetic field components at cell centers, necessary for performing the transformation from primitive to conserved variables and viceversa, are then found using interpolation from the cell faces. [100] showed that it is possible to find cell centred variants of CT schemes that go from the average field components at the cell center at a given time to those one (partial) time step ahead in a single step, without the need to compute magnetic fluxes at cell faces. The CT variant known as FCT is particularly well suited for finite volume conservative schemes as that employed by BHAC, as it calculates the electric fields necessary for the update as an average of the fluxes given by the Riemann solver. In this way, the time update for its cell centred version can be written using a form similar to (40). For example, for the update of theB 1 component, we obtain where the modified fluxes in the x 1 -direction are zero and the remaining fluxes are calculated as The derivation of equations (45) and (46) from the staggered version with magnetic fields located at cell faces is given in Appendix C. Since magnetic fields are stored at the cell center and not at the faces, the divergence conserved by the FCT method corresponds to a particular discretisation where ∆V * | i+1/2,j+1/2,k+1/2 = l1,l2,l3=0,1 Equation (47) is closely related to the integral over the surface of a volume containing eight cells in 3D (see Appendix D for the derivation), and it reduces to equation (27) from [100] in the special case of Cartesian coordinates. As mentioned before, this scheme can maintain ∇ · B = 0 to machine precision only if it was already zero at the initial condition. The corresponding curl operator used to setup initial conditions is derived in Appendix D.
In its current form, BHAC cannot handle both constrained transport and AMR. The reason is that special prolongation and restriction operators are required in order to avoid the creation of divergence when refining or coarsening. Due to the lack of information about the magnetic flux on cell faces, the problem of finding such divergence-preserving prolongation operators becomes underdetermined. However, storing the face-allocated (staggered) magnetic fluxes and applying the appropriate prolongation and restriction operators requires a large change in the code infrastructure on which we will report in an accompanying work.

Coordinates
Since one of the main motivations for the development of the BHAC code is to simulate BH accretion in arbitrary metric theories of gravity, the coordinates and metric data-structures have been designed to allow for maximum flexibility and can easily be extended. A list of the currently available coordinate systems is given in Table 2. In addition to the identifiers used in the code, the table lists whether numerical derivatives are used and whether the coordinates are initialised from the three-metric or the four-metric. The less well-known spacetimes and coordinates are described in the following subsection.

Modified Kerr-Schild coordinates
Modified KS coordinates were introduced by e.g., [17] with the purpose of stretching the grid radially and being able to concentrate resolution in the equatorial region.
The original coordinate transformation is equivalent to: where R 0 and h are parameters which control, respectively, how much resolution is concentrated near the horizon and near the equator. Unfortunately, the inverse of ϑ(θ) is a transcendental equation that has to be solved numerically. To avoid this complication and still capture the functionality of the modified coordinates, we instead use the following θ− transformation Now the solution to the cubic equation can be expressed in closed-form, and the only real root reads where This is compared with the original version (50) in Fig. 1 and shows a good match between the two versions of modified Kerr-Schild coordinates. The radial back-transformation follows trivially as and the derivatives for the diagonal Jacobian are With these transformations, we obtain the new metric g MKS = J T g KS J. Note that whenever the parameters h = 0 and R 0 = 0 are set, our MKS coordinates reduce to the standard logarithmic Kerr-Schild coordinates.

Rezzolla & Zhidenko parametrization
The Rezzolla-Zhidenko parameterisation [101] has been proposed to describe spherically-symmetric BH geometries in metric theories of gravity. In particular, using a continued-fraction expansion (Padé expansion) along the radial coordinate, deviations from general relativity can be expressed using a small number of coefficients. The line element reads with N (r) and B(r) being functions of the radial coordinate r. The radial position of the event horizon is fixed at r = r 0 > 0 which implies that N (r 0 ) = 0. Furthermore, the radial coordinate is compactified by means of the dimensionless coordinate in which x = 0 corresponds to the position of the event horizon, while x = 1 corresponds to spatial infinity. Through this dimensionless coordinate, the function N can be written as where A(x) > 0 for 0 ≤ x ≤ 1. Introducing additional coefficients , a n , and b n , the metric functions A and B are then expressed as follows Here A and B are functions describing the metric near the event horizon and at spatial infinity. In particular, A and B have rapid convergence properties, that is by Padé approximants where a 1 , a 2 , a 3 , . . . and b 1 , b 2 , b 3 , . . . are dimensionless coefficients that can, in principle, be constrained from observations. The dimensionless parameter is fixed by the ADM mass M and the coordinate of the horizon r 0 . It measures the deviation from the Schwarzschild case as It is easy to see that at spatial infinity (x = 1), all coefficients contribute to (62), while at event horizon only the first two terms remain, i.e.
Given a number of coefficients, any spherical spacetime can hence directly be simulated in BHAC. For example, the coefficients in the Rezzolla-Zhidenko parametrization for the Johannsen-Psaltis [103] metric and for Einstein-Dilaton BHs [104] have already been provided in [101]. Typically, expansion up to a 2 , b 2 yields sufficient numerical accuracy for the GRMHD simulations. The first simulations in the related horizon penetrating form of the Rezzolla-Zhidenko parametrization are discussed in [105].

Available reconstruction schemes
The second-order finite volume algorithm (40) requires numerical fluxes centered on the interface mid-point. As in any Godunov-type scheme [see e.g., 88,106], the fluxes are in fact computed by solving (approximate) Riemann problems at the interfaces (see Sect. 2.9). Hence for higher than first-order accuracy, the fluid variables need to be reconstructed at the interface by means of an appropriate spatial interpolation. Our reconstruction strategy is as follows. 1) Compute primitive variablesP from the averages of the conserved variablesŪ located at the cell barycenter. 2) Use the reconstruction formulae to obtain two representations for the state at the interface, one with a left-biased reconstruction stencil P L and the other with a right-biased stencil P R . 3) Convert the now point-wise values back to their conserved states U L and U R . The latter two states then serve as input for the approximate Riemann solver. A large variety of reconstruction schemes are available, which can be grouped into standard second-order total variation diminishing (TVD) schemes like "minmod", "vanLeer", "monotonized-central", "woodward" and "koren" [see 79, for details] and higher order methods like the third-order methods "PPM" [107], "LIMO3" [108] and the fifth-order monotonicity preserving scheme "MP5" due to [109]. While the overall order of the scheme will remain second-order, the higher accuracy of the spatial discretisation usually reduces the diffusion of the scheme and improves accuracy of the solution [see, e.g., 80]. For typical GRMHD simulations with nearevacuated funnel/atmosphere regions, we find the PPM reconstruction scheme to be a good compromise between high accuracy and robustness. For simple flows, e.g., the stationary toroidal field torus discussed in Sect. 3.4, the compact stencil LIMO3 method is recommended.

Characteristic speed and approximate Riemann solvers
The time-update of BHAC proceeds in a dimensionally unsplit manner, thus at each Runge-Kutta substep the interface-fluxes in all directions are computed based on the previous substep. The state is then advanced to the next substep with the combined fluxes of the cell. To compute these fluxes from the reconstructed conserved variables at the interface U L and U R , we provide two approximate Riemann solvers: 1) the Rusanov flux, also known as Total variation diminishing Lax-Friedrichs scheme (TVDLF) which is based on the largest absolute value of the characteristic waves normal to the interface c i , and 2) the HLL solver [110], which is based on the leftmost (c i − ) and rightmost (c i + ) waves of the characteristic fan with respect to the interface. The HLL upwind flux function for the conserved variable u ∈ U is calculated as and we set in accordance with [111]: with c i = max |c i − |, |c i + | . In addition to these two standard approximate Riemann solvers, we also provide a modified TVDLF solver that preserves positivity of the conserved density D. The algorithm was first described in the context of Newtonian hydrodynamics by [112] and was successfully applied in GRHD simulations by [113]. It takes advantage of the fact that the first-order Lax-Friedrichs flux F i,LO (u) is positivity preserving under a CFL condition CFL ≤ 1/2. Hence the fluxes can be constructed by combining the high order flux F i,HO (u) (obtained e.g., by PPM reconstruction) and F i,LO (u) such that the updated density does not fall below a certain threshold. [6] Specifically, the modified fluxes read where θ ∈ [0, 1] is chosen as a maximum value which ensures positivity of the cells adjacent to the interface (see [112] for details of its construction). Note that although we only stipulate the density be positive, the formula (68) must be applied to all conserved variables u ∈ U .
In relativistic MHD, the exact form of the characteristic wave speeds λ ± involves solution of a quartic equation [see, e.g., 86] which can add to the computational overhead. For simplicity, instead of calculating the exact characteristic velocities, we follow the strategy of [46] who propose a simplified dispersion relation for the fast MHD wave ω 2 = a 2 k 2 . As a trade-off, the simplification can overestimate the wavespeed in the fluid frame by up to a factor of 2, yielding a slightly more diffusive behaviour. The upper bound a for the fast wavespeed is given by which depends on the usual sound speed and Alfvén speed here given for an ideal EOS with adiabatic indexγ. As pointed out by [52], the 3+1 structure of the fluxes leads to characteristic waves of the form where λ i ± is the characteristic velocity in the corresponding special relativistic sys- For the simplified isotropic dispersion relation, the characteristics can then be obtained just like in special relativistic hydrodynamics [see, e.g., 89,114,115] . [6] In the general-relativistic hydrodynamic WhiskyTHC code [54,55], this desirable property allows to set floors on density close to the limit of floating point precision ∼ 10 −16 ρ ref .

Primitive variable recovery
It is well-known that the nonlinear inversion P (U ) is the Achilles heel of any relativistic (M)HD code and sophisticated schemes with multiple backup strategies have been developed over the years as a consequence (e.g., [116], [117], [77], [118], [119], [120]). Here we briefly describe the methods used throughout this work and refer to the previously mentioned references for a more detailed discussion.

Primary inversions
Two primary inversion strategies are available in BHAC. The first strategy, which we denote by "1D", is a straightforward generalisation of the one-dimensional strategy described in [121]. It involves a non-linear root finding algorithm which is implemented by means of a Newton-Raphson scheme on the auxiliary variable ξ. Once ξ is found, the velocity follows from (29) v and we calculate the second auxiliary variable Γ The thermal pressure p then follows from the particular EOS in use (Sect. 2.4). For example, for an ideal EOS we have For details of the consistency checks and bracketing, we refer the interested reader to [121].
In addition to the 1D scheme, we have implemented the "2DW" method of [52,116]. The 2DW inversion simultaneously solves the non-linear equations (25) and the square of the three-momentum S 2 , following (29) by means of a Newton-Raphson scheme on the two variables ξ and v 2 . Among all inversions tested by [116], the 2DW method was reported as the one with the smallest failure rate. We find the same trend, but also find that the lead of 2DW over 1D is rather minor in our tests.
With two distinct inversions that might fail under different circumstances, one can act as a backup strategy for the other. Typically we first attempt a 2DW inversion and switch to the 1D method when no convergence is found. The next layer of backup can be provided by the entropy method as described in the next section.

Entropy switch
To deal with highly magnetised regions, [60,77] introduced the advection of entropy to provide a backup strategy for the primitive variable recovery. Similar to [60,77], alongside the usual fluid equations, BHAC can be configured to solve an advection equation for the entropy S where we define given the adiabatic indexγ. This leads to the evolution equation for the conserved quantity ΓS. The primitive counterpart is the actual entropy κ = p/ργ, which can be recovered via κ = ΓS/D. In case of failure of the primary inversion scheme, using the advected entropy κ, we can attempt a recovery of primitive variables which does not depend on the conserved energy. Note that after the primitive variables are recovered from the entropy, we need to discard the conserved energy and set it to the value consistent with the entropy. On the other hand, after each successful recovery of primitive variables, the entropy is updated to κ = p/ργ, which is then advected to the next step. In addition, entropy-based inversion can be activated whenever β = 2p/b 2 ≤ 10 −2 since the primary inversion scheme is likely to fail in these highly magnetised regions. Tests of the dynamic switching of the evolutionary equations are described in Sect. 3.3. In GRMHD simulations of BH accretion, the "entropy region" is typically located in the BH magnetosphere, which is strongly magnetised and the error due to missing shock dissipation is thus expected to be small. In the rare instances where the entropy inversion also fails to converge to a physical solution, the code is normally stopped. To force a continuation of the simulation, last resort measures that depend on the physical scenario can be employed. Often the simulation can be continued when the faulty cell is replaced with averages of the primitive variables of the neighbouring healthy cells as described in [79]. In the GRMHD accretion simulations described below, failures could happen occasionally in the highly magnetised evacuated "funnel" region close to the outer horizon where the floors are frequently applied. We found that the best strategy is then to replace the faulty density and pressure values with the floor values and set the Eulerian velocity to zero. Note that in order to avoid generating spurious ∇ · B, the last resort measures should never modify the magnetic fields of the simulation.

Adaptive Mesh refinement
The computational grid employed in BHAC is provided by the MPI-AMRVAC toolkit and constitutes a fully adaptive block based (oct-) tree with a fixed refinement factor of two between successive levels. That is, the domain is first split into a number of blocks with equal amount of cells (e.g., 10 3 computational cells per block). Each block can be refined into two (1D), four (2D) or eight (3D) child-blocks with an again fixed number of cells. This process of refinement can be repeated ad libitum and the data-structure can be thought of a forest (collection of trees). All operations on the grid, for example time-update, IO and problem initialisation are scheduled via a loop over a space-filling curve. We adopt the Morton Z-order curve for ease of implementation via a simple recursive algorithm.
Currently, all cells are updated with the same global time-step and hence loadbalancing is achieved by cutting the space-filling curve into equal sections that are then distributed over the MPI-processes. The AMR strategy just described is applied in various astrophysical codes, for example codes employing the PARAMESH library [122][123][124], or the recent Athena++ framework [see, e.g., 125]. Compared to a patch-based approach [see, e.g., 126], the block based AMR has several advantages: 1) well-defined boundaries between neighbouring grids on different levels ,2) data is uniquely stored and updated, thus no unnecessary interpolations are performed, and 3) simple data-structure, e.g., straightforward integer arithmetic can be used to locate a particular computational block. For in-depth implementation details such as refinement/prolongation operations, indexing and ghost-cell exchange, we refer to [79]. Prolongation and restriction can be used on conservative variables or primitive variables. Typically primitive variables are chosen to avoid unphysical states which can otherwise result from the interpolations in conserved variables. The refinement criteria usually adapted is the Löhner's error estimator [127] on physical variables. It is a modified second derivative, normalised by the average of the gradient over one computational cell. The multidimensional generalization is given by The indices p, q run over all dimensions p, q = 1, ..., N D . The last term in the denominator acts as a filter to prevent refinement of small ripples, where f wave is typically chosen of order 10 −2 . This method is also used in other AMR codes such as FLASH [128], RAM [124], PLUTO [126] and ECHO [66].

Shock tube test with gauge effect
The first code test is considered in flat spacetime and therefore no metric source terms are involved. Herein we perform one-dimensional MHD shock tube tests with gauge effects by considering gauge transformations of the spacetime. Shock tube tests are well-known tests for code validation and emphasise the nonlinear behaviour of the equations, as well as the ability to resolve discontinuities in the solutions [see, e.g., 50,52].
The initial condition is given as and all other quantities are zero. In order to check whether the covariant fluxes are correctly implemented, we use different settings for the flat spacetime as detailed in Table 3. In the simulations, an ideal gas EOS is employed with an adiabatic index ofγ = 2. The 1D problem is run on a uniform grid in x−direction using 1024 cells spanning over x ∈ [−1/2, 1/2]. The simulations are terminated at t = 0.4. For the spatial reconstruction, we adopt the second order TVD limiter due to Koren [129]. Furthermore, RK3 timeintegration is used with Courant number set to 0.4. Case A is the reference solution without modification of fluxes due to the threemetric, lapse or shift [7] . By means of simple transformations of flat-spacetime, all other cases can be matched with the reference solution. Case B will coincide with solution A if B is viewed at t/2 = 0.2. Case C will agree with case A when it is shifted in positive x−direction by δx = β x t = 0. 16 In general, all cases agree very well with the rescaled solution. To give an example, Fig. 2 shows the rescaled simulation results of case F compared to the reference solution of case A. This test demonstrates the shock-capturing ability of the MHD code and enables us to conclude that the calculation of the covariant fluxes has been implemented correctly.

Boosted loop advection
In order to test the implementation of the GLM-GRMHD system, we perform the advection of a force-free flux-tube with poloidal and toroidal components of the magnetic field in a flat spacetime. [7] We note that for the reference solution we have relied here on the extensive set of tests performed in flat spacetime within the MPI-AMRVAC framework; however, we could also have employed as reference solution the "exact" solution as derived in Ref. [130].
The initial equilibrium configuration of a force-free flux-tube is given by a modified Lundquist tube [see e.g., 131], where we avoid sign changes of the vertical field component B z with the additive constant C = 0.01. Pressure and density are initialized as constant throughout the simulation domain. The initial pressure value is obtained from the central plasma-beta β 0 = B 2 (0)/2p, where B is the magnetic field in the co-moving system. The density is set to ρ = p/2 yielding a relativistic hot plasma. Consequently, an adiabatic index γ = 4/3 is used. We set β 0 = 0.01, which results in a high magnetisation σ 0 = B 2 (0)/(ρc 2 + 4p) 25. The equations for the magnetic field for r < 1 read and where t can be set to zero and where we assumed a vanishing electric field in the co-moving system. Therefore relativistic length contraction gives the loop a squeezed appearance. A simulation domain (x, y) ∈ [−1, 1] at a base resolution of N x × N y = 64 2 is initialised with an additional three levels of refinement. We advect the loop for one period (P = 2 √ 2/v) across the domain where periodic boundary conditions are used.
The advection over the coordinates can be counteracted by setting the shift vector appropriately, i.e. β = −v. This is an important consistency check of the implementation. Figure 3 shows the initial and final states of the force-free magnetic flux-tube advected for one period and for the case with spacetime shifted against the advection velocity. The advected and counter-shifted cases are in good agreement, with only the truly advected case being slightly more diffused, the effect of which is reflected in the activation of more blocks on the third AMR level.
To investigate the numerical accuracy the L 1 and L ∞ norms of the out-of-plane magnetic field component B z , as well as the divergence of magnetic field between the initial state and the simulation at a time after one advection period with different resolutions as seen in Fig.4 are checked. The error norms from analytically known solutions u * are defined as where the summation, respectively maximum operation, includes all cells in the domain and the integrals are performed over the volume of the cell ∆V i,j,k . In this sense, the reported errors correspond to the mean and maximal error in the computational domain. We should note that for the test of convergence, we use a uniform grid and choose v = 0.5 √ 0.5(1, 1, 0), β = √ 0.5(1, 0, 0) resulting in an advection in direction of the upper-left diagonal. A TVD "Koren" limiter is chosen. As expected, the convergence is second order for all cases.

Magnetised spherical accretion
A useful stress test for the conservative algorithm in a general-relativistic setting is spherical accretion onto a Schwarzschild BH with a strong radial magnetic field [46]. The steady-state solution is known as the Michel accretion solution [132] and represents the extension to general relativity of the corresponding Newtonian solution by [133]. The steady-state spherical accretion solution in general relativity is described in a number of works [see, e.g., 40,43]. It is easy to show that the solution is not affected when a radial magnetic field of the form B r ∝ γ −1/2 is added [45]. This test challenges the robustness of the code and of the inversion procedure P (U ) in particular. The calculation of the initial condition follows that outlined in [43]. Here, we parametrize the field strength via σ = b 2 /ρ at the inner edge of the domain (r = 1.9 M ). The simulation is setup in the equatorial plane using MKS coordinates corresponding to a domain of r KS ∈ [1.9, 20] M . The analytic solution remains fixed at the inner and outer radial boundaries. We run two cases, case 1 with magnetisation up to σ = 10 3 and case 2 with a very high magnetisation reaching up to σ = 10 4 . Since the problem is only 1D, the constraint ∇ · B = 0 has a unique solution which gets preserved via the FCT algorithm. Figure 5 illustrates the profiles for σ = 10 3 and two inversion strategies: 2DW (black +) and 2DW with entropy switching in regions of high magnetization b 2 /2p >

(red ×).
With the exception of the radial three-velocity near the BH horizon (r ≤ 5 M ), in both cases the simulations maintain well the steady-state solution. [8] Comparing theses results with and without entropy switching, the entropy strategy actually keeps the solution closer to the steady-state solution (solid black curves) even though the change of inversion strategy occurs in the middle of the domain, r 10. Figure 6 Error of density ρ in the highly magnetised Bondi flow: σ = 10 3 (left) and σ = 10 4 (right). The black data points are obtained with the standard 2D inversion and the red datapoints switch to the entropy inversion at β ≤ 10 −2 . One can see that both recipes are convergent with the expected order and that the error in the entropy strategy is decreased by roughly a factor of two.
The errors (L 1 and L ∞ norms) for the four cases are shown in Fig. 6. Again, the second-order accuracy of the algorithm is recovered. Using the entropy strategy increases the numerical accuracy by around a factor of two and we suggest its use in the highly magnetised regime of BH magnetospheres.

Magnetised equilibrium torus
As a final validation of the code in the GRMHD regime, we perform the simulation of a magnetised equilibrium torus around a spinning BH. A generalisation of the steady-state solution of the standard hydrodynamical equilibrium torus with constant angular momentum [see, e.g., 43,134,135] to MHD equilibria with toroidal magnetic fields was proposed by [136]. This steady-state solution is important since it constitutes a rare case of a non-trivial analytic solution in GRMHD [9] .
For the initial setup of the equilibrium torus, we adopt a particular relationship ω = ω(p), where ω = ρh is the fluid enthalpy andω =ω(p m ), wherep m = Lp m , ω = Lω, p m = b 2 /2 is the magnetic pressure, and L = g tφ g tφ − g tt g φφ . From these relationships, thermal and magnetic pressures are described as The analytical solutions can be constructed from [8] Note that the discrepancy in v r appears less dramatic when viewed in terms of the four-velocity u r . [9] We thank Chris Fragile for providing subroutines for this test case. for the introduced total potential W , where W = ln |u t |. The centre of the torus is located at (r c , π/2). At this point, we parametrize the magnetic field strength in terms of the pressure ratio β c = p g (r c , π/2)/p m (r c , π/2) .
The gas pressure and magnetic pressure at the centre of the torus are given by From these, the constants K and K m for barotropic fluids are obtained. The magnetic field distribution is given by the distribution of magnetic pressure p m . From the consideration of a purely toroidal magnetic field one obtains where A = g φφ + 2lg tφ + l 2 g tt and := −u φ /u t is the specific angular momentum. We perform 2D simulations using logarithmic KS coordinates with h = 0 and R 0 = 0. The simulation domain is θ ∈ [0, π], r KS ∈ [0.95 r h , 50 M ], where r h is the (outer) event horizon radius of the BH. The BH has the dimensionless spin parameter a = 0.9. For simplicity, we set the two indices to the same value of κ = η = 4/3 and also set the adiabatic index of the adopted ideal EOS to γ = 4/3. The remaining parameters are listed in the Table 4.
Initially, the velocity of the atmosphere outside of the torus is set to zero in the Eulerian frame, with density and gas pressure set to very small values of ρ = ρ min r −3/2 BL , p = p min r −5/2 BL with ρ min = 10 −5 and p min = 10 −7 . It is important to note that the atmosphere is free to evolve and only densities and pressures are floored according to the initial state. In the simulations we use the HLL approximate Riemann solver, third order LIMO3 reconstruction, two-step time update, and a CFL number of 0.5. We impose outflow conditions on the inner and outer boundaries of the radial direction and reflecting boundary conditions in the θ direction. As the magnetic field is purely toroidal, and will remain so during the time-evolution of this axisymmetric case, no particular ∇ · B = 0 treatment is used.
The top panels of Fig. 7 show the density distribution at the initial state and at t = 200 M , as well as the plasma beta distribution at t = 200 M . The rotational period of the disk centre is t r = 68 M . The initial torus configuration is well maintained after several rotation period. For a qualitative view of the simulations, the 1D radial and azimuthal distributions of the density are shown in the lower two panels in Fig. 7 with different grid resolutions. All but the low resolution case are visually indistinguishable from the initial condition in the bottom-left panel, showing ρ − r with a linear scale. Since the atmosphere is evolved freely, small density waves propagate in the ambient medium of the torus, as seen in the ρ − θ cut. This does  not adversely affect the equilibrium solution in the bulk of the torus however. Error quantification (L 1 and L ∞ ) is provided in Fig. 8. The second-order properties of the numerical scheme are well recovered.

Differences between FCT and GLM
Having implemented two methods for divergence control, we took the opportunity to compare the results of simulations using both methods. We analysed three tests: a relativistic Orszag-Tang vortex, magnetised Michel accretion, and magnetised accretion from a Fishbone-Moncrief torus. Although much less in-depth, this comparison is in the same spirit as those performed in previous works in non-relativistic MHD [100,137,138]. The well-known work by [100] compares seven divergencecontrol methods, including an early non-conservative divergence-cleaning method known as the eight-wave method [139], and three CT methods, finding that FCT is among the three most accurate schemes for the test problems studied. In [137], three divergence-cleaning schemes and one CT scheme were applied to the same test problem of supernova-induced MHD turbulence in the interstellar medium. It was found that the three divergence-cleaning methods studied suffer from, among other problems, spurious oscillations in the magnetic energy, which is attributed to the non-locality introduced by the loss of hyperbolicity in the equations. Finally, in [138], a non-staggered version of CT adapted to a moving mesh is compared to the divergence-cleaning Powell scheme [140], an improved version of the eight-wave method. They observe greater numerical stability and accuracy, and a better preservation of the magnetic field topology for the CT scheme. In their tests, the Powell scheme suffers from an artificial growth of the magnetic field. This is explained to be a result of the scheme being non-conservative.

Orszag-Tang Vortex
The Orszag-Tang vortex [141] is a common problem that can be used to test MHD codes for violations of ∇ · B. The relativistic version presented here was performed in 2D using Cartesian coordinates in a 128×128-resolution domain of 2π ×2π length units with periodic boundary conditions, and evolved for 10 time units (c = 1). The equation of state was chosen to be that of an ideal fluid withγ = 4/3 and the initial conditions were set to ρ = 1.0, p = 10.0, v x = −0.99 sin y, v y = 0.99 sin x, B x = − sin y and B y = sin 2x . Snapshots of the evolution are shown in Fig. 9.
As can be seen in Figs. 9 and 10, the general behaviour in both cases is quite similar qualitatively, with only slight differences at specific locations. For instance, when compared to GLM, FCT exhibited higher and sharper maxima for the magnitude of the magnetic field. In a similar fashion, some fine features in the Lorentz factor that can be seen in Fig. 9 for FCT appear to be smeared out when using GLM, giving a false impression of symmetry under 90 • rotations, while the actual symmetry of the problem is under 180 • rotations. This may be an evidence of FCT being less diffusive than GLM.

Spherical accretion
We tested the ability of both methods to preserve a stationary solution by evolving a magnetised Michel accretion in 2D, as shown in Fig. 11. We employed spherical The fluid obeyed an ideal equation of state withγ = 5/3 and the sonic radius was located at r c = 8, and the magnetic field was normalised so that the maximum magnetisation was σ = 10 3 . We repeated the numerical experiment of section Sect. 3.3, now in 2D. As shown in Fig. 11, numerical artefacts start to become noticeable at these later times. For instance, with these extreme magnetisations, for GLM we observe spurious features near the poles at θ = 0 and θ = π, as well as deviations in the velocity field near the outer boundary r = 10 M . The polar region is of special interest for jet simulations, where the divergence-control method must be robust enough to interplay with the axial boundary conditions. The bottom of Fig. 11 shows the profiles of several quantities at θ = π/2. Both divergence-control methods produce an excellent agreement between the solution at different times in the equatorial region. The rightmost column in the bottom of Fig. 11 shows the relative errors in the radial component of the magnetic field for each method. The errors for FCT are not only one order of magnitude lower than for GLM, but also behave differently, remaining at the same level near the more-magnetised inner region instead of growing as seen for GLM.

Accreting Torus
To compare both methods in a setting closer to our intended scientific applications, we simulated accretion from a magnetised perturbed Fishbone-Moncrief torus around a Kerr BH with M = 1 and a = 0.9375. We employed modified spherical MKS coordinates as described in Sect. 2.7.1 and a domain where r ∈ [1.29, 2500] and θ ∈ [0, π] with a resolution of N r ×N θ = 512×256, and evolved the system until t = 2000 M . At the radial boundaries, we imposed noinflow boundary conditions while at the boundaries with the polar axis we imposed symmetric boundary conditions for the scalar variables and the radial vector components and antisymmetric boundary conditions for the azimuthal and polar components. In the BHAC code, noinflow boundary conditions are implemented via continuous extrapolation of the primitive variables and by replacing the three-velocity with zero in case inflowing velocities are present in the solution. The fluid obeyed an ideal equation of state withγ = 4/3. The inner edge of the torus was located at r in = 6.0 and the maximum density was located at r max = 12.0. The initial magnetic field configuration consisted of a single loop with A φ ∝ (ρ/ρ max − ρ cut ) and zero for ρ < ρ cut = 0.2. To simulate vacuum, the region outside the torus was filled with a tenuous atmosphere as is customarily done in these types of simulation. In this case, the prescription for the atmosphere was ρ atm = ρ min r −3/2 and p atm = p min r −5/2 , where ρ min = 10 −5 and p min = 1/3 × 10 −7 . A qualitative difference can be seen even at early times of the simulation. The two upper panels of Fig. 12 show a snapshot of the simulation at t = 20 M using both GLM and FCT. For GLM some of the magnetic field has diffused out of the original torus, magnetising the atmosphere. This artefact is visible for GLM from almost the beginning of the simulation (t ≈ 20 M ), while for FCT it is minimal. Even though this particular artefact is not of crucial importance for the subsequent dynamics of the simulation, this points to a higher inclination of GLM to produce spurious magnetic field structures. At later times (bottom panels of Fig. 12), the most noticeable difference is the smaller amount of turbulent magnetic structures and the bigger plasma magnetisation inside the funnel in FCT, as compared to GLM. This latter difference indicates that the choice of technique to control ∇ · B may have an effect on the possibility of jet formation in GRMHD simulations, although this specific effect was not extensively studied.
To summarise this small section on the comparison between both divergencecontrol techniques, we found from the three tests performed that FCT seems to be less diffusive than GLM, is able to preserve for a longer time a stationary solution, and seems to create less spurious structures in the magnetic field. However, it still has the inconvenient property that it is not possible to implement a cell-entered version of it whilst fully incorporating AMR. As mentioned previously, we are currently working on a staggered implementation adapted to AMR, and this will be described in a separate work. Figure 11 Top: logarithmic density and streamlines in 2D magnetised Michel accretion at times t = 0 M (left) and t = 100 M using GLM (centre) and FCT (right). The horizon is marked by the black line at r = 2. Bottom: profiles at θ = π/2 of, from left to right, radial 3-velocity, density and radial magnetic field at t = 0 M (blue circles) and t = 100 M (red line) using GLM (upper) and FCT (lower). The last column shows the relative difference between the magnetic field at t = 100 M and at the initial condition.

Initial conditions
We consider a hydrodynamic equilibrium torus threaded by a weak magnetic field loop. The particular equilibrium torus solution with constant angular momentum was first presented by [134] and [142] and is now a standard test for GRMHD simulations [see, e.g., 40,50,125,135,143]. To facilitate cross-comparison, we set the initial conditions in the torus close to those adopted by [67,125]. Hence the spacetime is a Kerr BH with dimensionless spin parameter a = 0.9375. The inner radius of the torus is set to r in = 6 M and the density maximum is located at r max = 12 M , where radial and azimuthal positions refer to Boyer-Lindquist coordinates. With these choices, the orbital period of the torus at the density maximum becomes T = 247M . We adopt an ideal gas EOS with an adiabatic index ofγ = 4/3. A weak single magnetic field loop defined by the vector potential is added to the stationary solution. The field strength is set such that 2p max /b 2 max = 100, where global maxima of pressure p max and magnetic field strength b 2 max do not necessarily coincide. In order to excite the MRI inside the torus, the thermal pressure is perturbed by 4% white noise.
As with any fluid code, vacuum regions must be avoided and hence we apply floor values for the rest-mass density (ρ fl = 10 −5 r −3/2 ) and the gas pressure (p fl = 1/3 × 10 −7 r −5/2 ). In practice, for all cells which satisfy ρ ≤ ρ fl we set ρ = ρ fl , in addition if p ≤ p fl , we set p = p fl .
The simulations are performed using horizon penetrating logarithmic KS coordinates (corresponding to our set of modified KS coordinates with h = 0 and R 0 = 0). In the 2D cases, the simulation domain covers r KS ∈ [0.96r h , 2500 M ] and θ ∈ [0, π], where r h 1.35 M . In the 3D cases, we slightly excise the axial region θ ∈ [0.02π, 0.98π] and adopt φ ∈ [0, 2π]. We set the boundary conditions in the horizon and at r = 2500 M to zero gradient in primitive variables. The θ-boundary is handled as follows: when the domain extends all the way to the poles (as in our 2D cases), we adopt "hard" boundary conditions, thus setting the flux through the pole manually to zero. For the excised cone in the 3D cases, we use reflecting "soft" boundary conditions on primitive variables.
The time-update is performed with a two-step predictor corrector based on the TVDLF fluxes and PPM reconstruction. Furthermore, we set the CFL number to 0.4 and use the FCT algorithm. Typically, the runs are stopped after an evolution for t = 5000 M , ensuring that no signal from the outflow boundaries can disturb the inner regions. To check convergence, we adopt the following resolutions: N r × N θ ∈ {256 × 128, 512 × 256, 1024 × 512} in the 2D cases and N r × N θ × N φ ∈ {128 × 64 × 64, 192 × 96 × 96, 256 × 128 × 128, 384 × 192 × 192} in the 3D runs. In the following, the runs are identified via their resolution in θ-direction. For the purpose of validation, we ran the 2D cases also with the HARM3D code [77]. [10] To facilitate a quantitative comparison, we report radial profiles of disk-averaged quantities similar to [67,125,144]. For a quantity q(r, θ, φ, t), the shell average is defined as which is then further averaged over a given time interval to yield q(r) (note that we omit the weighting with the density as done by [67,125]). The limits θ min = π/3, θ max = 2π/3 ensure that atmosphere material is not taken into account in the averaging. The time-evolution is monitored with the accretion rateṀ and the magnetic flux threading the horizon φ Ḃ where both quantities are evaluated at the outer horizon r h . Figure 13 illustrates the qualitative time evolution of the torus by means of the restframe density ρ, plasma-β and the magnetisation σ = b 2 /ρ. After t 300 M , the MRI-driven turbulence leads to accretion onto the central BH. The accretion rate and magnetic flux threading the BH then quickly saturate into a quasi-stationary state (see also Fig. 14). The accreted magnetic flux fills the polar regions and gives rise to a strongly magnetised funnel with densities and pressures dropping to their floor values. For the adopted floor values we hence obtain values of plasma-β as low as 10 −5 and magnetisations peaking at σ ≈ 10 3 in the inner BH magnetosphere.

2D results
These extreme values pose a stringent test for the robustness of the code and, consequently, the funnel region must be handled with the auxiliary inversion based on the entropy switch (see Sect. 2.10.2).

Comparison to HARM3D
For validation purposes we simulated the same initial conditions with the HARM3D code. Wherever possible, we have made identical choices for the algorithm used in both codes, that is: PPM reconstruction, TVDLF Riemann solver and a two step time update. It is important to note that the outer radial boundary differs in both codes: while the HARM3D setup implements outflow boundary conditions at r = 50M , in the BHAC runs the domain and radial grid is doubled in the logarithmic Kerr-Schild coordinates, yielding identical resolution in the region of interest. This ensures that [10] The results were kindly provided by Monika Moscibrodzka. no boundary effects compromise the BHAC simulation. Next to the boundary conditions, also the initial random perturbation varies in both codes which can amount to a slightly different dynamical evolution. After verifying good agreement in the qualitative evolution, we quantify with both codesṀ and φ B according to equations (95) and (96). The results are shown in Fig. 14. Onset-time of accretion, magnitude and overall behaviour are in excellent agreement, despite the chaotic nature of the turbulent flow. We also find the same trend with respect to the resolution-dependence of the results: upon doubling the resolution, the accretion rate Ṁ , averaged over t ∈ [1000, 2000], increases significantly by a factor of 1.908 and 1.843 for BHAC and HARM , respectively. For φ B , the factors are 1.437 and 1.484. At a given resolution, the values for Ṁ and φ B agree between the two codes within their standard deviations. Furthermore, we have verified that these same resolution variations are within the run-to-run deviations due to a different random number seed for the initial perturbation. Further validation is provided in Fig. 15 where disk-averaged profiles for the two highest resolution 2D runs are shown according to equation (94). The quantities of interest are the rest-frame density ρ, the dimensionless temperature Θ := p/ρc 2 , the magnitude of the fluid-frame magnetic field |B| = √ b 2 , thermal and magnetic pressures P gas , P mag and the plasma-β. Again we set the averaging time t ∈ [1000, 2000] M with both codes. The agreement can be considered as very good, that is apart from a slightly higher magnetisation in HARM for r ∈ [20,30], the differences of which are well within the 1σ standard deviation over the averaging time. Small systematic departures at the outer edge of the HARM domain are likely attributable to boundary effects.

3D results
We now turn to the 3D runs performed with BHAC. The qualitative evolution of the high resolution run is illustrated in Fig. 16 showing rest-frame density and b 2 on the β(r) Figure 15 Disk-averaged quantities in the 2D validation runs. The blue curves are obtained with BHAC and the red curves with HARM3D in a two-dimensional setting. The shaded regions mark the 1σ standard deviation of the spatially-averaged snapshots (omitted for the highly fluctuating β ). Apart from a slightly higher magnetisation in HARM for r ∈ [20,30], we find excellent agreement between both codes. two slices z = 0 and y = 0. Overall, the evolution progresses in a similar manner to the 2D cases: MRI-driven accretion starts at t ≈ 300 M and enters saturation at around t 1000 M . Similar values for the magnetisation in the funnel region are also obtained. However, since the MRI cannot be sustained in axisymmetry as poloidal field cannot be re-generated via the ideal MHD induction equation [145], we expect to see qualitative differences between the 2D and 3D cases at late times.
Four different numerical resolutions were run which allows a first convergence analysis of the magnetised torus accretion scenario. Based on the convergence study, we can estimate which numerical resolutions are required for meaningful observational predictions derived from GRMHD simulations of this type.
Since we attempt to solve the set of dissipation-free ideal MHD equations, convergence in the strict sense cannot be achieved in the presence of a turbulent cascade [see also the discussion in 146,147]. [11] Instead, given sufficient scale separation, one might hope to find convergence in quantities of interest like the disk averages and accretion rates. The convergence of various indicators in similar GRMHD torus simulations was addressed for example by [67]. The authors found signs for convergence in most quantifications when adopting a resolution of 192 × 192 × 128, however no convergence was found in the correlation length of the magnetic field. Hence the question as to whether GRMHD torus simulations can be converged with the available computational power is still an open one.
From Figs. 17 and 18, it is clear that the resolution of the N θ = 64 run is insufficient: a peculiar mini-torus is apparent in the disk-averaged density which dimin- [11] Even when the dissipation length is well resolved, high-Reynolds number flows show indications for positive Lyapunov exponents and thus non-convergent chaotic behaviour [see, e.g., 148]. ishes with increasing resolution. Also the onset-time of accretion and the saturation values differ significantly between the N θ = 64 run and its high-resolution counterparts. These differences diminish between the high-resolution runs and we can see signs of convergence in the accretion rate: increasing resolution from N θ = 128 to N θ = 192 appears to not have a strong effect onṀ . Also the evolution of φ B agrees quite well between N θ = 128 and N θ = 192. Hence the systematic resolution dependence ofṀ and φ B in the (even higher resolution) 2D simulations appears to be an artefact of the axisymmetry. It is also noteworthy that the variability amplitude of the accretion rate is reduced in the 3D cases. It appears that the superposition of uncorrelated accretion events distributed over the φ-coordinate tends to smear out the sharp variability that results in the axisymmetric case. Although the simulations were run until t = 5000 M , in order to enable comparison with the 2D simulations, we deliberately set the averaging time to t ∈ [1000 M, 2000 M ]. Figure 18 shows that as the resolution is increased, the diskaveraged 3D data approaches the much higher resolution 2D results shown in Fig. 15, indicating that the dynamics are dominated by the axisymmetric MRI modes at early times. When the resolution is increased from N θ = 128 to N θ = 192, the disk-averaged profiles generally agree within their standard deviations, although we observe a continuing trend towards higher gas pressures and magnetic pressures in the outer regions r ∈ [30M, 50M ]. The overall computational cost quickly becomes significant: for the N θ = 128 simulation we spent 100 K CPU hours on the Iboga cluster equipped with Intel(R) Xeon(R) E5-2640 v4 processors. As the runtime scales with resolution according to N 4 θ , doubling resolution would cost a considerable 1.6 M CPU hours.

Effect of AMR
In order to investigate the effect of the AMR treatment, we have performed a 2D AMR-GRMHD simulation of the torus setup. It is clear that whether a simulation can benefit from adaptive mesh refinement is very much dependent on the physical scenario under investigation. For example, in the hydrodynamic simulations of recoiling BHs due to [37], refinement on the spiral shock was demonstrated to yield significant speedups at a comparable quality of solution. This is understandable as the numerical error is dominated by the shock hypersurface. In the turbulent accretion problem, whether automated mesh refinement yields any benefits is not clear.
The initial conditions for this test are the same as those used in Sect. 4.1. However, due to the limitation of current AMR treatment, we resort to the GLM divergence cleaning method. Three refinement levels are used and refinement is triggered by the error estimator due to [127] with the tolerance set to t = 0.1 (see Sect. 2.11). The numerical resolution in the base level is set to N r × N θ = 128 × 128. To test the validity and efficiency, we also perform the same simulation in a uniform grid with resolution of N r × N θ = 512 × 512 which corresponds to the resolution on the highest AMR level. Figure 19 shows the densities at t = 2000 M as well as the time-averaged density and plasma beta for the AMR and uniform cases. The averaged quantities are calculated in the time interval of t ∈ [1000 M, 2000 M ]. The overall behaviour is quite similar in both cases. Naturally, differences are seen in the turbulent structure in the torus and wind region for a single snapshot. However, in terms of averaged quantities, the difference becomes marginal. In order to better quantify the difference between the AMR and uniform runs, the mass accretion rate and horizon penetrating magnetic flux are shown in Fig. 20. These quantities exhibit a similar behaviour in both cases. In particular, the difference between the AMR run and the uniform run is smaller than the one from different resolution uniform runs and compatible with the run-to-run variation due to a different random number seed (cf. Sect. 4.2). This is unsurprising since the error estimator triggers refinement of the innermost  One of the important merits of using AMR is the possibility to resolve small and large scale dynamics simultaneously with lower computational cost than uniform grids. Figure 21 shows the large scale structure of the averaged magnetisation after 10000 M of simulation time. The averaged quantities are calculated in the time interval t ∈ [6000 M, 10000 M ]. In order to allow the large-scale magnetic field structure to settle down, we average over a later simulation time compared to the previous non-AMR cases. From the figure the collimation angle and magnetisation of the highly magnetised funnel in the AMR case are slightly wider than those in uniform case but the large-scale global structure is very similar in both cases.
A comparison of the computational time for a uniform resolution with 512 2 and the equivalent AMR run (three-level AMR) is shown in Table 5. It is encouraging that even in the naive three-level AMR simulation we obtain qualitatively similar results comparable to the high resolution uniform run, but with having spent only 64% of the computational time of the uniform run. [12] Figure 22 shows the evolution [12] Since we use the same Courant limited timestep for all grid-levels, the speedup is entirely due to saving in computational cells. The additional speedup that would be gained from [149]-type hierarchical timesteps can be estimated from the level population of the simulation: the expected additional gain is only ∼ 8% for this setup.  of the total number of cells during the simulations of AMR cases. Initially less than 2 16 cells are used even when we use three AMR levels, which is a similar number of cells as the uniform grid case with 256 × 256. When the simulation starts, the total cell number increases rapidly due to development of turbulence in the torus which is triggering higher refinement. We note that the total number of cells is still half of the total number of cells in the corresponding high-resolution uniform grid simulation (512 × 512), thus resulting in a direct reduction of computational cost. With increasing dynamic range, we expect the advantages of AMR to increase significantly, rendering it a useful tool for simulations involving structures spanning multiple scales. We leave a more detailed discussion on the effect of the AMR refinement strategy and various divergence-control methods to a future paper.

Radiation post-processing
In order to compute synthetic observable images of the BH shadow and surrounding accretion flow it is necessary to perform general-relativistic ray-tracing and GRRT post-processing [see, e.g., 73,[150][151][152][153][154][155][156]. In this article the GRRT code BHOSS (Black Hole Observations in Stationary Spacetimes) [78] is used to perform these calculations. From BHAC, GRMHD simulation data are produced which are subsequently used as input for BHOSS. Although BHAC has full AMR capabilities, for the GRRT it is most expedient to output GRMHD data that has been re-gridded to a uniform grid.
Since these calculations are performed in post-processing, the effects of radiation forces acting on the plasma during its magnetohydrodynamical evolution are not included. Additionally, the fast-light approximation has also been adopted in this study, i.e., it is assumed that the light-crossing timescale is shorter than the dynamical timescale of the GRMHD simulation and the dynamical evolution of the GRMHD simulation as light rays propagate through it is not considered. Such calculations are considered in an upcoming article [78].
Several different coordinate representations of the Kerr metric are implemented in BHOSS, including Boyer-Lindquist (BL), Logarithmic BL, Cartesian BL, Kerr-Schild (KS), Logarithmic KS, Modified KS and Cartesian KS. All GRMHD simulation data used in this study are specified in Logarithmic KS coordinates. Although BHOSS can switch between all coordinate systems on the fly, it is most straightforward to perform the GRRT calculations in the same coordinate system as the GRMHD data, only adaptively switching to e.g., Cartesian KS when near the polar region. This avoids the need to transform between coordinate systems at every point along every ray in the GRMHD data interpolation, saving computational time.

Radiative transfer equation
Electromagnetic radiation is described by null geodesics of the background spacetime (in this case Kerr), and these are calculated in BHOSS using a Runge-Kutta-Fehlberg integrator with fourth order adaptive step sizing and 5th order error control. Any spacetime metric may be considered in BHOSS, as long as the contravariant or covariant metric tensor components are specified, even if they are only tabulated on a grid. For the calculations presented in this article the Kerr spacetime is written algebraically and in closed-form.
The observer is calculated by constructing a local orthonormal tetrad using trial basis vectors. These basis vectors are then orthonormalized using the metric tensor through a modified Gram-Schmidt procedure. The initial conditions of each ray for the coordinate system under consideration are then calculated and the geodesics are integrated backwards in time from the observer, until they either: (i) escape to infinity (exit the computational domain), (ii) are captured by the BH, or (iii) are effectively absorbed by the accretion flow.
In order to perform these calculations the GRRT equation is integrated in parallel with the geodesic equations of motion of each ray. Written in covariant form, the (unpolarized) GRRT equation, in the absence of scattering, may be written [152] as where I := I ν /ν 3 is the Lorentz-invariant intensity, I ν is the specific intensity, ν is the frequency of radiation, α ν,0 is the specific absorption coefficient and j ν,0 is the specific emission coefficient. The subscript "ν" denotes evaluation of a quantity at a specific frequency, ν, and a subscript "0" denotes evaluation of a quantity in the local fluid rest frame. The terms k µ and u µ are the photon 4-momentum and the fluid 4-velocity of the emitting medium, respectively. The former is calculated from the geodesic integration and the latter is determined from the GRMHD simulation data. The affine parameter is denoted by λ. By introducing the optical depth along the ray together with the Lorentz-invariant emission coefficient η = j ν /ν 2 and Lorentzinvariant absorption coefficient (χ = να ν ), the GRRT equation (97) may be rewritten as Following [152], equation (99) may be reduced to two differential equations where is the relative energy shift between the observer ("obs") and the emitting fluid element. Integrating the GRRT equation in terms of the optical depth in the manner presented provides two major advantages. Firstly, the calculation of the geodesic and of the radiative transfer equation may be performed simultaneously, rather than having to calculate the entire geodesic, store it in memory, and then perform the radiative transfer afterwards. Secondly, by integrating in terms of the optical depth we may specify a threshold value (typically of order unity) whereby the geodesic integration is terminated when encountering optically thick media exceeding this threshold. The combination of these two methods saves significant computational time and expense.

BHOSS-simulated emission from Sgr A*
Having in mind the upcoming radio observations of the BH candidate Sgr A* at the Galactic Centre, the following discussion presents synthetic images of Sgr A*. The GRMHD simulations evolve a single fluid (of ions) and are scale-free in length and mass. Consequently a scaling must be applied before performing GRRT calculations.
Within BHOSS this means specifying the BH mass, which sets the length and time scales, and specifying either the mass accretion rate or an electron density scale, which scales the gas density, temperature and magnetic field strength to that of a radiating electron.
Since the GRMHD simulation is of a single fluid, it is necessary to adopt a prescription for the local electron temperature and rest-mass density. Several such prescriptions exist, some which scale using the mass accretion rate [see, e.g., 71,157,158], scale using density to determine the electron number density and physical accretion rate [see, e.g., 73,159], and some by employing a time-dependent smoothing model of the mass accretion rate [see, e.g., 67].
The dimensionless proton temperature, Θ p , is defined as where k B is the Boltzmann constant, T p is the geometrical (i.e., in physical units) proton temperature and m p is the proton mass. This is then calculated from the GRMHD simulation density (ρ) and pressure (p) as where the fact that the equation of state is ideal and thatγ = 4/3 has been assumed. The magnetic field strength in geometrical units, B geo , is readily obtained from the code magnetic field strength B = b µ b µ as What remains is to specify T e (or Θ e := k B T e /m e c 2 ) and ρ geo . For simplicity we adopt the prescription of [71], wherein T p /T e is assumed to be a fixed ratio. Whilst such an approximation is rather crude, to zeroth order the protons and electrons may be assumed to be coupled in this way. To scale the electron number density we adopt the method of [73], assuming a density scale typically of order 10 7 cm −3 . A somewhat more sophisticated approach is to employ a thresholding of the fluid plasma beta where, when the local plasma beta exceeds some threshold the electrons and protons are coupled as previously mentioned (disk region), but when not exceeded (typically in the funnel region) the electron temperature is assumed to be constant [73,76,157]. Since plasma beta is found to decrease with resolution [67] and in this paper we seek only to demonstrate the convergence of our simulated shadow images obtained from the GRMHD data in regions where the density is non-negligible, we adopt the former model. For the plasma emissivity we use the approximate formula for thermal magnetobremsstrahlung [160], which is given by where e is the electron charge, n e the electron number density, and and ϑ is the pitch angle of the photon with respect to the magnetic field. The absorption coefficient is readily obtained from Kirchoff's law. Each image is generated using a uniform grid of 1000 × 1000 rays, sampling 60 uniformly logarithmically spaced frequency bins between 10 9 Hz and 10 15 Hz. All panels depict the observed image as seen at an observational frequency of 230 GHz, i.e. the frequency at which the EHT will image Sgr A*. This resolution is chosen because the integrated flux over the entire ray-traced image is convergent: doubling the resolution from 500 × 500 to 1000 × 1000 yields an increase of 0.17%, and from 1000 × 1000 to 2000 × 2000 an increase of 0.09%.
In practical GRRT calculations only simulation data which has already reached a quasi-steady state, typically t > 2000 M , is used. In this study we focus on the observational appearance of the accretion flow and BH shadow image. The detailed discussion of the spectrum, variability and plasma models warrants a separate study.

Comparison of images
A natural and important question arises from GRRT calculations of BH shadows: do ray-traced images of GRMHD simulation data converge as the resolution of the GRMHD simulation is increased? The existence of an optimal resolution, beyond which differences in images are small, implies that one can save additional computational time and expense by running the simulation at this optimal resolution. It would also imply that the GRMHD data satisfactorily capture the small-scale structure, turbulence and variations of the accretion flow. As such, it is informative to investigate the convergence of BH shadow images obtained from GRMHD simulation data of differing resolutions, both quantitatively and qualitatively.
To address this question we first generate a series of four snapshot images at t = 2500 M of the the accretion flow and BH shadow from GRMHD simulation data. The resolution of these data are 2N × N × N in (r, θ, φ), i.e., twice as much resolution in the radial direction compared to the zenith and azimuthal directions. The images depicted in Fig. 23 correspond to N = 64, 96, 128 and 192 respectively. Here, the proton to electron temperature ratio was chosen as T p /T e = 3 (similar to [71,157]), the electron number density scaling as 5 × 10 7 cm −3 , the BH mass is set to 4.5 × 10 6 M , the source distance is 8.4 × 10 3 pc, the dimensionless BH spin parameter 0.9375 and the observer inclination angle with respect to the BH spin axis is 60 • .
A direct consequence of increasing the resolution of the GRMHD data is resolving the fine-scale turbulent structure of the accretion flow. The characteristic dark shadow delineating the BH shadow can be clearly seen in all images. As the resolution of the GRMHD data is increased, the images become less diffuse. It is difficult with the naked eye to draw firm physical conclusions, and so in the following we perform a quantitative pixel-by-pixel analysis.
With these snapshot images we may perform a quantitative measure of the difference between any two images through introducing the (normalised) crosscorrelation. For two given two-dimensional arrays f (x, y) and g(x, y) (i.e., 2D images), a measure of similarity or difference may be calculated through the crosscorrelation C, where C ∈ [−1, 1]. The normalised cross-correlation is defined as where µ f , σ f and µ g , σ g correspond to the mean and standard deviation of f and g respectively, and N is equal to the size of either f or g.
In the examples considered here the images are all of equal size and dimension, so N = N f = N g . Equation (109) may be interpreted as the inner product between two data arrays, with the value of C expressing the degree to which the data are aligned with respect to each other. When C = 1 the data are identical, save for a multiplicative constant, when C = 0 the data are completely uncorrelated, and when C < 0 the data are negatively correlated.
Each image pixel has an intensity value represented as a single greyscale value between zero and one. Given the relative intensity data of two different images, Equation (109) is then employed to determine the normalised cross-correlation between the two images. This procedure applied to the panels in Fig. (23) yields the following symmetric matrix of cross-correlation values between the images Indices i and j, where (i, j) = (1, 4), denote the images being cross-correlated. The rightmost column of equation (110) denotes the cross-correlation values, C i,4 , in descending order between images, i.e., the cross-correlation of image 4 with images 1, 2, 3 and 4 respectively. Since C i+1,4 > C i,4 it is clear that the similarity between images increases as the resolution of the GRMHD simulation is increased.
Similarly, for image 3 it is found that C i+1,3 > C i,3 . Finally, it also follows that C 3,4 > C 2,3 > C 1,2 , i.e., the correlation between successive pairs of images increases with increasing resolution, demonstrating the convergence of the GRMHD simulations with increasing grid resolution. Whilst the lowest resolution of 128 × 64 × 64 is certainly insufficient, both difference images and cross-correlation measures indicate that a resolution of 256 × 128 × 128 is sufficient and represents a good compromise. Figure 24 Matrix of image differences D i,j of the four panels in Fig. 23. Upper diagonal panels are greyscale differences. Lower diagonal panels are identical to corresponding upper diagonal panels but with differences illustrated with RGB pixel values. Black panels correspond to D i,i , i.e., trivially the difference between an image and itself.

Conclusions and outlook
We have described the capabilities of BHAC, a new multidimensional generalrelativistic magnetohydrodynamics code developed to perform hydrodynamical and MHD simulations of accretion flows onto compact objects in arbitrary stationary spacetimes exploiting the numerous advantages of AMR techniques. The code has been tested with several one-, two-and three-dimensional scenarios in specialrelativistic and general-relativistic MHD regimes. For validation, GRMHD simulations of MRI unstable tori have been compared with another well-known and tested GRMHD code, the HARM3D code. BHAC shows very good agreement with the HARM3D results, both qualitatively and quantitatively. As a first demonstration of the AMR capabilities in multi-scale simulations, we performed the magnetizedtorus accretion test with and without AMR. Despite the latter intrinsically implies an overhead of ∼ 10%, the AMR runtime amounted to 65% of that relative to the uniform grid, simply as a result of the more economical use of grid cells in the block based AMR. At the same time, the AMR results agree very well with the more expensive uniform-grid results. With increasing dynamic range, we expect the advantages of AMR to increase even more significantly, rendering it a useful tool for simulations involving structures of multiple physical scales.
Currently, two methods controlling the divergence of the magnetic field are available in BHAC and we compared them in three test problems. Although solutions obtained with the cell-centered flux-interpolated constrained transport (FCT) algorithm and the divergence cleaning scheme (GLM) yield the same (correct) physical behaviour in the case of weak magnetic fields, FCT performs considerably better in the presence of strong magnetic fields. In particular, FCT is less diffusive than GLM, is able to preserve a stationary solution, and it creates less spurious structures in the magnetic field. For example, the use of GLM in the case of accretion scenarios with strong magnetic fields leads to worrisome artefacts in the highly magnetised funnel region. The development of a constrained transport scheme compatible with AMR is ongoing and will be presented in a separate work [161].
The EHTC and its European contribution, the BlackHoleCam project [68], aim at obtaining horizon-scale radio images of the BH candidate at the Galactic Center. In anticipation of these results, we have used the 3D GRMHD simulations as input for GRRT calculations with the newly developed BHOSS code [78]. We found that the intensity maps resulting from different resolution GRMHD simulations agree very well, even when comparing snapshot data that was not time averaged. In particular, the normalised cross-correlation between images achieves up to 94.8% similarity between the two highest resolution runs. Furthermore, the agreement between two images converges as the resolution of the GRMHD simulation is increased. Based on this comparison, we find that moderate grid resolutions of 256 × 128 × 128 (corresponding to physical resolutions of ∆r KS × ∆θ KS × ∆φ KS = 0.04M × 0.024rad × 0.05rad at the horizon) yield sufficiently converged intensity maps. Given the large and likely degenerate parameter space and the uncertainty in modelling of the electron distribution, this result is encouraging, as it demonstrates that the predicted synthetic image is quite robust against the ever-present time variability, but also against the impact that the grid resolution of the GRMHD simulations might have. In addition, independent information on the spatial orientation and magnitude of the spin, such as the one that could be deduced from the dynamics of a pulsar near Sgr A* [162], would greatly reduce the space of degenerate solutions and further increase the robustness of the predictions that BHAC will provide in terms of synthetic images.
Finally, we have demonstrated the excellent flexibility of BHAC with a variety of different astrophysical scenarios that are ongoing and will be published shortly. These include: oscillating hydrodynamical equilibrium tori for the modelling of quasi-periodic oscillations [163], episodic jet formation [164] and magnetised tori orbiting non-rotating dilaton BHs [105].
Using again Eq. (114), the source term S := √ γαφg µν ∇ ν n µ can be rewritten as where the first term drops out due to the orthogonality n µ a µ = 0. For a symmetric tensor S µν , we have This follows from the relation Γ 0 ij = −K ij α −1 where Γ 0 ij are elements of the 4-Christoffel symbols [see e.g., (B.9) of 85]. Thus The fact that each of these integrals appear in the evolution equation of two magnetic fluxes guarantees the conservation of divergence, as will be explained in the next Section. On the other hand, the numerical fluxes corresponding to the magnetic field components that are returned by the Riemann solver are surface integrals of the electric field, for example, the flux in the x 2 -direction for B 1 is The innermost integral is the same as that of Eq. (130), so the average flux can be interpreted as ∆S 2F 2 i,j+1/2,k = −∆x iGi,j+1/2,k , whereG i,j+1/2,k is the mean value of the integral from Eq. (130). To second-order accuracy, this integral takes the valueG i,j+1/2,k at the middle of the cell, therefore G i+1/2,j+1/2,k can be found by interpolating the averaged fluxes from the four adjacent cell faces as Since we implemented a cell-centred version of FCT, we are interested in the evolution of the average magnetic field at the cell centres. To second order accuracy, the rate of change of the average value of the x 1 −component of the magnetic field is Now we substitute Eq. (133) into Eq. (128) and Eq. (128) into Eq. (134). After some algebra, we finally obtain eqs. 45 and 46.
Appendix D: Discretisation of ∇ · B and zero-divergence initial conditions CT schemes aim to maintain to zero at machine precision the discretisation of the divergence given by (∇ · B) i,j,k = 1 ∆V i,j,k Φ i+1/2,j,k − Φ i−1/2,j,k + Φ i,j+1/2,k − Φ i,j−1/2,k + Φ i,j,k+1/2 − Φ i,j,k−1/2 , which can be thought of as the volume average of the quantity ∂ a (γ 1/2 B a ) in the given cell. When calculating the evolution equation for (∇ · B) i,j,k , each of the integrals G appear with opposite signs in the expression for dΦ/dt (128) and cancel to machine precision. Therefore, if this discretisation of the divergence was originally zero, it will be zero to machine precision during the rest of the simulation.
However, in the cell-centred version of FCT employed here, we lack information concerning the magnetic flux at cell faces, so Equation (135) cannot be used to monitor the creation of divergence. We will therefore find a derived quantity that we can monitor based on the other available quantities.
Finally, summing over the three directions, we recover Eq. (47). Since the same second-order approximation is used both for the definition and for the time update ofB a , the definition of divergence given by Equation (47) is conserved to machine precision during each evolution step.