# A Fourth-Order Finite Volume Scheme for Resistive Relativistic Magnetohydrodynamics

Andrea Mignone, Vittoria Berta, Marco Rossazza, Matteo Bugli, Giancarlo Mattia, Luca del Zanna, Lorenzo Pareschi

► To cite this version:

Andrea Mignone, Vittoria Berta, Marco Rossazza, Matteo Bugli, Giancarlo Mattia, et al. A Fourth-Order Finite Volume Scheme for Resistive Relativistic Magnetohydrodynamics. Monthly Notices of the Royal Astronomical Society, 2024, 533 (2), pp.1670-1686. 10.1093/mnras/stae1729. hal-04661926

HAL Id: hal-04661926

https://hal.science/hal-04661926v1

Submitted on 21 Oct 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A fourth-order accurate finite volume scheme for resistive relativistic MHD

A. Mignone<sup>D, 1★</sup> V. Berta,<sup>1★</sup> M. Rossazza,<sup>1</sup> M. Bugli<sup>D, 1,2,3</sup> G. Mattia<sup>D, 4</sup> L. Del Zanna<sup>D 5,4,6</sup>
and L. Pareschi<sup>7,8</sup>

<sup>1</sup>Dipartimento di Fisica, Università degli Studi di Torino, Via Pietro Giuria 1, I-10125 Torino, Italy

<sup>2</sup>Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette, France

<sup>3</sup>INFN – Sezione di Torino, Via Pietro Giuria 1, I-10125 Torino, Italy

<sup>4</sup>INFN – Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino (FI), Italy

<sup>5</sup>Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino (FI), Italy

<sup>6</sup>INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

<sup>7</sup>Maxwell Institute for Mathematical Sciences and Department of Mathematics, Heriot-Watt University, Colin Maclaurin Building, Edinburgh EH14 4AS, UK

<sup>8</sup>Dipartimento di Matematica e Informatica, Università degli Studi di Ferrara, Via N. Machiavelli 30, I-44121 Ferrara, Italy

Accepted 2024 July 6. Received 2024 July 3; in original form 2024 March 19

## ABSTRACT

We present a finite-volume, genuinely fourth-order accurate numerical method for solving the equations of resistive relativistic magnetohydrodynamics in Cartesian coordinates. In our formulation, the magnetic field is evolved in time in terms of face-average values via the constrained-transport method, while the remaining variables (density, momentum, energy, and electric fields) are advanced as cell volume averages. Spatial accuracy employs fifth-order accurate WENO-Z reconstruction from point values (as described in a companion paper) to obtain left and right states at zone interfaces. Explicit flux evaluation is carried out by solving a Riemann problem at cell interfaces, using the Maxwell-Harten-Lax-van Leer with contact wave resolution. Time-stepping is based on the implicit-explicit Runge-Kutta (RK) methods, of which we consider both the third-order strong stability preserving SSP3(4,3,3) and a recent fourth-order additive RK scheme, to cope with the stiffness introduced by the source term in Ampere's law. Numerical benchmarks are presented in order to assess the accuracy and robustness of our implementation.

**Key words:** magnetic reconnection – MHD – plasmas – relativistic processes – methods: numerical – software: development.

## 1 INTRODUCTION

The investigation of relativistic plasma dynamics holds crucial importance in understanding the complex nature of high-energy astrophysical phenomena. Under typical astrophysical conditions, resistivity remains exceptionally low, and the ideal framework well describes processes occurring on rapid time-scales. Nevertheless, the evolving flow dynamics can cause the formation of local regions with steep gradients, such as current sheets, where the influence of resistivity can no longer be ignored.

It appears then clear the necessity of numerical simulations of astrophysical plasmas beyond the ideal magnetohydrodynamics (MHD) assumptions. Resistive relativistic magnetohydrodynamics (Res-RMHD) is indeed a framework that takes into account the contribution of finite plasma conductivity. In Res-RMHD, the comoving electric field no longer vanishes as in ideal MHD and it is directly proportional to the local current density, which is characterized by a scalar resistivity coefficient denoted as $\eta$.

Numerical methods for the solutions of the Res-RMHD equations have been now presented by a number of authors in the context of second-order conservative schemes (see e.g. Komissarov 2007; Palenzuela et al. 2009; Bucciantini & Del Zanna 2013;

Dionysopoulou et al. 2013; Mizuno 2013; Miranda-Aranguren, Aloy & Rembiasz 2018; Mignone et al. 2019; Cheong et al. 2022; Nakamura et al. 2023, and references therein) or higher order schemes (Dumbser & Zanotti 2009; Bugli, Del Zanna & Bucciantini 2014; Del Zanna, Bugli & Bucciantini 2014; Tomei et al. 2020).

A delicate aspect of their robust implementation is the stability of the method while dealing with the stiff relaxation term in Ampere's law, posing very strict constraints for time-explicit calculations. This makes the implementation of Res-RMHD numerical schemes more challenging than their ideal counterpart. In this sense, the employment of implicit-explicit (IMEX) Runge-Kutta (RK; see Pareschi & Russo 2005) represents an effective solution to the problem, allowing larger time-steps to overcome stiffness.

Another important aspect is related to the divergence-free condition and charge conservation that are not automatically fulfilled at the discrete level. To cope with this issue, here we follow an approach similar to Mignone et al. (2019, hereafter Paper I), by adopting the constrained transport (CT) approach to control the divergence of the magnetic field. In the CT formalism (see Evans & Hawley 1988; Balsara & Spicer 1999; Londrillo & del Zanna 2004; Mignone & Del Zanna 2021, and references therein), the magnetic field components are evolved as area-weighted quantities positioned on the faces orthogonal to the corresponding component direction. This ensures that the condition $\nabla \cdot \nabla \times \mathbf{E} = 0$ is satisfied also at a discrete level. In Paper I, the CT method has been extended to control

© 2024 The Author(s).

<footnote>* E-mail: andrea.mignone@unito.it (AM); vittoria.berta@unito.it (VB)</footnote>
the divergence of the electric field as well, thus ensuring a consistent approach for charge conservation. The resulting scheme has shown to be more robust when compared to cell-centred frameworks such as the generalized Lagrange multipliers (Munz et al. 2000; Dedner et al. 2002). In this work, however, we choose to constrain only the magnetic field and evolve the electric field as a zone-centred variable without explicitly enforcing its divergence, simply replacing the charge density with $\nabla \cdot E$ where needed. As demonstrated later and from extensive testing (not part of this work), we have not found any appreciable evidence of loss of accuracy or stability. Similar conclusions have been reported by Bucciantini & Del Zanna (2013), Del Zanna et al. (2016), and Tomei et al. (2020).

While second-order numerical schemes already provide an effective tool for solving such equations, higher order methods can play a crucial role in the context of Res-RMHD simulations. An increase in overall accuracy and convergence order of the numerical integration, in fact, can significantly decrease the numerical dissipation associated with the scheme's discretization errors. The benefits of such improvements are twofold. On the one hand, higher order schemes can generally lead to higher computational efficiencies, thus decreasing the computational time required to achieve a given accuracy (Berta et al. 2024). On the other hand, a decrease in numerical dissipation means that we can probe the dynamical impact of more realistic (i.e. lower) values of physical magnetic dissipation, which needs to dominate over the numerical one in any reliable Res-RMHD model (Mattia et al. 2023).

In this paper, we extend the fourth-order accurate finite volume (FV) method described in Berta et al. (2024, hereafter Paper II) to the Res-RMHD equations. Our scheme employs pointwise WENO-Z spatial reconstruction (see Paper II and the original paper by Borges et al. 2008) and either the strong stability preserving (SSP) SSP3(4,3,3) IMEX RK of Pareschi & Russo (2003, 2005) or the more recent IMEX additive Runge-Kutta scheme ARK4(3)7L[2]SA1 by Kennedy & Carpenter (2019) to achieve, respectively, third- or fourth-order accuracy in time. We refer to Carpenter et al. (2005), Conde et al. (2017), and Boscheri & Pareschi (2021) for other IMEX schemes of high order. Note that the use of pointwise spatial reconstructions is indeed crucial to achieve accuracy beyond the second order when evaluating the implicit stiff source terms within the IMEX formulation, as discussed recently in Boscarino, Russo & Semplice (2018).

The rest of the paper is structured as follows. In Section 2, we review the fundamental equations of the Res-RMHD framework. In Section 3, our fourth-order FV-CT method is illustrated, and in Section 4, we assess the robustness of the algorithm through extensive numerical benchmarking. Results are summarized in Section 5.

## 2 RESISTIVE RELATIVISTIC MHD EQUATIONS

When dealing with relativistic plasmas, it is convenient to start from the appropriate set of covariant equations. In the (single fluid) MHD approximation, these consist of the Maxwell equations:

$$ \nabla_{\mu} F^{\mu\nu} = -J^{\nu}, \quad \nabla_{\mu} F^{*\mu\nu} = 0, \qquad (1) $$

where $F^{\mu\nu}$, $F^{*\mu\nu}$, and $J^{\mu}$ are, respectively, the Faraday tensor, its dual, and the four-current density, and in the conservation laws for mass and total (gas and electromagnetic fields) energy–momentum:

$$ \nabla_{\mu} (\rho u^{\mu}) = 0, \quad \nabla_{\mu} T^{\mu\nu} = \nabla_{\mu} (T_{\text{gas}}^{\mu\nu} + T_{\text{em}}^{\mu\nu}) = 0, \qquad (2) $$

where $\rho$ is the rest-mass density, $u^{\mu}$ is the four-velocity of the fluid, and $T^{\mu\nu}$ is the total energy–momentum tensor. The two contributions

are, respectively,

$$ \begin{aligned} T_{\text{gas}}^{\mu\nu} &= \rho h u^{\mu} u^{\nu} + p g^{\mu\nu}, \\ T_{\text{em}}^{\mu\nu} &= F_{\lambda}^{\mu} F^{\nu\lambda} - \frac{1}{4} (F_{\lambda\kappa} F^{\lambda\kappa}) g^{\mu\nu}, \end{aligned} \qquad (3) $$

where $h$, $\varepsilon$, and $p$ are, respectively, the gas specific enthalpy, energy density, and pressure, as measured in the fluid rest frame and related by $\rho h = \varepsilon + p$. These thermodynamic quantities must be also connected by an equation of state (EoS), for instance, given in the form $h = h(\rho, p)$.

It is convenient to express the Faraday tensor and its dual in terms of the four-velocity and the magnetic $b^{\mu}$ and electric $e^{\mu}$ fields measured in the fluid rest frame:

$$ \begin{aligned} F^{\mu\nu} &= u^{\mu} e^{\nu} - u^{\nu} e^{\mu} + \epsilon^{\mu\nu\lambda\kappa} b_{\lambda} u_{\kappa}, \\ F^{*\mu\nu} &= u^{\mu} b^{\nu} - u^{\nu} b^{\mu} - \epsilon^{\mu\nu\lambda\kappa} e_{\lambda} u_{\kappa}, \end{aligned} \qquad (4) $$

where $\epsilon^{\mu\nu\lambda\kappa}$ represents the Levi-Civita pseudo-tensor. Ideal relativistic MHD is obtained by assuming a vanishing electric field in the frame comoving with the fluid, while in the resistive case we will assume a simple Ohmic relation

$$ e^{\mu} = \eta j^{\mu}, \qquad (5) $$

where $j^{\mu} = J^{\mu} + (J^{\nu}u_{\nu})u^{\mu}$ is the current density measured in the fluid frame and $\eta$ is the resistivity coefficient, supposed to be a scalar for the sake of simplicity.

Assuming from now on a Minkowski flat space-time and splitting time and space, it is convenient to write the fluid velocity as

