A Constrained Transport Method for the Solution of the Resistive Relativistic MHD Equations

TL;DR

This paper introduces a new Godunov-type numerical method for resistive relativistic MHD that ensures divergence-free magnetic fields and charge conservation, improving stability and accuracy over existing methods.

Contribution

The paper presents a novel constrained transport scheme with an efficient implicit-explicit Runge-Kutta method and a five-wave Riemann solver for resistive relativistic MHD.

Findings

Method is more stable and robust than divergence cleaning approaches.

Employs a less diffusive Riemann solver for better accuracy.

Numerical benchmarks demonstrate superior performance.

Abstract

We describe a novel Godunov-type numerical method for solving the equations of resistive relativistic magnetohydrodynamics. In the proposed approach, the spatial components of both magnetic and electric fields are located at zone interfaces and are evolved using the constrained transport formalism. Direct application of Stokes' theorem to Faraday's and Ampere's laws ensures that the resulting discretization is divergence-free for the magnetic field and charge-conserving for the electric field. Hydrodynamic variables retain, instead, the usual zone-centred representation commonly adopted in finite-volume schemes. Temporal discretization is based on Runge-Kutta implicit-explicit (IMEX) schemes in order to resolve the temporal scale disparity introduced by the stiff source term in Ampere's law. The implicit step is accomplished by means of an improved and more efficient Newton-Broyden multidimensional root-finding algorithm. The explicit step relies on a multidimensional Riemann solver to compute the line-averaged electric and magnetic fields at zone edges and it employs a one-dimensional Riemann solver at zone interfaces to update zone-centred hydrodynamic quantities. For the latter, we introduce a five-wave solver based on the frozen limit of the relaxation system whereby the solution to the Riemann problem can be decomposed into an outer Maxwell solver and an inner hydrodynamic solver. A number of numerical benchmarks demonstrate that our method is superior in stability and robustness to the more popular charge-conserving divergence cleaning approach where both primary electric and magnetic fields are zone-centered. In addition, the employment of a less diffusive Riemann solver noticeably improves the accuracy of the computations.

Tables

Table 1. Table 1: Growth rates for the tearing mode instability measured from the simulations. Twelve cases are shown corresponding to three different grid resolutions (left column) and to different combination of divergence control schemes (GLM and CT) and Riemann solvers (LF and MHLLC).

	GLM			CT
Resolution	LF	MHLLC		LF	MHLLC
$512\times 128$	0.49	0.31		0.47	0.32
$1024\times 256$	0.36	0.31		0.36	0.31
$2048\times 512$	0.28	0.30		0.28	0.30

Equations

$${\nabla }_{\mu }(\rho u^{\mu })=0,{\nabla }_{\mu }{T}_{\mathrm{tot}}^{\mu \nu }=0,$$

$${\nabla }_{\mu }{F}^{\mu \nu }=-J^{\nu },{\nabla }_{\mu }{F}^{\star \mu \nu }=0,$$

$${\nabla }_{\mu }{J}^{\mu }=0$$

$${\nabla }_{\mu }{T}_{\mathrm{g}}^{\mu \nu }=-{\nabla }_{\mu }{T}_{\mathrm{EM}}^{\mu \nu }=-J_{\mu }{F}^{\mu \nu },$$

$$T_{\mathrm{g}}^{\mu \nu }=(\epsilon +p)u^{\mu }{u}^{\nu }+p{g}^{\mu \nu },$$

$$F^{\mu \nu }$$

$$F^{\star \mu \nu }$$

$$e^{\mu }=F^{\mu \nu }{u}_{\nu },b^{\mu }=F^{\star \mu \nu }{u}_{\nu }=\frac{1}{2}{ϵ}^{\mu \nu \lambda \kappa }{F}_{\lambda \kappa }{u}_{\nu },$$

$$T_{\mathrm{EM}}^{\mu \nu }=(e^2+b^2)u^{\mu }{u}^{\nu }+\frac{1}{2}(e^2+b^2)g^{\mu \nu }-e^{\mu }{e}^{\nu }-b^{\mu }{b}^{\nu }.$$

$$J^{\mu }=q_0{u}^{\mu }+j^{\mu },$$

$$(\epsilon +p)u^{\nu }{\nabla }_{\nu }{u}_{\mu }+{\nabla }_{\mu }p+u_{\mu }{u}^{\nu }{\nabla }_{\nu }p=q_0{e}_{\mu }+{ϵ}_{\mu \nu \lambda \kappa }{j}^{\nu }{b}^{\lambda }{u}^{\kappa },$$

$$u^{\mu }{\nabla }_{\mu }\epsilon +(\epsilon +p){\nabla }_{\mu }{u}^{\mu }=j_{\mu }{e}^{\mu },$$

$$e^{\mu }=0,$$

$$e^{\mu }=\eta j^{\mu },$$

$$j^{\mu }={\sigma }_E{e}^{\mu }+{\sigma }_B{b}^{\mu },$$

$$u^{\mu }=(\gamma ,\mathbf{u}),$$

$$J^{\mu }=(q,\mathbf{J}),$$

$$e^{\mu }$$

$$b^{\mu }$$

$$\begin{array}[]{ll}&\displaystyle\frac{\partial D}{\partial t}+\mathbf{\nabla}\cdot(D\mathbf{v})=0,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\frac{\partial\mathbf{m}}{\partial t}+\mathbf{\nabla}\cdot(w\mathbf{u}\mathbf{u}+p\mathsf{I}+\mathsf{T})=0,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\frac{\partial\altmathcal{E}}{\partial t}+\mathbf{\nabla}\cdot\mathbf{m}=0,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\frac{\partial\mathbf{B}}{\partial t}+\mathbf{\nabla}\times\mathbf{E}=0,\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr&\displaystyle\frac{\partial\mathbf{E}}{\partial t}-\mathbf{\nabla}\times\mathbf{B}=-\mathbf{J},\end{array}$$

$$\mathcal{E}=w{\gamma }^2-p+\text{altmathcal}{P}_{\mathrm{EM}}.$$

$$\mathsf{T}=-\mathbf{EE}-\mathbf{BB}+\frac{1}{2}(E^2+B^2)\mathsf{I}$$

$$\nabla \cdot \mathbf{B}=0,\nabla \cdot \mathbf{E}=q,$$

$$\frac{\partial q}{\partial t}+\nabla \cdot \mathbf{J}=0,$$

$$\mathbf{J}=\frac{1}{\eta}\Big{[}\gamma\mathbf{E}+\mathbf{u}\times\mathbf{B}-(\mathbf{E}\cdot\mathbf{u})\mathbf{v}\Big{]}+q\mathbf{v},$$

$$\frac{\partial \text{altmathcal}U}{\partial t}=-\nabla \cdot \mathsf{F}(\text{altmathcal}U)+\text{altmathcal}{S}_e+\frac{1}{\eta }\text{altmathcal}S(\text{altmathcal}U)=\text{altmathcal}R(\text{altmathcal}U)+\frac{1}{\eta }\text{altmathcal}S(\text{altmathcal}U),$$

$$\mathsf{F}=\left(\begin{array}[]{c}\rho u_{j}\\ wu_{i}u_{j}+p\delta_{ij}+\mathsf{T}_{ij}\\ m_{j}\\ \varepsilon^{ijk}E_{k}\\ -\varepsilon^{ijk}B_{k}\end{array}\right)^{\intercal},$$

$$\altmathcal{S}_{e}=\left(\begin{array}[]{c}0_{\times 8}\\ q\mathbf{v}\end{array}\right)\,,\quad\frac{1}{\eta}\altmathcal{S}=\left(\begin{array}[]{c}0_{\times 8}\\ \tilde{\mathbf{J}}\end{array}\right).$$

$$\mathbf{B}_{f}\equiv\left(\begin{array}[]{l}B_{{\boldsymbol{x}_{f}}}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{{\boldsymbol{y}_{f}}}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{{\boldsymbol{z}_{f}}}\end{array}\right)=\left(\begin{array}[]{l}B_{x,i+\frac{1}{2},j,k}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{y,i,j+\frac{1}{2},k}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{z,i,j,k+\frac{1}{2}}\end{array}\right)\,,$$

$$\begin{array}[]{l}\Delta_{x}Q_{\boldsymbol{c}}\equiv Q_{\boldsymbol{c}}-Q_{{\boldsymbol{c}}-\hat{\boldsymbol{e}}_{x}}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\Delta_{y}Q_{\boldsymbol{c}}\equiv Q_{\boldsymbol{c}}-Q_{{\boldsymbol{c}}-\hat{\boldsymbol{e}}_{y}}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\Delta_{z}Q_{\boldsymbol{c}}\equiv Q_{\boldsymbol{c}}-Q_{{\boldsymbol{c}}-\hat{\boldsymbol{e}}_{z}}\,,\end{array}$$

$$\begin{array}[]{l}\displaystyle\left<{Q_{\boldsymbol{x}_{f}}}\right>_{x}\equiv\frac{Q_{\boldsymbol{x}_{f}}+Q_{{\boldsymbol{x}_{f}}-\hat{\boldsymbol{e}}_{x}}}{2}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\left<{Q_{\boldsymbol{y}_{f}}}\right>_{y}\equiv\frac{Q_{\boldsymbol{y}_{f}}+Q_{{\boldsymbol{y}_{f}}-\hat{\boldsymbol{e}}_{y}}}{2}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\left<{Q_{\boldsymbol{z}_{f}}}\right>_{z}\equiv\frac{Q_{\boldsymbol{z}_{f}}+Q_{{\boldsymbol{z}_{f}}-\hat{\boldsymbol{e}}_{z}}}{2}\,.\end{array}$$

Full text

††pubyear: 2019††pagerange: A Constrained Transport Method for the Solution of the Resistive Relativistic MHD Equations–B

1 Introduction

The study of the dynamics of relativistic plasmas is of utmost relevance for the interpretation of the phenomenology of high energy astrophysical sources. The ideal Magnetohydrodynamic (MHD) approximation, where dissipative processes are neglected, well describes the large scale dynamics of a plasma and it has been extended to the relativistic regime by Lichnerowicz (1967) and Anile (2005). Over the last decade, the ideal relativistic MHD (RMHD) approach has been used for describing the dynamics of objects like relativistic outflows and jets both from Active Galactic Nuclei (AGN) and gamma ray burst (GRB), accretion flows and pulsar wind nebulae (PWN) by means of numerical simulations (see e.g. McKinney & Blandford, 2009; Mimica et al., 2009; Mignone et al., 2010; McKinney et al., 2012; Mukherjee et al., 2013; Mizuno et al., 2015; Tchekhovskoy & Bromberg, 2016; Olmi et al., 2016; Rossi et al., 2017; Bromberg et al., 2018; Bugli et al., 2018) and important progresses have been reached in establishing robust and accurate numerical schemes for solving the ideal RMHD equations (see e.g. Komissarov, 1999; Balsara, 2001; Del Zanna et al., 2003; Mignone & Bodo, 2006; Gammie et al., 2003; Giacomazzo & Rezzolla, 2006; Del Zanna et al., 2007; Mignone et al., 2009). In typical astrophysical conditions, resistivity is very low and the ideal limit is very well suited for describing processes that occur on dynamical time scales.

However, the flow evolution may lead to the formation of localized region of large gradients (e.g. current sheets) where resistivity cannot be any longer neglected since its role becomes essential in the energy and momentum balance. Processes like magnetic reconnection can be very important for converting magnetic energy in other forms and may play a fundamental role for interpreting the phenomenology of high energy astrophysical objects and their description requires the treatment of resistive effects.

