General-relativistic Resistive Magnetohydrodynamics with Robust Primitive-variable Recovery for Accretion Disk Simulations

In the context of numerically solving the GRRMHD equations, it is useful to write the equations in the 3+1 form based on the Arnowitt–Deser–Misner (ADM) formalism (see e.g., Rezzolla & Zanotti 2013). We introduce the foliation of space-like hypersurfaces Σ_t, defined as isosurfaces of a scalar time function t, and a time-like unit vector that is normal to these hypersurfaces (Porth et al. 2017),

$$\begin{eqnarray}&&{n}_{\mu }:= -\alpha {{\rm{\nabla }}}_{\mu }t,\end{eqnarray} \tag{ 1 }$$

where α is the lapse function. The frame of the Eulerian observer is defined by the four-velocity n^μ, and the metric associated with each slice Σ_t can be written as

$$\begin{eqnarray}&&{\gamma }_{\mu \nu }:= {g}_{\mu \nu }+{n}_{\mu }{n}_{\nu }.\end{eqnarray} \tag{ 2 }$$

The spatial projection operator is then chosen,

$$\begin{eqnarray}&&{{\gamma }^{\mu }}_{\nu }:= {{\delta }^{\mu }}_{\nu }+{n}^{\mu }{n}_{\nu },\end{eqnarray} \tag{ 3 }$$

thus satisfying the constraint ${{\gamma }^{\mu }}_{\nu }{n}_{\mu }=0$. This can be used to project any four-vector or tensor into its spatial and temporal component. In this formulation, any metric can be written in the form

$$\begin{eqnarray}{g}_{\mu \nu }=\left(\begin{array}{cc}-{\alpha }^{2}+{\beta }_{k}{\beta }^{k} & {\beta }_{j}\\ {\beta }_{j} & {\gamma }_{{ij}}\end{array}\right),\end{eqnarray} \tag{ 4 }$$

where βⁱ is the shift three-vector and γ_ij is the three-metric representing the spatial part of g_μν, with determinant γ, for which (−g)^1/2 $:=$ α γ^1/2. The corresponding inverse metric reads

$$\begin{eqnarray}{g}^{\mu \nu }=\left(\begin{array}{cc}-1/{\alpha }^{2} & {\beta }^{j}/{\alpha }^{2}\\ {\beta }^{j}/{\alpha }^{2} & {\gamma }^{{ij}}-{\beta }^{i}{\beta }^{j}/{\alpha }^{2}\end{array}\right),\end{eqnarray} \tag{ 5 }$$

where γ^ij is the algebraic inverse of γ_ij, and βⁱ = γ^ijβ_j. It is generally straightforward to obtain the expressions of α, βⁱ, and γ_ij from the standard formulation of any general-relativistic metric (see Porth et al. 2017 for commonly used metrics in BHAC). Special relativity is trivially retrieved by setting α = 1, βⁱ = 0, and γ^ij = δ^ij. In the 3+1 formalism, the line element is written as

$$\begin{eqnarray}&&{{ds}}^{2}=-{\alpha }^{2}{{dt}}^{2}+{\gamma }_{{ij}}\left({{dx}}^{i}+{\beta }^{i}{dt}\right)\left({{dx}}^{j}+{\beta }^{j}{dt}\right),\end{eqnarray} \tag{ 6 }$$

describing the motion of coordinate lines as seen by an Eulerian observer

$$\begin{eqnarray}&&{x}_{t+{dt}}^{i}={x}_{t}^{i}-{\beta }^{i}(t,{x}^{j}){dt},\end{eqnarray} \tag{ 7 }$$

moving with four-velocity

$$\begin{eqnarray}&&{n}_{\mu }=(-\alpha ,0,0,0),\ {n}^{\mu }=(1/\alpha ,-{\beta }^{i}/\alpha ).\end{eqnarray} \tag{ 8 }$$

A fluid element with four-velocity u^μ has a Lorentz factor ${\rm{\Gamma }}:= -{u}^{\mu }{n}_{\mu }=\alpha {u}^{0}={(1-{v}^{2})}^{-1/2}$ with ${v}^{2}:= {v}_{i}{v}^{i}$. This defines the fluid three-velocity

$$\begin{eqnarray}&&{v}^{i}:= \displaystyle \frac{{{\gamma }^{i}}_{\mu }{u}^{\mu }}{{\rm{\Gamma }}}=\displaystyle \frac{{u}^{i}}{{\rm{\Gamma }}}+\displaystyle \frac{{\beta }^{i}}{\alpha },\ {v}_{i}:= {\gamma }_{{ij}}{v}^{j}=\displaystyle \frac{{u}_{i}}{{\rm{\Gamma }}}.\end{eqnarray} \tag{ 9 }$$

The fluid equations in general relativity are written as a set of conservation laws for mass,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }\left(\rho {u}^{\mu }\right)=0,\end{eqnarray} \tag{ 10 }$$

and energy and momentum,

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }=0,\end{eqnarray} \tag{ 11 }$$

where ρ is the rest-mass density. The stress–energy tensor for a magnetized perfect fluid is written as (Dionysopoulou et al. 2013; Qian et al. 2017)

$$\begin{eqnarray}&&{T}^{\mu \nu }\equiv {T}_{\mathrm{fluid}}^{\mu \nu }+{T}_{\mathrm{EM}}^{\mu \nu },\end{eqnarray} \tag{ 12 }$$

where the fluid part is expressed independently of the electromagnetic fields (see, e.g., Gammie et al. 2003),

$$\begin{eqnarray}&&{T}_{\mathrm{fluid}}^{\mu \nu }=\left[\rho \left(1+\epsilon \right)+p\right]{u}^{\mu }{u}^{\nu }+{{pg}}^{\mu \nu },\end{eqnarray} \tag{ 13 }$$

with fluid pressure p and specific internal energy . The electromagnetic part is generally given by

$$\begin{eqnarray}&&{T}_{\mathrm{EM}}^{\mu \nu }={F}^{\mu \alpha }{{F}^{\nu }}_{\alpha }-\displaystyle \frac{1}{4}{g}^{\mu \nu }{F}_{\alpha \beta }{F}^{\alpha \beta },\end{eqnarray} \tag{ 14 }$$

where F^μν is the Maxwell tensor with Hodge dual *F^μν, the Faraday tensor.

Applying the 3+1 split and assuming a stationary spacetime, the conservation Equations (10) and (11) can be written in the conservative form

$$\begin{eqnarray}&&{\partial }_{t}\left({\gamma }^{1/2}D\right)+{\partial }_{i}\left[{\gamma }^{1/2}\left(-{\beta }^{i}D+\alpha {v}^{i}D\right)\right]=0,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left({\gamma }^{1/2}{S}_{j}\right)+{\partial }_{i}\left[{\gamma }^{1/2}\left(-{\beta }^{i}{S}_{j}+\alpha {{W}^{i}}_{j}\right)\right]\\ \quad \ =\ {\gamma }^{1/2}\left(\displaystyle \frac{1}{2}\alpha {W}^{{ik}}{\partial }_{j}{\gamma }_{{ik}}+{S}_{i}{\partial }_{j}{\beta }^{i}-U{\partial }_{j}\alpha \right),\end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left({\gamma }^{1/2}\tau \right)+{\partial }_{i}\left[{\gamma }^{1/2}\left(-{\beta }^{i}\tau +\alpha \left({S}^{i}-{v}^{i}D\right)\right)\right]\\ \quad \ =\ {\gamma }^{1/2}\left(\displaystyle \frac{1}{2}{W}^{{ik}}{\beta }^{j}{\partial }_{j}{\gamma }_{{ik}}+{{W}^{j}}_{i}{\partial }_{j}{\beta }^{i}-{S}^{j}{\partial }_{j}\alpha \right),\end{array}\end{eqnarray} \tag{ 17 }$$

where we repeated the equations solved in Porth et al. (2017) for ideal-GRMHD. The purely spatial variant of the stress–energy tensor W^ij reads

$$\begin{eqnarray}\begin{array}{rcl}{W}^{{ij}} & := & {{\gamma }^{i}}_{\mu }{{\gamma }^{j}}_{\nu }{T}^{\mu \nu }\\ & = & \ \rho h{{\rm{\Gamma }}}^{2}{v}^{i}{v}^{j}-{E}^{i}{E}^{j}-{B}^{i}{B}^{j}+\left[p+\displaystyle \frac{1}{2}\left({E}^{2}+{B}^{2}\right)\right]{\gamma }^{{ij}},\end{array}\end{eqnarray} \tag{ 18 }$$

with h = h(ρ, p) as the specific enthalpy of the fluid, and E² $:=$ EⁱE_i, B² $:=$ BⁱB_i, and Eⁱ and Bⁱ as the three-vector spatial parts of the electric and magnetic field in the Eulerian frame, as defined in Equation (27).