$$ u^{\mu} = (\sqrt{1+u^2}, \mathbf{u}) \equiv (\gamma, \gamma \mathbf{v}), \qquad (6) $$

where $\mathbf{v}$ is the usual three-dimensional velocity of the fluid measured in the laboratory frame and $\gamma$ the associated Lorentz factor, while the electromagnetic fields can be split as

$$ \begin{aligned} e^{\mu} &= (\mathbf{u} \cdot \mathbf{E}, \gamma \mathbf{E} + \mathbf{u} \times \mathbf{B}), \\ b^{\mu} &= (\mathbf{u} \cdot \mathbf{B}, \gamma \mathbf{B} - \mathbf{u} \times \mathbf{E}), \end{aligned} \qquad (7) $$

where $\mathbf{E}$ and $\mathbf{B}$ are the usual electric and magnetic fields measured in the laboratory frame. The four-current can be also split

$$ J^{\mu} = (q, \mathbf{J}) = (q, q\mathbf{v} + \eta^{-1}\mathbf{E}), \qquad (8) $$

in which $q$ is the charge density measured in the laboratory frame and $\mathbf{J}$ the spatial current density, with

$$ \tilde{\mathbf{E}} = \gamma \mathbf{E} + \mathbf{u} \times \mathbf{B} - (\mathbf{E} \cdot \mathbf{u})\mathbf{v} \Rightarrow \tilde{\mathbf{E}} = \eta(\mathbf{J} - q\mathbf{v}), \qquad (9) $$

where equations (5) and (7) have been used. Given the above spatial vectors, Maxwell's system is written in the usual way, a couple of evolutionary equations

$$ \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0, \quad \frac{\partial \mathbf{E}}{\partial t} - \nabla \times \mathbf{B} = -\mathbf{J}, \qquad (10) $$

and a couple of non-evolutionary constraints

$$ \nabla \cdot \mathbf{B} = 0, \quad \nabla \cdot \mathbf{E} = q, \qquad (11) $$

and we use the latter, Gauss's law, to define the charge density $q$.

The final system of Res-RMHD is retrieved by splitting the covariant equations (2), and by the first couple of Maxwell's equations (10). We obtain 11 evolution equations that may be written in the form of a set of conservation laws with source terms

$$ \frac{\partial U}{\partial t} = -\nabla \cdot \mathbf{F}(U) + S_{\text{e}}(U) + S(U), \qquad (12) $$
where

$$U = \begin{pmatrix} D \\ m_i \\ \mathcal{E} \\ B_i \\ E_i \end{pmatrix}, \quad F_j = \begin{pmatrix} \rho u_j \\ \rho h u_i u_j - E_i E_j - B_i B_j + p_{\text{tot}} \delta_{ij} \\ m_j \\ \varepsilon^{ijk} E_k \\ -\varepsilon^{ijk} B_k \end{pmatrix}, \quad (13)$$

with $\varepsilon^{ijk}$ is the three-dimensional Levi-Civita symbol. Here the mass density $D = \rho\gamma$, the total momentum density $m = Dh\mathbf{u} + \mathbf{E} \times \mathbf{B}$, and the total energy density $\mathcal{E} = Dh\gamma - p + u_{\text{em}}$ (all measured in the laboratory frame) are the so-called conserved variables, where $u_{\text{em}} = (E^2 + B^2)/2$ is the electromagnetic energy density, and $p_{\text{tot}} = p + u_{\text{em}}$ is the total pressure. These variables are expressed in terms of the primitive variables $\rho$, $v$, and $p$. Notice that the electric and magnetic fields $\mathbf{E}$ and $\mathbf{B}$ act as both conserved and primitive variables.

The source term is non-zero only in the equation for $\mathbf{E}$ and it contains the spatial current density $\mathbf{J}$ in equation (8). This is split, for computational purposes, into standard and stiff ($\propto 1/\eta$) contributions, so that the source terms $S_e$ and $S$ in equation (12) are given by

$$S_e = \begin{pmatrix} 0 \times 8 \\ -qv_i \end{pmatrix}, \quad S = \frac{1}{\eta} \begin{pmatrix} 0 \times 8 \\ -\tilde{E}_i \end{pmatrix}, \quad (14)$$

where $\tilde{E}$ is defined in equation (9). Numerical stiffness is obviously given by the fact that the resistivity coefficient $\eta$ is negligible in real astrophysical plasmas, and finite but very small in computational applications, so the term $S$ can be very large and special care should be adopted for the evolution of the electric field (the terms $R$ and $S$ will be treated differently in the numerical time evolution of $U$).

We recall that in ideal MHD, i.e. assuming the limit $\eta \to 0$, this latter equation is not evolved at all, the current does not need to be computed, and the electric field is a derived quantity, simply given by $\mathbf{E} = -\mathbf{v} \times \mathbf{B}$, so that $\tilde{E} = 0$.

# 3 NUMERICAL METHOD

## 3.1 Basic discretization

We employ a Cartesian coordinate system with axes represented by the unit vectors $\hat{e}_x = (1, 0, 0)$, $\hat{e}_y = (0, 1, 0)$, and $\hat{e}_z = (0, 0, 1)$, uniformly discretized into a regular mesh with coordinate spacing $\Delta x$, $\Delta y$, and $\Delta z$. We will frequently use the notation $c = (i, j, k)$ as a shorthand to identify a single computational cell centred at $(x_i, y_j, z_k)$ and delimited by the six interfaces orthogonal to the coordinate axis centred, respectively, at $(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}})$. The grid interfaces location will be shorthanded with

$$\begin{aligned} \mathbf{x}_f &\equiv \left(i + \frac{1}{2}, j, k\right), \\ \mathbf{y}_f &\equiv \left(i, j + \frac{1}{2}, k\right), \\ \mathbf{z}_f &\equiv \left(i, j, k + \frac{1}{2}\right). \end{aligned} \quad (15)$$

Likewise, we define the location of the cell edges as

$$\begin{aligned} \mathbf{x}_e &\equiv \left(i, j + \frac{1}{2}, k + \frac{1}{2}\right), \\ \mathbf{y}_e &\equiv \left(i + \frac{1}{2}, j, k + \frac{1}{2}\right), \\ \mathbf{z}_e &\equiv \left(i + \frac{1}{2}, j + \frac{1}{2}, k\right). \end{aligned} \quad (16)$$

Integrating the original partial differential equation (PDE; equation (12)) over a cell volume for $U = \{D, m, \mathcal{E}, \mathbf{E}\}$ gives the semidiscrete form

$$\frac{d\langle U \rangle_e}{dt} = \hat{R}_e + \langle S \rangle_e, \quad (17)$$

where $\langle U \rangle_e$ denotes the volume average of the corresponding conserved quantity, $\langle S \rangle_e$ is the stiff part of the source term (second term in equation (14)), while

$$\begin{aligned} \hat{R}_e = & -\left(\frac{\hat{F}_{x,x_f} - \hat{F}_{x,x_f} - \hat{e}_x}{\Delta x}\right) - \left(\frac{\hat{F}_{y,y_f} - \hat{F}_{y,y_f} - \hat{e}_y}{\Delta y}\right) \\ & -\left(\frac{\hat{F}_{z,z_f} - \hat{F}_{z,z_f} - \hat{e}_z}{\Delta z}\right) + \langle S_e \rangle_e \end{aligned} \quad (18)$$

is the explicit right-hand side including surface flux contributions plus the volume average of the explicit source term $S_e$ (first term in equation (14)) whose only non-vanishing component is the $qv$ term in the electric field update.

In equation (18), $\hat{F}_{x,x_f}$ (and likewise $\hat{F}_{y,y_f}$ and $\hat{F}_{z,z_f}$) represents the surface-integrated numerical flux across an $x$-face, i.e.

$$\hat{F}_{x,x_f} = \frac{1}{\Delta y \Delta z} \int F_{x,x_f} \, dy \, dz, \quad (19)$$

which has to be computed with a high-order quadrature rule. Expressions for the $y$- and $z$-directions are obtained in a similar fashion.

Magnetic field variables are treated as staggered quantities and are evolved using the integral form of Faraday's law (the first part in equation (10)) and direct application of Stokes' theorem (see sections 2, 3, and 4 of Paper II for more details). Thus, if $\hat{B}_{x,x_f}$, $\hat{B}_{y,y_f}$, and $\hat{B}_{z,z_f}$ are the surface-averaged components of $\mathbf{B}$ lying, respectively, on the $x$-, $y$-, and $z$-interfaces, the CT (see Yee 1966; Brecht et al. 1981; Evans & Hawley 1988; Balsara & Spicer 1999) algorithm entails the following discrete update:

$$\begin{aligned} \frac{d\hat{B}_{x,x_f}}{dt} &= -\left(\frac{\bar{E}_{z,z_c} - \bar{E}_{z,z_c} - \hat{e}_y}{\Delta y} - \frac{\bar{E}_{y,y_c} - \bar{E}_{y,y_c} - \hat{e}_z}{\Delta z}\right), \\ \frac{d\hat{B}_{y,y_f}}{dt} &= -\left(\frac{\bar{E}_{x,x_c} - \bar{E}_{x,x_c} - \hat{e}_z}{\Delta z} - \frac{\bar{E}_{z,z_c} - \bar{E}_{z,z_c} - \hat{e}_x}{\Delta x}\right), \\ \frac{d\hat{B}_{z,z_f}}{dt} &= -\left(\frac{\bar{E}_{y,y_c} - \bar{E}_{y,y_c} - \hat{e}_x}{\Delta x} - \frac{\bar{E}_{x,x_c} - \bar{E}_{x,x_c} - \hat{e}_y}{\Delta y}\right). \end{aligned} \quad (20)$$

In the previous expressions, we have denoted the line-averaged electric field (electromotive force or emf) at zone edges with $\bar{E}_{x,x_c}$, $\bar{E}_{y,y_c}$, and $\bar{E}_{z,z_c}$. The actual computation of the emf is detailed in Section 3.7.

Finally, it is worth to remind that in the presented fourth-order scheme only magnetic field components are evolved as surface average, while electric field (unlike Paper I) is treated as a zone-centred variable rather than in a staggered fashion.

## 3.2 Temporal update

Equations (17) and (20) are evolved in time using IMEX RK schemes, for stiffly accurate (SA) computations. When applied to equation (17), an IMEX RK scheme can be expressed following the general formulation of Ascher, Ruuth & Spiteri (1997) and Pareschi & Russo
(2005), yielding

$$ 
\begin{aligned}
\langle U \rangle_c^{(k)} &= \langle U \rangle_c^n + \Delta t \sum_{j=1}^{k-1} \tilde{a}_{kj} \hat{R}_c + \Delta t \sum_{j=1}^{k} a_{kj} \langle S \rangle_c^{(j)}, \\
\langle U \rangle_c^{n+1} &= \langle U \rangle_c^n + \Delta t \sum_{j=1}^{v} \tilde{w}_j \hat{R}_c^{(j)} + \Delta t \sum_{j=1}^{v} w_j \langle S \rangle_c^{(j)},
\end{aligned} 
\quad (21) $$

where $v$ is the number of integration stages, $\tilde{A}$ is a $v \times v$ lower triangular matrix ($\tilde{a}_{ij} = 0, \forall j \ge i$) enclosing the coefficients of the explicit right-hand side $\hat{R}_c$ (equation 18 or the right-hand side of equation 20), while $A = (a_{ij})$ is also a $v \times v$ matrix including the coefficients for the implicit part of the right-hand side, $S^{(j)}$. We recall that the use of a diagonally implicit RK solver guarantees an explicit evaluation of the non-stiff terms and implicit evaluation only for the diagonally stiff terms.