The derivation of a consistent relativistic theory of non-ideal hydrodynamics and MHD has been achieved by several authors (Lichnerowicz, 1967; Israel, 1976; Stewart, 1977; Carter, 1991; Anile, 2005) and in particular the equations of resistive relativistic MHD have been derived in a relatively simple form. More recently several authors have discussed schemes for the numerical solution of such system of equations and have presented actual implementations (see, e.g. Komissarov, 2007; Dumbser & Zanotti, 2009; Palenzuela et al., 2009; Takamoto & Inoue, 2011; Bucciantini & Del Zanna, 2013; Dionysopoulou et al., 2013; Mizuno, 2013; Miranda-Aranguren et al., 2018).

From the numerical point of view, the solution of the resistive RMHD equations is more challenging than approaching their ideal counterpart as it draws on two main issues. First, the resistive RMHD equations are hyperbolic with a stiff relaxation term which accounts for the large difference between the dynamical and the diffusive time scales posing very strict constraints on the time step used for time-explicit calculations. One possibility for overcoming this problem is offered by implicit-explicit Runge-Kutta schemes (IMEX) (Pareschi & Russo, 2005). Such an approach represents a very effective solution to the problem, by combining the simplicity of an explicit treatment of the flux avoiding the time-step restrictions due to the stiffness. IMEX schemes for the resistive RMHD have been implemented and tested by several of the previously mentioned authors (e.g. Palenzuela et al., 2009; Bucciantini & Del Zanna, 2013; Dionysopoulou et al., 2013; Miranda-Aranguren et al., 2018), showing that they are very well suited for such system of equations. The second issue is related to the solenoidal condition and charge conservation. From an analytical point of view, both are direct consequences of Maxwell’s equations; at the discrete level, however, this does no longer hold because of the numerical errors introduced by the underlying algorithm. More specifically, the numerical counterparts of the divergence and the curl operator do not ensure that $\nabla\cdot(\nabla\times\mathbf{A})=0$, where $\mathbf{A}$ is any vector field. As a consequence, for the magnetic field, the constraint $\nabla\cdot\mathbf{B}=0$ may not be maintained during the evolution while, for the electric field, it is the condition $\nabla\cdot\mathbf{E}=q$ that is not respected (when $\mathbf{E}$ and $q$ are evolved through Ampere’s law and the charge conservation equation, respectively). This problem turns out to be particularly severe in shock-capturing schemes, where the stencils used during the reconstruction routines and the related accuracy may vary along the different directions, with the consequence that numerical partial derivatives do not commute (Londrillo & del Zanna, 2004).

Most of the previously mentioned investigators have devised numerical schemes based on a divergence cleaning approach whereby one solves a modified system of conservation laws where Faraday’s and Ampere’s law are coupled to generalized Lagrange multipliers (GLM, Munz et al., 2000; Dedner et al., 2002; Mignone & Tzeferacos, 2010). Divergence errors are convected out of the computational domain at the maximum characteristic speed and damped at the same time. The GLM approach offers ease of implementation since fluid variables and electromagnetic field retain a zone-centered representation. Conversely, in the constrained transport (CT) method originally introduced by Evans & Hawley (1988) (see also Balsara & Spicer, 1999; Londrillo & del Zanna, 2004, in the context of Godunov-type MHD schemes), the magnetic field has a staggered representation whereby the different components live on the face they are normal to. The numerical representation of the divergence and curl operators ensures that the condition $\nabla\cdot(\nabla\times\mathbf{E})=0$ is verified at discrete level, thus maintaining $\nabla\cdot\mathbf{B}=0$ to machine accuracy during the evolution. This approach has been used by Bucciantini & Del Zanna (2013) (see also Del Zanna et al., 2014, for a genuinely third-order scheme), in the context of resistive RMHD, to evolve the magnetic field while still keeping a zone-centered discretization for the electric field.

In the present work, we extend the constrained transport formalism also to the electric field and follow an approach similar to that outlined in Balsara et al. (2016), in the context of the two-fluid equations. There, an alternative staggered collocation for the electric field has been introduced that is both compatible with a Godunov scheme and gives an update of the Ampére’s law consistent with Gauss’s law. In the proposed method, the primary electric and magnetic field variables share the same staggered representation and are thus represented by their surface averages. Hydrodynamic quantities retain instead the usual zone-centered collocation and are interpreted as volume averages. The resistive RMHD equation solver has been implemented as part of the PLUTO code for astrophysical gas dynamics (Mignone et al., 2007) and includes both the standard GLM schemes as well as the newly proposed CT scheme.

The paper is structured as follows. In Section 2 we review the fundamental equations of resistive RMHD, starting from their covariant form. In Section 3 the new constrained transport formulation is presented while our five-wave resistive RMHD Riemann solver is derived in Section 4. Numerical benchmarks are presented in Section 5 and conclusions are finally drawn in Section 6.

2 Equations

In the present section we describe the equations for resistive relativistic MHD, first in general covariant form and later specialized to a Minkowski flat spacetime, separating time and space components and derivatives as needed for practical implementation in a numerical scheme. In the following we will adopt physical units where $c=4\pi=1$, a signature $(-1,+1,+1,+1)$, $g_{\mu\nu}$ will be the metric tensor, and $\nabla_{\mu}$ the covariant derivative associated to the metric. We will use Greek letters for covariant four-dimensional components and Latin letters for spatial three-dimensional ones.

2.1 Covariant formalism

The equations of relativistic MHD are composed by a first couple of conservation laws, one for baryon number (or equivalently mass, assuming a single fluid with particles of given rest mass) and one for total momentum-energy conservation:

[TABLE]

where $\rho$ is the rest mass density, $u^{\mu}$ the fluid four-velocity, and $T^{\mu\nu}_{\mathrm{tot}}$ the total (matter and fields) stress-energy tensor. The second couple is that of Maxwell’s equations

[TABLE]

where $F^{\mu\nu}$ is the Faraday electromagnetic (EM) tensor, $F^{\star\mu\nu}$ its dual, and $J^{\mu}$ the four-current density, which, due to the anti-symmetric property, satisfies the condition

[TABLE]

of electric charge conservation. Let us now split the total stress-energy tensor into the gas and electromagnetic field (EM) components, for which we have

[TABLE]

where the last term is the Lorentz force acting on the charged fluid.

We then decompose our quantities according to $u^{\mu}$. In the case of an ideal fluid, the matter contribution to the energy-momentum tensor is simply provided by ideal hydrodynamics as

[TABLE]

where $\varepsilon=T_{\rm g}^{\mu\nu}u_{\mu}u_{\nu}$ is the gas energy density, and $p$ is the kinetic pressure. The EM tensor and its dual can be expressed as

[TABLE]

where $\epsilon^{\mu\nu\lambda\kappa}$ is the Levi-Civita pseudo-tensor, and where the vectors

[TABLE]

are the electric and magnetic fields measured in the fluid rest frame. Given these definitions, the field component of the stress-energy tensor can be written as

[TABLE]

The four-current can be also split according to

[TABLE]

where $q_{0}=-J^{\mu}u_{\mu}$ is the proper electric charge density and $j^{\mu}$ the conduction current density. We show in Appendix A that the rest-frame charge density can be expressed in terms of the comoving fields and the kinematic vorticity (see Eq. 106). We also point out that even in the ideal limit ($\eta\to 0$) the rest-frame charge does not vanish but it can be expressed as $q_{0}=-b^{\mu}\omega_{\mu}$, where $\omega_{\mu}$ is the kinematic vorticity.

Projecting the momentum-energy conservation law Eq. (4) across the flow yields the equation of motion

[TABLE]

in which the right hand side is the Lorentz force, whereas the energy equation is obtained by projecting along the flow, so that

[TABLE]

where Joule heating acts as an energy source. We clearly need some sort of Ohm’s law to close the system and to specify the heating term.

In the literature three possibilities are most commonly adopted:

ideal plasma – the mobility of charge carriers is so high that the comoving electric field must vanish in order to prevent huge currents. Therefore Ohm’s law is

[TABLE]

and only the source-less Maxwell equation is needed to be solved (for $b^{\mu}$); 2. 2.

resistive plasma – we assume an isotropic tensor of electric conductivity, so that Ohm’s law is simply

[TABLE]

where $\eta$ is the (scalar) resistivity coefficient, and Joule heating retains the usual form $\eta j^{2}$; 3. 3.

dynamo-chiral resistive plasma – in addition to resistive dissipative effects, mean-field dynamo or chiral magnetic effects (CME) may be at work in the plasma, leading to an additional current component along the magnetic field and to magnetic field amplification. In Del Zanna & Bucciantini (2018) the following covariant form of Ohm’s law was proposed

[TABLE]

where subscripts have been added to distinguish the Ohmic and dynamo/chiral effects, coupling to $e^{\mu}$ and $b^{\mu}$, respectively. For $\sigma_{B}=0$ we retrieve the previous case, where $\sigma_{E}=1/\eta$ is the conduction coefficient.

From now on only the second option will be discussed, leaving the implementation of the dynamo/chiral case to a future work.

2.2 The equations for a flat spacetime

The equations are now rewritten for a flat Minkowski spacetime, that is for special resistive RMHD. We are going to separate time and space components and from now on we will use boldface notation and Latin indices for spatial vectors. The fluid four-velocity is now written as

[TABLE]

where $\gamma=\sqrt{1+u^{2}}$ is the Lorentz factor and $\mathbf{v}=\mathbf{u}/\gamma$ is the usual three-velocity. The four-current is now split as

[TABLE]

where $q$ is the charge density (now measured in the laboratory frame), and $\mathbf{J}$ the usual three-current. The EM fields are derived from the Faraday tensor and its dual as

[TABLE]

in which $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic field measured in the laboratory frame.

The set of resistive relativistic equations arising from the time and space split of the covariant Eqs. (1) and (2) are, in vectorial form,

[TABLE]

where $\mathsf{I}$ is the identity matrix and the fluid conserved variables are the density $D=\rho\gamma$ as measured in the laboratory frame, the total momentum density $\mathbf{m}=w\gamma\mathbf{u}+\mathbf{E}\times\mathbf{B}$, and the total energy density

[TABLE]

In the expressions above, $w=\varepsilon+p$ is the specific enthalpy and $\altmathcal{P}_{\rm EM}=(E^{2}+B^{2})/2$ denotes the EM energy density. Finally,

[TABLE]

is the Maxwell’s stress tensor. The remaining Maxwell’s equations give the constraints

[TABLE]

and the charge and current densities are also bound to satisfy the conservation equation

[TABLE]

which directly follows from Eq. (3).

These quantities both enter in Ohm’s law, Eq. (13), here rewritten in terms of spatial vectors alone. From its time component one can derive the $q_{0}\gamma$ term, so that the Ohm’s law for the spatial current becomes

[TABLE]

where $q=\mathbf{\nabla}\cdot\mathbf{E}$ from Gauss’ law, so that the current is determined once the fluid velocity and the electromagnetic fields are known, for a given value of the resistivity $\eta$. In the ideal MHD limit, $\eta\to 0$, the condition $\mathbf{E}+\mathbf{v}\times\mathbf{B}=0$ is retrieved (here $\mathbf{E}$ can be considered a secondary variable with respect to $\mathbf{v}$ and $\mathbf{B}$, to be derived by the above condition), whereas in the resistive limit, $\eta\to\infty$, we have $\mathbf{J}=q\mathbf{v}$, and the charge density $q$ satisfies a continuity equation.