Equations (15)–(17) describe the evolution of conserved quantities as measured from an Eulerian reference frame, namely the rest-mass density

$$\begin{eqnarray}&&D:=-\rho {u}^{\mu }{n}_{\mu }=\rho {\rm{\Gamma }},\end{eqnarray} \tag{ 19 }$$

the covariant 3-momentum density

$$\begin{eqnarray}&&{S}_{i}:={{\gamma }^{\mu }}_{i}{n}^{\alpha }{T}_{\alpha \mu }=\rho h{{\rm{\Gamma }}}^{2}{v}_{i}+{\gamma }^{1/2}{\eta }_{{ijk}}{E}^{j}{B}^{k},\end{eqnarray} \tag{ 20 }$$

and the (rescaled) conserved energy density τ $:=$ U − D, where

$$\begin{eqnarray}&&U:={T}^{\mu \nu }{n}_{\mu }{n}_{\nu }=\rho h{{\rm{\Gamma }}}^{2}-p+\displaystyle \frac{1}{2}\left({E}^{2}+{B}^{2}\right).\end{eqnarray} \tag{ 21 }$$

The electric and magnetic fields are evolved through Maxwell's equations, described in Section 2.3. In the absence of gravity, when α = 1, βⁱ = 0, γ^1/2 = 1, and ${\partial }_{t}\gamma =0$, these reduce to the special-relativistic conservation laws. The non-relativistic (Newtonian) limit is obtained by letting ${v}^{2}\ll 1$, p ≪ ρ and E² ≪ B² ≪ ρ, bearing in mind that c = 1.

Again, we emphasize that, unlike in ideal-GRMHD, the electric field in Equations (16) and (17) cannot be substituted as Eⁱ = −γ^−1/2η^ijkB_jv_k (with η_ijk the spatial Levi-Civita antisymmetric symbol). In the resistive-GRMHD limit, Eⁱ has to be obtained from Ampère's law (see next section), resulting in a larger system of equations and therefore implying a more complex solution procedure.

The covariant Maxwell equations in tensorial form are

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\nu }{F}^{\mu \nu }={{ \mathcal J }}^{\mu },\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\nu }{}^{* }{F}^{\mu \nu }=0,\end{eqnarray} \tag{ 23 }$$

where ${{ \mathcal J }}^{\mu }$ is the electric 4-current. Applying the 3+1 split, the tensors in the Maxwell Equations (22) and (23) can be decomposed in terms of the electromagnetic fields, as seen by an observer moving along the normal direction n^ν as

$$\begin{eqnarray}&&{F}^{\mu \nu }={n}^{\mu }{E}^{\nu }-{n}^{\nu }{E}^{\mu }-{\left(-g\right)}^{-1/2}{\eta }^{\mu \nu \lambda \kappa }{n}_{\lambda }{B}_{\kappa },\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{}^{* }{F}^{\mu \nu }={n}^{\mu }{B}^{\nu }-{n}^{\nu }{B}^{\mu }+{\left(-g\right)}^{-1/2}{\eta }^{\mu \nu \lambda \kappa }{n}_{\lambda }{E}_{\kappa },\end{eqnarray} \tag{ 25 }$$

with η^μνλκ the fully antisymmetric symbol (see e.g., Rezzolla & Zanotti 2013) and electric and magnetic field four-vectors,

$$\begin{eqnarray}&&{E}^{\mu }:= {F}^{\mu \nu }{n}_{\nu },\ {B}^{\mu }:= {}^{* }{F}^{\mu \nu }{n}_{\nu },\end{eqnarray} \tag{ 26 }$$

and their three-vector spatial parts in the Eulerian frame,

$$\begin{eqnarray}&&{E}^{i}={F}^{i\nu }{n}_{\nu }=\alpha {F}^{i0},\ {B}^{i}={}^{* }{F}^{i\nu }{n}_{\nu }=\alpha {}^{* }{F}^{i0}.\end{eqnarray} \tag{ 27 }$$

Equation (23) can be written in component form, resulting in Faraday's law,

$$\begin{eqnarray}&&{\partial }_{t}\left({\gamma }^{1/2}{B}^{j}\right)+{\partial }_{i}\left[{\gamma }^{1/2}\left({\beta }^{j}{B}^{i}-{\beta }^{i}{B}^{j}+{\gamma }^{-1/2}{\eta }^{{ijk}}\alpha {E}_{k}\right)\right]=0.\end{eqnarray} \tag{ 28 }$$

The temporal component of Equation (23) leads to the solenoidal constraint

$$\begin{eqnarray}&&{\gamma }^{-1/2}{\partial }_{i}\left({\gamma }^{1/2}{B}^{i}\right)=0.\end{eqnarray} \tag{ 29 }$$

In addition to Faraday's law, Equation (23) can be written in component form, resulting in Ampère's law for the electric field evolution in GRRMHD,

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left({\gamma }^{1/2}{E}^{j}\right)+{\partial }_{i}\left[{\gamma }^{1/2}\left({\beta }^{j}{E}^{i}-{\beta }^{i}{E}^{j}-{\gamma }^{-1/2}{\eta }^{{ijk}}\alpha {B}_{k}\right)\right]\\ \quad \ =\ -{\gamma }^{1/2}\left(\alpha {J}^{j}-q{\beta }^{j}\right).\end{array}\end{eqnarray} \tag{ 30 }$$

The current density ${{ \mathcal J }}^{\mu }$ is decomposed as

$$\begin{eqnarray}&&{{ \mathcal J }}^{\mu }={n}^{\mu }q+{J}^{\mu },\end{eqnarray} \tag{ 31 }$$

where J^μn_μ = 0, $q=-{{ \mathcal J }}^{\mu }{n}_{\mu }$ is the charge density, and J^μ the current density as measured by a Eulerian observer moving with four-velocity n^μ. The temporal component of Equation (23) then provides the charge density q in Equation (31),

$$\begin{eqnarray}&&{\gamma }^{-1/2}{\partial }_{i}\left({\gamma }^{1/2}{E}^{i}\right)=q.\end{eqnarray} \tag{ 32 }$$

The spatial current density is obtained from the resistive Ohm's law in 3+1 split formulation (see, e.g., Palenzuela et al. 2009; Bucciantini & Del Zanna 2013),

$$\begin{eqnarray}&&{J}^{i}={{qv}}^{i}+\displaystyle \frac{{\rm{\Gamma }}}{\eta }\left[{E}^{i}+{\gamma }^{-1/2}{\eta }^{{ijk}}{v}_{j}{B}_{k}-\left({v}_{k}{E}^{k}\right){v}^{i}\right],\end{eqnarray} \tag{ 33 }$$

with the resistivity η (not to be confused with the fully antisymmetric symbol η^μνλκ), as the reciprocal of the electrical conductivity (i.e., η = 1/σ). Substituting Equations (32) and (33) in (30), we obtain the final form of Ampère's law as

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}\left({\gamma }^{1/2}{E}^{j}\right)+{\partial }_{i}\left[{\gamma }^{1/2}\left({\beta }^{j}{E}^{i}-{\beta }^{i}{E}^{j}-{\gamma }^{-1/2}{\eta }^{{ijk}}\alpha {B}_{k}\right)\right]\\ \quad \ =\ -{\gamma }^{1/2}\displaystyle \frac{\alpha {\rm{\Gamma }}}{\eta }\left[{E}^{j}+{\gamma }^{-1/2}{\eta }^{{jik}}{v}_{i}{B}_{k}-\left({v}_{k}{E}^{k}\right){v}^{j}\right]\\ \qquad \ -\ \left(\alpha {v}^{j}-{\beta }^{j}\right){\partial }_{j}\left({\gamma }^{1/2}{E}^{j}\right).\end{array}\end{eqnarray} \tag{ 34 }$$

Note that the resistivity η can depend both on space and time and that this description of the current density is valid in any metric. Hall or dynamo terms can be added by extending Equation (33) to a generalized Ohm's law (e.g., Bucciantini & Del Zanna 2013; Palenzuela 2013).

Finally, it is useful to introduce the fluid-frame (comoving) electric and magnetic field (Bucciantini & Del Zanna 2013),

$$\begin{eqnarray}&&{e}^{\mu }={\rm{\Gamma }}({E}^{i}{v}_{i}){n}^{\mu }+{\rm{\Gamma }}({E}^{\mu }+{\gamma }^{-1/2}{\eta }^{\mu \nu \lambda }{v}_{\nu }{B}_{\lambda }),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{b}^{\mu }={\rm{\Gamma }}({B}^{i}{v}_{i}){n}^{\mu }+{\rm{\Gamma }}({B}^{\mu }-{\gamma }^{-1/2}{\eta }^{\mu \nu \lambda }{v}_{\nu }{E}_{\lambda }),\end{eqnarray} \tag{ 36 }$$