The final stage is always explicit and it is expressed in terms of the coefficient vectors $\tilde{w} = (\tilde{w}_1, \dots, \tilde{w}_v)^T$ and $w = (w_1, \dots, w_v)^T$. IMEX RK schemes are conveniently represented with a double tableau in the usual Butcher notation,

$$ \begin{array}{c|c} \tilde{c} & \tilde{c} \\ \hline \tilde{w}^T & w^T \end{array}, \quad (22) $$

where the coefficients $\tilde{c}$ and $c$ refer to the treatment of non-autonomous systems and are defined as

$$ \tilde{c}_i = \sum_{j=1}^{i-1} \tilde{a}_{ij}, \quad c_i = \sum_{j=1}^{i} a_{ij}. \quad (23) $$

In this paper, we consider two IMEX time-stepping methods: (i) the third-order accurate IMEX SSP3(4,3,3) scheme of Pareschi & Russo (2005) (with four implicit and three explicit stages) and (ii) the fourth-order IMEX ARK4(3)7L[2]SA1 of Kennedy & Carpenter (2019) (hereafter ARK4), which is L-stable, SA with seven stages. The explicit coefficients of both methods can be found in Appendix A.

From an algorithmic viewpoint, during the explicit-implicit update we found more convenient to separate the explicit contributions (acting on volume averages) from the implicit part, which, in our formulation, will be operated on pointwise values. For the generic $k$th stage in equation (21), we therefore act as follows:

$$ 
\begin{aligned}
\langle U \rangle_c^{(k*)} &= \langle U \rangle_c^n + \Delta t \sum_{j=1}^{k-1} (\tilde{a}_{kj} \hat{R}_c + a_{kj} \langle S \rangle_c^{(j)}), \\
U_c^{(k)} &= U_c^{(k*)} + \Delta t a_{kk} S_c^{(k)},
\end{aligned} 
\quad (24) $$

where the former step $(k*)$ is explicit, while the latter $(k)$ is implicit. The algorithm is summarized in Section 3.10.

### 3.3 Point value recovery of primitive variables

Our algorithm starts with the zone-averaged value of the solution $\langle U \rangle_c$ inside a cell and the face-averaged values of the staggered components of the magnetic field $\hat{B}$. At the very first step, we require the point value of the solution at the cell centre. Following Paper II, this operation may be carried out using the Laplacian operator,

$$ U_c = \left(1 - \frac{\Delta}{24}\right) \langle U \rangle_c + O(h^4), \quad (25) $$

where $\Delta$ is the Laplacian operator first introduced in the context of high-order methods by McCorquodale & Colella (2011), i.e.

$$ \Delta \langle U \rangle_c \equiv (\Delta^x + \Delta^y + \Delta^z) \langle U \rangle_c, \quad (26) $$

where e.g.

$$ \Delta^x \langle U \rangle_c = (\langle U \rangle_{c-\hat{e}_x} - 2 \langle U \rangle_c + \langle U \rangle_{c+\hat{e}_x}), \quad (27) $$

and similarly of the $y$- and $z$-directions. Note that this step applies only to cell-centred quantities, i.e. lab density, momentum, energy, and electric fields. For the magnetic field, instead, one has to first recover its point values at face centres using only the transverse Laplacian, so that

$$ 
\begin{aligned}
B_{x,x_f} &= \left(1 - \frac{\Delta^y + \Delta^z}{24}\right) \hat{B}_{x,x_f}, \\
B_{y,y_f} &= \left(1 - \frac{\Delta^z + \Delta^x}{24}\right) \hat{B}_{y,y_f}, \\
B_{z,z_f} &= \left(1 - \frac{\Delta^x + \Delta^y}{24}\right) \hat{B}_{z,z_f}.
\end{aligned} 
\quad (28) $$

Then, following Felker & Stone (2018), the pointwise value of the magnetic field at the cell centre is obtained using an unlimited high-order interpolation, i.e.

$$ 
\begin{aligned}
B_{x,c} &= -\frac{1}{16}(B_{x,x_f+\hat{e}_x} + B_{x,x_f-2\hat{e}_x}) + \frac{9}{16}(B_{x,x_f} + B_{x,x_f-\hat{e}_x}), \\
B_{y,c} &= -\frac{1}{16}(B_{y,y_f+\hat{e}_y} + B_{y,y_f-2\hat{e}_y}) + \frac{9}{16}(B_{y,y_f} + B_{y,y_f-\hat{e}_y}), \\
B_{z,c} &= -\frac{1}{16}(B_{z,z_f+\hat{e}_z} + B_{z,z_f-2\hat{e}_z}) + \frac{9}{16}(B_{z,z_f} + B_{z,z_f-\hat{e}_z}).
\end{aligned} 
\quad (29) $$

Now that all conservative variables are available at zone centres as point values, we convert them to primitive variables, i.e. $V_c = V(U_c)$. The conversion between conservative and primitive variables cannot be written in closed analytical form as it requires the inversion of a non-linear function. Differently from Paper I, here we prefer to follow an approach similar to that of Palenzuela et al. (2009), which 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. In our implementation, we employ the root solver of Mignone, Plewa & Bodo (2005, see also section 2 in Mignone & Bodo 2005) that can be extended to different EoSs.

Note that equation (25) may produce unphysical states in the presence of sharp gradients and in regions of rapidly changing profiles. A limiting strategy to identify troubled cells is indeed necessary before operating any conversion operation to maintain numerical stability. We discuss our fallback approach in Section 3.9 of this paper, while the reader may consult section 3.4 of Paper II for a full description.

### 3.4 Implicit update of pointwise quantities

During the implicit part of our update (see the second part in equation 24), we operate directly on the point value of the electric field vector:

$$ E = E^* - \frac{1}{\tilde{\eta}} [\gamma E + u \times B - (E \cdot u)v], \quad (30) $$

where, by analogy with equation (24), we identify $E^* \equiv E^{(k*)}$ (available from the last explicit stage) and $\tilde{\eta} = \eta/(a_{kk}\Delta t)$, which is now a non-dimensional quantity. The previous equation provides a relation between the electric field $E$ and the four-velocity $u$. Using

$$ E \cdot u = \frac{\tilde{\eta}\gamma}{1 + \tilde{\eta}\gamma} (E^* \cdot u), \quad (31) $$
we can rewrite equation (30) as an explicit function $E \equiv E(u)$:

$$ (\tilde{\eta} + \gamma)E(u) = \tilde{\eta} E^* - u \times B + \frac{\tilde{\eta}}{1 + \tilde{\eta}\gamma} (E^* \cdot u)u, \quad (32) $$

$$ \text{where } \gamma = \sqrt{1 + u^2}. $$

Note that the four-velocity is not a conservative variable and it must be derived from the set of conserved variables $U$. Since momentum does not change during the implicit step, we use the implicit relation

$$ f(u) = m - \left( Dh(u)u + E(u) \times B \right) = 0 \quad (33) $$

to derive $u$, as in Paper I. Note also that $D$ and $B$ are held constant during this step. Following Bucciantini & Del Zanna (2013), we adopt a Newton–Broyden root-finding method to derive the four-velocity (and the other primitive variables) from the set of conserved variables. Additionally, we refer the reader to Mattia et al. (2023) for a generalization of the Taub EoS (Mignone et al. 2005) and to Tomei et al. (2020) for a more comprehensive treatment of the dynamo term in Ohm's law along with for any space–time metric in general relativity. The iterative multidimensional Newton–Broyden algorithm can be sketched as follows:

(i) At the beginning of any step $\kappa = 0, 1, 2, \dots$, use equation (32) to obtain $E = E(u)$. When $\kappa = 0$, an initial guess for the four-velocity is needed: if $\tilde{\eta} > 1$, employ the value provided at the current time level, otherwise reckon that obtained from the ideal RMHD equations.

(ii) Given $u^{(\kappa)}$ and $E^{(\kappa)}$, compute $f^{(\kappa)} \equiv f(u^{(\kappa)})$ from the conservation of the total momentum density:

$$ f^{(\kappa)} = m - \left[ Dh(u)u + E(u) \times B \right]^{(\kappa)}. \quad (34) $$

Notice that the laboratory density, magnetic field, total energy, and conserved momentum remain constant during this iteration scheme. Equation (34) represents a set of three non-linear equations to be solved for $u^{(\kappa)}$. Thermodynamics quantities are also written as functions of the four-velocity; for an ideal gas law one finds

$$ h(u) = 1 + \Gamma_1 \frac{p\gamma}{D}, \quad (35) $$

with $\Gamma_1 = \Gamma/(\Gamma - 1)$ ($\Gamma = 4/3$ for relativistically hot plasmas), and, from the conservation of the total energy, the pressure is derived as

$$ p(u) = \frac{\mathcal{E} - D\gamma - (E^2 + B^2)/2}{\Gamma_1 \gamma^2 - 1}. \quad (36) $$

(iii) Using the Jacobian $J = \partial f / \partial u$, obtain an improved guess of the four-velocity through

$$ u^{(\kappa+1)} = u^{(\kappa)} - (J^{(\kappa)})^{-1} f^{(\kappa)}, \quad (37) $$

with the explicit form of the Jacobian being

$$ J_{ij} \equiv \frac{\partial f_i}{\partial u_j} = -Dh \delta_{ij} - D u_i \frac{\partial h}{\partial u_j} - \epsilon_{ilm} \frac{\partial E_l}{\partial u_j} B_m. \quad (38) $$

The derivatives of $h$ and $E$ are given by

$$ D \frac{\partial h}{\partial u^j} = - \frac{\Gamma_1}{\gamma^2 \Gamma_1 - 1} \left[ (\gamma h D + p) v_j + \gamma E_i \frac{\partial E_i}{\partial u_j} \right], \quad (39) $$

given the relation $\partial\gamma/\partial u_j = u_j/\gamma = v_j$, and, from equation (32), by

$$ (\tilde{\eta} + \gamma) \frac{\partial E_i}{\partial u_j} = -E_i v_j - \epsilon_{ijk} B_k + \frac{\tilde{\eta}}{1 + \tilde{\eta}\gamma} \left[ u_i E_j^* + (E^* \cdot u) \left( \delta_{ij} - \frac{\tilde{\eta}}{1 + \tilde{\eta}\gamma} u_i v_j \right) \right], \quad (40) $$

which can be also compared to the Jacobian presented in Tomei et al. (2020) (for a vanishing dynamo coefficient $\xi = 0$).

(iv) Exit the iteration cycle if the error is less than some prescribed accuracy $|f^{(\kappa)}| < \epsilon$, otherwise go back to step (i) and let $\kappa \to \kappa + 1$. For standard applications, our algorithm converges in 2–5 iterations given a relative tolerance of $10^{-11}$.

The above procedure is to be implemented at the beginning of all IMEX substeps at cell centres using pointwise quantities. Furthermore, it automatically operates the conversion from conservative to primitive variables, which are then readily available for reconstruction at cell interfaces and flux computation. For this reason, this procedure must be extended also to a layer of ghost zones (e.g. three for a five-point stencil reconstruction) in order to avoid extra communication at the end of the implicit step.

## 3.5 Point value reconstruction of primitive variables

Reconstruction at interfaces is operated directly on the function point values rather than on the one-dimensional average quantities as first introduced by Paper II. In the following paragraphs, we summarize the pointwise fifth-order weighted essentially non-oscillatory scheme of Borges et al. (WENO-Z; see Borges et al. 2008).

The interface value presented in equation (23) of Paper II at $x = x_{i+\frac{1}{2}}$ is computed as the convex combination of third-order accurate interface values built on the three possible three-point substencils $\{i-2, i-1, i\}$, $\{i-1, i, i+1\}$, and $\{i, i+1, i+2\}$, i.e.