For numerical purposes, we use more compact notations and rewrite the system in quasi-conservative form as

[TABLE]

where $\altmathcal{U}=(D,m_{i},{\cal E},B_{i},E_{i})^{\intercal}$ is the full array of conserved quantities, $\mathsf{F}$ is the flux tensor

[TABLE]

where $\varepsilon^{ijk}$ is the three-dimensional Levi-Civita symbol. The source term is non-zero only in the Ampere’s law and it contains the current density $\mathbf{J}$ (Eq. 24) which we split, for computational purposes, into a stiff ($\tilde{\mathbf{J}}$) and non-stiff ($q\mathbf{v}$) contribution. The source terms $\altmathcal{S}$ and $\altmathcal{S}_{e}$ in Eq. (25) take thus the form

[TABLE]

In Mignone et al. (2018) it has been shown that the system of hyperbolic PDE given by (25) admits 10 propagating modes which are easily recognized in the limits of small or large conductivities. In the $\eta\to\infty$ limit, matter and electromagnetic fields decouple and solution modes approach pairs of light and acoustic waves as well as a number of purely damped (non-propagating) modes. In the $\eta\to 0$ (ideal) limit, modes of propagation coincide with a pair of fast magnetosonic, a pair of slow and Alfvén modes, as expected. The contact mode is always present and it is unaffected by the conductivity.

3 Description of the CT-IMEX Scheme

3.1 Notations and General Formalism

We adopt a Cartesian coordinate system with unit vectors $\hat{\boldsymbol{e}}_{x}=(1,0,0)$, $\hat{\boldsymbol{e}}_{y}=(0,1,0)$ and $\hat{\boldsymbol{e}}_{z}=(0,0,1)$ uniformly discretized into a regular mesh with coordinate spacing $\Delta x$, $\Delta y$ and $\Delta z$. Computational zones are centered at $(x_{i},\,y_{j},\,z_{k})$ and delimited by the six interfaces orthogonal to the coordinate axis aligned, respectively, with $(x_{i\pm\frac{1}{2}},\,y_{j},\,z_{k})$, $(x_{i},\,y_{j\pm\frac{1}{2}},\,z_{k})$ and $(x_{i},\,y_{j},\,z_{k\pm\frac{1}{2}})$.

Primary zone-centered flow variables include density, momentum and energy and are stored by their volume averages inside the zone and labeled as $\altmathcal{U}_{{\boldsymbol{c}}}=(D,\,\mathbf{m},\,\altmathcal{E})_{{\boldsymbol{c}}}$ where the ${\boldsymbol{c}}$ subscript is a shorthand notation for $(i,j,k)$.

Conversely, electric and magnetic fields share the same staggered representation and are understood as surface-averaged quantities located at cell interfaces, as shown in Fig. 1. We denote them as $\mathbf{E}_{f},\mathbf{B}_{f}$ where the ${f}$ subscript tags the different face-centered electric and magnetic field components, e.g.,

[TABLE]

and likewise for the face-centered electric field $\mathbf{E}_{f}$. The subscripts ${\boldsymbol{x}_{f}}$, ${\boldsymbol{y}_{f}}$ and ${\boldsymbol{z}_{f}}$ identify the different component as well as their staggered location inside the control volume, i.e., ${\boldsymbol{x}_{f}}\equiv(x,i+\frac{1}{2},j,k)$, ${\boldsymbol{y}_{f}}\equiv(y,i,j+\frac{1}{2},k)$ and ${\boldsymbol{z}_{f}}\equiv(z,i,j,k+\frac{1}{2})$. Within the CT approach, it also convenient to define the location of cell edges by introducing the “e” subscript, i.e., ${\boldsymbol{x}_{e}}\equiv(x,i,j+\frac{1}{2},k+\frac{1}{2})$, ${\boldsymbol{y}_{e}}\equiv(y,i+\frac{1}{2},j,k+\frac{1}{2})$ and ${\boldsymbol{z}_{e}}\equiv(z,i+\frac{1}{2},j+\frac{1}{2},k)$. Our formulation draws on the two fluid approach of Balsara et al. (2016) and assumes a staggered representation for both the electric and magnetic field. Zone-centered variables are updated using the standard finite volume approach. Staggered quantities are updated using a discrete version of Stokes’ theorem.

Choosing a staggered representation for $\mathbf{E}$ may seem unusual and somehow not consistent with the ideal limit, where $\mathbf{E}=-\mathbf{v}\times\mathbf{B}$ is not a continuous quantity across an interface. However, this choice is consistent with a frozen Riemann solver where (in absence of source terms) the Ampere’s law retains the same form as the induction equation thus leading, for a 1D problem, to the $E_{x}=const$ assumption.

In the following we will make frequent use of the backward difference operators $\Delta_{x}$, $\Delta_{y}$ and $\Delta_{z}$ defined as

[TABLE]

where $Q$ can be any flow quantity. The $\Delta$ operators can be equivalently applied to face-centered or edge-centered values.

In a similar fashion, we also introduce the face-to-center average operators which, to second-order, read

[TABLE]

3.2 IMEX Runge-Kutta Time Stepping

A major challenge when dealing with the numerical solution of Eq. (25) is the presence of a stiff source term in Ampere’s law. Owing to small resistivities typical of astrophysical plasma, this equation may easily become stiff and, for an almost ideal plasma (very small $\eta$), it cannot be solved through an explicit method in a computationally efficient way as the stiff time scale would impose prohibitively small time steps. In order to overcome the time step restriction, we rely on the strong-stability-preserving (SSP) IMplicit-EXplicit (IMEX) Runge-Kutta method, introduced by Pareschi & Russo (2005) and already employed in the context of the resistive RMHD equations by several authors, as cited in the introduction. Applying the IMEX formalism to the original conservation law (25), the resulting time-stepping scheme is explicit in $\altmathcal{R}(\altmathcal{U})$ and implicit in the stiff source term $\altmathcal{S}(\altmathcal{U})/\eta$.

We employ the $2^{\rm nd}$ order IMEX-SSP2(2,2,2) scheme which, when applied to the system (25), consists of the following three stages,

[TABLE]

where $a=1-1/\sqrt{2}$ and the array $\altmathcal{U}=(D,\,\mathbf{m},\,\altmathcal{E},\,\mathbf{B},\,\mathbf{E})$ contains zone-centered as well as face-centered conserved quantities. The arrays $\altmathcal{R}(\altmathcal{U})=-\nabla\cdot\mathsf{F}(\altmathcal{U})+\altmathcal{S}_{e}(\mathbf{\altmathcal{U}})$ and $\altmathcal{S}(\altmathcal{U})$ defined in the previous section embed, respectively, the explicit and implicit contributions. It is important to realize that the first and second stage of the previous time-marching scheme are (locally) implicit as the source term is evaluated at the same intermediate stage as the left hand side. This step involves therefore an implicit update of the electric field which, in our case, is located at zone faces. This is discussed in Section 3.5.

Without loss of generality, we write the single integration stage $s$ of a SSP Runge-Kutta-IMEX scheme by separately working out the zone-centered variables $\altmathcal{U}_{\boldsymbol{c}}=(D,\,\mathbf{m},\,\altmathcal{E})_{\boldsymbol{c}}$ from the staggered fields $\mathbf{B}_{f},\,\mathbf{E}_{f}$ as

[TABLE]

where $\tilde{a}_{sp}$ and $a_{sp}$ are the IMEX coefficients relative to the explicit and implicit temporal discretization (see the Butcher tableau’s available from tables II, III and IV in Pareschi & Russo, 2005), respectively. The tensor $\mathsf{F}$ incorporates the flux components of the hydrodynamical variables only (columns 1-5 of Eq. 26).

At the implementation level, a generic IMEX stage can be decomposed into a sequence of explicit and implicit steps, that is, by first evolving Eqns. (32)–(34) without the last source term in Ampere’s law to some intermediate value and then by solving for the implicit source term alone. The explicit and implicit steps are described in the next two sections.

3.3 Recovery of Primitive Variables

Although the primary set to be evolved in time is that of conserved variables ($\altmathcal{U}_{\boldsymbol{c}}$, $\mathbf{E}_{f}$ and $\mathbf{B}_{f}$), primitive variables defined by $\altmathcal{V}=(\rho,\mathbf{u},p,\mathbf{E},\mathbf{B})$ are required at any stage of the IMEX time stepping. The conversion from primitive to conservative variables poses no difficulty, but the inverse transformation cannot be written in closed analytical form and must be recovered numerically. Here we adopt the approach of Dumbser & Zanotti (2009) in which the problem is reduced to the solution of a quartic function in the Lorentz factor,

[TABLE]

where $C_{1}=(\mathbf{m}-\mathbf{E}\times\mathbf{B})^{2}$ and $C_{2}=(\altmathcal{E}-\altmathcal{P}_{\rm EM})$ are, respectively, the square of the hydrodynamical momentum and the energy. In the expressions above we have assumed an ideal equation of state,

[TABLE]

where $\Gamma_{1}=\Gamma/(\Gamma-1)$ and $\Gamma$ is the (constant) specific heat ratio.

As shown in the appendix of Zenitani et al. (2009), physically consistent solutions correspond to the larger of the two real roots of the previous quartic function. We solve the quartic function using a combination of brackets, bisection and Newton-Raphson method. Once the Lorentz factor is found, the remaining primitive variables can be computed as:

[TABLE]

where the second one is obtained from the energy equation (20) while $\altmathcal{P}_{\rm EM}$ - the EM energy density - has been defined after Eq. (20).

A distinct inversion scheme (used, for instance, by Dionysopoulou et al., 2013) consists of subtracting the electromagnetic contributions from momentum and energy densities and then resorting to a standard relativistic hydro inversion scheme to find the pressure. This can easily achieved using, e.g., the approach outlined by Mignone et al. (2005) which is also valid for different equations of state.

3.4 Explicit Step

We now describe the spatial discretization adopted for the evolution of the zone-centered and staggered variables during the explicit stage of an IMEX stage (Eqns 32-34 without the last term in Ampere’s law).

3.4.1 Explicit Update of Zone-Centered Quantities

In the finite volume approach, Eq. (32) is naturally interpreted as an integral relation relating the change of a volume-averaged conserved quantity to its surface-averaged flux integral across the cell boundary. The discrete form of the divergence operator is thus computed using fluid and electromagnetic quantities available at the $p$-th stage of the time-marching scheme:

[TABLE]

where the $\Delta$’s are the backward difference operators defined in Eq. (29). Here the different $\altmathcal{F}^{*}$ are the hydrodynamic components of the flux computed with a one-dimensional Riemann solver applied at cell interfaces between left and right states at the $(p)$-th stage. To second-order accuracy, a midpoint quadrature rule suffices so that, e.g.,

[TABLE]

where $\altmathcal{V}^{L}_{\boldsymbol{x}_{f}}$ and $\altmathcal{V}^{R}_{\boldsymbol{x}_{f}}$ are the one-sided limit values of the piecewise polynomial reconstruction from within the two neighbor zones adjacent to the interface (higher than second-order finite-volume schemes requires more quadrature points, see, for instance, in Balsara et al., 2016). The reconstruction is carried out on the set of primitive variables $\altmathcal{V}$ (rather than $\altmathcal{U}$) as it is known to produce less oscillatory results. For a second-order reconstruction one has, e.g.,

[TABLE]

where $\delta_{x}\altmathcal{V}_{\boldsymbol{c}}$ are limited slopes in the $x$-direction. A widespread choice, which will also be used by default in the present work, is the van Leer (or harmonic mean) limiter:

[TABLE]

A slightly more compressive option is given by

[TABLE]

where $\rm{mm}(\cdot,\cdot)$ is the minmod limiter. Eq. (42) is known as the monotonized central (MC) limiter.

A popular choice for solving the Riemann problem at cell interfaces in the context of Res-RMHD is the Lax-Friedrichs (LF) scheme, where the fastest characteristic speed is chosen to be the speed of light,

[TABLE]

One-dimensional flux functions in the $y$- and $z$-directions are computed in a similar fashion. Despite its simplicity, the LF scheme (43) can lead to excessive smearing of discontinuous waves as only the two outermost (light) waves are retained in the Riemann fan. For this reason we derive, in section 4, an improved five-wave solver designed to capture the light waves as well as the intermediate acoustic modes, including the contact discontinuity. Our approach is similar to the Harten-Lax-van Leer contact wave (HLLC) Riemann solver recently proposed in the appendix of Miranda-Aranguren et al. (2018).

3.4.2 Explicit Update of Face-Centered Quantities

In the constrained-transport method, Eq. (33) and (34) are understood as surface-averaged relations so that a discrete version of Stokes’ theorem can be applied. This leads to the appearance of edge-centered values (or line-averaged integrals) of the electric and magnetic fields, that is, $\mathbf{E}_{e}\equiv(E_{\boldsymbol{x}_{e}},E_{\boldsymbol{y}_{e}},E_{\boldsymbol{z}_{e}})$ and $\mathbf{B}_{e}\equiv(B_{\boldsymbol{x}_{e}},B_{\boldsymbol{y}_{e}},B_{\boldsymbol{z}_{e}})$. Exact integration over the control volume surfaces yields

[TABLE]

while, for the Ampere’s law we have similarly:

[TABLE]

The $\Delta$ operators are again defined in Eq. (29).

The edge-centered electric and magnetic fields are tagged with a star and are obtained by using a two-dimensional Riemann solver at the $(p)$-th stage of integration. These upwind flux functions result indeed from a four-state representation since two surfaces of discontinuity intersect at a zone edge so that modes of Riemann fans coming from different directions overlap (see, for instance, Del Zanna et al., 2003; Londrillo & del Zanna, 2004; Balsara et al., 2016). Since the Maxwell’s equations are linear and involve propagation of light waves (in the frozen limit approach used here), a proper combination of upwind numerical fluxes along the corresponding orthogonal coordinates leads to the single-valued numerical flux out of these four-states. For the $z$-components one has, for instance,

[TABLE]

and likewise for the $z$-component of the edge-centered magnetic field,

[TABLE]

The four terms with a double superscript in square bracket denote reconstructed values of the $z$-components of the electric and magnetic fields from the face centers (where they are primarily defined) to the zone edge ${\boldsymbol{z}_{e}}$. This is schematically illustrated in Fig. 2 where a top view of the four zones intersecting at an edge (represented by the central point common to all squares) is shown. The first and second superscript refer to the state lying to the left (L) or to the right (R) of the edge with respect to the $x$- and $y$-direction, respectively. For a second-order reconstruction one may use, e.g.,

[TABLE]

where the $\left<{\cdot}\right>_{z}$ is the arithmetic average operator in the $z$-direction, see Eq. (30).

The $x$- and $y$-components of the staggered magnetic field, on the other hand, are continuous at $x$- and $y$-faces and the dissipative terms in Eq. (46) come from two independent wave fans involving jumps in the transverse directions. The same arguments hold for the $x$- and $y$-components of the staggered electric field when evaluating the dissipative terms in Eq. (47). Flux functions located at $x$- and $y$-edges are obtained similarly by cyclic index permutation.

The last term in square brackets in Eq. (45) is a non-stiff source responsible for passive charge advection and it is discretized, for reasons that will be clear in the following section, using an upwind selection rule. We have explored two different options, giving essentially very similar results. In the first one, slightly more diffusive, a simple Rusanov Lax-Friedrichs solver may be employed:

[TABLE]

and similarly for at the $y$- and $z$-interfaces. Alternatively, one could select the upwind direction by looking at the sign of the density flux $F^{*,[\rho]}_{{\boldsymbol{x}_{f}}}$, computed during the evolution of zone-centered variables (Section 3.4.1):

[TABLE]

where $D=\gamma\rho$ is the laboratory density.

The choice of an upwind rule ensures that the charge, computed as the discrete divergence of the electric field, remains oscillation-free (see the discussion at the end of Section 3.5 and 3.6).

3.5 Implicit Step

During the implicit step one has to deal with

[TABLE]

where $\delta t=a_{ss}\Delta t$, $\mathbf{R}^{(s-1)}_{e}$ contains only explicit terms, $\mathbf{B}^{(s)}$ is also computed explicitly, while the electric field $\mathbf{E}$ and four-velocity $\mathbf{u}=\gamma\mathbf{v}$ must be determined. Eq. (51) holds at any location inside the computational zone and has to be solved iteratively by adjusting primitive variables ($\rho$, $\mathbf{v}$ and $p$) while leaving laboratory density, momentum, energy and magnetic field unchanged. In the proposed CT formalism, a major complication now arises since the primary electric field variable is located at zone faces while conserved variables retain a zone-centered representation. An implicit update would have now to be carried at the three different zone faces, thus requiring simultaneous reconstructions from neighboring cells with the unwanted effect of making the algorithm non-local anymore.

This complication could be overcome if one notices that the implicit relation (51) is linear in the electric field and could be easily inverted if the velocity field $\mathbf{u}$ at the next level is known. This suggests that one could first solve the implicit step using the cell-centered electric field and then reconstruct the resulting velocity field at the faces. Then Eq. (51) could be readily solved for the staggered electric field yielding

[TABLE]

where now the velocity has been reconstructed using the values obtained during the cell-center implicit step and we have dropped the superscript $(s)$ for ease of notations. The $(\mathbf{E}\cdot\mathbf{u})$ term (which contains also the transverse components of $\mathbf{E}$) can be rewritten by taking the scalar product of Eq. (51) with $\mathbf{u}$:

[TABLE]

Equation (52) together with (53) now directly gives the staggered components of $\mathbf{E}$ at the chosen interface once the values of velocity and magnetic fields are available at the same location. The interface values of $\mathbf{u}$ and $\mathbf{B}$ can be obtained by interpolating the cell-centered four-velocity and staggered magnetic fields at the same location. However, both $\mathbf{u}$ and the transverse components of $\mathbf{B}$ will, in general, be discontinuous at a given interface.

In order to overcome this problem, we have examined two different options. In the first one, a single-value of velocity and magnetic field can be simply obtained by taking, e.g., the arithmetic average of the left and right states adjoining the interface. Then Eq. (52) can be readily solved for the staggered electric field. In the second option, Eq. (52) could be solved separately for the left and right states thus yielding two different values of the staggered electric field, $\mathbf{E}^{L}_{f}$ and $\mathbf{E}^{R}_{f}$. A single-valued electric field can then be obtained by some averaging procedure and, although in principle several options are possible to produce a single-valued electric field, the rationale for choosing one average versus another is better understood by inspecting at the resulting discretization of the charge (Section 3.6). Here we choose the simple arithmetic average

[TABLE]

which from Eq. (51) allows the (stiff part of the) current to be expressed as

[TABLE]

As it will be later shown in Section 3.6, this choice is equivalent to a conservative update of the charge in which the interface fluxes are computed by means of a local Rusanov Lax-Friedrichs Riemann solver. The amount of dissipation is introduced by the upwind discretization of the non-stiff source term $q\mathbf{v}$ (see the discussion at the end of Section 3.4.2). Although both options have been implemented in the code and no substantial difference has been found, we use the second options throughout this work.

The implicit solver for the cell-centered electric field is the subject of next sub-section.

3.5.1 A Newton-Broyden Scheme for Cell-Center Implicit Update

We propose a novel method to solve Eq. (51) at the cell center based on a multidimensional Newton-Broyden root-finder scheme. Our method exploits conservation of momentum, energy and density through the iterative cycle so that, instead of solving directly Eq. (51), we search for the roots of

[TABLE]

where the specific enthalpy $h=w/\rho$ and the electric field are functions of the four-velocity while $\mathbf{m}\equiv\mathbf{m}^{(s)}$, $D\equiv D^{(s)}$ and $\mathbf{B}\equiv\mathbf{B}^{(s)}$ are known quantities at the beginning of the implicit step. Equation (56) gives a nonlinear system of equations in the four-velocity since $\mathbf{E}=\mathbf{E}(\mathbf{u})$ through Eq. (52) while $h=h(\mathbf{u})$ through the energy equation (20). It simply states that the momentum must not change during the implicit update. We use the four velocity (rather than the electric field) as the independent variable as it offers the advantage that conservative to primitive inversion can be avoided during the cycle.

The Newton-Brodyen method can be sketched as follows:

At the beginning of the step, compute the electric field at zone centers by simple arithmetic average of the face value,

[TABLE]

where $d=x,y,z$. Also, provide a suitable guess of the four-velocity111 We employ the four velocity $\mathbf{u}_{0}$ at the current time level if $\eta/\Delta t>1$; otherwise we use the four-velocity obtained from the ideal RMHD equations. The iteration counter is set to $\kappa=0$; 2. 2.

Using the current value of $\mathbf{u}^{(\kappa)}$, obtain the improved values of the cell-centered electric field using Eq. (52) as well as pressure and enthalpy using Eq. (37); 3. 3.

Compute $\mathbf{f}^{(\kappa)}\equiv\mathbf{f}(\mathbf{u}^{(\kappa)})$:

[TABLE]

Note that the laboratory density, magnetic field and total energy cannot change during this step. 4. 4.

Using the Jacobian, $\mathsf{J}=\partial\mathbf{f}/\partial\mathbf{u}$, obtain an improved guess of the four-velocity through

[TABLE]

where

[TABLE]

and an explicit derivation is presented in Appendix B. 5. 5.

Exit from the iteration cycle if the error is less than some prescribed accuracy $|\mathbf{f}^{(\kappa)}|<\epsilon$, otherwise we go back to step 2 and let $\kappa\to\kappa+1$.

Note that conversion from conservative to primitive variables is not required at each cycle since the primitive variables are automatically updated. For standard applications, our algorithm converges in $2-5$ iterations with a relative tolerance of $10^{-11}$.

3.6 Charge Conservation

In the CT approach, the charge is not an independent variable (as it is the case for the GLM), but it is directly obtained from the divergence of the electric field,

[TABLE]

The charge density is thus collocated at the cell center and it is conserved by construction. In order to see this, let us take the discrete divergence of Ampere’s law (Eq. 34) at any integration stage. Since CT schemes automatically fulfill the condition $\nabla\cdot\nabla\times=0$ at the discrete level, one finds

[TABLE]

The previous relation shows that the charge density obeys a conservative equation by construction.

We now show that when the choices given by Eq. (55) and Eq. (50) are applied to Ampere’s law, one obtains a stable and oscillation-free equivalent discretization for the charge evolution. For simplicity we restrict our attention to the $1^{\rm st}$ order IMEX scheme which, using our standard notations, reads

[TABLE]

By taking the divergence of the last expression and using Eq. (55) and Eq. (50) one finds

[TABLE]

and using Eq. (49), for instance, one finds that

[TABLE]

is a Rusanov Lax-Friedrichs numerical flux.

3.7 Summary of the CT Method

Hereafter we summarize the main steps followed by our IMEX SSP2-CT scheme in order to ease up the implementation of the different computational tasks.