allowing to rewrite the electromagnetic part ${T}_{\mathrm{EM}}^{\mu \nu }$ of Equation (12) as (Qian et al. 2017),

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{EM}}^{\mu \nu } & = & ({b}^{2}+{e}^{2})\left({u}^{\mu }{u}^{\nu }+\displaystyle \frac{1}{2}{g}^{\mu \nu }\right)-{b}^{\mu }{b}^{\nu }-{e}^{\mu }{e}^{\nu }\\ & & -\ {u}_{\lambda }{e}_{\beta }{b}_{\kappa }\left({u}^{\mu }{\gamma }^{-1/2}{\eta }^{\nu \lambda \beta \kappa }+{u}^{\nu }{\gamma }^{-1/2}{\eta }^{\mu \lambda \beta \kappa }\right).\end{array}\end{eqnarray} \tag{ 37 }$$

The comoving electric and magnetic field strength, ${e}^{2}:={e}^{\mu }{e}_{\mu }$, b² $:=$ b^μb_μ, are also employed in the definition of useful dimensionless plasma quantities (e.g., the magnetization σ_mag $:=$ b²/ρ and the gas-to-magnetic pressure ratio, or plasma-β_th $:=$ p_gas/p_mag = 2p/b²).

To adopt a conservative scheme, we write the full system of equations treated in BHAC in the form

$$\begin{eqnarray}&&{\partial }_{t}({\gamma }^{1/2}\,{\boldsymbol{U}})+{\partial }_{i}({\gamma }^{1/2}\,{{\boldsymbol{F}}}^{i})={\gamma }^{1/2}\,{\boldsymbol{S}},\end{eqnarray} \tag{ 38 }$$

where ${\boldsymbol{U}}$ represents conserved variables and ${{\boldsymbol{F}}}^{i}$ are the fluxes,

$$\begin{eqnarray}&&{\boldsymbol{U}}=\left[\begin{array}{c}D\\ {S}_{j}\\ \tau \\ {B}^{j}\\ {E}^{j}\end{array}\right],\ {{\boldsymbol{F}}}^{i}=\left[\begin{array}{c}{{ \mathcal V }}^{i}D\\ \alpha {W}_{j}^{i}-{\beta }^{i}{S}_{j}\\ \alpha ({S}^{i}-{v}^{i}D)-{\beta }^{i}\tau \\ {\beta }^{j}{B}^{i}-{\beta }^{i}{B}^{j}+{\gamma }^{-1/2}{\eta }^{{ijk}}\alpha {E}_{k}\\ {\beta }^{j}{E}^{i}-{\beta }^{i}{E}^{j}-{\gamma }^{-1/2}{\eta }^{{ijk}}\alpha {B}_{k}\end{array}\right],\end{eqnarray} \tag{ 39 }$$

with the transport velocity ${{ \mathcal V }}^{i}:= \alpha {v}^{i}-{\beta }^{i}$. The sources read

$$\begin{eqnarray}&&{\boldsymbol{S}}=\left[\begin{array}{c}0\\ \displaystyle \frac{1}{2}\alpha {W}^{{ik}}{\partial }_{j}{\gamma }_{{ik}}+{S}_{i}{\partial }_{j}{\beta }^{i}-U{\partial }_{j}\alpha \\ \displaystyle \frac{1}{2}{W}^{{ik}}{\beta }^{j}{\partial }_{j}{\gamma }_{{ik}}+{W}_{i}^{j}{\partial }_{j}{\beta }^{i}-{S}^{j}{\partial }_{j}\alpha \\ 0\\ -\alpha {J}^{j}+{\beta }^{j}q\end{array}\right].\end{eqnarray} \tag{ 40 }$$

The form of the GRRMHD equations as evolved in BHAC allows for a temporally and spatially dependent scalar resistivity η(xⁱ, t). In our implementation of GRRMHD the resistivity can depend on any dynamic or static quantity (e.g., rest-mass density, current density, or the position explicitly; see Ripperda et al. 2019 for an application of non-uniform current-dependent resistivity).

The characteristic velocities are required by the Riemann solver and the Courant–Friedrichs–Lewy (CFL) condition that limits the time step. Given the 3+1 structure of the fluxes, we obtain characteristic waves of the form

$$\begin{eqnarray}&&{\lambda }_{\pm }^{i}=\alpha {\lambda }_{\pm }^{{\prime} i}-{\beta }^{i},\end{eqnarray} \tag{ 41 }$$

with ${\lambda }_{\pm }^{{\prime} i}$ the characteristic velocity in the ith direction in the locally flat frame α → 1, β^j → 0 (Anile 1989; Del Zanna et al. 2007). For simplicity we assume the characteristics to be in the limit of maximum diffusivity—that is, the fastest waves locally travel with the speed of light, which after transforming to the Eulerian frame (Pons et al. 1998; White et al. 2016) yields for each component (Del Zanna et al. 2007; Bucciantini & Del Zanna 2013),

$$\begin{eqnarray}&&{\lambda }_{\pm }^{{\prime} i}=\pm \sqrt{{\gamma }^{{ii}}}.\end{eqnarray} \tag{ 42 }$$

Note that a multidimensional Riemann solver for the SRRMHD equations was recently presented by Mignone et al. (2018, 2019) and Miranda-Aranguren et al. (2018).

Our implementation of the GRRMHD equations enforces Equation (29) to roundoff error by means of the staggered CT scheme of Balsara & Spicer (1999), whose implementation in BHAC has been presented in detail by Olivares et al. (2019). The charge density is obtained by numerically taking the divergence of the evolved electric field as in Equation (32).

When the resistivity of the plasma is very small yet finite, the system of Equations (38) becomes stiff. An explicit integration, which is commonly used in ideal-GRMHD codes, then requires time-steps that essentially scale with the resistivity, resulting in prohibitive computational costs. In Komissarov (2007), a Strang-splitting technique is applied in SRRMHD simulations such that the stiff resistive terms can be explicitly computed for the electric field evolution. However, the procedure relies on the assumption that magnetic field and the fluid velocity field remains constant during the (faster) evolution of the resistive electric field Eⁱ. For small values of η, the solution becomes inaccurate or requires extremely small time-steps. An alternative solution is to split off the stiff part of the system of Equations (38) and treat it with an implicit step. With this approach no assumptions are necessary, and in principle all resistivity regimes can be treated without time-step restrictions other than a standard CFL condition. This method was first proposed by Palenzuela et al. (2009) for SRRMHD and later extended by Bucciantini & Del Zanna (2013) and Dionysopoulou et al. (2013) to GRRMHD. Several improvements of the IMEX method for SRRMHD have been proposed recently by Mignone et al. (2019).

Here, we adopt the first–second order IMEX scheme as proposed by Bucciantini & Del Zanna (2013). In particular, we split the system (38) into non-stiff

$$\begin{eqnarray}&&{\partial }_{t}({\gamma }^{1/2}\,{\boldsymbol{X}})={{\boldsymbol{Q}}}_{{\boldsymbol{X}}}\left({\gamma }^{1/2}{\boldsymbol{X}},{\gamma }^{1/2}{\boldsymbol{Y}}\right),\end{eqnarray} \tag{ 43 }$$

and stiff equations

$$\begin{eqnarray}&&{\partial }_{t}({\gamma }^{1/2}\,{\boldsymbol{Y}})={{\boldsymbol{Q}}}_{{\boldsymbol{Y}}}\left({\gamma }^{1/2}{\boldsymbol{X}},{\gamma }^{1/2}{\boldsymbol{Y}}\right)+\displaystyle \frac{1}{\eta }{{\boldsymbol{R}}}_{{\boldsymbol{Y}}}\left({\gamma }^{1/2}{\boldsymbol{X}},{\gamma }^{1/2}{\boldsymbol{Y}}\right),\end{eqnarray} \tag{ 44 }$$

with the conserved quantities ${\boldsymbol{U}}$ split into two subsets $\left\{{\boldsymbol{X}},{\boldsymbol{Y}}\right\}$

$$\begin{eqnarray}&&{\boldsymbol{X}}:= \left[\begin{array}{c}D\\ {S}_{i}\\ \tau \\ {B}^{i}\end{array}\right],\ {\boldsymbol{Y}}:= \left[\begin{array}{c}{E}^{i}\end{array}\right],\end{eqnarray} \tag{ 45 }$$

containing the non-stiff and the stiff variables, respectively.