$$ V_{i+\frac{1}{2}}^{L} = \omega_0 \frac{3V_{i-2} - 10V_{i-1} + 15V_i}{8} + \omega_1 \frac{-V_{i-1} + 6V_i + 3V_{i+1}}{8} + \omega_2 \frac{3V_i + 6V_{i+1} - V_{i+2}}{8}, \quad (41) $$

with $\{V_{i\pm 2}, V_{i\pm 1}, V_i\}$ point values. The weights $\omega_l$, for $l = \{0, 1, 2\}$, are defined as

$$ \omega_l = \frac{\alpha_l}{\sum_m \alpha_m}, \quad \alpha_l = d_l \left( 1 + \frac{|\beta_0 - \beta_2|}{\beta_l + \epsilon} \right). \quad (42) $$

Here $d_l = \{d_0 = 1/16, d_1 = 5/8, d_2 = 5/16\}$ denotes the optimal weights that reckon the fifth-order accurate approximation of equation (23) of Paper II, $\epsilon = 10^{-40}$ is a small number avoiding division by zero, and $\beta_l$ are the smoothness indicators (equation 26 of Paper II). A detailed description of pointwise reconstructions is presented in section 3 of Paper II.

## 3.6 Riemann solver and flux computation

The integrand in equation (19) is approximated by means of a Riemann solver that, following the traditional formalism of FV Godunov schemes, returns a stable and consistent numerical flux as a function of left/right (L/R) interface states obtained in the previous section. We will employ the Maxwell–Harten–Lax–van Leer contact (MHLLC) solver described in Paper I, which considers Maxwell’s and gas dynamics equations decoupled during this step. This leads to (i) a pair of outermost electromagnetic waves leading to jumps in $E$ and $B$ (and therefore total momentum and energy), and (ii) an inner wave fan including two hydrodynamics shocks enclosing a contact wave.

At the practical level, we first obtain the transverse components of $E$ and $B$ using Maxwell’s jump conditions. Specializing at an
x-interface:

$$
\begin{align}
\tilde{B}_y &= \frac{B_{y,L} + B_{y,R}}{2} + \frac{E_{z,R} - E_{z,L}}{2}, \nonumber \\
\tilde{B}_z &= \frac{B_{z,L} + B_{z,R}}{2} - \frac{E_{y,R} - E_{y,L}}{2}, \nonumber \\
\tilde{E}_y &= \frac{E_{y,L} + E_{y,R}}{2} - \frac{B_{z,R} - B_{z,L}}{2}, \nonumber \\
\tilde{E}_z &= \frac{E_{z,L} + E_{z,R}}{2} + \frac{B_{y,R} - B_{y,L}}{2}. \tag{43}
\end{align}
$$

Normal components ($B_x$ and $E_x$) remain continuous.

Next, we compute the hydrodynamics flux $\mathcal{H}^*$ using the Harten-Lax-van Leer contact (HLLC) relativistic solver of Mignone & Bodo (2005):

$$
\mathcal{H}^* = \begin{cases} \mathcal{H}_L & \text{if } \lambda_L > 0, \\ \mathcal{H}_L + \lambda_L(\mathcal{W}_L^* - \mathcal{W}_L) & \text{if } \lambda_L < 0 < \lambda^*, \\ \mathcal{H}_R + \lambda_R(\mathcal{W}_R^* - \mathcal{W}_R) & \text{if } \lambda^* < 0 < \lambda_R, \\ \mathcal{H}_R & \text{if } \lambda_R < 0, \end{cases} \quad (44)
$$

where $\mathcal{H}_S = (D, Q_k v_x + p \delta_{kx}, Q_y)_S$, $S = L, R$ denotes the left or right states, $k = x, y, z$ is the spatial component, $\mathcal{W}_S = (D, Q_k, \mathcal{E}_h)_S$, $\mathcal{E}_h = w\gamma^2 - p$ is the gas energy, and $Q_k = w\gamma^2 v_k$ is the $k$th component of the gas momentum. The wave speeds $\lambda_{L/R}$ are estimated from the fastest and the slowest speeds as in Mignone & Bodo (2005). The states in the star region enclosing the contact mode are obtained as

$$
\mathcal{W}_S^* = \frac{1}{\lambda_S - \lambda^*} \begin{pmatrix} D(\lambda - v_x) \\ Q_k(\lambda - v_x) + (p^* - p)\delta_{kx} \\ E(\lambda - v_x) + p^*\lambda^* - pv_x \end{pmatrix}_S \quad (45)
$$

for $S = L, R$ and the speed of the contact mode $\lambda_L^* = \lambda_R^* = \lambda^*$ is found from the negative branch of the quadratic equation

$$
\mathcal{H}_{[\mathcal{E}_h]}^{\text{hl}} (\lambda^*)^2 - (\mathcal{W}_{[\mathcal{E}_h]}^{\text{hl}} + \mathcal{H}_{[Q_x]}^{\text{hl}}) \lambda^* + \mathcal{W}_{[Q_x]}^{\text{hl}} = 0. \quad (46)
$$

Here the subscript [.] picks a specific component of the array, while $p^* = \mathcal{H}_{[Q_x]}^{\text{hl}} - \lambda^* \mathcal{H}_{[\mathcal{E}_h]}^{\text{hl}}$ is the pressure in the star region, and the superscript hll indicates the HLL-average state (see equations 73 and 74 in Paper I). Putting it all together, we retrieve the final flux:

$$
\mathcal{F}_{x,x_f} = \begin{pmatrix} \mathcal{H}^* \\ (0_{\times 6}) \end{pmatrix} + \begin{pmatrix} 0 \\ \hat{e}_x \cdot \tilde{T} \\ \hat{e}_x \cdot (\tilde{E} \times \tilde{B}) \\ \hat{e}_x \cdot \tilde{E} \\ -\hat{e}_x \times \tilde{B} \end{pmatrix} . \quad (47)
$$

A similar procedure is followed to calculate the y- and z-interface fluxes.

Finally, fourth-order accuracy at zone interfaces is ensured by computing the surface-averaged flux through

$$
\hat{F}_{x,x_f} = \left(1 + \frac{\Delta^y + \Delta^z}{24}\right) \mathcal{F}_{x,x_f}, \quad (48)
$$

where $\Delta^{y,z}$ have been defined in equations (26) and (27). Note that the fluxes are used to update density, momentum, energy, and electric fields. Magnetic fields, instead, are advanced in a separate step using the CT formalism (see Section 3.7).

3.7 Constrained transport update

During the staggered update, we evolve the face-averaged compo-
nents of magnetic fields using equation (20). In order to retain fourth-
order accuracy we follow the steps already illustrated in Paper II. In
what follows, we recap the procedure.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top view of the intersection between four neighbour zones in the xy plane. N, S, E, and W indicate the four cardinal directions with respect to the zone edge (here represented by the intersection between four neighbour zones), $\mathcal{R}_x(F_{y_f})$ and $\mathcal{R}_y(F_{x_f})$ are 1D reconstruction operators applied to each zone face, see equations (49) and (51).</figcaption></figure>

First, one needs to reconstruct the L/R transverse electric field
components from the four faces sharing the same edge to that edge.
With reference to Fig. 1, which shows a top view of a z-edge in the
xy plane, this entails a reconstruction of the four adjacent values
of $E_z$ (available at the x- and y-face centres) along the transverse
directions (y and x, respectively) to obtain four discontinuous values
coming up at $z_e \equiv (i + \frac{1}{2}, j + \frac{1}{2}, k)$:

$$
\begin{align}
E_{z,z_e}^S &= \mathcal{R}_y^+(E_{z,x_f}), & E_{z,z_e}^N &= \mathcal{R}_y^-(E_{z,x_f+\hat{e}_y}), \nonumber \\
E_{z,z_e}^W &= \mathcal{R}_x^+(E_{z,y_f}), & E_{z,z_e}^E &= \mathcal{R}_x^-(E_{z,y_f+\hat{e}_x}), \tag{49}
\end{align}
$$

where the four cardinal directions in parenthesis are referred to the
edge location $z_e$. For convenience, only one component per face has
to be reconstructed:<fnref n="1" />

$$
E_{z,x_f} = \frac{E_{z,x_f}^L + E_{z,x_f}^R}{2}. \quad (50)
$$

In 3D, one has to repeat the reconstruction (equation 49) for $E_{x,y_f}$ and $E_{x,z_f}$, respectively, along the z- and y-direction to $\mathbf{x}_e$ and also for $E_{y,x_f}$ and $E_{y,z_f}$ along the z- and x-direction to $\mathbf{y}_e$.

Next, we also reconstruct the face-centred point values of the
normal components of the magnetic field towards the same edge:

$$
\begin{align}
B_{x,z_e}^S &= \mathcal{R}_y^+(B_{x,x_f}), & B_{x,z_e}^N &= \mathcal{R}_y^-(B_{x,x_f+\hat{e}_y}), \nonumber \\
B_{y,z_e}^W &= \mathcal{R}_y^+(B_{y,y_f}), & B_{y,z_e}^E &= \mathcal{R}_y^-(B_{y,y_f+\hat{e}_x}). \tag{51}
\end{align}
$$

Equations (49) and (51) allow us to obtain the edge-centred point
value of the upwind electric field, using the CT-Maxwell averaging
scheme (see section 3.4.2 of Paper I):

$$
E_{z,z_e} = \left[ \frac{E_z^N + E_z^S + E_z^W + E_z^E}{4} + \frac{B_y^E - B_y^W}{2} - \frac{B_x^N - B_y^S}{2} \right]_{z_e} \quad (52)
$$

Similar expressions hold for $E_{x,x_e}$ and $E_{y,y_e}$ by suitable index permutation.

Finally, we obtain the emf (i.e. the averaged line integral of the
electric field needed in the discrete version of Faraday’s law equation

<footnote n="1">¹The approach followed here differs from the one employed in Paper I, where the reconstruction is performed from the cell centre using linear interpolation. Also, a typo is present in equation (48) of Paper I, since the arithmetic average operator (.)z is not needed.</footnote>
20) as

$$ \bar{E}_{z,z_c} = \left(1 + \frac{\Delta^z}{24}\right) E_{z,z_c}. \qquad (53) $$

Equations (51)–(53) are extended to the 3D case by suitable index permutation.

## 3.8 Explicit source term computation

In addition to the flux difference operator, the right-hand side (equation 18) requires an explicit contribution from the $qv$ source term during the electric field update. To fourth order, this term should be included as a volume average over the zone, i.e.

$$ \langle S \rangle_c = \left(1 + \frac{\Delta}{24}\right) q_c \mathbf{v}_c, \qquad (54) $$

where $q_c$ and $v_c$ are, respectively, the point values of the charge and the fluid velocity. The former can be obtained from a cell-centred discretization of $\nabla \cdot E$ yielding

$$ q_c = \delta_x E_x + \delta_y E_y + \delta_z E_z. \qquad (55) $$

Using the point value of the electric field already at disposal during the explicit step, we approximate (to fourth order) the $x$-contribution to the divergence as

$$ \delta_x E_x = \frac{8(E_{x,c+\hat{e}_x} - E_{x,c-\hat{e}_x}) - (E_{x,c+2\hat{e}_x} - E_{x,c-2\hat{e}_x})}{12\Delta x}. \qquad (56) $$

Contributions along the $y$- and $z$-directions are obtained in a similar way.

We point out that the employment of unlimited differences in equation (56) may lead to the appearance of spurious oscillations in the presence of strong gradients. This does not represent an issue for the tests presented below, but it could be further improved by including information from the reconstructed interface values (limited). We leave this investigation to a forthcoming work.

## 3.9 Order reduction at discontinuities