At the beginning of integration, we start with $\altmathcal{U}^{n}_{\boldsymbol{c}}$, $\mathbf{E}^{n}_{f}$ and $\mathbf{B}^{n}_{f}$ as our primary fields. Set $s=1$. 2. 2.

Average the primary staggered fields ($\mathbf{E}$ and $\mathbf{B}$) to zone centers, using Eq. (30). Perform the IMEX implicit update at zone centers to obtain $\mathbf{E}^{(s)}_{\boldsymbol{c}}$ and primitive hydrodynamic variables $\altmathcal{V}^{(s)}$, as described in Section 3.5.1. 3. 3.

Using $\mathbf{E}^{(s)}_{\boldsymbol{c}}$ and $\mathbf{u}^{(s)}_{\boldsymbol{c}}$, achieve the implicit step on the staggered electric field as well, as detailed in Section 3.5. This yields $\mathbf{E}^{(s)}_{f}$. 4. 4.

Compute interface Godunov fluxes and the multidimensional line-averaged electric and magnetic fields to obtain the divergence and curl operators contained in the predictor step ($s=2$) of the IMEX scheme, that is, Eqns (32, 33, 34) for $s=2$. Also, compute the explicit source terms needed for this stage and add all terms up to obtain the explicit contributions needed for $\altmathcal{U}^{(2)}_{\boldsymbol{c}}$, $\mathbf{B}^{(2)}_{f}$ and $\mathbf{E}^{(2)}_{f}$. 5. 5.

Repeat the implicit steps 2 and 3 with $s=2$ to achieve $\altmathcal{U}^{(2)}_{\boldsymbol{c}}$ $\mathbf{E}^{(2)}_{f}$. 6. 6.

Perform the final corrector stage as in step 4 and obtain the solution at the next time level $n+1$.

4 The MHLLC Riemann Solver

Within the IMEX formalism, the solution to the Riemann problem can be obtained under the frozen limit condition (infinite resistivity) by ignoring the effect of the stiff source term in Ampere’s law on the characteristic wave propagation. In this limit, the current-density is $\mathbf{J}=q\mathbf{v}$ and the momentum and energy equations in (19) can be rewritten as

[TABLE]

Thus the coupling between EM and hydrodynamic fields results only from a non-zero charge density. At the linear level Mignone et al. (2018) have shown that, in this regime, the characteristic structure of the resistive RMHD equations entails a pair of acoustic modes and a pair of light waves. The contribution of the charge, being a second-order (nonlinear) term, is neglected in the dispersion relation. In addition, assuming that the presence of the source term on the right hand side of Eq. (66) does not change the Rankine-Hugoniot conditions (a similar inference has been done by Miranda-Aranguren et al., 2018), we can consider Maxwell’s and hydrodynamic equations to be decoupled during the solution of the Riemann problem.

Our approach therefore is based on the direct combination of two Riemann solvers: one for the outermost EM waves across which only transverse components of electric and magnetic fields can change, and a second one across the sound waves where only hydrodynamical variables have non-trivial jumps. We shall refer to the former and the latter as the outer and the inner Riemann solver, respectively. Across the outermost EM waves, jump conditions for $\mathbf{E}$ and $\mathbf{B}$ follow directly from Maxwell’s equations. The presence of the outermost waves modifies total momentum and energy so that conservative variables behind the discontinuities serve as input left and right states to the inner Riemann problem for which any any relativistic hydrodynamic solver may be employed. Here we choose to rely on relativistic HLLC scheme of Mignone & Bodo (2005).

Across each wave $\lambda$ the jump conditions

[TABLE]

must be satisfied for any pair of states $\altmathcal{U}_{r}$ and $\altmathcal{U}_{l}$ and corresponding flux functions $F_{l}$ and $F_{r}$. The Riemann fan comprises five waves and its structure is schematically depicted in Fig. 3. For a discontinuity propagating in the $n$ ($=1,2,3$) direction, our state and flux functions are more conveniently written by separating hydrodynamical from EM terms as

[TABLE]

where $\mathsf{T}$ is the Maxwell stress tensor (Eq. 21) while

[TABLE]

are the hydrodynamic momentum, energy and the EM Poynting vector, respectively.

Owing to its structure, we label our Riemann solver as “MHLLC”, where “M” denotes the outer Maxwell solver while “HLLC” stands for the inner the Harten-Lax-van Leer-Contact approximate Riemann solver originally introduced by Toro et al. (1994) in the context of the classical equations of gas-dynamics. The HLLC scheme improves over the traditional HLL method (Harten et al., 1983; Einfeldt et al., 1991) by restoring the contact wave in the solution. We point out that, although devised from different assumptions, our method of solution ends up being similar to the version of the HLLC solver reported in the appendix of Miranda-Aranguren et al. (2018). The outer and inner solvers are described in the next two sections.

4.1 Outer Solver

The jump conditions (67) must hold across the leftmost wave ($\lambda_{-}=-1$) between states $\altmathcal{U}_{L}$ and $\tilde{\altmathcal{U}}_{L}$ and, likewise, across the rightmost wave ($\lambda_{+}=1$) between $\altmathcal{U}_{R}$ and $\tilde{\altmathcal{U}}_{R}$. In the frozen limit, EM fields can be discontinuous across the outermost waves but do not experience further jumps across the inner waves $\lambda_{L}$ and $\lambda_{R}$. For this reason we set $\tilde{\mathbf{E}}_{L}=\mathbf{E}^{*}_{L}=\mathbf{E}^{*}_{R}=\tilde{\mathbf{E}}_{R}$ (and similarly for $\mathbf{B}$). Specializing to the $x$-direction and solving Maxwell’s jump conditions for the transverse components yields

[TABLE]

Normal components $B_{x}$ and $E_{x}$ are continuous at the interface and do not experience any jump.

With the above definitions, the continuity equation is trivially fulfilled since both density and fluid velocity are continuous across $\lambda_{\pm}$. Likewise, it can also be verified that the jump conditions (67) of the momentum and energy equations are also automatically satisfied. In particular, by combining the jump conditions of the $x$-component of the momentum and energy, one consistently finds

[TABLE]

where $\tilde{\altmathcal{P}}_{\rm em}$ is the electromagnetic pressure.

4.2 Inner Solver

As discussed before, the frozen limit allows us to employ any hydrodynamic relativistic solver inside the inner fan. In order to capture the three-wave pattern characterizing the actual solution, our method of choice is based on the approximate HLLC solver of Mignone & Bodo (2005). Writing explicitly the jump conditions across $\lambda_{L}$ and $\lambda_{R}$ in the $x$-direction (similar results are obtained by index permutation) and taking advantage of the fact that EM fields do not experience any jump inside the inner fan, one finds

[TABLE]

where $S=L,R$ and $\altmathcal{W}=(D,\,Q_{x},\,Q_{y},\,Q_{z},\altmathcal{E}_{h})$, $\altmathcal{H}=(Dv_{x},\,Q_{x}v_{x}+p,\,Q_{y}v_{x},\,Q_{z}v_{x},\,Q_{x})$ are the subset of $\altmathcal{U}$ and $\altmathcal{F}$ containing only hydrodynamical quantities. Notice that, by virtue of Eq. (71), EM quantities have disappeared from the jump conditions and that the tilde sign can be dropped since since hydrodynamical variables are continuous across the outer waves, that is, $\tilde{\altmathcal{W}}_{L,R}=\altmathcal{W}_{L,R}$ and $\tilde{\altmathcal{H}}_{S}=\altmathcal{H}_{S}$ for $S=L,R$.

The system of equations (72) is thus identical to Eq. [16] of Mignone & Bodo (2005) and can be solved likewise by imposing continuity of pressure and normal velocity across the contact wave, $p^{*}=p^{*}_{L}=p^{*}_{R}$ and $\lambda^{*}=v^{*}_{x,L}=v^{*}_{x,R}$. To this purpose, we introduce the average state and flux

[TABLE]

and

[TABLE]

which allow us to rewrite the jump conditions (72) across $\lambda_{L}$ and $\lambda_{R}$ equivalently as

[TABLE]

where, again, $S=L,R$. Imposing the momentum-energy relation $Q_{x}=(\altmathcal{E}_{h}+p)v_{x}$ and combining the energy and $x-$component of the momentum jump conditions (75) leads to a quadratic equation for $\lambda^{*}$:

[TABLE]

where $[.]$ denote a specific component of the array. Eq. (76) has one physical admissible solution given by the negative branch. Once $\lambda^{*}$ is obtained, gas pressure in the star region is obtained from the momentum-energy equations as

[TABLE]

The hydrodyamic contribution to the interface flux is then given by

[TABLE]

The total flux is finally given by the sum of the hydro and EM contributions, the later being computed with the outer solver:

[TABLE]

where $n$ labels the direction of propagation of the discontinuity ($x$ in our example). Finally, the wave speed estimate for the fastest and slowest speeds $\lambda_{L}$ and $\lambda_{R}$ is done as in Mignone & Bodo (2005).

5 Numerical Benchmarks

We now present a number of standard numerical benchmarks aiming at verifying and assessing the robustness, accuracy and computational performance of our resistive CT relativistic scheme. Unless otherwise stated, test calculations will be carried out using the ideal equation of state (36) with $\Gamma=4/3$, the MHLLC Riemann solver and the van Leer limiter (41).

Since the maximum characteristic velocity is equal to the speed of light, the time step is easily determined by the relation

[TABLE]

where $C_{a}$ is the Courant number, $N_{\rm dim}$ is the number of spatial dimensions and $\alpha_{m}=1$ if $N_{\rm dim}\geq m$ and $\alpha_{m}=0$ otherwise. We set the default Courant number to $C_{a}=0.4$ for two-dimensional problems and $C_{a}=0.25$ for three-dimensional ones.

5.1 Telegraph Equation

In the first test we consider the propagation of light waves in presence of finite value of conductivity. Specifically, we want to solve the Maxwell’s equation in the fluid rest frame,

[TABLE]

Solutions of Eq. (81) also satisfy the telegraph equation which is obtained upon differentiating the previous system with respect to time,

[TABLE]

The system of Maxwell’s equation (81) admits plane wave solutions with wavenumber $k$ and frequency $\omega$ tied by the dispersion relation

[TABLE]

which, as noted in Mignone et al. (2018), is also the dispersion relation for the telegraph equation. Eq. (83) admits propagating modes when $\sigma<2|k|$. From the eigenvector of the system (81), we obtain the exact solution for a single mode,

[TABLE]

where $\phi(x,t)=kx-\mu t$, while $B_{1}$ is the initial perturbation amplitude. For $\sigma\neq 0$, Eq. (84) describes the propagation of monochromatic damped light waves with phase speed $v_{\phi}=\mu/k$ and thus with dispersive character. Damping is suppressed in the limit $\sigma\to 0$ (and thus $v_{\phi}\to 1$) where we recover the classical non-dispersive light wave propagation.

For numerical purposes we consider here an inclined wave front by rotating the 1D solution around the $z$-axis by an angle $\alpha$. The wavevector has orientation $\mathbf{k}=k_{x}(1,\tan\alpha,0)$, where $k_{x}=2\pi/L_{x}$ while $\tan\alpha=L_{x}/L_{y}$. The computational domain is $L_{x}=1$, $L_{y}=1/2$ so that $\tan\alpha=2$. Eq. (84) with $B_{1}=1$ and $\phi(x,0)=\mathbf{k}\cdot\mathbf{x}$ is used to initialize the electric and magnetic field inside the computational box. Three dimensional vectors are then rotated according to