The time stepping involves a second-order time discretization for the non-stiff variables in ${\boldsymbol{X}}$, which are evolved explicitly as in Porth et al. (2017), and a first-order scheme for the stiff variables in ${\boldsymbol{Y}}$, evolved implicitly. The overall solution step from time level n to n + 1 is written

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\boldsymbol{X}}}}^{(1)} & = & {\tilde{{\boldsymbol{X}}}}^{n}+\displaystyle \frac{{\rm{\Delta }}t}{2}{{\boldsymbol{Q}}}_{{\boldsymbol{X}}}({\tilde{{\boldsymbol{X}}}}^{n},{\tilde{{\boldsymbol{Y}}}}^{n}),\\ {\tilde{{\boldsymbol{Y}}}}^{(1)} & = & {\tilde{{\boldsymbol{Y}}}}^{n}+\displaystyle \frac{{\rm{\Delta }}t}{2}{{\boldsymbol{Q}}}_{{\boldsymbol{Y}}}({\tilde{{\boldsymbol{X}}}}^{n},{\tilde{{\boldsymbol{Y}}}}^{n})+\displaystyle \frac{{\rm{\Delta }}t}{2\eta }{{\boldsymbol{R}}}_{{\boldsymbol{Y}}}({\tilde{{\boldsymbol{X}}}}^{(1)},{\tilde{{\boldsymbol{Y}}}}^{(1)}),\\ {\tilde{{\boldsymbol{X}}}}^{n+1} & = & {\tilde{{\boldsymbol{X}}}}^{n}+{\rm{\Delta }}t{{\boldsymbol{Q}}}_{{\boldsymbol{X}}}({\tilde{{\boldsymbol{X}}}}^{(1)},{\tilde{{\boldsymbol{Y}}}}^{(1)}),\\ {\tilde{{\boldsymbol{Y}}}}^{n+1} & = & {\tilde{{\boldsymbol{Y}}}}^{n}+{\rm{\Delta }}t{{\boldsymbol{Q}}}_{{\boldsymbol{Y}}}({\tilde{{\boldsymbol{X}}}}^{(1)},{\tilde{{\boldsymbol{Y}}}}^{(1)})+\displaystyle \frac{{\rm{\Delta }}t}{\eta }{{\boldsymbol{R}}}_{{\boldsymbol{Y}}}({\tilde{{\boldsymbol{X}}}}^{n+1},{\tilde{{\boldsymbol{Y}}}}^{n+1}),\end{array}\end{eqnarray} \tag{ 46 }$$

where we have incorporated the γ^1/2 factors in $\tilde{{\boldsymbol{X}}}\,:= {\gamma }^{1/2}{\boldsymbol{X}}$, $\tilde{{\boldsymbol{Y}}}:= {\gamma }^{1/2}{\boldsymbol{Y}}$.

The implicit step represented by the ${{\boldsymbol{R}}}_{{\boldsymbol{Y}}}$ terms can be treated analytically due to the linearity (in the electric field) of the resistive Ohm's law (33). For simplicity, but without loss of generality, consider the last electric field update step in the algorithm provided, from intermediate level (1) to the next time level n + 1. Using Equation (33), this can be written explicitly as

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{E}}^{i,n+1} & = & {\tilde{E}}^{i,* }+\displaystyle \frac{\alpha {\rm{\Delta }}t{\rm{\Gamma }}}{\eta }\left[{\tilde{E}}^{i,n+1}+{\gamma }^{-1/2}{\eta }^{{ijk}}{v}_{j}^{n+1}{\tilde{B}}_{k}^{n+1}\right.\\ & & \left.-\ \left({\tilde{E}}^{k,n+1}{v}_{k}^{n+1}\right){v}^{i,n+1}\right],\end{array}\end{eqnarray} \tag{ 47 }$$

where ${\tilde{E}}^{i}:= {\gamma }^{1/2}{E}^{i}$ and ${\tilde{B}}^{i}:= {\gamma }^{1/2}{B}^{i}$. Here, the explicitly updated electric field is ${\tilde{E}}^{i,* }={\tilde{E}}^{i,n}+{\rm{\Delta }}t{{\boldsymbol{Q}}}_{{\boldsymbol{Y}}}({\tilde{{\boldsymbol{X}}}}^{(1)},{\tilde{{\boldsymbol{E}}}}^{(1)})$. Equation (47) only involves local operations (no spatial derivatives needed); hence its inversion is straightforward if the terms on the right-hand side are known, and leads to an explicit expression for the new electric field,

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{E}}^{i,n+1} & = & \displaystyle \frac{{\tilde{E}}^{i,* }}{1+{\sigma }_{{\rm{H}}}{{\rm{\Gamma }}}^{n+1}}-\displaystyle \frac{{{\rm{\Gamma }}}^{n+1}}{1+{\sigma }_{{\rm{H}}}{{\rm{\Gamma }}}^{n+1}}\left[{\gamma }^{-1/2}{\eta }^{{ijk}}{v}_{j}^{n+1}{\tilde{B}}_{k}^{n+1}\right.\\ & & \left.-\ {{\rm{\Gamma }}}^{n+1}\displaystyle \frac{{\tilde{E}}^{k,* }{v}_{k}^{n+1}}{{{\rm{\Gamma }}}^{n+1}+{\sigma }_{{\rm{H}}}}{v}^{i,n+1}\right],\end{array}\end{eqnarray} \tag{ 48 }$$

where σ_H $:=$ αΔt/η. In order to avoid singularities in the ideal-MHD limit $\eta \to 0$, the equation above can be recast as

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{E}}^{i,n+1} & = & \displaystyle \frac{\eta {\tilde{E}}^{i,* }}{\eta +{\sigma }_{{\rm{L}}}{{\rm{\Gamma }}}^{n+1}}-\displaystyle \frac{{\sigma }_{{\rm{L}}}{{\rm{\Gamma }}}^{n+1}}{\eta +{\sigma }_{{\rm{L}}}{{\rm{\Gamma }}}^{n+1}}\left[{\gamma }^{-1/2}{\eta }^{{ijk}}{v}_{j}^{n+1}{\tilde{B}}_{k}^{n+1}\right.\\ & & \left.-\ \eta {{\rm{\Gamma }}}^{n+1}\displaystyle \frac{{\tilde{E}}^{k,* }{v}_{k}^{n+1}}{\eta {{\rm{\Gamma }}}^{n+1}+{\sigma }_{{\rm{L}}}}{v}^{i,n+1}\right],\end{array}\end{eqnarray} \tag{ 49 }$$

where σ_L $:=$ αΔt = η σ_H. Note that Bucciantini & Del Zanna (2013) have a typo in their Equations (33), (35), and (36) for the formulation of the implicit and explicit updates. The analytic inversion is applied in the same way at each substep of the time-stepping algorithm, hence making Equations (48) and (49) completely general by adjusting σ_H and σ_L with the coefficients from a Butcher tableau corresponding to the current substep (Pareschi & Russo 2005; Palenzuela et al. 2009). Therefore, the simple first–second order IMEX algorithm (46) can be extended to arbitrary high-order accuracy while keeping the update equations for Eⁱ unchanged. However, contrary to higher-order schemes, the first–second algorithm (46) naturally includes the ideal-MHD limit, η = 0, without suffering from numerical singularities (Bucciantini & Del Zanna 2013).

Note that the electric field update in Equations (48) and (49) involves the three-velocity vⁱ to be known at the same time level of Eⁱ. However, vⁱ is a primitive quantity that depends nonlinearly on Eⁱ. This dependence makes the update equations intrinsically implicit, and implies that the electric field update in the IMEX algorithm (46) must be carried out concurrently to a conserved-to-primitive-variable inversion. The inversion recovery strategy is of key importance for GR(R)MHD simulations and high sensitivity to the physical parameters makes it a particularly challenging part of the algorithm.

Throughout the solution process of the GRRMHD Equation (38), a transformation of the conserved variables D, S_i, τ into the primitive variables ρ, v_i, and p is necessary. This is a local operation that requires to solve the system of nonlinear Equations (19)–(21). The solution of such a system cannot be written in closed form, requiring a root-finding algorithm that constitutes one of the most expensive and sensitive parts of the whole solution procedure of relativistic MHD codes.⁷

For most of the operations during the GRMHD evolution (i.e., as long as the electric field does not depend on primitive variables), we carry out the conserved-to-primitive inversion by solving the system of equations

$$\begin{eqnarray}\begin{array}{rcl}D & := & \rho {\rm{\Gamma }},\\ {S}_{i} & := & \xi {v}_{i}+{\gamma }^{1/2}{\eta }_{{ijk}}{E}^{j}{B}^{k},\\ \tau & := & \xi -p-D+\displaystyle \frac{1}{2}\left({E}^{2}+{B}^{2}\right),\end{array}\end{eqnarray} \tag{ 50 }$$

where ξ $:=$ ρhΓ². Provided that Eⁱ and Bⁱ are known, such a system can usually be reduced to one single equation (in ξ, p, or other scalar variables; see e.g., Noble et al. 2006 or Siegel et al. 2018) and solved with a one-dimensional (1D) Newton–Raphson (NR) iteration (where 1D refers to the single scalar equation that has to be solved and not to a spatial dimension). This is the standard approach in BHAC for the solution of the ideal-GRMHD equations (van der Holst et al. 2008; Keppens et al. 2012; Porth et al. 2017).