As pointed out in Paper II, fourth-order accuracy cannot be maintained in proximity of discontinuous solutions or steep gradients where the function may not be differentiable. In order to detect such features we compute, at the beginning of the RK step, the derivative ratio sensor as outlined in section 3.4 of Paper II. The detector flags computational cells where the order of the scheme can be locally reduced (this strategy is also known as ‘fallback’ process) by (i) not using the Laplacian during point value recovery and flux integration (equations 25 and 48) and (ii) switching to either third-order WENO or second-order linear reconstruction.

The detector employs a five-zone stencil and does not require additional ghost zones. A similar approach is also followed by Verma et al. (2019), albeit the authors make use of ‘global smoothness indicators’, defined in terms of weighted averages computed from the density and the three magnetic field components.

## 3.10 Algorithm summary

To better outline our algorithm, let $\mathcal{D}$ be the set of zones comprising the active computational domain (e.g. no ghost zones). Then $\mathcal{D} + 1$ denotes the set of active computational zones augmented by one layer of ghost zones. Since the last step of an IMEX RK scheme is always explicit, our implementation can be best summarized starting

from the end of the explicit stage ($k$*) (first part of equation 24) and moving on to the next, $k \to (k+1)^*$:

(i) Start with the volume average of conservative variable, $\langle U \rangle_c^{(k*)}$, defined on $\mathcal{D}$ at the end of the more recent explicit stage.

(ii) Assign boundary conditions on $\langle U \rangle_c$ in the first layer of ghost zones only so that $\langle U \rangle_c$ is correctly specified on $\mathcal{D} + 1$.

(iii) Recover the point value of conservative ($U_c$) and primitive variables, $V_c = \mathcal{V}(U_c)$ (as shown in Section 3.3), on $\mathcal{D}$ only.

(iv) Assign boundary conditions on $V_c$. This will define primitive variables on $\mathcal{D} + n_g$, where $n_g = 3$.

(v) Carry out the implicit step by updating $V_c$ on $\mathcal{D} + n_g$ as outlined in Section 3.4. Notice that conservative variables $U_c$ do not change during this step, except for the electric field components.

(vi) Using the point value of primitive variables, compute the implicit source term point value in $\mathcal{D} + 1$, i.e. $S_c = S(V_c)$.

(vii) Compute the zone-averaged implicit source term $\langle S \rangle_c$ on $\mathcal{D}$ using Simpson quadrature rule.

(viii) Achieve the explicit step by constructing the right-hand side $\hat{R}_c$ (equation 18) and its staggered version (right-hand side of equation 20):

(a) perform spatial reconstruction on primitive point value variables (Section 3.5);

(b) solve a Riemann problems with L/R states at zone interfaces, see Section 3.6;

(c) compute the surface-averaged flux using equation (48);

(d) compute the point value of the electric field at zone edges (equation 52) following the procedure outlined in Section 3.7;

(e) obtain the line-averaged emf through equation (53);

(f) compute the volume average of the explicit source term (Section 3.8).

(ix) Complete the explicit step by obtaining the volume-average conserved variables at the next RK stage, i.e. $\langle U \rangle_c^{(k+1)*}$.

Note that our formulation requires two boundary calls per stage: the first one is operated on $\langle U \rangle_c$ but it extends only to one layer of ghost zones, while the second one specifies boundary values on the point value or primitive variables $V_c$ and it extends to three ghost zones. Unlike previous works, our formulation does not require more than three ghost zones.

## 4 NUMERICAL BENCHMARKS

We now verify and assess the accuracy and robustness of our implementation through standard 2D and 3D reference solutions. Unless otherwise stated, we will employ the pointwise WENO-Z reconstruction (Paper II) and the MHLLC Riemann solver coupled with the CT-Maxwell emf average in all test problems. By default, order reduction (see Section 3.9) is disabled unless otherwise stated.

Initial conditions are integrated with a four-point Gaussian quadra-

For quantitative purposes, we evaluate $L_1$ norm errors for a generic quantity $Q$ against a reference solution $Q_{ref}$ as $\epsilon_1(Q) = \sum_i |Q(x_i) - Q_{ref}(x_i)|/N$, where $x_i$ is a generic point of the domain and $N$ is the total number of grid zones.

All computations have been carried out by means of the PLUTO astrophysical code (Mignone et al. 2007, 2012), where the fourth-order resistive method has been implemented. The Courant number is set to 0.4 (in 2D) and 0.3 for 3D calculations.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Numerical results of the telegraph equation test, for $\theta = 30^{\circ}$ (top panels) and $\theta = 90^{\circ}$ (bottom panels). Left column: 1D horizontal profiles of $B'_z$ using different values of the $\sigma = 1, 10,$ and $20$ for $N_x = 32$. Central and right columns: L1 norm errors of $B^* = B \cdot \hat{n}$ and $E^* = E \cdot (\hat{n} \times \hat{e}_x)$ as functions of the resolution. The dashed lines give the ideal convergence slope for the fourth order.</figcaption></figure>

## 4.1 Telegraph equation

In the fluid rest frame, the Res-RMHD equations can be employed to study the propagation of light waves in the presence of a finite conductivity, as originally presented in Mignone, Mattia & Bodo (2018) and Paper I. For a large plasma inertia ($\rho \to \infty$), only the Maxwell equations may be evolved in time:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0, \\ \\ \displaystyle \frac{\partial \mathbf{E}}{\partial t} - \nabla \times \mathbf{B} = -\sigma \mathbf{E}. \end{array} \right. \qquad (57) $$

Equation (57) admits damped propagating plane wave solutions when $\sigma < 2|k|$, where $k$ is the wavenumber. The exact mode in one dimension is

$$ \begin{aligned} \mathbf{B} &= B_1 e^{-\sigma t/2} \cos \phi(x, t) \hat{n}, \\ \mathbf{E} &= B_1 e^{-\sigma t/2} \left[ \frac{\mu}{k} \cos \phi(x, t) + \frac{\sigma}{2k} \sin \phi(x, t) \right] \hat{n} \times \hat{e}_x, \end{aligned} \qquad (58) $$

where $\phi(x, t) = kx - \mu t$, $\mu = (k^2 - \sigma^2/4)^{1/2}$, while $B_1$ is the initial perturbation amplitude. Here, in addition to Paper I, we have introduced $\hat{n} = (0, \cos \Theta, \sin \Theta)$ as the oscillation direction of the magnetic field in the yz plane.

Following Paper I, we rotate the 1D solution by an angle $\alpha$ around the z-axis and consider the 2D computational domain $x' \in [0, L_x]$, $y' \in [0, L_y]$ with $L_x = 1$, $L_y = 1/2$ discretized with $N_x \times N_x/2$ zones. Solving in a 2D domain allows us to assess the correct discretization of multidimensional terms. Electric and magnetic field vectors in this frame ($E'$ and $B'$) are obtained through the rotation

matrix

$$ \{E, B\}' = \begin{pmatrix} \cos \alpha - \sin \alpha & 0 \\ \sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{pmatrix} \{E, B\}. \qquad (59) $$

The wave vector has orientation $k = k_x(1, \tan \alpha, 0)$, where $k_x = 2\pi/L_x$ while $\tan \alpha = L_x/L_y = 2$. Equation (58) with $B_1 = 1$ and $\phi(x', 0) = k \cdot x'$ is used to initialize the electric and magnetic field. We also set density to a very large value ($\rho = 10^{12}$), while velocity are set to zero. Periodic boundary conditions are imposed everywhere.

We employ the fourth-order ARK4 IMEX scheme and run computations at different grid resolutions for exactly one wave period $T = 2\pi/\mu$, where $\mu$ is defined after equation (58). In Fig. 2, we plot, from left to right, the profiles of $B'_z$ (at the resolution $N_x = 32$), the L1-norm errors of $B \cdot \hat{n}$ and $E \cdot (\hat{n} \times \hat{e}_x)$ as a function of the grid resolution (vectors are transformed back to the 1D frame), for $\theta = 30^{\circ}$ (top panels) and $\theta = 90^{\circ}$ (bottom panels). Green, blue, and red colours refer to computations obtained with $\sigma = 1, 10,$ and $20$, respectively. In all cases, we recover the expected order of accuracy (4).

## 4.2 Stationary charged vortex

We now consider an exact equilibrium solution of the full Res-RMHD equations, originally introduced in Paper I. It consists of a rotating charged flow embedded in a vertical magnetic field. Assuming purely radial dependence, hydrodynamic and electromagnetic variables are best expressed in a cylindrical coordinate system
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Numerical results of the charged vortex problem. Top-left panel: We show eight equally spaced contour levels of $E$ chosen as $(q_0/2)r_k/(r_k^2 + 1)$ (where $r_k = 1, \dots, 8$) overlaid on the coloured maps of pressure at $t = 5$. Resistivity is set to $\eta = 10^{-3}$ and the grid resolution is $256^2$ zones. Top-right panel: 1D horizontal cuts, at $y = 0$, of charge at $t = 3$ obtained with CT. Bottom-left panel: $L_1$ norm errors of the charge as functions of the resistivity using CT with $256^2$ zones. Bottom-right panel: $L_1$ norm errors of the gas pressure as functions of the resolution for the SSP(4,3,3), ARK4, and the second-order scheme. The dotted lines give the expected second-order scaling.</figcaption></figure>

$(r, \phi, z)$ as

$$ egin{cases} E_r = \frac{q_0}{2} \frac{r}{r^2 + 1}, \\ B_z = \frac{\sqrt{(r^2 + 1)^2 - q_0^2/4}}{r^2 + 1}, \\ v_\phi = -\frac{q_0}{2} \frac{r}{\sqrt{(r^2 + 1)^2 - q_0^2/4}}, \\ p = -\frac{\rho}{\Gamma_1} + \left[ \frac{4r^2 + 4 - q_0^2}{(r^2 + 1)(4 - q_0^2)} \right]^{\Gamma_1/2} \left( p_0 + \frac{\rho}{\Gamma_1} \right), \\ q = \frac{q_0}{(r^2 + 1)^2}, \end{cases} \quad (60) $$

where $\Gamma_1 = \Gamma/(\Gamma - 1)$, $p_0 = 0.1$, $q_0 = 0.7$ is the charge at $r = 0$. Density is set to unity. The initial equilibrium condition given by equation (60) is mapped to a 2D Cartesian domain $(x, y \in [-10, 10])$, and the system is evolved until $t = 5$. Boundary conditions replicate the equilibrium solution throughout the integration. Since the equilibrium condition does not depend on the resistivity, numerical computations with different values of the resistivity $\eta$ depend solely on the stability of the algorithm.

For the present purpose, we have tested both the third-order SSP3(4,3,3) RK IMEX and the fourth-order ARK4 schemes. Results

are shown in the four panels of Fig. 3. The integrity of the solution is effectively maintained. This is demonstrated by both the pressure colour map (top left panel) showing also isocontour levels of the electric field magnitude) and the 1D horizontal cuts of the charge (top right).

In the bottom left panel of Fig. 3, we plot the charge error as a function of the resistivity and compare our results with those of Paper I using a second-order method. Our fourth-order scheme is able to reduce the error by two orders of magnitude and we have found that both the SSP(4,3,3) and the ARK4 yield nearly identical results. However, the ARK4 method becomes unstable for $\eta \lesssim 10^{-4}$. For this reason, we plot only the results obtained with the SSP method. Overall, the charge error exhibits minimal variation over time, underscoring the stability of the algorithm for any value of the resistivity parameter in the chosen range.