[TABLE]

where a prime indicates quantities in the rotated frame. Instead of solving just Maxwell’s equation, we solve the full set of the resistive RMHD equations by prescribing a large fluid inertia ($\rho=10^{12}$) so that the fluid velocities become negligible and Maxwell’s equations decouple from the remaining conservation laws. Periodic boundary conditions are imposed everywhere.

Numerical results are shown in Fig. 4 for the CT (circles) and GLM (plus signs) schemes using different values of $\sigma=1,10,20$, corresponding to green, blue and red coloured lines. Computed profiles for $B_{z}$ are compared against the exact solution after one period $T=2\pi/\mu$ in the top left panel of Fig. 4 at the resolution of $64\times 32$ grid zones. Corresponding errors in L1 norm are plotted, as a function of the resolution $N_{x}=16,32,...,256$ ($N_{y}=N_{x}/2$), in the top right panel indicating that CT and GLM yields very similar results.

Although similar errors are produced in the $z$ component of magnetic field, results show significant differences by inspecting the normal component of the electric field. Note that no charge should be produced during the evolution since $E_{x}$ - the normal component of the electric field in the unrotated frame - should vanish identically. Nevertheless, propagation along an oblique direction (which is not the main diagonal) does not lead to perfect cancellation of multidimensional terms so that a non-solenoidal component of the current is generated at the truncation level of the scheme. This best illustrated in the bottom panels of Fig. 4 where we plot the L1 norm errors of $E_{x}$ (bottom left) and the maximum value of the charge (bottom right). Our CT method yields, overall, more uniform convergence with resolution when compared to GLM and the discrepancy between the two schemes worsen for larger values of $\sigma$. While the charge remains below $10^{-6}$ with the CT scheme, the GLM appears to produce spurious divergence with poor convergence rates. At the largest conductivity ($\sigma=20$, red lines), for instance, the two methods differ by over 3 orders of magnitude.

5.2 Rotated Shock-Tube

The shock-tube problem is a standard numerical benchmark consisting of an initial discontinuity separating two constant states. Here we adopt a configuration similar to the one presented in Palenzuela et al. (2009); Dionysopoulou et al. (2013); Miranda-Aranguren et al. (2018) and assign 1D left and right states as

[TABLE]

while $B_{x}=B_{y}=0$ and the electric field is computed from the ideal condition, $\mathbf{E}=-\mathbf{v}\times\mathbf{B}$. The ideal EoS (36) with $\Gamma=2$ is used.

We consider two variants of the problem. In the first one, the standard configuration with zero initial velocity is adopted: the electric field evolves by remaining perpendicular to the both the velocity and the magnetic field and the current density has a non-vanishing component only in the $y$-direction. No charge is therefore produced during the evolution since $E_{x}=0$ holds at any time. In the second variant, we prescribe an initial velocity $v_{y}=0.2$ everywhere so that the initial electric field presents a discontinuity at the origin and the charge is therefore $q(x)=-v_{y}(B_{zR}-B_{zL})\,\delta(x)$.

We rotate the initial condition by an angle $\alpha={\rm atan}(1/2)$ around the $z$ axis. The computational domain is defined by a rectangular box of width $[-1/2,1/2]$ in the $x$-direction while in the transverse direction we have $y\in[-r/2,r/2]$ where $r=N_{y}/N_{x}$. Zero-gradient boundary conditions hold at the rightmost and leftmost sides of the box, whereas for any flow variables $Q$ at the lower and upper boundaries we impose the translational invariance $Q(i,j)=Q(i\mp 1,j\pm 2)$. Vectors are rotated accordingly as in the previous section.

In the first variant of this test, we set $N_{x}=400,\,N_{y}=8$ and stop computations at $t=0.4/\sqrt{1+\tan^{2}\alpha}$ using different values of the resistivity in the range $\eta=10^{3}$ to $\eta=10^{-9}$ one decade apart. The $z$-component of the magnetic field (unchanged under rotation) is plotted in the upper left panel of Fig. 5 for the same values of $\eta$ also considered by Palenzuela et al. (2009). In the top right panel of the same figure we plot the error (in L1 norm) of the normal component of the electric field $E_{n}=E_{x}\cos\alpha+E_{y}\sin\alpha$ obtained with our CT schemes as well as with the GLM scheme for the entire range of values of $\eta$. While the errors are similar for the two schemes, integration with the GLM scheme was not possible for values of $\eta$ below $10^{-5}$ due to numerical instabilities. Errors increase sharply at around $\eta\approx 10^{-4}$ as the solution becomes progressively steeper and the scheme accuracy asymptotically reduces by one order of magnitude. For smaller values of $\eta$, numerical resistivity becomes comparable and no difference can be noticed at this resolution.

In the second variant of this test, the resistivity has been fixed to a large value ($\eta=10^{9}$) and a jump in the transverse velocity component is also present. A non-zero delta-spike charge appears since the normal component of $\mathbf{E}$ is now discontinuous. Both $E_{n}$ and the charge are advected at the local fluid velocity as shown in the bottom left panel of Fig. 5 where a closeup view in the range $x\in[0.1,0.2]$ is plotted. Our CT method yields a sharper representation of the discontinuity and the value of the charge is twice the value obtained with GLM (the exact value being infinitely large). Last, in the bottom-right panel, we show the maximum value of the charge density for different grid resolutions $N_{x}=100,200,400,800,1600$.

5.3 Magnetized Blast Wave

The magnetized blast wave problem consists of a cylindrical (in 2D) or spherical (in 3D) explosion propagating in a uniform magnetized medium. It is mainly used to test the ability of the scheme in handling oblique shock waves propagating in strongly magnetized environments as well as the fidelity in preserving the symmetric/antisymmetric properties of the solution. Potential flaws in the numerical scheme may easily lead to unphysical densities or pressures if the divergence-free condition is not properly controlled and the scheme does not introduce adequate dissipation across oblique discontinuous features (see, for instance, Mignone & Bodo, 2006; Mignone & Tzeferacos, 2010, and references therein).

5.3.1 Cylindrical Blast Wave

Two-dimensional versions of this problem have been considered by previous investigators, e.g., Komissarov (2007); Palenzuela et al. (2009); Mizuno (2013); Dionysopoulou et al. (2013); Miranda-Aranguren et al. (2018). Computations are carried out inside the Cartesian square $x,y\in[-6,6]$ filled with uniform density and pressure $\rho=p=10^{-3}$ for the exception of a small circular region, for $r\leq 0.8$, where density and pressure take the values $\rho=10^{-2}$ and $p=1$. Here $r=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius. For $0.8\leq r\leq 1$ continuity with the external environment is imposed through an exponential decay. Velocity and electric field are initialized to zero, while the magnetic field is oriented along the $x$-direction $\mathbf{B}=(0.1,0,0)$ (this is the configuration adopted by Komissarov, 2007; Miranda-Aranguren et al., 2018). Computations have been performed using $N_{x}\times N_{y}=200\times 200$ grid zones and the system is evolved until $t=4$. Zero-gradient boundary conditions have been imposed on all sides.

Results obtained with the CT algorithm are shown in the top panels of Fig. 6 where we display two-dimensional maps of the normal component of magnetic field, gas pressure and Lorentz factor. The explosion produces a fast forward shock that propagates (nearly) radially leaving behind a reverse shock that delimits the inner region where expansion takes place radially. Magnetic field lines pile up in the $y$-direction building up a shell of higher magnetic pressure. The gas moves preferably in the $x$-direction where it achieves a higher Lorentz factor ($\gamma_{\max}\approx 3.62$). Electric field and current have a non-vanishing component only in the $z$-direction and no charge is produced as $\nabla\cdot\mathbf{J}=0$ trivially holds. We have checked our results with the GLM scheme and found no appreciable differences, as confirmed by the one-dimensional profiles along the $x$- and $y$-directions in the corresponding bottom panels.

The computational overhead brought by the CT scheme, which is intrinsically multidimensional, is partially balanced out by the lower number of variables to be evolved ($14$ vs $11$). For this particular 2D configuration, for instance, we found that our CT method is approximately $\sim 5\%$ more expensive than the GLM scheme.

5.3.2 Spherical Blast Wave

The three-dimensional version of this problem has been formerly examined by Komissarov (2007); Dionysopoulou et al. (2013) by extending the previous configuration to a square Cartesian box $x,y,z\in[-6,6]$ with initial conditions identical to the 2D case but with $r=\sqrt{x^{2}+y^{2}+z^{2}}$ now being the spherical radius. Parameters are the same as for the cylindrical explosion. In particular, note that our initial magnetic field is $\mathbf{B}=0.1\hat{\boldsymbol{e}}_{x}$ and thus twice the value used by Dionysopoulou et al. (2013). The same configuration is used in the original work by Komissarov (2007) although we adopt here a much smaller value of resistivity, $\eta=10^{-6}$.

Fig. 7 shows 2D slices of the solution in the $xy$ plane at $t=4$ obtained with the CT scheme at the resolution of $200^{3}$ grid zones. The solution is qualitatively similar to the 2D case although a few differences are worth noticing. The inner region is delimited by a stronger reverse shock ($\gamma_{\max}\approx 5.5$) and encloses a stronger rarefaction wave when compared to the 2D case where gas pressure reaches much smaller values, $\min(p)\approx 4\times 10^{-6}$. The plasma is accelerated mostly in the $x$-direction to larger Lorentz factor, $\gamma\approx 5.5$, when compared to the 2D case ($\gamma\approx 3.5$). Another crucial difference is the local production of non-zero charge which was absent from the 2D case: this reveals an important difference between the stability and robustness of the two methods. Although both CT and GLM conserve charge up to machine accuracy ($\sim 10^{-15}$ is the total integrated charge for the two methods), our CT method runs without any problem whereas GLM failed already for values of the resistivity smaller than $10^{-2}$ (unless the time-step is lowered) owing to large-amplitude oscillations in the charge.

The relative computational cost between CT and GLM is larger in 3D than in 2D and, in our implementation, it reaches approximately $15\%$ for this particular case. This owes to the increased operation count which, in the 3D staggered method, accounts not only for the spatial dimensionality but also for the additional spatial reconstructions required by the multidimensional Riemann solver.

5.4 Stationary Charged Vortex

We propose, for the first time to the extent of our knowledge, a new exact equilibrium solution of the fully Res-RMHD equations. The solution is best described by adopting cylindrical coordinates $(r,\,\phi,\,z)$ and consists of a rotating flow with uniform density embedded in a vertical magnetic field $\mathbf{B}=(0,0,B_{z})$. In the ideal limit, this gives rise to a purely radial electric field, $\mathbf{E}=(E_{r},0,0)$ with $E_{r}=-v_{\phi}B_{z}$. Assuming purely radial dependence, the only non-trivial equations are the radial component of the momentum equation together with the $\phi$-component of Ampere’s law:

[TABLE]

where $q=\partial_{r}(rE_{r})/r$ is the charge and we have used the fact $\mathbf{J}=(0,qv_{\phi},0)$. The appearance of a charge is a consequence of the fact that we now have a radial electric field generated by a rotating flow (see, e.g., the review by Spruit, 2013). Using the ideal electric field condition, $E_{r}+v_{\phi}B_{z}=0$, the right hand side of Eq. (87) vanishes identically and the previous system of equations is rearranged as two independent ordinary differential equations, that is,

[TABLE]