However, the IMEX scheme presented in the previous Section for the GRRMHD equations involves an implicit update of Eⁱ, where both the new electric field and the new three-velocity are unknown. In this case, the system of nonlinear Equations (50) provided cannot be inverted, as Eⁱ is not known a priori but rather an additional variable determined by Equations (48) or (49). Therefore, during the implicit step the conserved to primitive transformation must be carried out concurrently to the implicit electric field update. The system of Equation (50) is thus augmented with Equation (48) or (49) for the electric field, forming again a closed set in the variables ρ, v_i, ξ, and Eⁱ. The new system requires a robust and accurate nonlinear solution method, typically an iterative algorithm. This combined update-transformation step is a crucial operation, which heavily influences the overall performance and accuracy of the IMEX algorithm. If the iteration fails to converge, the electric field cannot be updated and the GRRMHD solution becomes inaccurate. The high failure rate in this step is an issue reported in several relativistic resistive MHD implementations, particularly in low plasma-${\beta }_{\mathrm{th}}\lesssim 0.5$ regimes (Palenzuela et al. 2009; Dionysopoulou et al. 2013; Del Zanna et al. 2016; Qian et al. 2017). A robust inversion method that is reliable in particularly demanding physical regimes (e.g., low-β_th, high-σ_mag) is essential to model accretion flows onto black holes.

Here we present a set of strategies for the inversion-update step. Based on the performance of each strategy, we design a robust approach that yields a minimal amount of failures, allowing for a wide range of simulation parameters that are unattainable with currently available methods.

3.5.1. "1D" Fixed-point Strategy

The most commonly used approach in GRRMHD consists of reducing the system of nonlinear equations to one (Palenzuela et al. 2009; Dionysopoulou et al. 2013), or two (Bucciantini & Del Zanna 2013; Qian et al. 2017) scalar equation(s). A usual choice is to solve the energy equation

$$\begin{eqnarray}&&\xi =p+D+\tau -\displaystyle \frac{1}{2}\left({E}^{2}({v}_{i})+{B}^{2}\right)\end{eqnarray} \tag{ 51 }$$

for the scalar variable ξ—hence the "1D" fixed-point notation or alternatively "2D" fixed-point for two scalar variables (see, e.g., Noble et al. 2006 and Qian et al. 2017 for an iteration on Γ and $W:= (\rho +p\hat{\gamma }/(\hat{\gamma }-1)){{\rm{\Gamma }}}^{2}$). The pressure p(ρ, ξ) is determined by eliminating the dependence on ρ by substitution with ρ = D/Γ, and reduced to a function of ξ only via the relation

$$\begin{eqnarray}&&{\rm{\Gamma }}=\sqrt{1+\displaystyle \frac{S{{\prime} }^{2}}{{\xi }^{2}-S{{\prime} }^{2}}},\end{eqnarray} \tag{ 52 }$$

which follows from Equation (50). Here, $S{{\prime} }^{2}:= S{{\prime} }^{i}{S}_{i}^{{\prime} }$, with ${S}_{i}^{{\prime} }:= {S}_{i}-{\gamma }^{1/2}{\eta }_{{ijk}}{E}^{j}{B}^{k}$. The dependence of Eⁱ on v_i (Equations (48) or (49)), however, is nonlinear and cannot be recast into an explicit relation Eⁱ(ξ). As a consequence, the usual solution approach consists of a hybrid NR iteration, where the dependence of Eⁱ on ξ is not taken into account and the electric field is obtained with a fixed-point iteration. Starting from an initial guess for ξ and Eⁱ, each nonlinear iteration is composed of the following steps:

• 1.  
  Compute the velocity as

  $$\begin{eqnarray}&&{v}_{i}=\displaystyle \frac{{S}_{i}-{\gamma }^{1/2}{\eta }_{{ijk}}{E}^{j}{B}^{k}}{{\xi }^{(m)}}.\end{eqnarray} \tag{ 53 }$$
• 2.  
  Compute Γ and p, and the electric field from Equation (48) or (49).
• 3.  
  Compute the residual,

  $$\begin{eqnarray}&&f(\xi )=\xi -p-D-\tau +\displaystyle \frac{1}{2}({E}^{2}+{B}^{2}),\end{eqnarray} \tag{ 54 }$$

  and its derivative neglecting the dependence of Eⁱ on ξ,

  $$\begin{eqnarray}&&\displaystyle \frac{{df}}{d\xi }=1-\displaystyle \frac{{dp}}{d\xi }.\end{eqnarray} \tag{ 55 }$$
• 4.  
  Update the value of ξ at the mth iteration with a NR step,

  $$\begin{eqnarray}&&{\xi }^{(m+1)}={\xi }^{(m)}-f({\xi }^{(m)}){\left(\displaystyle \frac{{df}({\xi }^{(m)})}{d\xi }\right)}^{-1}.\end{eqnarray} \tag{ 56 }$$
• 5.  
  Track the absolute change in the iteration variables, $| {\xi }^{(m+1)}-{\xi }^{(m)}|$ and $| {E}^{i,(m+1)}-{E}^{i,(m)}|$. The iteration is stopped if this difference falls below a prescribed tolerance, which we normally take to be 10⁻¹⁴.

Step 3 is where a crucial assumption is introduced. Computing the electric field with a fixed-point strategy of this type and neglecting the dependence Eⁱ(ξ) is equivalent to assuming that the electric field only varies slightly between successive Newton steps. This is not always true, especially when the system is very stiff. The stiffness of the nonlinear equations can originate from a parameter choice (e.g., for low values of the resistivity η), or from the physical regime described by the conserved quantities (e.g., large electromagnetic energy density compared to the rest-mass density or pressure resulting in low β_th and high σ_mag). In such cases, the electric field becomes dynamically important, and its variation with respect to other quantities cannot be neglected.

Most implementations of IMEX schemes employing the 1D (or 2D) fixed-point scheme above report numerical issues related to combinations of low-η, high-σ_mag, and low-β_th regimes (Palenzuela et al. 2009; Qian et al. 2017). Failures in the inversion-update step can sometimes be mitigated by reducing the time step (Palenzuela et al. 2009), which effectively reduces the stiffness in the nonlinear system of equations to invert. However, this is an undesirable constraint especially for production runs of accretion flows, where large regions of low-β_th (i.e., the ambient surrounding the accretion disk) or high-σ_mag (i.e., the jet) can rapidly determine a degradation of computational performance.

3.5.2. "3D/4D" Fully Consistent Strategies

3.5.3. Entropy Inversion

In highly magnetized regions of accretion flows (e.g., the relativistic jet in accretion flow simulations), the evolution equation for the conserved energy τ can sometimes become too numerically inaccurate, resulting in unphysical solutions to the conserved to primitive inversion problem. In such situations, BHAC relies on an additional backup strategy for the conserved to primitive inversion, based on the entropy κ. For most of the code operations that involve a standard inversion step (i.e., when Eⁱ is known at the current time step), this entropy "switch" consists of finding the root of a single nonlinear equation, typically in the unknown Γ. The energy τ is discarded in the process, and replaced with a value consistent with the newly recovered primitives. For details on this standard procedure, we refer the reader to the corresponding Section in Porth et al. (2017).

When the entropy switch is needed during the more complicated inversion-update step in the implicit part of our IMEX scheme, the system cannot be reduced to one single equation. We approach the problem by applying the same 3D/4D strategies presented previously, with slight modifications. For the 3D scheme in u_i and the 4D scheme in (z, Eⁱ), we replace any closure relation with an ideal-gas law for the enthalpy (Rezzolla & Zanotti 2013)

$$\begin{eqnarray}&&h=1+\displaystyle \frac{\hat{\gamma }}{\hat{\gamma }-1}\displaystyle \frac{p}{\rho },\end{eqnarray} \tag{ 61 }$$

augmented with the polytropic (isentropic) EOS $p=\kappa {\rho }^{\hat{\gamma }}$. The pressure can still be computed explicitly from the iteration variables, therefore leaving the scheme essentially unchanged. For the 4D scheme in (ξ,  u_i), we replace the equation for ξ such that the residuals read

$$\begin{eqnarray}&&{\boldsymbol{f}}(\xi ,{\boldsymbol{u}})=\left[\begin{array}{c}\xi -\rho h{{\rm{\Gamma }}}^{2}\\ {u}_{i}-{\rm{\Gamma }}\left({S}_{i}-{\gamma }^{1/2}{\eta }_{{ijk}}{E}^{j}{B}^{k}\right)/\xi \end{array}\right],\end{eqnarray} \tag{ 62 }$$