In the lower right quadrant of Fig. 3, we present the $L_1$ norm errors of gas pressure against resolution, showcasing the achievement of fourth-order accuracy for both methods. Notably, in comparison with the results presented in Paper I, both charge and pressure errors are markedly reduced, with a two-order improvement in convergence.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Snapshots of the gas pressure *p* (upper panels) and profiles of $E_z$ (lower panels) for the resistive rotor test at its final time *t* = 0.3 at 400 × 400 grid points using the SSP3(4,3,3) IMEX time-stepping. From left to right, results correspond to computations with different resistivities $\eta = 10^{-6}$, $\eta = 10^{-3}$, and $\eta = 10^{-1}$.</figcaption></figure>

## 4.3 2D relativistic rotor

We propose a resistive variant of the ideal special relativistic MHD test proposed by Del Zanna, Bucciantini & Londrillo (2003), as also presented by Dumbser & Zanotti (2009), Bucciantini & Del Zanna (2013), Miranda-Aranguren et al. (2018), and Nakamura et al. (2023). An overdense clump ($\rho = 10$) concentrated in a circular region ($r = 0.1$) at the centre of the computational square $[-0.5, 0.5]^2$ is rapidly spinning ($\omega = 8.5$) in a uniform medium ($\rho = 1, p = 1, \Gamma = 4/3$) at rest. The vortex is embedded in a constant magnetic field $\mathbf{B} = (1, 0, 0)$, which is progressively wrapping around the spinning vortex, launching torsional Alfvén waves in the ambient fluid. At the final time ($t = 0.3$), the torsional Alfvén waves have approximately been rotated at a 90° angle, with the magnetic pressure deforming the initially circular vortex into its characteristic oval shape.

We perform a set of 400 × 400 grid points simulations with increasing resistivity, from a quasi-ideal regime ($\eta = 10^{-6}$ and $10^{-3}$) to a fully resistive one ($\eta = 10^{-1}$). Both the SSP(4,3,3) and the ARK4 were used for this test using order reduction to WENO3. As the resistivity increases, the vortex angular momentum as well as the torsional Alfvén waves that cause the well-known deformation are dissipated, leading to a more circular shape. Fig. 4 presents the snapshots of thermal pressure *p* (top panels) and the electric field $E_z$ (bottom panels) at the final time *t* = 0.3. While results for SSP(4,3,3) appear essentially identical to those obtained with ARK4, only the former is able to successfully complete the integration with $\eta = 10^{-6}$, again confirming the reduced stability of the fourth-order method, which is not SSP.

It is relevant to notice how the pressure distribution is affected by the onset of spurious modes ($m = 4$) in the most resistive regime ($\eta = 10^{-1}$). This artefact results from the employment of Cartesian coordinates and it becomes more evident when using a MHLLC solver rather than a more diffusive one (with which the effect is substantially mitigated).

## 4.4 3D blast wave

The magnetized blast wave problem is used here to assess the robustness of the numerical method in handling the propagation of shock waves through magnetic field of varying strength as well as the capacity of preserving the symmetric properties of the solution. In fact, unphysical features could emerge if the divergence-free condition is not properly controlled (see e.g. Mignone & Bodo 2006; Del Zanna et al. 2007).

Following Paper I, we consider the computational domain $x, y, z \in [-6, 6]$ threaded by a uniform magnetic field $\mathbf{B} = B_0\hat{e}_z$ with $B_0 = 0.1$. The plasma is initially at rest while density and pressure experience a sharp transition across a sphere of radius $r_c = 0.8$:

$$ \begin{pmatrix} \rho \\ p \end{pmatrix} = f_i(r) \begin{pmatrix} \rho_{\text{in}} \\ p_{\text{in}} \end{pmatrix} + (1 \cdots f_i(r)) \begin{pmatrix} \rho_{\text{out}} \\ p_{\text{out}} \end{pmatrix}, \quad (61) $$

where $f_i(r)$ is a taper (smoothing) function,

$$ f_i(r) = \max \left[ \min \left( \frac{e^x - e}{1 - e}, 1 \right), 0 \right], \quad (62) $$

with $\chi = (r - r_c)/(1 - r_c)$. Similar 3D configurations have been used by Komissarov (2007) and Dionysopoulou et al. (2013), although with larger resistivity or lower magnetic fields. The system is evolved until $t = 4$ using the SSP(4,3,3) method. We also enable order reduction to be linear for this test. Zero-gradient boundary conditions apply everywhere.

Fig. 5 shows the outcomes of the simulation, presenting, in the top panel, 2D maps of the *x*-component of magnetic field (left), gas pressure (middle) while, in the bottom panel, electromagnetic pressure (left) and Lorentz factor (middle). In the rightmost column, we show one-dimensional profiles along the *z*-axis of thermal pressure (top) and Lorentz factor (bottom). As a consequence of the explosion, a swiftly advancing forward shock moves in an (almost) radial manner, while concurrently, a trailing reverse shock marks the boundary of the internal region where radial expansion
<figure><img src="image_6.png" /><figcaption>**Figure 5.** 3D blast wave at $t = 4$ using $\eta = 10^{-6}$ and the CT scheme. In the top panels, we show 2D slices in the $xy$ plane of $B_x$ (left), gas pressure $p_g$ (middle), and its 1D profile along the $z$-axis (right). In the bottom panels, we show 2D slices of the electromagnetic pressure (left), Lorentz factor (middle), and its cuts along the $z$-axis (right).</figcaption></figure>

occurs. Electromagnetic fields pile up in the $y$-direction, forming a shell characterized by elevated magnetic pressure. Gas exhibits a pronounced preference for motion along the $x$-direction, attaining a heightened Lorentz factor ($\gamma_{max} \approx 5.4$).

## 4.5 Tearing mode

In what follows we validate our method on the linear growth of the ideal tearing mode instability with a Harris current sheet profile, following the results presented by Del Zanna et al. (2016, see also Miranda-Aranguren et al. 2018; Paper I).

In this problem, a static ($v = 0$) plasma is uniformly distributed in the computational domain with constant initial density and pressure $\rho_0$ and $p_0$. The initial magnetic field profile satisfies the force-free condition and it is given by

$$ \mathbf{B} = B_0 \left[ \tanh\left(\frac{x}{a}\right) \hat{e}_y + \operatorname{sech}\left(\frac{x}{a}\right) \hat{e}_z \right], \quad (63) $$

while the initial electric field is null everywhere. The physical parameters that set the proper conditions for the instability to grow are the current sheet thickness $a$, 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}$ obtained by assuming an ideal gas law with relativistic adiabatic index $\Gamma = 4/3$, and the Lundquist number $S = v_a L/\eta$, where $\eta$ is the resistivity and $L$ is the reference spatial scale of variation. In order to trigger the instability, $S \gg 1$.

Generalizing the MHD results obtained by Pucci & Velli (2014) and Landi et al. (2015) to this regime, provided a sufficiently thin current sheet, the reconnection process is expected to occur on the ideal Alfvénic time-scale $\tau_a = L/v_a$. The critical (inverse) aspect ratio determining this regime is

$$ a = S^{-1/3} L. \quad (64) $$

In the limit of large $S$, the growth rate $\gamma_{TM}$ of the instability becomes $\gamma_{TM} \approx 0.6v_a/L$, which, in a relativistic regime ($v_a \to c = 1$), results in a very efficient process commonly referred to as ideal tearing instability.

The initial equilibrium is perturbed by a small amplitude ($\epsilon = 10^{-4}$) variation in the magnetic field employing only one wavenumber $k$:

$$ \delta \mathbf{B} = \epsilon B_0 \operatorname{sech} \left( \frac{x}{a} \right) \left( \frac{\cos(ky)}{ka} \tanh \left( \frac{x}{a} \right) \right). \quad (65) $$

Following section 5.5 of Paper I, we initialize the magnetic field in the $xy$ plane with the vector potential

$$ A_z = -B_0 \left[ a \log \left( \cosh \frac{x}{a} \right) - \frac{\epsilon}{k} \sin(ky) \operatorname{sech} \left( \frac{x}{a} \right) \right]. \quad (66) $$

In accordance with the work developed by Del Zanna et al. (2016) and later resumed in Paper I, the linear theory predicts that the wavenumber associated with the fastest growing mode of the ideal tearing instability is $k_{max} = 1.4S^{1/6} = 14$. The corresponding growth rate is $\gamma_{TM} \approx 0.3$. None the less, due to the neglect of the compressible terms in the linear approximation as presented in Del Zanna et al. (2016), the maximum value of the dispersion relation is shifted to $k = 12$ with the consequent modification of the expected growth rate to $\gamma_{TM} \approx 0.27$. We therefore align with the works of Del Zanna et al. (2016), Miranda-Aranguren et al. (2018), and Paper I and choose $k = 12$ for our simulations. Likewise, we set $L = B_0 = \Sigma = \beta = 1$, deriving therefore $\rho_0 = 1$, $p_0 = 0.5$, $v_a = 0.5$, and $S = 10^6$. In this way $a = 0.01$ and $\eta = 5 \times 10^{-7}$. In order to recover the predicted growth rate, we perform a set of simulations with increasing resolution $N_x \times N_x/4$ ($N_x = \{256, 512, 1024, 2048\}$), on a rectangular computational domain with $x \in [-20a, 20a] = [-0.2, 0.2]$ and outflow boundary conditions, and $y \in [0, 0.523598]$ imposing
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Horizontal cuts of the x-component of the magnetic field (left panel), and of the velocity (middle panel) for the tearing mode instability at *t* = 10. The profiles are normalized to unity and show the trends along *y* = 0 for the magnetic field (left panel), and along (3/4)L_y for the velocity (middle panel). Each resolution (*N*<sub>x</sub> = {256, 512, 1024, 2048}) is demarked with a different colour. Measured growth rates (right panel) are plotted as solid (dashed) lines for the fourth-order (second-order) data and are compared to the ideal trend (dashed-dotted line). Crosses mark the time intervals for which each growth rate has been calculated.</figcaption></figure>

<table><caption>**Table 1.** Growth rates of the linear phase of the tearing instability as measured from the second- and fourth-order simulation data, listed by resolution. The measured growth rates refer to the onset of the first linear phase of the instability. The time intervals chosen to evaluate each slope are marked as crosses in the right panel of Fig. 6.</caption><thead><tr><th>Resolution</th><th>Second order</th><th>Fourth order</th></tr></thead><tbody><tr><td>256 × 64</td><td>1.19</td><td>0.75</td></tr><tr><td>512 × 128</td><td>0.89</td><td>0.44</td></tr><tr><td>1024 × 256</td><td>0.58</td><td>0.30</td></tr><tr><td>2048 × 512</td><td>0.38</td><td>0.27</td></tr></tbody></table>

periodicity along *y*. Simulations have been carried out until *t* = 20 using the IMEX SSP(4,3,3) time-stepping.

In the left and middle panels of Fig. 6, we plot the horizontal cuts at *y* = 0 of the x-components of the magnetic field (left panel), velocity (middle panel) at *t* = 10 for all resolutions. The magnetic field at all resolutions closely resembles the eigenfunction solutions presented by Del Zanna et al. (2016) (see fig. 1 of that paper). The same holds for the velocity profile, although a more pronounced dependency on resolution is observed along with some clipping. Following Miranda-Aranguren et al. (2018) and Paper I, we measure the growth rate by first computing

$$ \gamma(t) = \frac{1}{2} \log \left( \int B_x^2(t) \, dS \right), \qquad (67) $$

which we plot, as a function of time, in the right panel of Fig. 6 (solid lines) along with the second-order data from Paper I (dotted lines) for each resolution. The growth rate is then computed from a linear fit of equation (67), while the red dashed-dotted line represents the ideal value ($\gamma_{TM} \simeq 0.27$). Table 1 compares the results obtained with the second- and fourth-order scheme (middle and right columns, respectively). For each slope, a different time window has been chosen in order to properly fit the linear growth of the fastest growing mode (crosses in Fig. 6) during the first linear phase. While both the second- and fourth-order schemes show systematic convergence towards the predicted value, the latter requires approximately half of the resolution (or less) to reproduce the results of the former. The best agreement is reached with the fourth-order scheme using 2048 zones.