where, using the same formalism already presented by Bodo et al. (2013, 2016), we have introduced $H^{2}(r)=B_{z}^{2}(r)-E_{r}^{2}(r)$. The $H^{2}(r)$ function can be chosen arbitrarily provided the following conditions are met: i) $H^{2}(r)>0$, which guarantees that velocity remain subluminal and ii) $dH^{2}(r)/dr(r=0)\geq 0$ which ensures that $E_{r}^{2}\geq 0$. The equality sign holds at $r=0$ where the electric field must vanish. A simple solution that satisfies the previous requisites is

[TABLE]

where $q_{0}$ is the charge at $r=0$. Using Eq. (91), $E_{r}$ can be found by differentiating $H^{2}(r)$ with respect to $r$ (equation 90), $B_{z}$ and $v_{\phi}$ follow from Eq. (91) and the ideal condition while gas pressure can be obtained by solving the differential equation (89). The final result is:

[TABLE]

where $\Gamma_{1}=\Gamma/(\Gamma-1)$. We set density to unity while $p_{0}=0.1$ gives the pressure at infinity. Charge is set to $q_{0}=0.7$. Notice that the previous solution is also an exact solution of the ideal RMHD equations and that, using Eq. (109), the rest-frame charge can be written as $q_{0}=-(\nabla\times\mathbf{v})\cdot\mathbf{B}$ at $r=0$.

We carry out computations on the two-dimensional Cartesian square $x,y\in[-10,10]$ using a uniform resolution of $N_{x}\times N_{y}$ zones and evolve the system until $t=5$. Boundary values are held fixed to the equilibrium solution throughout the integration. Note that the equilibrium condition (92) does not depend on the resistivity and numerical solutions carried out with different values of $\eta$ depend solely on the stability of the algorithm used for this particular problem. This has been verified for a wide range of the resistivity parameter, namely, $\eta\in[10^{3},10^{-8}]$ using a grid resolution $N_{x}=N_{y}=256$. Errors in L1 norm for the charge are plotted in the top left panel of Fig. 8 as a function of $\eta$. Our results confirm that the CT scheme remains stable for any value of the resistivity parameter in the chosen range. In contrast, results obtained with GLM scheme give good agreement only for large values of $\eta$ while numerical instabilities is exhibited for $\eta\lesssim 10^{-3}$. In the bottom left panel, we plot the errors (in $L1$ norm) of gas pressure as function of the resolution ($N_{x}=N_{y}=32,64,128,256,512$) showing second-order convergence for both CT and GLM schemes. Here $\eta=1$ and $\eta=10^{-2}$ have been used for the computations.

Numerical results for $\eta=10^{-3}$, which at the resolution of $256$ zones sets the verge of stability for GLM, can be analyzed in the central panels of Fig. 8 where we show coloured maps of pressure overlaid with iso-contour levels of the electric field magnitude (top and bottom). Large oscillations are found in proximity of the coordinates origin with the GLM scheme (bottom central panel). These large overshoots can also be recognized in the 1D horizontal cuts of charge $q$ and $E_{y}$ shown, respectively, in the top and bottom right panels using both CT (blue circles) and GLM (red X symbol) at $t=3$. Note that the exact solution for the $y$-component of the electric field should vanish for $y=0$ but large-amplitude oscillations are clearly visible with GLM. We point out that an increase in the resolution - which is accompanied by a reduction of the time step - extends the stability of the GLM method to lower values of the resistivity. This result is consistent with the large errors introduced by the stiffness of the charge evolution equation.

5.5 Tearing Mode

In the next example, we investigate the linear growth of the ideal tearing mode in a Harris-like current sheet, following the work of Del Zanna et al. (2016) and Miranda-Aranguren et al. (2018). The initial configuration consists of an initially static ($\mathbf{v}=0$) and uniform plasma with constant density and pressure value, $\rho_{0}$ and $p_{0}$. The initial magnetic field satisfies the force-free condition

[TABLE]

where $a$ is the current sheet thickness. Useful parameters are the magnetization $\Sigma=B_{0}^{2}/\rho_{0}$, the plasma beta $\beta=2p_{0}/B_{0}^{2}$, the Alfvén velocity $v_{a}=B_{0}/\sqrt{B_{0}^{2}+w_{0}}=B_{0}/\sqrt{B_{0}^{2}+\rho_{0}+4p_{0}}$ (for an ideal gas law with adiabatic index $\Gamma=4/3$), and the Lundquist number $S=v_{a}L/\eta\gg 1$, where $L$ is a typical spatial scale and $\eta$ the resistivity. The electric field is initially zero everywhere. According to the MHD works by Pucci & Velli (2014); Landi et al. (2015), when extremely thin current sheets are considered the tearing mode ceases to be a slow process (with growth rate $\gamma_{\mathrm{TM}}\propto S^{-1/2}$) and reconnection occurs instead on the ideal Alfvénic time $\tau_{a}=L/v_{a}$. The threshold is provided by the critical (inverse) aspect ratio

[TABLE]

for which, in the asymptotic limit of large $S$ but independently on the actual value of $S$, $\gamma_{\rm TM}\simeq 0.6v_{a}/L$. In the relativistic regime, where $v_{a}\to c=1$, this so-called ideal tearing instability thus becomes a very efficient process.

In order to trigger the instability we perturb the initial configuration with a single-mode magnetic field equal to

[TABLE]

where $\epsilon=10^{-4}$ is the initial perturbation amplitude and $k$ the wave number. For computational purposes, we initialize the magnetic field in the $x-y$ plane using the $z$-component of the vector potential

[TABLE]

which ensures an initial divergence-free discretization of $\mathbf{B}$ at the beginning. Following Del Zanna et al. (2016), we use for this test $L=B_{0}=\Sigma=\beta=1$, hence $\rho_{0}=1$, $p_{0}=0.5$, $v_{a}=0.5$ and $S=10^{6}$, so that $a=0.01$ and $\eta=5\times 10^{-7}$. The theory predicts, for the ideal tearing mode, a wave number for the fastest growing mode of $k_{\mathrm{max}}=1.4\,S^{1/6}=14$, providing $\gamma_{\mathrm{TM}}\simeq 0.3$. However, probably due to the diffusion of the initial force-free field or to compressible effects, Del Zanna et al. (2016) found a maximum in the dispersion relation curve for $k=12$ and $\gamma_{\mathrm{TM}}\simeq 0.27$, hence, also following Miranda-Aranguren et al. (2018), we will adopt $k=12$ as the standard test value for the wave vector.

We perform computations on the rectangular domain $x\in[-20a,+20a]=[-0.2,+0.2]$ and $y\in[0,2\pi/k]=[0,0.5236]$, using free outflow conditions at the $x$ boundaries and periodicity along $y$. For the sake of comparison, we have repeated calculations using our CT implementation as well as the GLM scheme with different grid resolutions corresponding to $N_{x}\times N_{x}/4$ zones with $N_{x}=512,1024,2048$. The monotonized central limiter was used for these computations. We note that GLM requires, for stability, a smaller CFL number $C_{a}=0.2$, whereas computations with CT are carried at twice this value, $C_{a}=0.4$.

In the left panels of Fig. 9 we show horizontal cuts of the $x$-components of magnetic field and velocity at $t=10$ using CT (solid lines) and GLM (dashed lines) at the chosen grid size using different colors. At the resolutions of $Nx=1024,2048$ the profiles well approximate the eigenfunctions given by Del Zanna et al. (2016) (see their Figure 1). Overall, the CT scheme performs similarly to the GLM method albeit with slightly reduced numerical dissipation and better convergence with resolution.

Perturbations are expected to grow exponentially as $Q_{1}\propto e^{\gamma_{TM}t}$ and we have measured the numerical growth rate by considering, as suggested by Miranda-Aranguren et al. (2018), the integral of the $x$-component of magnetic field (squared),

[TABLE]

Plots of $f(t)$ are shown in the central panels of Fig. (9) for the LF (top) and the MHLLC (bottom) Riemann solvers. Our results indicate that increasing the resolution leads to smaller growth rates, in agreement with the previous findings. However, two distinct phases can be discerned using the MHLLC solver: a steeper growth for $t\lesssim 6-8$ followed by a softer one for $t\gtrsim 6-8$, the actual value depending on the resolution. The behavior remains unaltered when switching from CT to GLM and it is not observed by Del Zanna et al. (2016) and Miranda-Aranguren et al. (2018) who used $5$-th or higher-order reconstructions. For our second-order scheme, instead, we attribute this behavior to compressible effects enhanced by the resolution of density jumps, probably triggering spurious modes that grow faster in the early stage of evolution. For this reason, the growth rates are computed as the slope obtained by the linear fit of Eq. (97) versus time over the interval $6\leq t\leq 12$ and reported in Table 1. Results with the MHLLC solver show better convergence rates with resolution when compared to the more diffusive LF scheme (only few percent with resolution doubling). At the largest resolution, the MHLLC Riemann solver yields larger growth rates than LF and converges to the actual value ($\gamma_{\rm TM}\approx 0.30$).

It is also instructive to compare the charge evolution obtained with the GLM and CT methods, as shown in the right panels of Fig. 9 where we plot the maximum value of the charge versus time. With the GLM scheme, a systematic excess of charge is produced which is noticeably reduced by doubling the resolution (approximately one order of magnitude with the LF solver). Fluctuations with the CT scheme are much less affected by the grid size and are restrained within a factor of 2. Finally we show, in Fig. 10, the charge density distribution obtained with the GLM and CT schemes at $t=10$ by narrowing the view down to the $x$-axis. With the CT method, charge density is mostly concentrated around the current sheet and the solution appears to be well-behaved and oscillation-free. On the contrary, with the GLM scheme, charge distribution spreads out on the sides and the solution develops large overshoots as well as spurious oscillations which appear to have numerical origin.

5.6 Kelvin-Helmholtz Flow

As a test application, we consider the evolution of a double shear layer as already presented, in the context of resistive relativistic flows by Mizuno (2013) and in the ideal case by Mignone et al. (2009); Beckwith & Stone (2011). The initial condition consists of a background (double) shear layer with non-uniform density distribution,

[TABLE]

where $\phi(y)=(|y|-1/2)/a$. The background flow is then perturbed by

[TABLE]

Following Mizuno (2013) we take the velocity of the shear layer to be $v_{\rm sh}=1/2$ with thickness $a=10^{-2}$ and the density contrast $\rho_{h}=1$, $\rho_{l}=10^{-2}$. The parameters in the last equations, $\alpha=10^{-1}$ and $A_{0}=0.1$, characterize the perturbation cutoff height and amplitude. The magnetic field is initially constant and uniform with a poloidal ($B_{x}$) and toroidal ($B_{z}$, out of the plane) components parameterized by

[TABLE]

where $\mu_{p}=0.01$ and $\mu_{t}=1$ are magnetization parameters.

The Cartesian box used for the integration of the resistive RMHD equations is defined by $x\in[-1/2,1/2]$, $y\in[-1,1]$ using different values of the conductivity, $\sigma=1/\eta=0,10,10^{2},10^{3}$ and $10^{5}$. We perform computations using the MHLLC solver and the MC limiter (42) until $t=15$ with the nominal resolution of $256\times 512$ grid zones. Computations with the CT scheme remained stable at the nominal Courant number ($C_{a}=0.4$) for any value of the conductivity while numerical instabilities occurred at large $\sigma$’s with the GLM scheme unless the CFL number was lowered to $0.1$.