where the enthalpy is still given by the polytropic EOS as previously.

If the entropy-switch strategy is activated, it is applied by default whenever the primary inversion fails. Particularly for accretion flow simulations, the entropy switch is also applied upon successful primary inversion in regions where a low β_th < 10⁻² is detected. Switching between different systems of nonlinear equations, which is required for the entropy-based inversion, is easily handled in BHAC with the NK subroutines, which do not require the full Jacobian for each different inversion strategy. In case of failure of the entropy-based inversion, the last-resort solutions include the replacement of the primitive variables in the faulty cells with averages from nearby converged zones. Alternatively, floor values in pressure, velocity, and rest-mass density can be set if convergence in the inversion step is not reached. Note that the electric field is not floored to an arbitrary value, but rather recalculated with Equation (48) or (49) by using the floor value of v_i.

In order to test the inversion-update strategies presented in Section 3.5, we explore the performance of each algorithm in a large parameter space. We choose to study vastly different GRMHD regimes with σ_mag ∈ [10⁻², 10²] and ${\beta }_{\mathrm{th}}\,\in [{10}^{-10},{10}^{5}]$, by selecting appropriate sets of primitive variables. By picking primitive quantities Γ − 1 ∈ [10⁻², 10³] and η ∈ [10⁻¹⁴, 10⁶], we can construct complete sets of corresponding conserved variables. These are then run through the inversion-update algorithm, whose performance can be evaluated by comparing the retrieved primitives with the manufactured initial sets. Our calculations are performed in flat spacetime, hence without considering the effect of spacetime curvature onto the inversion-update process.

The 3D–u_i, 4D–(ξ, u_i), and 4D–(z, Eⁱ) inversion schemes are applied with and without the entropy-switch strategy discussed in Section 3.5.3. These results are compared to the performance of the 1D-ξ scheme with fixed-point calculation of Eⁱ. Our comparison includes the results of a "backup" strategy designed to yield the maximum rate of successful inversions. This is constructed by applying one of the three multi-D strategies (typically the 3D–u_i) and then switching to another scheme upon failure of the inversion process. As a final backup measure, the entropy-switch strategy in u_i is called in case the three primary strategies fail in retrieving a valid set of primitives (for these tests, we do not apply the β_th threshold described in Section 3.5.3). This backup combined strategy shows dramatic improvements over the standard 1D-ξ scheme, and is the default choice in BHAC for production runs. In the unit tests below, the order of the strategies used in the "backup" approach is 3D–u_i, 4D–(ξ, u_i), 4D–(z, Eⁱ), and finally the entropy switch. In all cases, we use a NR approach with hardcoded Jacobian (we find no significant differences in the convergence rates when applying a NK scheme instead).

As a first test, we explore the (η, σ_mag) parameter space. Considering a wide range of values for the resistivity allows for investigating the importance of the dynamics of Eⁱ compared to the ideal-GRMHD limit for both high and low magnetization. The corresponding sets of primitives are constructed by choosing B² = 1, Γ = 2, and β_th = 0.1 as fixed parameters, as applicable for a relativistic magnetized plasma. For the electric field update, we choose an electric field strength (E*)² = 0.1 (i.e., normally obtained from an explicit update) and a time step Δt = 0.01. The polytropic index is fixed to $\hat{\gamma }=2$.

Figure 1 shows the results of 10⁶ conserved-to-primitive inversions expressed in terms of the number of iterations needed to reach a solution of the nonlinear system with an absolute accuracy of 10⁻¹⁴ in the computed primitives. The maximum amount of iterations allowed is 100, after which the algorithm is stopped and the inversion is considered as having failed (denoted by dark red dots in the plots). Note that, in production runs, additional checks are applied on the iteration error when the maximum iteration number (we typically allow for 100 iterations) is reached. In this case, if the final error is only slightly larger than the prescribed tolerance, the solution can still be considered valid via a larger, user-defined acceptance tolerance. The two 4D strategies (middle columns) show rather complementary regions of failure, with the (ξ, u_i) scheme being more reliable for high-σ_mag zones and the (z, Eⁱ) scheme converging more easily for low-σ_mag zones. The 3D scheme in u_i (leftmost column) shows the highest rate of successful inversions, with no specific regions of failed recovery. The entropy switch applied to the three strategies as backup options (bottom row) shows slightly larger success rate, but does not change the convergence regions qualitatively. All strategies show superior performance compared to the standard 1D strategy in ξ with fixed-point calculation of Eⁱ (top-right panel), both in terms of convergence rate and number of required iterations. The combined "backup" strategy (our default choice for calculations in BHAC) therefore provides a dramatic improvement over the often applied fixed-point strategies, with zero failures in the explored parameter space, compared to a ∼13% failure rate for the 1D scheme in ξ, although the set of equations are admittedly slightly different if the entropy switch is applied.

Zoom In Zoom Out Reset image size

As a second test, we explore the (η, β_th) parameter space. Considering the variation of β_th relates the resistive dynamics of Eⁱ (defined by η) with the case of magnetically dominated (i.e., low-β_th) or thermally dominated energy (i.e., high-β_th) plasma. The primitive sets are constructed upon fixing B² = 1, Γ = 2, σ_mag = 10, (E*)² = 0.1, and ${\rm{\Delta }}t=0.01$. The results are shown in Figure 2, which illustrates how the performance of the three new schemes is similar to the previous case, with almost complementary convergence zones for the 4D schemes, and seemingly scattered failures for the 3D scheme in u_i. For this case, the entropy switch greatly increases the success rate of the inversion procedure. Overall, the new backup strategy combined with entropy switch (bottom-right panel) yields a ∼0.005% failure rate, a major improvement over the standard 1D scheme in ξ (top-right panel). The latter shows a large region of no convergence, with an overall ∼76% failure rate. Failures in the fixed-point strategy seem to be mostly driven by simultaneous conditions of low resistivity and low-β_th, which could preclude modeling large zones of accretion flows which are magnetically dominated and nearly ideal ($\eta \to 0$).

Zoom In Zoom Out Reset image size

As a final test, we consider the (Γ, σ_mag) parameter space. High-σ_mag, high-Γ regions are common in accretion flow simulations (e.g., in the highly magnetized jet emerging from compact objects, where the fluid can be accelerated to high Lorentz factors). Here we manufacture the primitive sets by fixing B² = 1, η = 0.1, β_th = 0.1, (E*)² = 0.1, and ${\rm{\Delta }}t=0.01$. The results in Figure 3 show a generally large region of failure at high-Γ for all strategies. The entropy switch significantly improves the convergence rate for all the strategies. Compared to the standard 1D strategy in ξ (showing a ∼48% failure rate), our combined backup strategy shows failures only in a restricted ∼0.007% of the considered parameter space.

Zoom In Zoom Out Reset image size

In summary, our tests clearly indicate that the new inversion-update strategies presented in Section 3.5 are necessary to properly handle a wide range of regimes typically encountered in GR(R)MHD simulations of accretion flows and of compact-binary mergers. We observe dramatic improvements over the standard approach of a 1D inversion scheme on a scalar with fixed-point calculation of Eⁱ, with recorded failures only in an extremely limited number of cases (e.g., 49 failures over a total of 10⁶ points considered in the [η, β_th] space, or zero failures in the considered [η, σ_mag] space). However, one should be careful when considering the entropy switch as a reliable backup strategy: using the entropy results in a different physical system, where a lower temperature is assumed as applicable in the isentropic limit $d(p/{\rho }^{\hat{\gamma }})/{dt}=0$. In all cases we detect a final error on the computed primitives of the same order of the iteration error (hence below 10⁻¹⁴ in case of convergence).

As a second test actually solving the set of GRRMHD equations, we have considered a one-dimensional shock tube in flat spacetime. Such tests are very restrictive for code validation and show strong nonlinear behavior and steep discontinuities. The ability of the code to handle a range of resistivities for such a problem is essential for astrophysical applications where shocks are ubiquitous. We use the shock-tube setup as proposed by Brio & Wu (1988) to test the code performance, and compare to the results in the ideal-GRMHD limit from Porth et al. (2017). For nonzero resistivity we compare the efficiency and performance of the different inversion methods, including the benchmark Strang-split scheme of Komissarov (2007). Considering the lack of exact solutions for shock tubes with nonzero resistivity, we postpone testing convergence properties to Sections 4.3, 4.4, and 4.5.

The initial conditions are given by