## 4.6 3D Kelvin-Helmholtz instability

As a final application benchmark, we consider a 3D version of Kelvin-Helmholtz instability in relativistic MHD (see e.g. Buciantini & Del Zanna 2006) as presented in Paper I and earlier, in 2D, also by Beckwith & Stone (2011) and Mizuno (2013). At *t* = 0, we prescribe a double shear layer configuration with non-uniform density and constant magnetic field:

$$ \begin{align} \rho &= \frac{1 + \tanh \phi(y)}{2} \rho_h + \frac{1 - \tanh \phi(y)}{2} \rho_h, \notag \\ \mathbf{v} &= v_{sh} \tanh \phi(y) \hat{e}_x + A_0 \left[ \operatorname{sign}(y) \sin(2\pi x) \hat{e}_y + \zeta \hat{e}_z \right] f(y), \notag \\ \mathbf{B} &= \left( \sqrt{2p\sigma\mu_p}, 0, \sqrt{2p\sigma\mu_t} \right), \tag{68} \end{align} $$

where $v_{sh} = 1/2$ is the shear velocity, $a = 0.01$ is the shear thickness, $\phi(y) = (|y| - 1/2)/a$ while $A_0 = 0.1$ is the perturbation amplitude, and $\zeta$ is a random number in the range [-1, 1]. The function

$$ f(y) = e^{-\phi^2 a^2 / \alpha^2} \qquad (69) $$

(with $\alpha = 0.1$) is used to dampen the perturbations as we move away from the two interfaces located at $y = \pm 1/2$. The parameters $\mu_p$ and $\mu_t$ control the in plane ($\mu_p = 0.01$) and out of plane ($\mu_t = 1$) magnetic field strength. We set resistivity to $\eta = 1/\sigma = 10^{-5}$ and evolve the system until $t = 30$ using the MHLLC Riemann solver. Following Paper I, we perform simulations on the Cartesian domain $x, z \in [-1/2, 1/2]$ and $y \in [-1, 1]$ with uniform grid resolution $N \times 2N \times N$, with $N = 128$ (low resolution) and $N = 256$ (high resolution).

Simulations have been repeated with both the high-order scheme – employing IMEX SSP3(4,3,3) time-stepping – as well as using also the second-order scheme presented in Paper I, based on piecewise linear reconstruction and the second-order IMEX time-stepping. For the fourth-order method, we employ order reduction (linear).

The time evolution of the instability can be observed in Fig. 7 where we show 2D slices at $z = 0$ for density (first and second row from top) and poloidal to toroidal magnetic field (third and fourth) with $N = 256$. The large-scale motion remains confined along the initial shear direction and the turbulent state is characterized by sheet-like thin structures, where most of the magnetic field energy is trapped. These rounded slabs remain roughly parallel to the $xz$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 2D slice cut of the density and poloidal to toroidal magnetic field at $z = 0$ for the 3D Kelvin-Helmholtz instability at $t = 10$ and 30 at the resolution of 256 × 512 × 256 zones. Here $B_p = \sqrt{B_x^2 + B_y^2}$.</figcaption></figure>

plane, although significant medium-scale structures develop in this direction for $t \gtrsim 20$ around the shear layer.

A direct measure of the growth rate, quantified by $\Delta v_y = (\max(v_y) - \min(v_y))/2$, is shown in the top-left panel of Fig. 8 for both the second- and fourth-order simulations at low and high resolutions (blue and red colours, respectively). As noticed in Paper I (see fig. 11 of that paper), with the MHLLC Riemann solver, an increase in resolution (or scheme accuracy) is accompanied by a gradual reduction of the growth rate, until convergence is finally reached. This is also the case here and clearly visible with the second-order method (red lines) but only marginally with the fourth-order scheme that shows better convergence properties at low resolution. A quantitative measure of the growth rate, in the range given by the circle symbols in the figure, yields indeed $\approx 2.94$ and $\approx 2.14$ for the second-order scheme (low and high resolution, respectively) and $\approx 2.17$ and $\approx 2.31$ for the fourth-order scheme.

The enhanced dissipation of electromagnetic and kinetic energies (top middle and right panels in Fig. 8) using our high-order scheme stems from the increased turbulence that comes from having an effective dissipation at scales smaller than the low-order case. In this sense, the performance of our high-order scheme at low resolutions ($N_x = 128$) already outperforms the one obtained by the second-order scheme at twice the resolution.

The amount of small-scale features is an indirect measurement of the numerical dissipation, which we quantify through the spectral

energy density,

$$P(k_x) = \int |\mathcal{F}(k_x, y, z)|^2 dy dz, \qquad (70)$$

where $\mathcal{F}(k_x, y, z)$ is the 1D fast Fourier transform of the density taken across the $x$ dimension. $P(k_y)$ and $P(k_z)$ are constructed similarly, by corresponding index permutations. Results, shown in the bottom row of Fig. 8 ($P(k_y)$ and $P(k_z)$), indicate that large-scale modes are better resolved with the fourth-order scheme.

## 5 SUMMARY

We have presented a genuinely fourth-order FV algorithm for the solution of the Res-RMHD equations in multiple dimensions. To our knowledge, this is one of the few contributions addressing the design of the fourth-order methods based on an IMEX approach for multidimensional PDEs with stiff terms.

High-order quadrature rules to obtain volume, surface, or line averages are evaluated using appropriate Simpson rules, as outlined in the seminal paper by McCorquodale & Colella (2011). Point values are recovered from volume (or surface) averages following the opposite procedure. The fourth-order spatial accuracy is achieved by reconstructing point values (rather than 1D volume averages) using the modified WENO-Z reconstruction of Paper II, while for temporal integration, we have implemented and compared the third-order RK IMEX SSP3(4,3,3) scheme of Pareschi & Russo (2003, 2005) with
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panels: Growth rate (left) and volume-integrated electromagnetic $P_{em} = (B^2 + E^2)/2$ (middle) and kinetic $E_{kin} = \rho\gamma(\gamma - 1)$ (right) as a function of time for the 3D Kelvin-Helmholtz instability. Circles in the top-left plot mark the time windows used to compute the growth rate. Bottom panels: From left to right, spectral energy densities in the x-, y-, and z-direction at $t = 15$. Solid and dashed lines are used for the fourth- and second-order schemes, while different colors are used for the resolution of 128 × 256 × 128 and 256 × 512 × 256 zones, respectively.</figcaption></figure>

the fourth-order IMEX-ARK4(3)7L[2]SA_1 scheme of Kennedy & Carpenter (2019) (ARK). The divergence-free condition is ensured by evolving magnetic fields according to the CT formalism, while the electric field is treated as a volume average quantity. Equations are evolved in conservation form, requiring the solution of a Riemann problem between discontinuous left and right states at zone interfaces. Our scheme employs only three ghost zones, and thus a more compact stencil when compared to its previous version (Paper II).

Numerical benchmarks in two and three dimensions demonstrate the accuracy of the method along with reduced numerical dissipation. For smooth problems, our results show superior accuracy and faster convergence when compared to the standard second-order schemes. While the ARK4 retains the temporal fourth-order accuracy for smooth problems, we found in our experience that the third-order SSP scheme is more robust for lower values of the resistivity, $\eta \lesssim 10^{-5}$, and thus preferable for highly stiff problems.

The scheme can also cope with complex flow features including steep gradients and discontinuous waves, although the spatial order of the scheme may have to be locally reduced ('fallback' approach) so as to enhance robustness in critical regions. To identify such regions, we employ the discontinuity detector proposed in Paper II, based on the ratio of derivative, specifically designed to detect high-frequency oscillations typically arising in proximity of a discontinuous front.

The employment of genuinely fourth-order methods may represent a major upgrade in the field of astrophysical plasma computations, allowing more accurate computations at a reduced cost.

## ACKNOWLEDGEMENTS

Computational resources were provided by the Centro di Compe-
tenza sul Calcolo Scientifico (C3S) of the University of Torino
(c3s.unito.it).

This project has received funding from the European Union's Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101064953 (GR-PLUTO).

This work has received funding from the European High Per-
formance Computing Joint Undertaking (JU) and Belgium, Czech
Republic, France, Germany, Greece, Italy, Norway, and Spain under
grant agreement no. 101093441 (SPACE).

This paper was supported by the Fondazione ICSC, Spoke 3 Astrophysics and Cosmos Observations. National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) Project ID CN_00000013 ‘Italian Research Center on High-Performance Computing, Big Data and Quantum Computing’ funded by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e creazione di ‘campioni nazionali di R&S (M4C2-19)’ – Next Generation EU (NGEU).

The partial support from Gruppo Nazionale Calcolo Scientifico-Istituto Nazionale di Alta Matematica (GNCS INdAM) and MIUR-PRIN Project 2022, no. 2022KKJP4X ‘Advanced numerical methods for time-dependent parametric partial differential equations with applications’ is also acknowledged.
DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

Ascher U. M., Ruuth S. J., Spiteri R. J., 1997, *Appl. Numer. Math.*, **25**, 151

Balsara D. S., Spicer D. S., 1999, *J. Comput. Phys.*, **149**, 270

Beckwith K., Stone J. M., 2011, *ApJS*, **193**, 6

Berta V., Mignone A., Bugli M., Mattia G., 2024, *J. Comput. Phys.*, **499**, 112701 (Paper II)

Borges R., Carmona M., Costa B., Don W. S., 2008, *J. Comput. Phys.*, **227**, 3191

Boscarino S., Russo G., Semplice M., 2018, *Comput. Fluids*, **169**, 155

Boscheri W., Pareschi L., 2021, *J. Comput. Phys.*, **434**, 110206

Brecht S. H., Lyon J., Fedder J. A., Hain K., 1981, *Geophys. Res. Lett.*, **8**, 397

Bucciantini N., Del Zanna L., 2006, *A&A*, **454**, 393

Bucciantini N., Del Zanna L., 2013, *MNRAS*, **428**, 71

Bugli M., Del Zanna L., Bucciantini N., 2014, *MNRAS*, **440**, L41

Carpenter M. H., Kennedy C. A., Bijl H., Viken S. A., Vatsa V. N., 2005, *J. Sci. Comput.*, **25**, 157

Cheong P. C.-K., Pong D. Y. T., Yip A. K. L., Li T. G. F., 2022, *ApJS*, **261**, 22

Conde S., Gottlieb S., Grant Z. J., Shadid J. N., 2017, *J. Sci. Comput.*, **73**, 667

Dedner A., Kemm F., Kröner D., Munz C.-D., Schnitzer T., Wesenberg M., 2002, *J. Comput. Phys.*, **175**, 645

Del Zanna L., Bucciantini N., Londrillo P., 2003, *A&A*, **400**, 397

Del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, *A&A*, **473**, 11

Del Zanna L., Bugli M., Bucciantini N., 2014, in Pogorelov N. V., Audit E., Zank G. P., eds, ASP Conf. Ser. Vol. 488, 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013). Astron. Soc. Pac., San Francisco, p. 217

Del Zanna L., Papini E., Landi S., Bugli M., Bucciantini N., 2016, *MNRAS*, **460**, 3753

Dionysopoulou K., Alic D., Palenzuela C., Rezzolla L., Giacomazzo B., 2013, *Phys. Rev. D*, **88**, 044020

Dumbser M., Zanotti O., 2009, *J. Comput. Phys.*, **228**, 6991