The growth rates, computed as $\Delta v_{y}=(\max(v_{y})-\min(v_{y}))/2$, are shown in the top left panel of Fig. 11 for selected values of the conductivity. Our measured growth rates favourably compare to those of Mizuno (2013), indicating that different values of the conductivity have a negligible impact on the growth of the instability. The $\sigma=0$ case (purple dashed-triple dotted line), which in the work by Mizuno (2013) yielded a smaller growth rate, did not make particular difference in our case. In all likelihood, this discrepancy can be attributed to the choice of the Riemann solver as it can be inferred from the bottom left panel of Fig. 11, where GLM and CT schemes are compared using the LF (red) and MHLLC (blue) solvers. Switching from the former to the latter leads to a larger growth rate which become closer to the high-conductivity case (black dotted line). In the top right panel, we plot the poloidal field amplification which is enhanced with increasing conductivity in accordance with Mizuno (2013). However, our five-wave solver leads to a steeper poloidal field amplification when compared to the scheme of Mizuno (2013) and to an earlier start of the saturation phase ($t\approx 4$ instead of $t\approx 5$, see Figure 11 of Mizuno, 2013). A resolution study, shown in the bottom right panel of Fig. 11 for $\sigma=10^{5}$, confirms that the MHLLC solver yields larger growth rates and faster convergence when compared to the simpler LF scheme (similar results were found with the HLLD Riemann solver in the Kelvin-Helmholtz test presented by Mignone et al., 2009).

The time evolution is shown in Fig. 12 for the CT scheme with $\sigma=10^{5}$ at three different snapshots, $t=3,7$ (left and middle panels) and $t=10$ (right panel). Top and bottom panels show, respectively, coloured distributions of density and the poloidal to toroidal magnetic field ratio ($\sqrt{B_{x}^{2}+B_{y}^{2}}/B_{z}$). The linear phase is followed by vortex formation and the transition to the nonlinear regime during which the mixing layer enlarges and magnetic field becomes amplified and stretched into filamentary structures. Note also the formation of the intermediate vortex which does not appear when using the LF Riemann solver.

Computations with the CT scheme required approximately $5\%$ more time than the GLM case. The same result was established for the 2D Blast wave test problem, see section 5.3.1.

5.6.1 Extension to Three-Dimensions

We extend the previous configuration to three-dimensions by choosing (as in Beckwith & Stone, 2011) the computational domain $x,z\in[-1/2,1/2]$ and $y\in[-1,1]$, using a uniform grid of $256\times 512\times 256$ zones. We employ the same initial condition with the exception of the $z$-component of velocity which is now prescribed as in Eq. (99) with the $\sin()$ function replaced by a random number distribution in the range $[-1,1]$. In order to assess the robustness of our algorithm we employ the most stringent value of $\sigma=10^{5}$ using the CT scheme with the MHLLC and the LF solvers. The system is evolved until $t=30$.

The growth rates, computed again as $\Delta v_{y}=(\max(v_{y})-\min(v_{y})/2$, are plotted in the top panel of Fig. 13 for the two Riemann solvers. In analogy to the 2D case, we observe a faster growth and an earlier transition to the nonlinear regime when the MHLLC solver is employed.

The nonlinear evolution is illustrated in Fig. 14 through a series of volume rendering of density (top) and poloidal to toroidal magnetic field (bottom) at $t=10,20$ and $t=30$. The large-scale motion remains confined along the initial shear direction and sheet-like thin structures, where most of the magnetic field energy is trapped, characterize the turbulent state. These rounded slabs remain roughly parallel to the $z$-direction although significant medium-scale structures develop in this direction for $t\gtrsim 20$ around the shear layer.

In order to quantify the numerical diffusion of the two different solvers (which affects the amount of turbulent structure), we evaluate the spectral energy density as

[TABLE]

where $\mathscr{F}_{x}(k_{x},y,z)$ is the one-dimensional fast Fourier transform of the density taken across the $x$ dimension. Likewise, we construct $P(k_{y})$ and $P(k_{z})$ by index permutation of (101). The spectral densities across the three directions are plotted in the bottom panel of Fig. 13 for the MHLLC (blue lines) and the LF (red lines) solvers using solid, dashed and dashed-dotted lines corresponding, respectively, to $P(k_{x})$, $P(k_{y})$ and $P(k_{z})$. Overall we see that power spectra in the MHLLC case are systematically larger by $\sim 2$ orders of magnitudes (at large wavenumbers) in the $x$- and $y$-directions when compared to the LF case. While most small-scale power resides in these two directions, the same trend is also observed in turbulent structures across the $z$-direction which make a significant contribution at large-moderate wavelengths. Our results favourably compare to those of Beckwith & Stone (2011) where the HLLD and HLL solvers have been compared on the same test problem.

While simulations with the GLM scheme could not be completed for such a small value of the resistivity (if not by considerably lowering the Courant number), the results presented in this section demonstrates the our CT scheme is remarkably more stable and robust also in the three-dimensional case. The computational performance of our CT method yielded a $\sim 13\%$ additional overhead when compared to the same calculation with GLM method (with lower value of $\sigma$). A comparable performance was obtained in the 3D Blast wave problem, §5.3.2.

6 Summary

In this work we have presented a new second-order Godunov-type constrained transport method for the solution of the resistive relativistic MHD equations in the context of IMplicit-EXplicit (IMEX) Runge-Kutta methods. Our method follows the work of Balsara et al. (2016) and sets the primary components of electric and magnetic fields at zone interfaces. A discrete version of Stoke’s theorem is employed for the evolution of the electromagnetic fields while hydrodynamic variables are instead located at the zone-center and are treated in the usual finite-volume sense.

During the explicit stages of the IMEX time-stepping, numerical fluxes needed for the update of electromagnetic fields are obtained by applying a two-dimensional Maxwell solver at zone edges while a standard one-dimensional Riemann solver is used at zone interfaces to advance zone-centered hydrodynamic variables. This introduces proper upwinding and it ensures that Faraday’s law for the the magnetic field is advanced in a divergence-free fashion while the ensuing discretization for the electric field conserves charge to machine precision. When dealing with the stiff source term, a solution approach for the implicit update of staggered electric field at zone faces that retains the point-local character has been developed. This has been shown to be formally equivalent to a conservative scheme in which charge is upwinded using a local Lax-Friedrichs flux. The proposed method of discretization is consistent with Ampere’s law from which charge conservation directly follows at the continuous level.

In addition, we have also introduced a new Riemann solver based on the frozen condition of the underlying hyperbolic system of conservation laws with stiff relaxation source terms. Owing to the weak coupling between Maxwell’s and hydrodynamics equations inherent in this limit, the solution to the Riemann problem can be approached by the combination of an outer solver for the EM waves and an inner solver for resolving hydrodynamic waves. Our composite MHLLC Riemann solver (where “M” denotes the outer Maxwell while HLLC is the Harten-Lax-van Leer of Mignone & Bodo (2005) applied to the hydrodynamic equations) has reduced numerical diffusion when compared to the traditional Lax-Friedrichs (or HLL) solver and it shares a solution procedure analogous to the one outlined in the appendix of Miranda-Aranguren et al. (2018).

An extensive suite of two- and three-dimensional numerical benchmarks with some new analytical solutions has been used to assess the performance of the newly proposed CT method. A direct comparison with the widespread GLM scheme of Palenzuela et al. (2009), that employs a conservative cell-centered discretization, reveals that our CT scheme gives comparable (albeit less diffusive) results for moderate or large resistivities although its benefits are more evident in the ideal limit - small value of the resistivity - in problems where a net charge is produced. In this regime, we have found that our CT scheme is markedly more robust than the GLM method which instead fails when the time step becomes smaller than the resistivity owing to large spurious oscillations in the charge. We argue that this may result from the (unstable) explicit discretization of the charge equation used in the standard GLM formalism, where the divergence of the current introduces stiffness. In support of this, we point out that variants of the GLM scheme in which the charge is not an evolutionary equation and it is computed directly from the divergence of the electric fields (e.g. Dionysopoulou et al., 2013) do not seem to suffer from this behavior although more investigation is certainly needed. Similar conclusions can be drawn for the CT scheme of Bucciantini & Del Zanna (2013), in which the electric field retain a zone-centered representation.

Finally, we point out that spurious (local) charge production may occur at the numerical level in both schemes, owing to discretization errors introduced when taking the divergence of the current. Most likely, these issues can be ameliorated by introducing schemes with spatial order of accuracy greater than two. Higher-order reconstruction (such as WENO or PPM) can be easily accommodated for in our formulation although we postpone genuinely third (or higher) finite volume schemes to forthcoming works.

Acknowledgements

The authors acknowledge support from the PRIN-MIUR project Multi-scale Simulations of High-Energy Astrophysical Plasmas (Prot. 2015L5EE2Y). We would also like to thank the referee for his/her constructive comments which helped to improve the quality of this manuscript.

Appendix A On the Rest Frame Charge Density

We illustrate here the relationship between the vorticity of the flow and the charge density, using the covariant formalism. The kinematic vorticity four-vector is defined as (Rezzolla & Zanotti, 2013)

[TABLE]

where clearly $\omega^{\mu}u_{\mu}=0$. Thanks to this definition, the covariant derivative of the fluid velocity can be split as

[TABLE]

where $a^{\mu}=(u^{\nu}\nabla_{\nu})u^{\mu}$ is the acceleration, normal to $u^{\mu}$ too. Using the definitions of the comoving electromagnetic fields in Eq. (7) we derive the following relation

[TABLE]

Let us now take the divergence of the comoving electric field. Using Maxwell’s equations we write

[TABLE]

hence, recalling the definition of the comoving charge density $q_{0}$, we find

[TABLE]

which differs from the usual Gauss’ law.

The above equation provides a link between the evolution of the comoving charge and the electromagnetic fields. Notice that in the ideal limit (or for small values of the resistivity) the terms with the (comoving) electric field $e^{\mu}=\eta j^{\mu}$ are negligible compared to the one with the magnetic field. Hence in this case we find the simple relation

[TABLE]

providing $q_{0}$ directly and involving the kinematic vorticity and the magnetic field alone. The above scalar product is a relativistic invariant, therefore it is convenient to calculate it in the comoving frame of the fluid, that is, for a flat spacetime metric

[TABLE]

where $\boldsymbol{\omega}=\mathbf{\nabla}\times\mathbf{u}$ is the vorticity three-vector (recall that the velocity vanishes but not its spatial derivatives). The comoving charge density then becomes, in this case

[TABLE]

Appendix B Jacobian of the IMEX-Newton Method

In order to apply the Newton-Broyden scheme during the implicit stage of our SSP-IMEX scheme, the Jacobian (60) must be computed:

[TABLE]

where $\varepsilon_{ijk}$ is the Levi-Civita symbol.

For the ideal equation of state (36) the gradient of the specific enthalpy $h=w/\rho$ is obtained as:

[TABLE]

while the gradient of the pressure is obtained by differentiating the second in (37) while keeping ${\cal E}$ and $\mathbf{B}$ constant:

[TABLE]

The second term on the right hand side of Eq. (110) involves derivatives of the electric field. Using Eqns. (52) and (53), we split them as the sum of three terms:

[TABLE]

where

[TABLE]

We proceed by first computing the derivatives of ${\cal A}$ which is a scalar function:

[TABLE]

Next we differentiate $\mathbf{H}$:

[TABLE]

To obtain the derivatives of the last term ($\mathbf{K}$) in Eq. (114) we first define the quantity

[TABLE]

together with its derivatives,

[TABLE]

so that the gradient of $\mathbf{K}$ with respect to velocity is calculated as

[TABLE]

Putting all toegther, the derivatives of the electric field are finally given by

[TABLE]