$$\begin{eqnarray}&&\begin{array}{l}(\rho ,p,{v}^{x},{v}^{y},{v}^{z},{B}^{x},{B}^{y},{B}^{z})\\ \ =\ \left\{\begin{array}{cc}(1.0,1.0,0.0,0.0,0.0,0.5,1.0,0.0) & x\lt 0\\ (0.125,0.1,0.0,0.0,0.0,0.5,-1.0,0.0) & x\gt 0\end{array}\right.,\end{array}\end{eqnarray} \tag{ 63 }$$

with an adiabatic index of $\hat{\gamma }=2$. These settings result in β_th = 1.6 for x < 0 and β_th = 0.16 for x > 0. We use a uniform grid with 1024 points spanning $x\in [-1/2,1/2]$. We adopt a second-order TVD limiter (Koren 1993) for spatial reconstruction with CFL number of 0.4.

All tests have been performed with both the Strang-split scheme of Komissarov (2007) and with all inversion-update methods for the IMEX scheme. Note that for these tests we do not activate the entropy fix, nor do we replace faulty cells or apply the floor models since for the multi-D inversion strategies no failures are encountered.

In the right panel of Figure 4, the results for the B^y-component of the magnetic field are shown for all resistivities η ∈ [0, 10⁴] considered, at t = 0.2. The results obtained with the IMEX and Strang schemes cannot be distinguished visually for $\eta \geqslant {10}^{-5}$ cases. For η ≥ 10 the results correspond to the zero conductivity case and for η ≤ 10⁻⁵ no visual differences are observed between resistive and ideal-MHD.

Zoom In Zoom Out Reset image size

In the left panel we compare the runtime for the different primitive-recovery methods, normalized to the runtime of ideal-GRMHD. The Strang-split method of Komissarov (2007) performs best for high resistivity η ≥ 10⁻³, since no additional iterations on the electric field are included in the conserved to primitive transformation. However, for η < 10⁻³ the proportionality of the time step to the resistivity rapidly decreases the performance to prohibitively long runtimes.

For the IMEX scheme with a 1D ξ inversion method we observe a similar trend, and without reducing the CFL condition, no convergence is reached for η < 10⁻³. For example, the 1D ξ inversion scheme needs a CFL condition of 0.06 for η = 10⁻⁴, a CFL condition of 0.006 for η = 10⁻⁵, and 0.0006 for η = 10⁻⁶ (these results are not shown in Figure 4, where we keep the CFL fixed). Palenzuela et al. (2009) reached a similar conclusion, stating that the 1D primitive-variable recovery for more demanding Riemann problems (such as this shock-tube case) lacks robustness for ratios of β_th ≲ 0.4. Note that the shock tube as tested in Dumbser & Zanotti (2009), Bucciantini & Del Zanna (2013), and Qian et al. (2017) is less restrictive for the inversion scheme, due to the smaller discontinuity in plasma-β_th ∈ [0.45, 0.4] between the left and right state in the initial conditions. Our tests confirm that this less demanding setup can be correctly modeled for any η without reducing the time step with all IMEX inversion schemes (including the 1D ξ method) considered in this work.

The 3D and 4D inversion schemes presented in Section 3.5 perform similarly and always produce the correct solution, with a runtime that is comparable to the 1D method. The runtime always remains within a factor ∼2 of the runtime of the ideal-GRMHD solver in BHAC. If using a hardcoded Jacobian for the primitive recovery, such that a NR iterative method can be applied, these schemes become even faster and always produce a correct solution for the shock tube within ∼1.5 times the ideal-GRMHD runtime, for any value of the resistivity.

The third test case is the evolution of a thin current sheet first considered by Komissarov (2007). Once the layer has expanded over several times its initial width, a self-similar evolution ensues. The analytic solution at time t is described by

$$\begin{eqnarray}&&{B}^{y}(x,t)=\mathrm{erf}\left(\displaystyle \frac{x}{2\sqrt{\eta t}}\right),\end{eqnarray} \tag{ 64 }$$

for the magnetic field, while the electric field evolves as

$$\begin{eqnarray}&&{E}^{z}(x,t)=\sqrt{\displaystyle \frac{\eta }{\pi t}}\exp \left(-\displaystyle \frac{{x}^{2}}{4\eta t}\right),\end{eqnarray} \tag{ 65 }$$

and we set t = 1 as initial condition, to start with a resolved state in the self-similar phase. Rest-mass density and pressure are homogeneous and set to ρ = 1 and p = 5000, while all remaining GRRMHD variables are set to zero. The dynamics take place in the x-direction, which is resolved between x ∈ [−1.5, 1.5] by 256 grid points. Here we fix the resistivity to η = 0.01. The test is reproduced with all 3D and 4D inversion methods. Note that for these tests we do not activate the entropy fix, nor do we replace faulty cells or apply the floor models. In all cases we apply an HLL reconstruction scheme with a Koren-type limiter (Koren 1993), and we keep a CFL ratio of 0.5.

The analytic solution for B^y and E^z is shown in Figure 5 at time t = 10 (black line). The numerical results (red line) cannot be distinguished visually. In order to assess the accuracy of the evolution, we study the order of convergence of the numerical solution. We measure the L₁ and ${L}_{\infty }$ norm of the error in the numerical solution by progressively increasing the number of grid points and comparing to a high-resolution run with 8192 grid points. We choose not to compare the numerical results with the analytic solution above, which is only valid in the limit of infinite pressure (as pointed out by Bucciantini & Del Zanna 2013). The error trend thus obtained is reported in Figure 6. We observe that, for low-resolution runs, the accuracy of the scheme is above second order (as expected by the use of a Koren limiter for a smooth solution), a sign that spatial errors dominate over temporal inaccuracies. For high resolutions, where spatial errors become progressively less important, the scheme tends to first-order accuracy, as is expected from the application of the first–second order IMEX scheme from Bucciantini & Del Zanna (2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, in order to test the implementation of the fluxes in 3+1 split formulation, we run the setup under different gauges. These are summarized in Table 1. The results including gauge effects are indistinguishable from Figure 5 once the coordinate-transformations have been accounted for.

Mignone et al. (2019) has recently proposed the first exact two-dimensional equilibrium solution of the SRRMHD equations, which describes a rotating flow with a uniform rest-mass density in a vertical magnetic field and a radial electric field. Adopting a set of cylindrical coordinates (r, ϕ, z), the solution is given by

$$\begin{eqnarray}\begin{array}{rcl}{E}_{r} & = & \displaystyle \frac{{q}_{0}}{2}\displaystyle \frac{r}{{r}^{2}+1},\\ {B}_{z} & = & \displaystyle \frac{\sqrt{{\left({r}^{2}+1\right)}^{2}-{q}_{0}^{2}/4}}{{r}^{2}+1},\\ {v}_{\phi } & = & -\displaystyle \frac{{q}_{0}}{2}\displaystyle \frac{r}{\sqrt{{\left({r}^{2}+1\right)}^{2}-{q}_{0}^{2}/4}},\\ p & = & -\displaystyle \frac{\rho (\hat{\gamma }-1)}{\hat{\gamma }}+\left({p}_{0}+\displaystyle \frac{\rho (\hat{\gamma }-1)}{\hat{\gamma }}\right)\\ & & \times \ {\left(\displaystyle \frac{4{r}^{2}+4-{q}_{0}^{2}}{({r}^{2}+1)(4-{q}_{0}^{2})}\right)}^{\tfrac{\hat{\gamma }}{2(\hat{\gamma }-1)}},\end{array}\end{eqnarray} \tag{ 66 }$$

with radial coordinate r $:=$ x² + y², and on the axis r = 0, we choose the charge density q₀ = 0.7, pressure p₀ = 0.1, and uniform rest-mass density ρ = 1 in the whole domain, in accordance with Mignone et al. (2019). The adiabatic index is set as $\hat{\gamma }=4/3$ and the resistivity as η = 10⁻³. Note that Mignone et al. (2019) evolve the charge density with a separate evolution equation and set it initially as q = q₀/(r² + 1)², whereas in BHAC it is obtained as the divergence of the evolved electric field (see Equation (32)).

The simulation is carried out on a two-dimensional Cartesian grid with x, y ∈ [−10, 10] with a uniform resolution of [N_x × N_y] and (N_x = N_y = 32, 64, 128, 256, 512) cells until time t = 5. We apply continuous extrapolation of all quantities at the boundaries. We use the 3D–u_i primitive-recovery method from Section 3.5 (similar to the inversion method used by Mignone et al. 2019), and we do not activate the entropy fix, nor do we replace faulty cells or apply the floor models since no inversion failures are encountered.

In Figure 7 we show a horizontal cut at y = 0 of the charge density q at t = 0 and t = 2 (left panel), and the second-order convergence of the L₁ and L_∞ norms on the difference in pressure p, between the initial and final time, as a function of resolution (right panel). We point out that our solution is in accordance with the results obtained by the application of the CT method to control the divergence of both the electric and magnetic field variables of Mignone et al. (2019). The large-amplitude oscillations observed by Mignone et al. (2019) in q and E_y and attributed to a general Lagrange-multiplier method for the divergence cleaning are absent in our evolution, where we only control the magnetic field divergence by a CT method instead of both the electric and magnetic field divergences. Second-order convergence is obtained by maintaining the equilibrium solution, showing that spatial errors dominate over temporal errors.

Zoom In Zoom Out Reset image size

As a first test in general relativity, we consider the problem of spherical accretion onto a Schwarzschild black hole with a strong radial magnetic field (Gammie et al. 2003; Villiers & Hawley 2003) and compare to the steady-state solution (Bondi 1952; Michel 1972). This test is a particular challenge for the primitive-recovery method due to the low β_th regions close to the event horizon. We set the mass of the black hole M = 1 and dimensionless spin a = 0, such that distances and times are measured in terms of M. We follow the initial conditions given by Hawley et al. (1984), and we parameterize the field strength through the magnetization σ_mag = 10³ at the inner edge of the domain at r = 1.9 M. To test the inversion method, the resistivity is set to η = 0, while the full set of GRRMHD equations is solved, allowed by the first–second IMEX scheme. We employ two-dimensional modified spherical Kerr-Schild (MKS) coordinates as described in Porth et al. (2017) with r ∈ [1.9 M, 10 M], θ ∈ [0, π] and a uniform radial resolution N_r = 200, and angular resolution N_θ = 100. The steady-state effectively reduces to a one-dimensional problem due to the purely radial dependence of the equilibrium solution. The analytic solution is fixed at the radial boundaries. The test has been reproduced with all 3D and 4D inversion methods. The entropy switch is activated if β_th ≤ 10⁻² (for r ≲ 8 M, see Figure 8) or if the primary inversion procedure fails.

Zoom In Zoom Out Reset image size

Figure 8 shows the radial profiles of rest-mass density ρ, radial three-velocity v^r, β_th, and σ_mag as found with primitive recovery with entropy switch (green dashed line) and without (red dashed line) compared to the analytic solution (black solid line). For the radial three-velocity (top right panel), it is clear that the entropy switch improves the solution close to the event horizon r ≲ 6 M. The improvement provided by the entropy switch is also visible in the L₁ and L_∞ norm of the rest-mass density (Figure 9), increasing the numerical accuracy by approximately a factor four. Second-order convergence is obtained both with and without the entropy switch.

Zoom In Zoom Out Reset image size

Finally, we simulate accretion from a magnetized test-fluid torus (Fishbone & Moncrief 1976) around a Kerr black hole. We again set the mass of the black hole M = 1 and the dimensionless spin a = 0.9375. We employ MKS coordinates on a two-dimensional domain where $r\in [1.29,2500]$ and θ ∈ [0, π] with a uniform resolution of N_r × N_θ = 512 × 256 cells. At the initial state, the inner edge of the torus is located at r = 6 and the maximum rest-mass density is localized at r = 12. To simulate the vacuum region outside the torus we set the rest-mass density and the pressure in the atmosphere as ρ_atm = ρ_minr^−3/2 and p_atm = p_minr^−5/2. The rest-mass density and the pressure are reset whenever they fall below these floor values. The normalization of the power-law floor model is set to ρ_min = 10⁻⁴, p_min = (1/3) × 10⁻⁶.

The initial magnetic field configuration consists of a weak single loop given by the vector potential

$$\begin{eqnarray}&&{A}_{\phi }\propto \max (\rho /{\rho }_{\max }-0.2,0),\end{eqnarray} \tag{ 67 }$$

where ρ_max is the global maximum rest-mass density in the torus. The field strength is determined such that $2{p}_{\max }/{b}_{\max }^{2}=100$, where the spatial locations where p_max and ${b}_{\max }^{2}$ are found do not necessarily coincide. Note that this configuration does not result in an exact MHD equilibrium.

At the polar axis, we impose symmetric boundary conditions for all scalar variables, the radial and poloidal vector components v^r, B^r, v^ϕ, B^ϕ, and the azimuthal component E^θ; antisymmetric boundary conditions are imposed for the azimuthal vector components v^θ and B^θ, and for the radial and poloidal components E^r, E^ϕ (see Porth et al. 2017 and Porth et al. 2019 for a discussion on the boundary conditions in GRMHD simulations of magnetized accretion flows). At the inner and outer radial boundaries, we impose zero-inflow boundary conditions. We choose an ideal-gas EOS with $\hat{\gamma }=4/3$.

The initial equilibrium configuration is perturbed with random, low-amplitude pressure oscillations. This triggers the MRI during the accretion of gas from the torus onto the central object. The instability amplifies the initial magnetic field and drives the disruption of the equilibrium toward a quasi-steady state around t ∼ 500 M. During the process, the resistivity determines the development of diffusive processes. Here, we choose a range of values for $\eta \in [{10}^{-14},{10}^{-2}]$, and we simulate the development of the MRI until $t=2000M$. We also study the exact ideal-MHD limit η = 0 (allowed by the first–second IMEX scheme). Additionally, we perform the same simulation with the ideal-GRMHD version of BHAC. The latter differs from the η = 0 case simulated with our resistive algorithm in many aspects, most notably by solving the induction equation for the magnetic field by assuming that the electric field is a purely dependent quantity Eⁱ = −γ^−1/2η^ijkv_jB_k (as presented in Porth et al. 2017). Comparing the results of the resistive scheme in the ideal-MHD limit with the purely ideal-MHD implementation is therefore a strict and important benchmark. The resistive runs were all completed within a runtime of a factor ∼2–3 longer compared to the ideal-MHD case.

In Figure 10 we show the spatial distribution of the characteristic quantities β_th and σ_mag for progressively decreasing η, and for the simulation run with the ideal-GRMHD version of BHAC. The results are averaged in time between t = 500 M and t = 1000 M to account for statistical fluctuations in the quasi-steady-state accretion. The η = 10⁻² case most evidently shows the diffusive effects of resistivity, quenching the MRI-induced turbulence. Turbulent features become progressively more apparent as η decreases, until the point where numerical resistivity dominates over the explicit physical resistivity. For the resolution considered here, this threshold can be identified around η ∼ 10⁻⁴, at which point the results become visually indistinguishable from simulations at lower resistivity values (including the η = 0 limit). Minor visual differences between the ideal-GRMHD result of BHAC and the η ≤ 10⁻⁴ cases are attributed to the differences in the numerical scheme in the two cases (e.g., the different characteristic speed employed; see Section 3.2).

Zoom In Zoom Out Reset image size

To remove the smoothing introduced by the time averaging, Figure 11 shows a close-up view of the accretion region at time t = 1600 M during the evolution of the system. The η = 10⁻² run (left) shows no sign of turbulence, which is almost completely suppressed by the diffusive processes introduced by the high resistivity. The η = 10⁻¹⁴ case (right), on the contrary, clearly shows the formation of characteristically turbulent structures, with steep gradients both in β_th and σ_mag.

Zoom In Zoom Out Reset image size

Finally, for a more quantitative comparison, we monitor the accretion rate $\dot{M}$ and magnetic flux through the black hole event horizon Φ_B for all runs, defined as

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

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{B}:= \displaystyle \frac{1}{2}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }| {B}^{r}| \sqrt{-g}d\theta d\phi .\end{eqnarray} \tag{ 69 }$$

In Figure 12 the evolution in time of both quantities is shown for the η = 10⁻⁴, 10⁻³, 10⁻² runs, together with the results from the η = 0 run and the ideal-MHD run. The plots show how large resistivity, which is above the numerical resistivity threshold η > 10⁻⁴, can affect the evolution of the system, delaying the instability in time and decreasing the final semi-steady-state values. These results are consistent with the findings by Qian et al. (2017), establishing that a high resistivity significantly quenches the MRI. Additionally, due to the robust conserved to primitive strategies presented in Section 3.5, we find no difficulty in simulating the demanding cases where $\eta \to 0$. As shown in Figure 12, the development time of the MRI and the final steady-state values for both $\dot{M}$ and Φ_B are in good agreement between the η = 0 run and the ideal-GRMHD simulation. These are also consistent with the η = 10⁻⁴ results, confirming that the numerical resistivity is of the order of η < 10⁻³ for the resolution considered here. Identifying this threshold is of major importance, since dissipative length scales that need to be resolved in resistive simulations are proportional to the resistivity. Hence, the necessary resolution depends on the resistivity as $N\propto {\eta }^{-1}$. If the resolution is lower than the necessary threshold to capture the resistive dynamics, the numerical resistivity is prevailing. With explicit resistivity, we can explore new physical regimes that are unattainable in ideal-GRMHD, and explore the effect of dissipative length scales on the development of the MRI. Explicit treatment of viscosity in GRMHD (e.g., Fragile et al. 2018; Fujibayashi et al. 2018) in combination with resistivity will soon allow for the investigation of turbulent black hole accretion without relying on numerical dissipation.

Zoom In Zoom Out Reset image size