Evans C. R., Hawley J. F., 1988, *ApJ*, **332**, 659

Felker K. G., Stone J. M., 2018, *J. Comput. Phys.*, **375**, 1365

Kennedy C. A., Carpenter M. H., 2019, *Appl. Numer. Math.*, **136**, 183

Komissarov S. S., 2007, *MNRAS*, **382**, 995

Landi S., Del Zanna L., Papini E., Pucci F., Velli M., 2015, *ApJ*, **806**, 131

Londrillo P., del Zanna L., 2004, *J. Comput. Phys.*, **195**, 17

McCorquodale P., Colella P., 2011, *Commun. Appl. Math. Comput. Sci.*, **6**, 1

Mattia G., Del Zanna L., Bugli M., Pavan A., Ciolfi R., Bodo G., Mignone A., 2023, *A&A*, **679**, A49

Mignone A., Bodo G., 2005, *MNRAS*, **364**, 126

Mignone A., Bodo G., 2006, *MNRAS*, **368**, 1040

Mignone A., Del Zanna L., 2021, *J. Comput. Phys.*, **424**, 109748

Mignone A., Plewa T., Bodo G., 2005, *ApJS*, **160**, 199

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, *ApJS*, **170**, 228

Mignone A., Zanni C., Tzeferacos P., van Straalen B., Colella P., Bodo G., 2012, *ApJS*, **198**, 7

Mignone A., Mattia G., Bodo G., 2018, *Phys. Plasmas*, **25**, 092114

Mignone A., Mattia G., Bodo G., Del Zanna L., 2019, *MNRAS*, **486**, 4252 (Paper I)

Miranda-Aranguren S., Aloy M. A., Rembiasz T., 2018, *MNRAS*, **476**, 3837

Mizuno Y., 2013, *ApJS*, **205**, 7

Munz C. D., Omnes P., Schneider R., Sonnendrücker E., Voß U., 2000, *J. Comput. Phys.*, **161**, 484

Nakamura K., Miyoshi T., Nonaka C., Takahashi H. R., 2023, *Eur. Phys. J. C*, **83**, 229

Palenzuela C., Lehner L., Reula O., Rezzolla L., 2009, *MNRAS*, **394**, 1727

Pareschi L., Russo G., 2003, in Hou T. Y., Tadmor E., eds, Hyperbolic Problems: Theory, Numerics, Applications. Springer, Berlin, p. 241

Pareschi L., Russo G., 2005, *J. Sci. Comput.*, **25**, 129

Pucci F., Velli M., 2014, *ApJ*, **780**, L19

Tomei N., Del Zanna L., Bugli M., Bucciantini N., 2020, *MNRAS*, **491**, 2346

Verma P. S., Teissier J.-M., Henze O., Müller W.-C., 2019, *MNRAS*, **482**, 416

Yee K., 1966, *IEEE Trans. Antennas Propagation*, **14**, 302
APPENDIX A: EXPLICIT BUTCHER TABLEAUX

The following tables list the tableaux for the explicit (left) and implicit (right) coefficients for the IMEX-SSP3(4,3,3) L-stable scheme (see also Pareschi & Russo 2003, 2005, table 6), and for the IMEX-ARK4(3)7L[2]SA₁ of Kennedy & Carpenter (2019). More specifically, the explicit form of Butcher's tableaux for the IMEX-SSP3(4,3,3) consists in

$$
\begin{array}{c|cccc}
  0 & 0 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 \\
  1 & 0 & 1 & 0 & 0 \\
  1/2 & 0 & 1/4 & 1/4 & 0 \\
  \hline
  0 & 1/6 & 1/6 & 2/3 & 0 \\
\end{array}
\quad
\begin{array}{c|cccc}
  \alpha & \alpha & 0 & 0 & 0 \\
  0 & -\alpha & \alpha & 0 & 0 \\
  1 & 0 & 1-\alpha & \alpha & 0 \\
  1/2 & \beta & \eta & 1/2-\beta-\eta-\alpha & \alpha \\
  \hline
  0 & 1/6 & 1/6 & 2/3 & 0 \\
\end{array}
\qquad \text{(A1)}
$$

where $\alpha = 0.24169426078821$, $\beta = 0.06042356519705$, and $\eta = 0.12915286960590$. Similarly, the full form of the Butcher's tableaux for the IMEX-ARK4(3)7L[2]SA₁ consists in

$$
\begin{array}{c|cccccc}
  0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  \tilde{c}_2 & \tilde{a}_{21} & 0 & 0 & 0 & 0 & 0 \\
  \tilde{c}_3 & \tilde{a}_{31} & \tilde{a}_{32} & 0 & 0 & 0 & 0 \\
  \tilde{c}_4 & \tilde{a}_{41} & \tilde{a}_{42} & \tilde{a}_{43} & 0 & 0 & 0 \\
  \tilde{c}_5 & \tilde{a}_{51} & \tilde{a}_{52} & \tilde{a}_{53} & \tilde{a}_{54} & 0 & 0 \\
  \tilde{c}_6 & \tilde{a}_{61} & \tilde{a}_{62} & \tilde{a}_{63} & \tilde{a}_{64} & \tilde{a}_{65} & 0 \\
  \tilde{c}_7 & \tilde{a}_{71} & \tilde{a}_{72} & \tilde{a}_{73} & \tilde{a}_{74} & \tilde{a}_{75} & \tilde{a}_{76} \\
  \hline
  0 & 0 & \tilde{w}_3 & \tilde{w}_4 & \tilde{w}_5 & \tilde{w}_6 & \kappa \\
\end{array}
\quad
\begin{array}{c|cccccc}
  0 & 0 & 0 & 0 & 0 & 0 & 0 \\
  c_2 & \kappa & \kappa & 0 & 0 & 0 & 0 \\
  c_3 & a_{31} & a_{32} & \kappa & 0 & 0 & 0 \\
  c_4 & a_{41} & a_{42} & a_{43} & \kappa & 0 & 0 \\
  c_5 & a_{51} & a_{52} & a_{53} & a_{54} & \kappa & 0 \\
  c_6 & a_{61} & a_{62} & a_{63} & a_{64} & a_{65} & \kappa \\
  c_7 & a_{71} & a_{72} & a_{73} & a_{74} & a_{75} & a_{76} \\
  \hline
  0 & 0 & w_3 & w_4 & w_5 & w_6 & \kappa \\
\end{array}
\qquad \text{(A2)}
$$

where $\kappa = 0.1235$. The remaining 21 explicit ($\tilde{a}$) and the 20 implicit coefficients (a) are, respectively,

$$
\tilde{a}_{21} = \frac{247}{1000} = 0.247 \qquad \text{(A3)}
$$

$$
\tilde{a}_{31} = \frac{247}{4000} = 0.06175 \qquad \tilde{a}_{32} = \frac{2694949928731}{7487940209513} = 0.35990537495
$$

$$
\tilde{a}_{41} = \frac{464650059369}{8764239774964} = 0.05301658458 \qquad \tilde{a}_{42} = \frac{878889893998}{2444806327765} = 0.35949264529
$$

$$
\tilde{a}_{43} = -\frac{952945855348}{12294611323341} = -0.07750922988
$$

$$
\tilde{a}_{51} = \frac{476636172619}{8159180917465} = 0.05841715944 \qquad \tilde{a}_{52} = -\frac{1271469283451}{7793814740893} = -0.16313824817
$$

$$
\tilde{a}_{53} = -\frac{859560642026}{4356155882851} = -0.19732090979 \qquad \tilde{a}_{54} = \frac{1723805262919}{4571918432560} = 0.37704199852
$$

$$
\tilde{a}_{61} = \frac{6338158500785}{11769362343261} = 0.5385303227 \qquad \tilde{a}_{62} = -\frac{4970555480458}{10924838743837} = -0.45497746895
$$

$$
\tilde{a}_{63} = \frac{3326578051521}{2647936831840} = 1.25629056234 \qquad \tilde{a}_{64} = -\frac{880713585975}{1841400956686} = -0.47828452721
$$

$$
\tilde{a}_{65} = -\frac{1428733748635}{8843423958496} = -0.16155888888
$$

$$
\tilde{a}_{71} = \frac{760814592956}{3276306540349} = 0.23221715782 \qquad \tilde{a}_{72} = \frac{760814592956}{3276306540349} = 0.23221715782
$$

$$
\tilde{a}_{73} = -\frac{47223648122716}{6934462133451} = -6.80999437504 \qquad \tilde{a}_{74} = \frac{71187472546993}{9669769126921} = 7.36185855242
$$

$$
\tilde{a}_{75} = -\frac{13330509492149}{9695768672337} = -1.37487907794 \qquad \tilde{a}_{76} = \frac{11565764226357}{8513123442827} = 1.35858058491
$$
and

$$a_{31} = \frac{624185399699}{4186980696204} = 0.14907768747 \qquad a_{32} = \frac{624185399699}{4186980696204} = 0.14907768747 \tag{A4}$$

$$a_{41} = \frac{1258591069120}{10082082980243} = 0.12483442871 \qquad a_{42} = \frac{1258591069120}{10082082980243} = 0.12483442871$$

$$a_{43} = -\frac{322722984531}{8455138723562} = -0.03816885743$$

$$a_{51} = -\frac{436103496990}{5971407786587} = -0.0730319403 \qquad a_{52} = -\frac{436103496990}{5971407786587} = -0.0730319403$$

$$a_{53} = -\frac{2689175662187}{11046760208243} = -0.24343568716 \qquad a_{54} = \frac{4431412449334}{12995360898505} = 0.34099956776$$

$$a_{61} = -\frac{2207373168298}{14430576638973} = -0.15296500088 \qquad a_{62} = -\frac{2207373168298}{14430576638973} = -0.15296500088$$

$$a_{63} = \frac{242511121179}{3358618340039} = 0.07220562047 \qquad a_{64} = \frac{3145666661981}{7780404714551} = 0.40430630248$$

$$a_{65} = \frac{5882073923981}{14490790706663} = 0.4059180788$$

$$a_{71} = 0 \qquad a_{72} = 0$$

$$a_{73} = \frac{9164257142617}{17756377923965} = 0.51611072831 \qquad a_{74} = -\frac{10812980402763}{74029279521829} = -0.14606356393$$

$$a_{75} = \frac{1335994250573}{5691609445217} = 0.23473048589 \qquad a_{76} = \frac{2273837961795}{8368240463276} = 0.27172234973$$

where, for the sake of simplicity, the horizontal lines are used to highlight the different stages of the algorithm. From equation (23), the explicit form of the coefficients $\tilde{c} = c$ becomes

$$\begin{align}
\tilde{c}_2 &= \frac{247}{1000} = 0.247, & \tilde{c}_3 &= \frac{4276536705230}{10142255878289} = 0.42165537495, \tag{A5} \\
\tilde{c}_4 &= \frac{67}{200} = 0.335, & \tilde{c}_5 &= \frac{7}{40} = 0.075,
\end{align}$$

$$\tilde{c}_6 = \frac{7}{10} = 0.7, \qquad \tilde{c}_7 = 1. \qquad \text{The Butcher's tableaux are finally completed by the coefficients } w = \tilde{w}:$$

$$\begin{align}
\tilde{w}_3 &= \frac{9164257142617}{17756377923965} = 0.51611072831, & \tilde{w}_4 &= -\frac{10812980402763}{74029279521829} = -0.14606356393, \notag \\
\tilde{w}_5 &= \frac{1335994250573}{5691609445217} = 0.23473048589, & \tilde{w}_6 &= \frac{2273837961795}{8368240463276} = 0.27172234973. \tag{A6}
\end{align}$$

This paper has been typeset from a TeX/LA TeX file prepared by the author.