# A 4^th^-order accurate finite volume method for ideal classical and special relativistic MHD based on pointwise reconstructions

Vittoria Berta, Andrea Mignone, Matteo Bugli, Giancarlo Mattia

► To cite this version:

Vittoria Berta, Andrea Mignone, Matteo Bugli, Giancarlo Mattia. A 4^th^-order accurate finite volume method for ideal classical and special relativistic MHD based on pointwise reconstructions. J.Comput.Phys., 2024, 499, pp.112701. 10.1016/j.jcp.2023.112701. hal-04267752

HAL Id: hal-04267752

https://hal.science/hal-04267752

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
<figure><img src="image_1.png" /></figure>

# A 4^th^-order accurate finite volume method for ideal classical and special relativistic MHD based on pointwise reconstructions

V. Berta<sup>a,*</sup>, A. Mignone<sup>a</sup>, M. Bugli<sup>a,b</sup>, G. Mattia<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Magnetohydrodynamics (MHD)

High-order finite volume methods

Constrained transport

Riemann solvers

## ABSTRACT

We present a novel implementation of a genuinely 4^th^-order accurate finite volume scheme for multidimensional classical and special relativistic magnetohydrodynamics (MHD) based on the constrained transport (CT) formalism. The scheme introduces several novel aspects when compared to its predecessors yielding a more efficient computational tool. Among the most relevant ones, our scheme exploits pointwise to pointwise reconstructions (rather than one- dimensional finite volume ones), employs the generic upwind constrained transport averaging and sophisticated limiting strategies that include both a discontinuity detector and an order reduction procedure. Selected numerical benchmarks demonstrate the accuracy and robustness of the method.

## 1. Introduction

The modeling of astrophysical plasmas is nowadays demanding for more efficient and accurate numerical schemes for solving the equations of classical and special relativistic magnetohydrodynamics (MHD, RMHD). Physical interest is increasingly directed towards the investigation of highly nonlinear flows featuring both smooth and discontinuous solutions, with a great deal of attention being devoted to enhance the accuracy of the underlying numerical methods while retaining stability and computational efficiency.

Traditionally, in the context of finite volume (FV) or finite difference (FD) schemes, a great variety of 2^nd^-order methods have been developed to solve the MHD and RMHD equations (e.g., [1–18]). However, in the last decade several efforts have been made to design higher than 2^nd^-order numerical methods both for FD (e.g., [19,20,15,21–24]) as well as for FV (e.g., [25–32]) approaches, leading to the dawn of high-order numerical codes (e.g. [33–38]) in plasma physics and computational fluid-dynamics. Most 4^th^- order schemes, in fact, can reach unprecedented performances overstepping the limits of traditional 2^nd^-order frameworks because of their intrinsic lower dissipation properties. In addition, smooth solutions are improved at much faster rates, yielding enhanced convergence and henceforth a computational efficiency that increases with the dimensionality of the problem and with resolution.

In this paper we present a genuine 4^th^-order accurate finite volume method for the numerical solution of the ideal classical and special relativistic MHD equations. Such method draws on the work originally introduced by McCorquodale & Colella [28] and then refined by Felker & Stone [36], but with a number of important improvements and differences. In fact, our numerical scheme is introducing, for the first time in the literature of FV schemes, a version of the MP5 [39] and the WENOZ [40–42] spatial

<footnote><sup>a</sup> Dipartimento di Fisica, Università di Torino, via Pietro Giuria 1, 10125 Torino, Italy</footnote>

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

<footnote><sup>c</sup> INFN, Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: vittoria.berta@unito.it, vittoria.berta@edu.unito.it (V. Berta).</footnote>
<figure><img src="image_2.png" /></figure>

# A 4^th^-order accurate finite volume method for ideal classical and special relativistic MHD based on pointwise reconstructions

V. Berta<sup>a,*</sup>, A. Mignone<sup>a</sup>, M. Bugli<sup>a,b</sup>, G. Mattia<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Magnetohydrodynamics (MHD)

High-order finite volume methods

Constrained transport

Riemann solvers

## ABSTRACT

We present a novel implementation of a genuinely 4^th^-order accurate finite volume scheme for multidimensional classical and special relativistic magnetohydrodynamics (MHD) based on the constrained transport (CT) formalism. The scheme introduces several novel aspects when compared to its predecessors yielding a more efficient computational tool. Among the most relevant ones, our scheme exploits pointwise to pointwise reconstructions (rather than one- dimensional finite volume ones), employs the generic upwind constrained transport averaging and sophisticated limiting strategies that include both a discontinuity detector and an order reduction procedure. Selected numerical benchmarks demonstrate the accuracy and robustness of the method.

## 1. Introduction

The modeling of astrophysical plasmas is nowadays demanding for more efficient and accurate numerical schemes for solving the equations of classical and special relativistic magnetohydrodynamics (MHD, RMHD). Physical interest is increasingly directed towards the investigation of highly nonlinear flows featuring both smooth and discontinuous solutions, with a great deal of attention being devoted to enhance the accuracy of the underlying numerical methods while retaining stability and computational efficiency.

Traditionally, in the context of finite volume (FV) or finite difference (FD) schemes, a great variety of 2^nd^-order methods have been developed to solve the MHD and RMHD equations (e.g., [1–18]). However, in the last decade several efforts have been made to design higher than 2^nd^-order numerical methods both for FD (e.g., [19,20,15,21–24]) as well as for FV (e.g., [25–32]) approaches, leading to the dawn of high-order numerical codes (e.g. [33–38]) in plasma physics and computational fluid-dynamics. Most 4^th^- order schemes, in fact, can reach unprecedented performances overstepping the limits of traditional 2^nd^-order frameworks because of their intrinsic lower dissipation properties. In addition, smooth solutions are improved at much faster rates, yielding enhanced convergence and henceforth a computational efficiency that increases with the dimensionality of the problem and with resolution.

In this paper we present a genuine 4^th^-order accurate finite volume method for the numerical solution of the ideal classical and special relativistic MHD equations. Such method draws on the work originally introduced by McCorquodale & Colella [28] and then refined by Felker & Stone [36], but with a number of important improvements and differences. In fact, our numerical scheme is introducing, for the first time in the literature of FV schemes, a version of the MP5 [39] and the WENOZ [40–42] spatial

<footnote><sup>a</sup> Dipartimento di Fisica, Università di Torino, via Pietro Giuria 1, 10125 Torino, Italy</footnote>

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

<footnote><sup>c</sup> INFN, Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: vittoria.berta@unito.it, vittoria.berta@edu.unito.it (V. Berta).</footnote>
reconstruction algorithms based on pointwise values, which are instead traditionally used in FD schemes. As we shall demonstrate, pointwise reconstructions lower the algorithmic complexity of the overall scheme by retaining the 4th-order accuracy with fewer operations per step and without even enlarging the stencil size. Moreover, this new feature reduces the number of Riemann problems to be solved per direction, leading to a more efficient and cost-effective scheme. Concurrently, robustness is ensured by means of a local discontinuity detector to distinguish smooth solutions from discontinuous ones when determining the function local point value from its cell average. Regarding temporal integration, we employ a semi-discrete approach based on a five-stage 4th-order explicit strong stability preserving Runge-Kutta method (eSSPRK(5,4), [43,44]). The solenoidal condition of the magnetic field is fulfilled by means of a high-order formulation of the generalized upwind constrained transport (UCT, [45,3,5,46,47]) algorithm, which ensures a divergence-free magnetic field up to machine accuracy, and works well with any generic upwind average [18].

The paper is structured as follows: in §2 we will briefly describe the MHD and RMHD equations and how they are discretized in a general FV framework together with the notation used throughout this work. In §3 an extensive description of the novel numerical method will be provided. In §4 we will illustrate our method to control the solenoidal condition. Numerical tests and conclusions will be given, respectively, in §5 and §6.

## 2. The ideal MHD and RMHD equations

The ideal MHD equations can be split out in two coupled sub-systems [19,47,18], the first one being a time-dependent hyperbolic system for the evolution of the conservative flow variables (the mass density $ρ$, the momentum density $ρv$, and the energy density $ε$). This sub-system can be written in compact form by introducing a state vector of conservative quantities, $U$, and a rank-2 tensor, $F$, whose rows represent the fluxes relative to each component of $U$

$$ \frac{\partial U}{\partial t} + \nabla \cdot F = 0, \qquad (1) $$

where

$$ U = \begin{pmatrix} \rho \\ \rho v \\ \epsilon \end{pmatrix}, \quad F = \begin{pmatrix} \rho v \\ \rho v v - BB + I p_t \\ (\epsilon + p_t)v - (v \cdot B)B \end{pmatrix}^T. \qquad (2) $$

In Eq. (2), $v = (v_x, v_y, v_z)$ is the fluid velocity, $B = (B_x, B_y, B_z)$ is the magnetic field, $p_t = p + B^2/2$ is the total pressure expressed as the sum of the thermal and magnetic pressure, and $I$ is the identity matrix. The total energy density is composed by a kinetic, a thermal and a magnetic term

$$ \epsilon = \frac{1}{2} \rho v^2 + \frac{p}{\Gamma - 1} + \frac{B^2}{2}, \qquad (3) $$

while $Γ$ is the specific heat ratio for an adiabatic equation of state. On the other hand, the second sub-system consists of the induction equation for the evolution of the magnetic field

$$ \frac{\partial \mathbf{B}}{\partial t} + c \nabla \times \mathbf{E} = 0. \qquad (4) $$

Due to the presence of the curl operator instead of the divergence, the latter equation cannot be considered a simple extension of the HD subsystem [19]. Additionally, in the ideal MHD framework, the electric field $E$ is a function of the fluid velocity $v$ and the magnetic field $B$

$$ c\mathbf{E} = -\mathbf{v} \times \mathbf{B}, \qquad (5) $$

where $c$ is the speed of light.

Analytically, the commutativity of spatial derivatives endows the induction equation with the solenoidal condition for the magnetic field

$$ \nabla \cdot \mathbf{B} = 0. \qquad (6) $$

This stationary condition, once satisfied at $t=0$, is analytically preserved at $t>0$. On the other hand, since numerical derivatives do not commute, the condition expressed in Eq. (6) is not automatically fulfilled by numerical schemes, leading to the rising of magnetic monopoles which cause unphysical plasma transport and dissipation of momentum and energy in ideal frameworks [48]. For these reasons, HD numerical schemes cannot be straightforwardly extended to MHD systems without introducing some corrections to numerically constraint Eq. (6) to be satisfied.

Similarly to the MHD case, the ideal RMHD equations can be expressed in the form given by Eq. (1) with

$$ U = \begin{pmatrix} D \\ \mathbf{m} \\ \epsilon \end{pmatrix}, \quad F = \begin{pmatrix} Dv \\ w_t \gamma^2 \mathbf{v} \mathbf{v} - \mathbf{b} \mathbf{b} + I p_t \\ \mathbf{m} \end{pmatrix}^T, \qquad (7) $$

where $D = \gamma \rho$, $m = (\rho h \gamma^2 + |B|^2)v - (v \cdot B)B$, and $\epsilon = \rho h \gamma^2 - p + |B|^2/2 + (|v|^2|B|^2 - (v \cdot B)^2)/2$ are respectively the relativistic mass density, momentum density, and energy density, while $w_t = \rho h + |B|^2/\gamma^2 + (v \cdot B)^2$ and $b = B/\gamma + \gamma(v \cdot B)v$. In the previous
reconstruction algorithms based on pointwise values, which are instead traditionally used in FD schemes. As we shall demonstrate, pointwise reconstructions lower the algorithmic complexity of the overall scheme by retaining the 4th-order accuracy with fewer operations per step and without even enlarging the stencil size. Moreover, this new feature reduces the number of Riemann problems to be solved per direction, leading to a more efficient and cost-effective scheme. Concurrently, robustness is ensured by means of a local discontinuity detector to distinguish smooth solutions from discontinuous ones when determining the function local point value from its cell average. Regarding temporal integration, we employ a semi-discrete approach based on a five-stage 4th-order explicit strong stability preserving Runge-Kutta method (eSSPRK(5,4), [43,44]). The solenoidal condition of the magnetic field is fulfilled by means of a high-order formulation of the generalized upwind constrained transport (UCT, [45,3,5,46,47]) algorithm, which ensures a divergence-free magnetic field up to machine accuracy, and works well with any generic upwind average [18].

The paper is structured as follows: in §2 we will briefly describe the MHD and RMHD equations and how they are discretized in a general FV framework together with the notation used throughout this work. In §3 an extensive description of the novel numerical method will be provided. In §4 we will illustrate our method to control the solenoidal condition. Numerical tests and conclusions will be given, respectively, in §5 and §6.

## 2. The ideal MHD and RMHD equations

The ideal MHD equations can be split out in two coupled sub-systems [19,47,18], the first one being a time-dependent hyperbolic system for the evolution of the conservative flow variables (the mass density $ρ$, the momentum density $ρv$, and the energy density $ε$). This sub-system can be written in compact form by introducing a state vector of conservative quantities, $U$, and a rank-2 tensor, $F$, whose rows represent the fluxes relative to each component of $U$

$$ \frac{\partial U}{\partial t} + \nabla \cdot F = 0, \qquad (1) $$

where

$$ U = \begin{pmatrix} \rho \\ \rho v \\ \epsilon \end{pmatrix}, \quad F = \begin{pmatrix} \rho v \\ \rho v v - BB + I p_t \\ (\epsilon + p_t)v - (v \cdot B)B \end{pmatrix}^T. \qquad (2) $$

In Eq. (2), $v = (v_x, v_y, v_z)$ is the fluid velocity, $B = (B_x, B_y, B_z)$ is the magnetic field, $p_t = p + B^2/2$ is the total pressure expressed as the sum of the thermal and magnetic pressure, and $I$ is the identity matrix. The total energy density is composed by a kinetic, a thermal and a magnetic term

$$ \epsilon = \frac{1}{2} \rho v^2 + \frac{p}{\Gamma - 1} + \frac{B^2}{2}, \qquad (3) $$

while $Γ$ is the specific heat ratio for an adiabatic equation of state. On the other hand, the second sub-system consists of the induction equation for the evolution of the magnetic field

$$ \frac{\partial \mathbf{B}}{\partial t} + c \nabla \times \mathbf{E} = 0. \qquad (4) $$

Due to the presence of the curl operator instead of the divergence, the latter equation cannot be considered a simple extension of the HD subsystem [19]. Additionally, in the ideal MHD framework, the electric field $E$ is a function of the fluid velocity $v$ and the magnetic field $B$

$$ c\mathbf{E} = -\mathbf{v} \times \mathbf{B}, \qquad (5) $$

where $c$ is the speed of light.

Analytically, the commutativity of spatial derivatives endows the induction equation with the solenoidal condition for the magnetic field

$$ \nabla \cdot \mathbf{B} = 0. \qquad (6) $$

This stationary condition, once satisfied at $t=0$, is analytically preserved at $t>0$. On the other hand, since numerical derivatives do not commute, the condition expressed in Eq. (6) is not automatically fulfilled by numerical schemes, leading to the rising of magnetic monopoles which cause unphysical plasma transport and dissipation of momentum and energy in ideal frameworks [48]. For these reasons, HD numerical schemes cannot be straightforwardly extended to MHD systems without introducing some corrections to numerically constraint Eq. (6) to be satisfied.

Similarly to the MHD case, the ideal RMHD equations can be expressed in the form given by Eq. (1) with

$$ U = \begin{pmatrix} D \\ \mathbf{m} \\ \epsilon \end{pmatrix}, \quad F = \begin{pmatrix} Dv \\ w_t \gamma^2 \mathbf{v} \mathbf{v} - \mathbf{b} \mathbf{b} + I p_t \\ \mathbf{m} \end{pmatrix}^T, \qquad (7) $$

where $D = \gamma \rho$, $m = (\rho h \gamma^2 + |B|^2)v - (v \cdot B)B$, and $\epsilon = \rho h \gamma^2 - p + |B|^2/2 + (|v|^2|B|^2 - (v \cdot B)^2)/2$ are respectively the relativistic mass density, momentum density, and energy density, while $w_t = \rho h + |B|^2/\gamma^2 + (v \cdot B)^2$ and $b = B/\gamma + \gamma(v \cdot B)v$. In the previous
expressions we introduced the specific enthalpy $h$, the Lorentz factor $\gamma = (1 - |v|^2)^{-1/2}$ and we set the speed of light $c$ to unity. In addition, the same induction equation as in the MHD regime needs to be added to Eq. (7), thus requiring the same numerical scheme to control the divergence-free condition shown in Eq. (6).

## 2.1. Notations

We employ a Cartesian coordinate system with unit vectors $\hat{e}_x = (1, 0, 0)$, $\hat{e}_y = (0, 1, 0)$, and $\hat{e}_z = (0, 0, 1)$ and maintain the same notation formalism adopted by Mignone & Del Zanna [18]. In each direction the coordinate spacing $\Delta x$, $\Delta y$, and $\Delta z$ is uniform. Computational cells are centered at $(x_i, y_j, z_k)$ and delimited by the six interfaces orthogonal to the coordinate axis centered, 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}})$.

According to Gauss' theorem, the conserved quantities are evolved as volume averages $\langle U \rangle_c$ over the cell volume

$$ \langle U \rangle_c = \frac{1}{\Delta x \Delta y \Delta z} \int U(x, y, z, t) dxdydz, \qquad (8) $$

where the $c$ subscript is a shorthand notation for $(i, j, k)$, by means of the neat difference of the surface-averaged representation of the fluxes at zone interfaces:

$$ \begin{aligned} \hat{F}_{x,x_f} &= \frac{1}{\Delta y \Delta z} \int \hat{e}_x \cdot F(U(x_{i+\frac{1}{2}}, y, z, t)) dy dz, \\ \hat{F}_{y,y_f} &= \frac{1}{\Delta x \Delta z} \int \hat{e}_y \cdot F(U(x, y_{j+\frac{1}{2}}, z, t)) dz dx, \\ \hat{F}_{z,z_f} &= \frac{1}{\Delta x \Delta y} \int \hat{e}_z \cdot F(U(x, y, z_{k+\frac{1}{2}}, t)) dx dy. \end{aligned} \qquad (9) $$

In the finite-volume fashion, a semi-discrete method of lines approximates the partial differential equation (PDE) associated with Eq. (1), yielding an ordinary differential equation (ODE) in the time variable

$$ \frac{d \langle U \rangle_c}{dt} = - \left( \frac{\hat{F}_{x,x_f} - \hat{F}_{x,x_f} - \hat{e}_x}{\Delta x} + \frac{\hat{F}_{y,y_f} - \hat{F}_{y,y_f} - \hat{e}_y}{\Delta y} + \frac{\hat{F}_{z,z_f} - \hat{F}_{z,z_f} - \hat{e}_z}{\Delta z} \right). \qquad (10) $$

Eqns. (8)-(10) frame the discretization for zone-centered variables. In constrained transport MHD instead, the magnetic field update relies on a discrete version of Stoke's theorem [49,50,45,5], where the magnetic field is treated as a staggered primary variable. In this representation, each component of the magnetic field

$$ \hat{\mathbf{B}}_f \equiv \begin{pmatrix} \hat{B}_{x,i+\frac{1}{2},j,k} \\ \hat{B}_{y,j,i+\frac{1}{2},k} \\ \hat{B}_{z,i,j,k+\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} \hat{B}_{x,x_f} \\ \hat{B}_{y,y_f} \\ \hat{B}_{z,z_f} \end{pmatrix} \qquad (11) $$

is evolved as an area-weighted average on the face orthogonal to the field component. The subscripts $x_f$, $y_f$ and $z_f$ identify the location of the faces, i.e., $x_f \equiv (i + \frac{1}{2}, j, k)$, $y_f \equiv (i, j + \frac{1}{2}, k)$, and $z_f \equiv (i, j, k + \frac{1}{2})$.

Furthermore, Eq. (4) prescribes that the variation in time of the magnetic flux is equal to the line integral of the electric field $\bar{E}_e$ over zone edges, namely, the electromotive force (EMF)

$$ \bar{E}_e \equiv \begin{pmatrix} \bar{E}_{x,i,j+\frac{1}{2},k+\frac{1}{2}} \\ \bar{E}_{y,i+\frac{1}{2},j,k+\frac{1}{2}} \\ \bar{E}_{z,i+\frac{1}{2},j+\frac{1}{2},k} \end{pmatrix} = \begin{pmatrix} \bar{E}_{x,x_e} \\ \bar{E}_{y,y_e} \\ \bar{E}_{z,z_e} \end{pmatrix}, \qquad (12) $$

where:

$$ \begin{aligned} \bar{E}_{x,x_e} &\equiv \frac{1}{\Delta x} \int E_x(x, y_{j+\frac{1}{2}}, z_{k+\frac{1}{2}}, t) dx, \\ \bar{E}_{y,y_e} &\equiv \frac{1}{\Delta y} \int E_y(x_{i+\frac{1}{2}}, y, z_{k+\frac{1}{2}}, t) dy, \\ \bar{E}_{z,z_e} &\equiv \frac{1}{\Delta z} \int E_z(x_{i+\frac{1}{2}}, y_{j+\frac{1}{2}}, z, t) dz. \end{aligned} \qquad (13) $$

Similarly to the magnetic field, the edge-centered electric field positions are labeled as $x_e \equiv (i, j + \frac{1}{2}, k + \frac{1}{2})$, $y_e \equiv (i + \frac{1}{2}, j, k + \frac{1}{2})$, and $z_e \equiv (i + \frac{1}{2}, j + \frac{1}{2}, k)$. Given all these elements, it is possible to retrieve the discretized version of Stokes theorem
expressions we introduced the specific enthalpy h, the Lorentz factor $\gamma = (1 - |v|^2)^{-1/2}$ and we set the speed of light c to unity. In addition, the same induction equation as in the MHD regime needs to be added to Eq. (7), thus requiring the same numerical scheme to control the divergence-free condition shown in Eq. (6).

## 2.1. Notations

We employ a Cartesian coordinate system with unit vectors $\hat{e}_x = (1, 0, 0)$, $\hat{e}_y = (0, 1, 0)$, and $\hat{e}_z = (0, 0, 1)$ and maintain the same notation formalism adopted by Mignone & Del Zanna [18]. In each direction the coordinate spacing $\Delta x$, $\Delta y$, and $\Delta z$ is uniform. Computational cells are centered at $(x_i, y_j, z_k)$ and delimited by the six interfaces orthogonal to the coordinate axis centered, respectively, at $(x_{i\pm\frac{1}{2}}, y_{j\pm\frac{1}{2}}, z_{k\pm\frac{1}{2}})$.

According to Gauss' theorem, the conserved quantities are evolved as volume averages $\langle U \rangle_c$ over the cell volume

$$ \langle U \rangle_c = \frac{1}{\Delta x \Delta y \Delta z} \int U(x, y, z, t) \, dxdydz, \qquad (8) $$

where the c subscript is a shorthand notation for $(i, j, k)$, by means of the neat difference of the surface-averaged representation of the fluxes at zone interfaces:

$$ \begin{aligned} \hat{F}_{x,x_f} &\equiv \frac{1}{\Delta y \Delta z} \int \hat{e}_x \cdot F(U(x_{i+\frac{1}{2}}, y, z, t)) \, dy \, dz, \\ \hat{F}_{y,y_f} &\equiv \frac{1}{\Delta x \Delta z} \int \hat{e}_y \cdot F(U(x, y_{j+\frac{1}{2}}, z, t)) \, dz \, dx, \\ \hat{F}_{z,z_f} &\equiv \frac{1}{\Delta x \Delta y} \int \hat{e}_z \cdot F(U(x, y, z_{k+\frac{1}{2}}, t)) \, dx \, dy. \end{aligned} \qquad (9) $$

In the finite-volume fashion, a semi-discrete method of lines approximates the partial differential equation (PDE) associated with Eq. (1), yielding an ordinary differential equation (ODE) in the time variable

$$ \frac{d \langle U \rangle_c}{dt} = - \left( \frac{\hat{F}_{x,x_f} - \hat{F}_{x,x_f} - \hat{e}_x}{\Delta x} + \frac{\hat{F}_{y,y_f} - \hat{F}_{y,y_f} - \hat{e}_y}{\Delta y} + \frac{\hat{F}_{z,z_f} - \hat{F}_{z,z_f} - \hat{e}_z}{\Delta z} \right). \qquad (10) $$

Eqns. (8)-(10) frame the discretization for zone-centered variables. In constrained transport MHD instead, the magnetic field update relies on a discrete version of Stoke's theorem [49,50,45,5], where the magnetic field is treated as a staggered primary variable. In this representation, each component of the magnetic field

$$ \hat{\mathbf{B}}_f \equiv \begin{pmatrix} \hat{B}_{x,i+\frac{1}{2},j,k} \\ \hat{B}_{y,j,i+\frac{1}{2},k} \\ \hat{B}_{z,i,j,k+\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} \hat{B}_{x,x_f} \\ \hat{B}_{y,y_f} \\ \hat{B}_{z,z_f} \end{pmatrix} \qquad (11) $$

is evolved as an area-weighted average on the face orthogonal to the field component. The subscripts $x_f$, $y_f$ and $z_f$ identify the location of the faces, i.e., $x_f \equiv (i + \frac{1}{2}, j, k)$, $y_f \equiv (i, j + \frac{1}{2}, k)$, and $z_f \equiv (i, j, k + \frac{1}{2})$.

Furthermore, Eq. (4) prescribes that the variation in time of the magnetic flux is equal to the line integral of the electric field $\bar{E}_e$ over zone edges, namely, the electromotive force (EMF)

$$ \bar{E}_e \equiv \begin{pmatrix} \bar{E}_{x,i,j+\frac{1}{2},k+\frac{1}{2}} \\ \bar{E}_{y,i+\frac{1}{2},j,k+\frac{1}{2}} \\ \bar{E}_{z,i+\frac{1}{2},j+\frac{1}{2},k} \end{pmatrix} = \begin{pmatrix} \bar{E}_{x,x_e} \\ \bar{E}_{y,y_e} \\ \bar{E}_{z,z_e} \end{pmatrix}, \qquad (12) $$

where:

$$ \begin{aligned} \bar{E}_{x,x_e} &\equiv \frac{1}{\Delta x} \int E_x(x, y_{j+\frac{1}{2}}, z_{k+\frac{1}{2}}, t) \, dx, \\ \bar{E}_{y,y_e} &\equiv \frac{1}{\Delta y} \int E_y(x_{i+\frac{1}{2}}, y, z_{k+\frac{1}{2}}, t) \, dy, \\ \bar{E}_{z,z_e} &\equiv \frac{1}{\Delta z} \int E_z(x_{i+\frac{1}{2}}, y_{j+\frac{1}{2}}, z, t) \, dz. \end{aligned} \qquad (13) $$

Similarly to the magnetic field, the edge-centered electric field positions are labeled as $x_e \equiv (i, j + \frac{1}{2}, k + \frac{1}{2})$, $y_e \equiv (i + \frac{1}{2}, j, k + \frac{1}{2})$, and $z_e \equiv (i + \frac{1}{2}, j + \frac{1}{2}, k)$. Given all these elements, it is possible to retrieve the discretized version of Stokes theorem
$$ 
\begin{aligned}
\frac{d \hat{B}_{x,x_f}}{dt} &= - \left( \frac{\hat{E}_{z,x_e} - \hat{E}_{z,x_e} - \hat{\epsilon}_y}{\Delta y} - \frac{\hat{E}_{y,y_e} - \hat{E}_{y,y_e} - \hat{\epsilon}_z}{\Delta z} \right), \\
\frac{d \hat{B}_{y,y_f}}{dt} &= - \left( \frac{\hat{E}_{x,x_e} - \hat{E}_{x,x_e} - \hat{\epsilon}_z}{\Delta z} - \frac{\hat{E}_{z,x_e} - \hat{E}_{z,x_e} - \hat{\epsilon}_x}{\Delta x} \right), \\
\frac{d \hat{B}_{z,z_f}}{dt} &= - \left( \frac{\hat{E}_{y,y_e} - \hat{E}_{y,y_e} - \hat{\epsilon}_x}{\Delta x} - \frac{\hat{E}_{x,x_e} - \hat{E}_{x,x_e} - \hat{\epsilon}_y}{\Delta y} \right),
\end{aligned}
\qquad (14) $$

where we set $c = 1$, thus verifying

$$ \frac{d}{dt} \left( \frac{\hat{B}_{x,x_f} - \hat{B}_{x,x_f} - \hat{\epsilon}_x}{\Delta x} + \frac{\hat{B}_{y,y_f} - \hat{B}_{y,y_f} - \hat{\epsilon}_y}{\Delta y} + \frac{\hat{B}_{z,z_f} - \hat{B}_{z,z_f} - \hat{\epsilon}_z}{\Delta z} \right) = 0, \qquad (15) $$

which is valid exactly and guarantees the fulfillment of the solenoidal condition for an initially divergence-free magnetic field. Notice that no approximation has been made in the treatment so far. To summarize our notation convention, given the generic physical quantity $Q$ we indicate with $\langle Q \rangle_c$ its volume average. We use $\hat{Q}_{x_f, y_f, z_f}$ and $\hat{Q}_{x_e, y_e, z_e}$ to refer to the face and edge averages, respectively. Similarly, the notations $Q_c$, $Q_{x_f, y_f, z_f}$, and $Q_{x_e, y_e, z_e}$ refer to the pointwise values at the center of the cell, face, and edge, respectively.

### 3. The 4th-order accurate method

In this section we introduce the main features of our 4th-order accurate numerical scheme. Finite volume Godunov-type methods evolve the cell-averaged conservative variables (e.g., $\langle U \rangle_c = \{\langle \rho \rangle_c, \langle \rho v \rangle_c, \langle \varepsilon \rangle_c\}$, while the representation of B depends on the choice of the algorithm chosen to control the solenoidal condition and may be either a cell-centered or a staggered quantity) and typically consist of three main stages: a reconstruction step at cell faces, flux computation at zone interfaces via Riemann solvers, and time integration. At 2nd-order accuracy, point values can be interchanged with volume or surface averages, so that the transformation between conservative and primitive variables (e.g., $V_c = \{\rho, v, p\}$) same arguments as before hold for B) can be simply defined as $V_c \approx \langle V \rangle_c \approx \mathcal{V}(\langle U \rangle_c)$. This property, however, does no longer hold for higher order schemes, where one has to distinguish the cell-centered point value from its volume average and, likewise, the face-centered flux from its surface average. On top of that, one-dimensional FV reconstruction schemes have to be reformulated in order to obtain 4th-order accurate left and right states at zone interfaces. These states are best obtained using primitive or characteristic variables in order to avoid unwanted numerical oscillations. These steps are described in the following sections.

#### 3.1. Point value recovery of primitive variables

Before reconstruction of primitive variables can take place, their values must become available at cell centers. Since conservative schemes evolve the volume average of conservative quantities $\langle Q \rangle_c$, their local point value $Q_c$ has to be obtained first. To 4th-order accuracy, point values can be retrieved from cell averages using the relation

$$ Q_c = \langle Q \rangle_c - \frac{\Delta \langle Q \rangle_c}{24} + O(h^4), \qquad (16) $$

and conversely

$$ \langle Q \rangle_c = Q_c + \frac{\Delta Q_c}{24} + O(h^4), \qquad (17) $$

where Eq. (17) is effectively a high-order Simpson-like quadrature rule. In Eqns. (16) and (17) Q is a component of $U_c$, while $\Delta$ is the Laplacian operator first introduced in the context of high-order methods by McCorquodale & Colella [28], i.e.

$$ \Delta Q_c \equiv \Delta^x Q_c + \Delta^y Q_c + \Delta^z Q_c, \qquad (18) $$

with components given by:

$$ \Delta^x Q_c \equiv (Q_{c-\hat{\epsilon}_x} - 2Q_c + Q_{c+\hat{\epsilon}_x}), \quad \Delta^y Q_c \equiv (Q_{c-\hat{\epsilon}_y} - 2Q_c + Q_{c+\hat{\epsilon}_y}), \quad \Delta^z Q_c \equiv (Q_{c-\hat{\epsilon}_z} - 2Q_c + Q_{c+\hat{\epsilon}_z}). \qquad (19) $$

The interested reader may refer to Appendix A for a full derivation of Eqns. (16) and (17).

Similarly, high-order transformations need to be applied also in the conversion between point values and face- or edge-averaged quantities, by means of the corresponding 2D or 1D operators. At zone faces, for instance, we use transverse Laplacian operators defined by:

$$ \Delta^x Q_{x_f} = \Delta^y Q_{x_f} + \Delta^z Q_{x_f}, \quad \Delta^y Q_{y_f} = \Delta^x Q_{y_f} + \Delta^z Q_{y_f}, \quad \Delta^z Q_{z_f} = \Delta^x Q_{z_f} + \Delta^y Q_{z_f}. \qquad (20) $$

This is employed, for example, to recover the face-centered point value of the magnetic field from its surface average:
$$ \begin{aligned} \frac{d \hat{B}_{x,x_f}}{dt} &= - \left( \frac{\hat{E}_{z,x_e} - \hat{E}_{z,x_e} - \hat{\epsilon}_y}{\Delta y} - \frac{\hat{E}_{y,y_e} - \hat{E}_{y,y_e} - \hat{\epsilon}_z}{\Delta z} \right), \\ \frac{d \hat{B}_{y,y_f}}{dt} &= - \left( \frac{\hat{E}_{x,x_e} - \hat{E}_{x,x_e} - \hat{\epsilon}_z}{\Delta z} - \frac{\hat{E}_{z,x_e} - \hat{E}_{z,x_e} - \hat{\epsilon}_x}{\Delta x} \right), \\ \frac{d \hat{B}_{z,z_f}}{dt} &= - \left( \frac{\hat{E}_{y,y_e} - \hat{E}_{y,y_e} - \hat{\epsilon}_x}{\Delta x} - \frac{\hat{E}_{x,x_e} - \hat{E}_{x,x_e} - \hat{\epsilon}_y}{\Delta y} \right), \end{aligned} \quad (14) $$

where we set $c = 1$, thus verifying

$$ \frac{d}{dt} \left( \frac{\hat{B}_{x,x_f} - \hat{B}_{x,x_f} - \hat{\epsilon}_x}{\Delta x} + \frac{\hat{B}_{y,y_f} - \hat{B}_{y,y_f} - \hat{\epsilon}_y}{\Delta y} + \frac{\hat{B}_{z,z_f} - \hat{B}_{z,z_f} - \hat{\epsilon}_z}{\Delta z} \right) = 0, \quad (15) $$

which is valid exactly and guarantees the fulfillment of the solenoidal condition for an initially divergence-free magnetic field. Notice that no approximation has been made in the treatment so far. To summarize our notation convention, given the generic physical quantity $Q$ we indicate with $\langle Q \rangle_c$ its volume average. We use $\hat{Q}_{x_f,y_f,z_f}$ and $\hat{Q}_{x_e,y_e,z_e}$ to refer to the face and edge averages, respectively. Similarly, the notations $Q_c$, $Q_{x_f,y_f,z_f}$, and $Q_{x_e,y_e,z_e}$ refer to the pointwise values at the center of the cell, face, and edge, respectively.

### 3. The 4th-order accurate method

In this section we introduce the main features of our 4th-order accurate numerical scheme. Finite volume Godunov-type methods evolve the cell-averaged conservative variables (e.g., $\langle U \rangle_c = \{\langle \rho \rangle_c, \langle \rho v \rangle_c, \langle \varepsilon \rangle_c\}$, while the representation of $\mathbf{B}$ depends on the choice of the algorithm chosen to control the solenoidal condition and may be either a cell-centered or a staggered quantity) and typically consist of three main stages: a reconstruction step at cell faces, flux computation at zone interfaces via Riemann solvers, and time integration. At 2nd-order accuracy, point values can be interchanged with volume or surface averages, so that the transformation between conservative and primitive variables (e.g., $V_c = \{\rho, v, p\}$) same arguments as before hold for $\mathbf{B}$) can be simply defined as $V_c \approx \langle V \rangle_c \approx \mathcal{V}(\langle U \rangle_c)$. This property, however, does no longer hold for higher order schemes, where one has to distinguish the cell-centered point value from its volume average and, likewise, the face-centered flux from its surface average. On top of that, one-dimensional FV reconstruction schemes have to be reformulated in order to obtain 4th-order accurate left and right states at zone interfaces. These states are best obtained using primitive or characteristic variables in order to avoid unwanted numerical oscillations. These steps are described in the following sections.

#### 3.1. Point value recovery of primitive variables

Before reconstruction of primitive variables can take place, their values must become available at cell centers. Since conservative schemes evolve the volume average of conservative quantities $\langle Q \rangle_c$, their local point value $Q_c$ has to be obtained first. To 4th-order accuracy, point values can be retrieved from cell averages using the relation

$$ Q_c = \langle Q \rangle_c - \frac{\Delta \langle Q \rangle_c}{24} + O(h^4), \quad (16) $$

and conversely

$$ \langle Q \rangle_c = Q_c + \frac{\Delta Q_c}{24} + O(h^4), \quad (17) $$

where Eq. (17) is effectively a high-order Simpson-like quadrature rule. In Eqns. (16) and (17) $Q$ is a component of $U_c$, while $\Delta$ is the Laplacian operator first introduced in the context of high-order methods by McCorquodale & Colella [28], i.e.

$$ \Delta Q_c \equiv \Delta^x Q_c + \Delta^y Q_c + \Delta^z Q_c, \quad (18) $$

with components given by:

$$ \Delta^x Q_c \equiv (Q_{c-\hat{\epsilon}_x} - 2Q_c + Q_{c+\hat{\epsilon}_x}), \quad \Delta^y Q_c \equiv (Q_{c-\hat{\epsilon}_y} - 2Q_c + Q_{c+\hat{\epsilon}_y}), \quad \Delta^z Q_c \equiv (Q_{c-\hat{\epsilon}_z} - 2Q_c + Q_{c+\hat{\epsilon}_z}). \quad (19) $$

The interested reader may refer to Appendix A for a full derivation of Eqns. (16) and (17).

Similarly, high-order transformations need to be applied also in the conversion between point values and face- or edge-averaged quantities, by means of the corresponding 2D or 1D operators. At zone faces, for instance, we use transverse Laplacian operators defined by:

$$ \Delta^x Q_{x_f} = \Delta^y Q_{x_f} + \Delta^z Q_{x_f}, \quad \Delta^y Q_{y_f} = \Delta^x Q_{y_f} + \Delta^z Q_{y_f}, \quad \Delta^z Q_{z_f} = \Delta^x Q_{z_f} + \Delta^y Q_{z_f}. \quad (20) $$

This is employed, for example, to recover the face-centered point value of the magnetic field from its surface average:
$$B_{x,x_f} = \hat{B}_{x,x_f} - \frac{\Delta_{\perp}^{x} \hat{B}_{x,x_f}}{24}, \quad B_{y,y_f} = \hat{B}_{y,y_f} - \frac{\Delta_{\perp}^{y} \hat{B}_{y,y_f}}{24}, \quad B_{z,z_f} = \hat{B}_{z,z_f} - \frac{\Delta_{\perp}^{z} \hat{B}_{z,z_f}}{24}. \qquad (21)$$

Following Felker & Stone [36], the pointwise value of the magnetic field at the cell center 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} \qquad (22)$$

Once the pointwise values of the conservative variables are obtained by means of Eqns. (16) and (22), the conversion to primitive variables takes place.

## 3.2. The pointwise MP5 and WENOZ spatial reconstruction algorithms

Here we introduce for the first time in a FV framework, to the extent of our knowledge, reconstruction schemes that operate directly on the function point values rather than on one-dimensional average quantities as traditionally done by previous investigations (see, e.g., [40,39,51,28,31,37]). Compared to other FV 4th-order methods, such as the one proposed by Felker & Stone [36], we do not need to convert the average value $\langle U \rangle_c$ into the mapped value $\mathcal{V}(\langle U_c \rangle)$, nor to compute from it the average value $\langle V \rangle_c$ before reconstructing, thus resulting in a more efficient and cost-effective scheme. In addition, the advantages of a point value reconstruction include that left and right states after reconstruction can be directly fed into the Riemann solver, avoiding the computation of interface-averaged Riemann problems. A similar feature can also be found in the scheme presented by Núñez-de La Rosa & Munz [34], as they reconstruct point values of the conservative quantities at the interfaces starting from cell-averaged values via a multi-dimensional Gaussian quadrature rule [51]. However, the approach from Núñez-de La Rosa & Munz is meant to reconstruct either conservative or characteristic quantities (whose volume-averaged values are the only available at the beginning of the integration step), while in our implementation we can either reconstruct primitive or characteristic variables (which in both approaches increase the stability of the scheme). More importantly, the high-order reconstruction based on the approach by [51] leads to a number of Riemann problems at each cell's interface equal to the number of Gaussian points used during the integration (e.g. 4 in the case of a three-dimensional domain when using 2 Gaussian points). Our scheme, on the other hand, always produces one Riemann problem per cell's interface.

The monotonicity preserving (MP5) scheme of Suresh & Huynh [39] provides a 5th-order accurate polynomial reconstruction at cell faces and, when necessary, limits the reconstructed value to preserve monotonicity near discontinuities while ensuring accuracy in smooth regions. While in the original formulation the unlimited interface value is retrieved by one-dimensional cell averages, our scheme replaces its definition by using point values in the same fashion as FD schemes (see, e.g., Appendix A.2 of Del Zanna et al. [20] for a brief summary). The unlimited value at $i + \frac{1}{2}$ is now given by

$$P(x)|_{i+\frac{1}{2}} = V_{i+\frac{1}{2}}^{L} = \frac{1}{128} (3V_{i-2} - 20V_{i-1} + 90V_{i} + 60V_{i+1} - 5V_{i+2}). \qquad (23)$$

where $\{V_{i\pm 2}, V_{i\pm 1}, V_i\}$ are now primitive quantities evaluated at cell-centers. The scheme employs, as in the original case, a five-point stencil in order to distinguish between local extrema and a genuine O(1) discontinuity. The MP5 limiting procedure remains unchanged as monotonicity-preserving bounds are equally well expressed in terms of point values. We refer the reader to the original work of Suresh & Huynh [39] or to Appendix 2 of Del Zanna et al. [20].

Likewise, we also introduce a pointwise version of the 5th-order weighted essentially non-oscillatory (WENO, [40]) scheme in the version proposed by Borges et al. (WENOZ, [41]), which attains better spatial resolution with reduced numerical dissipation at a modest computational cost. The same interface value of Eq. (23) at $x = x_{i+\frac{1}{2}}$ is computed as the convex combination of 3rd-order accurate interface values built on the three possible three-point sub-stencils $\{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}, \qquad (24)$$

where, again, $\{V_{i\pm 2}, V_{i\pm 1}, V_i\}$ are point values. The weights $\omega_i$, 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). \qquad (25)$$

Here $d_l$ denotes the optimal weights that reckon the 5th-order accurate approximation of Eq. (23), while $\epsilon = 10^{-40}$ is a small number which avoids division by zero. In our new framework, the optimal weights $\{d_0 = 1/16, d_1 = 5/8, d_2 = 5/16\}$ can be determined from the uniqueness of the interpolating polynomial with $d_0 + d_1 + d_2 = 1$. On the other hand, the smoothness indicators $\beta_l$ estimate the regularity of the polynomial approximation and, for a reconstruction based on point values, have been found to be the same as in Borges [41]:
$$B_{x,x_f} = \hat{B}_{x,x_f} - \frac{\Delta_{\perp}^{x} \hat{B}_{x,x_f}}{24}, \quad B_{y,y_f} = \hat{B}_{y,y_f} - \frac{\Delta_{\perp}^{y} \hat{B}_{y,y_f}}{24}, \quad B_{z,z_f} = \hat{B}_{z,z_f} - \frac{\Delta_{\perp}^{z} \hat{B}_{z,z_f}}{24}. \qquad (21)$$

Following Felker & Stone [36], the pointwise value of the magnetic field at the cell center 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} \qquad (22)$$

Once the pointwise values of the conservative variables are obtained by means of Eqns. (16) and (22), the conversion to primitive variables takes place.

## 3.2. The pointwise MP5 and WENOZ spatial reconstruction algorithms

Here we introduce for the first time in a FV framework, to the extent of our knowledge, reconstruction schemes that operate directly on the function point values rather than on one-dimensional average quantities as traditionally done by previous investigations (see, e.g., [40,39,51,28,31,37]). Compared to other FV 4th-order methods, such as the one proposed by Felker & Stone [36], we do not need to convert the average value $\langle U \rangle_c$ into the mapped value $\mathcal{V}(\langle U_c \rangle)$, nor to compute from it the average value $\langle V \rangle_c$ before reconstructing, thus resulting in a more efficient and cost-effective scheme. In addition, the advantages of a point value reconstruction include that left and right states after reconstruction can be directly fed into the Riemann solver, avoiding the computation of interface-averaged Riemann problems. A similar feature can also be found in the scheme presented by Núñez-de La Rosa & Munz [34], as they reconstruct point values of the conservative quantities at the interfaces starting from cell-averaged values via a multi-dimensional Gaussian quadrature rule [51]. However, the approach from Núñez-de La Rosa & Munz is meant to reconstruct either conservative or characteristic quantities (whose volume-averaged values are the only available at the beginning of the integration step), while in our implementation we can either reconstruct primitive or characteristic variables (which in both approaches increase the stability of the scheme). More importantly, the high-order reconstruction based on the approach by [51] leads to a number of Riemann problems at each cell's interface equal to the number of Gaussian points used during the integration (e.g. 4 in the case of a three-dimensional domain when using 2 Gaussian points). Our scheme, on the other hand, always produces one Riemann problem per cell's interface.

The monotonicity preserving (MP5) scheme of Suresh & Huynh [39] provides a 5th-order accurate polynomial reconstruction at cell faces and, when necessary, limits the reconstructed value to preserve monotonicity near discontinuities while ensuring accuracy in smooth regions. While in the original formulation the unlimited interface value is retrieved by one-dimensional cell averages, our scheme replaces its definition by using point values in the same fashion as FD schemes (see, e.g., Appendix A.2 of Del Zanna et al. [20] for a brief summary). The unlimited value at $i + \frac{1}{2}$ is now given by

$$P(x)|_{i+\frac{1}{2}} = V_{i+\frac{1}{2}}^{L} = \frac{1}{128} (3V_{i-2} - 20V_{i-1} + 90V_{i} + 60V_{i+1} - 5V_{i+2}). \qquad (23)$$

where $\{V_{i\pm 2}, V_{i\pm 1}, V_i\}$ are now primitive quantities evaluated at cell-centers. The scheme employs, as in the original case, a five-point stencil in order to distinguish between local extrema and a genuine O(1) discontinuity. The MP5 limiting procedure remains unchanged as monotonicity-preserving bounds are equally well expressed in terms of point values. We refer the reader to the original work of Suresh & Huynh [39] or to Appendix 2 of Del Zanna et al. [20].

Likewise, we also introduce a pointwise version of the 5th-order weighted essentially non-oscillatory (WENO, [40]) scheme in the version proposed by Borges et al. (WENOZ, [41]), which attains better spatial resolution with reduced numerical dissipation at a modest computational cost. The same interface value of Eq. (23) at $x = x_{i+\frac{1}{2}}$ is computed as the convex combination of 3rd-order accurate interface values built on the three possible three-point sub-stencils $\{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}, \qquad (24)$$

where, again, $\{V_{i\pm 2}, V_{i\pm 1}, V_i\}$ are point values. The weights $\omega_i$, 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). \qquad (25)$$

Here $d_l$ denotes the optimal weights that reckon the 5th-order accurate approximation of Eq. (23), while $\epsilon = 10^{-40}$ is a small number which avoids division by zero. In our new framework, the optimal weights $\{d_0 = 1/16, d_1 = 5/8, d_2 = 5/16\}$ can be determined from the uniqueness of the interpolating polynomial with $d_0 + d_1 + d_2 = 1$. On the other hand, the smoothness indicators $\beta_l$ estimate the regularity of the polynomial approximation and, for a reconstruction based on point values, have been found to be the same as in Borges [41]:
$$
\begin{align}
\beta_0 &= \frac{13}{12}(V_{i-2} - 2V_{i-1} + V_i)^2 + \frac{1}{4}(V_{i-2} - 4V_{i-1} + 3V_i)^2, \nonumber \\
\beta_1 &= \frac{13}{12}(V_{i-1} - 2V_i + V_{i+1})^2 + \frac{1}{4}(V_{i-1} - V_{i+1})^2, \tag{26} \\
\beta_2 &= \frac{13}{12}(V_i - 2V_{i+1} + V_{i+2})^2 + \frac{1}{4}(3V_i - 4V_{i+1} + V_{i+2})^2. \nonumber
\end{align}
$$

Note that the right state, $V_{i-\frac{1}{2}}^R$, is obtained similarly by reversing the index order (e.g. $i-1 \to i+1$, etc).

We will also occasionally employ the 3rd-order WENO scheme (here denoted as WENO3) [40] in the improved version of Yamaleev & Carpenter [52]. The point value approach leads to the following expression for the left interface value

$$ V_{i+\frac{1}{2}}^{L} = \omega_0 \frac{V_{i+1} + V_i}{2} + \omega_1 \frac{-V_{i+1} + 3V_i}{2}, \qquad (27) $$

where the weights are given by:

$$ \omega_l = \frac{\alpha_l}{\alpha_0 + \alpha_1}, \quad \alpha_l = d_l \left( 1 + \frac{|\Delta_{i+\frac{1}{2}} - \Delta_{i-\frac{1}{2}}|^2}{\beta_l + \epsilon} \right), \quad \beta_0 = \Delta_{i-\frac{1}{2}}^2, \quad \beta_1 = \Delta_{i+\frac{1}{2}}^2, \qquad (28) $$

with $l = \{0, 1\}$. In Eq. (28), $\Delta_{i+\frac{1}{2}} = (V_{i+1} - V_i)$ whereas $\{d_0 = 3/4, d_1 = 1/4\}$ (note that in the standard finite volume formulation, starting from 1D cell averages, one has $d_0 = 2/3, d_1 = 1/3$). Finally, in order to avoid loss of accuracy at critical points (see [52] and [15]) we adopt here $\epsilon = \Delta x^2$.

### 3.3. Flux computation and temporal integration

Spatial reconstruction provides the 4th-order accurate pointwise values for left/right states at a given face center (e.g. $V_{x_f}^L, V_{x_f}^R$). These serve as input states to the Riemann solver for the computation of the high-order interface fluxes. Our implementation comes with a variety of Riemann solvers, such as the linearized Roe solver [53], HLL-type Riemann solvers [54], including the HLLC [55,56,11] and HLLD [57,58], the Lax-Friedrichs [59], and the GFORCE [60] solver. Denoting with $\mathcal{F}_{x,y,z}(\cdot, \cdot)$ a generic 1-D Riemann solver along the $\{x, y, z\}$ direction, the interface upwind punctual fluxes are obtained as:

$$ F_{x,x_f} = \mathcal{F}_x(V_{x_f}^L, V_{x_f}^R), \quad F_{y,x_f} = \mathcal{F}_y(V_{y_f}^L, V_{y_f}^R), \quad F_{z,x_f} = \mathcal{F}_z(V_{z_f}^L, V_{z_f}^R). \qquad (29) $$

4th-order accurate surface averages are then obtained by the quadrature rule in Eq. (17) using the transverse Laplacian operators:

$$ \hat{F}_{x,x_f} = F_{x,x_f} + \frac{\Delta_{\perp}^x F_{x,x_f}}{24}, \quad \hat{F}_{y,y_f} = F_{y,y_f} + \frac{\Delta_{\perp}^y F_{y,y_f}}{24}, \quad \hat{F}_{z,z_f} = F_{z,z_f} + \frac{\Delta_{\perp}^z F_{z,z_f}}{24}. \qquad (30) $$

Note that this requires an additional integration in the first layer of ghost zones in order to retrieve the flux information. Moreover, thanks to the deployment of point value reconstructions, we do not need to solve a Riemann problem with interface-averaged primitive states as it is instead done by Felker & Stone (Eq. 17 in [36]).

These fluxes are used to build the right-hand side of Eq. (10) which will be hereafter denoted with $\mathcal{L}(U)$. Our numerical scheme employs the five stage 4th-order explicit Strong Stability Preserving Runge-Kutta method (eSSPRK(5,4), [43,44]) to discretize Eq. (10)

$$
\begin{align}
U^{(1)} &= U^n + 0.391752226571890\Delta t \mathcal{L}(U^n) \nonumber \\
U^{(2)} &= 0.444370493651235U^n + 0.555629506348765U^{(1)} + 0.368410593050371\Delta t \mathcal{L}(U^{(1)}) \nonumber \\
U^{(3)} &= 0.620101851488403U^n + 0.379898148511597U^{(2)} + 0.251891774271694\Delta t \mathcal{L}(U^{(2)}) \nonumber \\
U^{(4)} &= 0.178079954393132U^n + 0.821920045606868U^{(3)} + 0.544974750228521\Delta t \mathcal{L}(U^{(3)}) \tag{31} \\
U^{n+1} &= 0.517231671970585U^{(2)} + 0.096059710526147U^{(3)} + 0.063692468666290\Delta t \mathcal{L}(U^{(3)}) \nonumber \\
&\quad +0.386708617503268U^{(4)} + 0.226007483236906\Delta t \mathcal{L}(U^{(4)}). \nonumber
\end{align}
$$

The eSSPRK(5,4) method is strongly stable (namely, no bounded temporal growth is allowed: $\|U^{n+1}\| \le \|U^n\|$) under the assumption that the forward Euler method employed is strongly stable under the Courant-Friedrichs-Levy (CFL) restriction ($\|U^n + \Delta t \mathcal{L}(U^n)\| \le \|U^n\|$, [61,62]). In our implementation, the CFL condition determines the time step according to

$$ \Delta t = C \min_{ijk} \left[ \min_d \left( \frac{\Delta x^d}{|\lambda_{\max}^d|} \right) \right] \qquad (32) $$

where C is the CFL number, $\lambda_{\max}^d$ is the outermost fast magnetosonic wave and $\Delta x^d$ is the cell length through the different directions d.

For a detailed description of the stability properties of the method we refer to [62,44,63–65].
$$
\begin{align}
\beta_0 &= \frac{13}{12}(V_{i-2} - 2V_{i-1} + V_i)^2 + \frac{1}{4}(V_{i-2} - 4V_{i-1} + 3V_i)^2, \nonumber \\
\beta_1 &= \frac{13}{12}(V_{i-1} - 2V_i + V_{i+1})^2 + \frac{1}{4}(V_{i-1} - V_{i+1})^2, \tag{26} \\
\beta_2 &= \frac{13}{12}(V_i - 2V_{i+1} + V_{i+2})^2 + \frac{1}{4}(3V_i - 4V_{i+1} + V_{i+2})^2. \nonumber
\end{align}
$$

Note that the right state, $V_{i-\frac{1}{2}}^R$, is obtained similarly by reversing the index order (e.g. $i-1 \to i+1$, etc).

We will also occasionally employ the 3rd-order WENO scheme (here denoted as WENO3) [40] in the improved version of Yamaleev & Carpenter [52]. The point value approach leads to the following expression for the left interface value

$$ V_{i+\frac{1}{2}}^{L} = \omega_0 \frac{V_{i+1} + V_i}{2} + \omega_1 \frac{-V_{i+1} + 3V_i}{2}, \qquad (27) $$

where the weights are given by:

$$ \omega_l = \frac{\alpha_l}{\alpha_0 + \alpha_1}, \quad \alpha_l = d_l \left( 1 + \frac{|\Delta_{i+\frac{1}{2}} - \Delta_{i-\frac{1}{2}}|^2}{\beta_l + \epsilon} \right), \quad \beta_0 = \Delta_{i-\frac{1}{2}}^2, \quad \beta_1 = \Delta_{i+\frac{1}{2}}^2, \qquad (28) $$

with $l = \{0, 1\}$. In Eq. (28), $\Delta_{i+\frac{1}{2}} = (V_{i+1} - V_i)$ whereas $\{d_0 = 3/4, d_1 = 1/4\}$ (note that in the standard finite volume formulation, starting from 1D cell averages, one has $d_0 = 2/3, d_1 = 1/3$). Finally, in order to avoid loss of accuracy at critical points (see [52] and [15]) we adopt here $\epsilon = \Delta x^2$.

### 3.3. Flux computation and temporal integration

Spatial reconstruction provides the 4th-order accurate pointwise values for left/right states at a given face center (e.g. $V_{x_f}^L, V_{x_f}^R$). These serve as input states to the Riemann solver for the computation of the high-order interface fluxes. Our implementation comes with a variety of Riemann solvers, such as the linearized Roe solver [53], HLL-type Riemann solvers [54], including the HLLC [55,56,11] and HLLD [57,58], the Lax-Friedrichs [59], and the GFORCE [60] solver. Denoting with $\mathcal{F}_{x,y,z}(\cdot, \cdot)$ a generic 1-D Riemann solver along the $\{x, y, z\}$ direction, the interface upwind punctual fluxes are obtained as:

$$ F_{x,x_f} = \mathcal{F}_x(V_{x_f}^L, V_{x_f}^R), \quad F_{y,x_f} = \mathcal{F}_y(V_{y_f}^L, V_{y_f}^R), \quad F_{z,x_f} = \mathcal{F}_z(V_{z_f}^L, V_{z_f}^R). \qquad (29) $$

4th-order accurate surface averages are then obtained by the quadrature rule in Eq. (17) using the transverse Laplacian operators:

$$ \hat{F}_{x,x_f} = F_{x,x_f} + \frac{\Delta_{\perp}^x F_{x,x_f}}{24}, \quad \hat{F}_{y,x_f} = F_{y,x_f} + \frac{\Delta_{\perp}^y F_{y,x_f}}{24}, \quad \hat{F}_{z,x_f} = F_{z,x_f} + \frac{\Delta_{\perp}^z F_{z,x_f}}{24}. \qquad (30) $$

Note that this requires an additional integration in the first layer of ghost zones in order to retrieve the flux information. Moreover, thanks to the deployment of point value reconstructions, we do not need to solve a Riemann problem with interface-averaged primitive states as it is instead done by Felker & Stone (Eq. 17 in [36]).

These fluxes are used to build the right-hand side of Eq. (10) which will be hereafter denoted with $\mathcal{L}(U)$. Our numerical scheme employs the five stage 4th-order explicit Strong Stability Preserving Runge-Kutta method (eSSPRK(5,4), [43,44]) to discretize Eq. (10)

$$
\begin{align}
U^{(1)} &= U^n + 0.391752226571890\Delta t \mathcal{L}(U^n) \nonumber \\
U^{(2)} &= 0.444370493651235U^n + 0.555629506348765U^{(1)} + 0.368410593050371\Delta t \mathcal{L}(U^{(1)}) \nonumber \\
U^{(3)} &= 0.620101851488403U^n + 0.379898148511597U^{(2)} + 0.251891774271694\Delta t \mathcal{L}(U^{(2)}) \nonumber \\
U^{(4)} &= 0.178079954393132U^n + 0.821920045606868U^{(3)} + 0.544974750228521\Delta t \mathcal{L}(U^{(3)}) \tag{31} \\
U^{n+1} &= 0.517231671970585U^{(2)} + 0.096059710526147U^{(3)} + 0.063692468666290\Delta t \mathcal{L}(U^{(3)}) \nonumber \\
&\quad +0.386708617503268U^{(4)} + 0.226007483236906\Delta t \mathcal{L}(U^{(4)}). \nonumber
\end{align}
$$

The eSSPRK(5,4) method is strongly stable (namely, no bounded temporal growth is allowed: $\|U^{n+1}\| \le \|U^n\|$) under the assumption that the forward Euler method employed is strongly stable under the Courant-Friedrichs-Levy (CFL) restriction ($\|U^n + \Delta t \mathcal{L}(U^n)\| \le \|U^n\|$, [61,62]). In our implementation, the CFL condition determines the time step according to

$$ \Delta t = C \min_{ijk} \left[ \min_d \left( \frac{\Delta x^d}{|\lambda_{\max}^d|} \right) \right] \qquad (32) $$

where C is the CFL number, $\lambda_{\max}^d$ is the outermost fast magnetosonic wave and $\Delta x^d$ is the cell length through the different directions d.

For a detailed description of the stability properties of the method we refer to [62,44,63–65].
### 3.4. Order reduction in the presence of discontinuities

The recovery of point values from volume averages follows the assumption that the function $Q$ is locally continuous up to the 3rd derivative in its expansion. However, this assumption manifestly breaks down in the immediate neighbor of a discontinuity, where using Eq. (16) may easily produce unphysical values such as negative energies or densities. In order to address this shortcoming, we modify the recovery of pointwise conservative variables at the beginning of any integration stage by modifying Eq. (16) as

$$U_c = \langle U \rangle_c - \theta_c \frac{\Delta \langle U \rangle_c}{24} + O(h^4), \qquad (33)$$

where $\theta_c$ has been introduced to detect the presence of a discontinuity

$$\theta_c = \begin{cases} 1 & \text{if } \eta_c < \eta_d, \\ 0 & \text{otherwise.} \end{cases} \qquad (34)$$

The parameter $\eta_c$ is defined as $\eta_c = \sqrt{\eta_{c,x}^2 + \eta_{c,y}^2 + \eta_{c,z}^2}$, while $\eta_d$ is a threshold parameter. When $\theta_c = 0$ in Eq. (33), the local order of the scheme is reduced in proximity of discontinuous fronts, including avoiding the integration of fluxes (Eq. (30)). In these regions, we further increase the robustness of the scheme by also lowering the order of the reconstruction scheme by switching to the 3rd-order WENO scheme (see end of §3.2) or to piecewise linear reconstruction. This approach (known as the “fallback”) has been employed, in a similar fashion, also by [34] in the context of high-order schemes. We now propose two different discontinuity detectors in order to evaluate $\eta_c$.

The first one is based on Jameson's pressure-based shock sensor [66]

$$\eta_{c,x} = \frac{\langle Q \rangle_{c+\hat{e}_x} - 2 \langle Q \rangle_c + \langle Q \rangle_{c-\hat{e}_x}}{|\langle Q \rangle_{c+\hat{e}_x}| + 2|\langle Q \rangle_c| + |\langle Q \rangle_{c-\hat{e}_x}| + \epsilon}, \qquad (35)$$

(similar expressions hold for $\eta_{c,y}$ and $\eta_{c,z}$) where $\epsilon$ is a small number to prevent negligible variations from triggering the switch. Eq. (35) is first evaluated with $\langle Q \rangle$ being the volume-averaged density, pressure, and magnetic pressure and the maximum over these three values is finally considered. Jameson's switch has the advantage of keeping the stencil compact (no additional ghost zones are needed), albeit three zones may not be sufficient to distinguish a discontinuity from a smooth extremum [39].

The second discontinuity sensor that we propose employs a derivatives ratio that attempts to identify high-frequency oscillations which typically arise near a discontinuous front. Since high wave-number Fourier modes will not be dumped around these critical points, we expect the contribution of higher order derivatives to be non-negligible. To this end, we first evaluate the normalized contribution of odd and even modes separately by evaluating the parameters $\eta_c^o$ and $\eta_c^e$

$$\eta_c^o = \frac{|\delta^{(3)} \langle Q \rangle_c|}{|\langle Q \rangle_{c,ref}| + |\delta^{(1)} \langle Q \rangle_c| + |\delta^{(3)} \langle Q \rangle_c| + \epsilon}, \qquad \eta_c^e = \frac{|\delta^{(4)} \langle Q \rangle_c|}{|\langle Q \rangle_{c,ref}| + |\delta^{(2)} \langle Q \rangle_c| + |\delta^{(4)} \langle Q \rangle_c| + \epsilon}, \qquad (36)$$

where $\delta^{(m)}Q_c$ (with $m = 1, .., 4$) are undivided 2nd-order accurate approximations to the $m$-th derivative. In the $x$-direction, for instance:

$$\begin{aligned} \delta^{(1)} \langle Q \rangle_c &= \frac{1}{2} (\langle Q \rangle_{c+\hat{e}_x} - \langle Q \rangle_{c-\hat{e}_x}), \\ \delta^{(2)} \langle Q \rangle_c &= \langle Q \rangle_{c+\hat{e}_x} - 2 \langle Q \rangle_c + \langle Q \rangle_{c-\hat{e}_x}, \\ \delta^{(3)} \langle Q \rangle_c &= \frac{1}{2} (\langle Q \rangle_{i+2\hat{e}_x} - 2(\langle Q \rangle_{c+\hat{e}_x} - 2\langle Q \rangle_{c-\hat{e}_x}) - \langle Q \rangle_{c-2\hat{e}_x}), \\ \delta^{(4)} \langle Q \rangle_c &= \langle Q \rangle_{c+2\hat{e}_x} - 4\langle Q \rangle_{c+\hat{e}_x} + 6\langle Q \rangle_c - 4\langle Q \rangle_{c-\hat{e}_x} + \langle Q \rangle_{c-2\hat{e}_x}. \end{aligned} \qquad (37)$$

In our experience, using primitive variables has shown to yield a more robust identification method. For this reason, we take $\langle Q \rangle_c = V(\langle U \rangle_c)$ at the beginning of the Runge-Kutta stage. The reference value $\langle Q \rangle_{c,ref}$ in Eq. (36) is equal to $\langle Q \rangle_c$ for density and pressure while we employ $\sqrt{\rho_c/\rho_c}$ (isothermal sound speed) for velocity components and $|B|_c$ for magnetic field components. Once $\eta_c^o$ and $\eta_c^e$ have been computed, we simply take the maximum of the two: $\eta_c = \max(\eta_c^o, \eta_c^e)$.

The derivative ratio detector needs a larger stencil thereby asking for one additional ghost zone. Along the same arguments, a more sophisticated approach that requires an even larger stencil, is shown by Bambozzi & Pires [67].

### 4. The 4th-order constrained transport

As we outlined in §2, the HD sub-system cannot be straightforwardly extended to the MHD system without introducing numerical strategies to control the rising of spurious magnetic monopoles. Failure to fulfill this requisite may lead to unphysical effects such as plasma acceleration parallel to the field, incorrect jump conditions, wrong propagation speed of discontinuities and odd-even decoupling [48,6,5]. In what follows we employ the upwind constrained transport method [45,3,5,46,47,18] in order to ensure that the divergence-free condition is maintained to machine accuracy.

In the FV-CT formalism the field's surface integrals $\hat{B}_{x,x_f}$, $\hat{B}_{y,y_f}$, $\hat{B}_{z,z_f}$ are the primary representation of the magnetic field evolved in time. Eq. (4) regulates the time variation of the magnetic field fluxes, implying that the formal correspondence
### 3.4. Order reduction in the presence of discontinuities

The recovery of point values from volume averages follows the assumption that the function $Q$ is locally continuous up to the 3rd derivative in its expansion. However, this assumption manifestly breaks down in the immediate neighbor of a discontinuity, where using Eq. (16) may easily produce unphysical values such as negative energies or densities. In order to address this shortcoming, we modify the recovery of pointwise conservative variables at the beginning of any integration stage by modifying Eq. (16) as

$$U_c = \langle U \rangle_c - \theta_c \frac{\Delta \langle U \rangle_c}{24} + O(h^4), \quad (33)$$

where $\theta_c$ has been introduced to detect the presence of a discontinuity

$$\theta_c = \begin{cases} 1 & \text{if } \eta_c < \eta_d, \\ 0 & \text{otherwise.} \end{cases} \quad (34)$$

The parameter $\eta_c$ is defined as $\eta_c = \sqrt{\eta_{c,x}^2 + \eta_{c,y}^2 + \eta_{c,z}^2}$, while $\eta_d$ is a threshold parameter. When $\theta_c = 0$ in Eq. (33), the local order of the scheme is reduced in proximity of discontinuous fronts, including avoiding the integration of fluxes (Eq. (30)). In these regions, we further increase the robustness of the scheme by also lowering the order of the reconstruction scheme by switching to the 3rd-order WENO scheme (see end of §3.2) or to piecewise linear reconstruction. This approach (known as the “fallback”) has been employed, in a similar fashion, also by [34] in the context of high-order schemes. We now propose two different discontinuity detectors in order to evaluate $\eta_c$.

The first one is based on Jameson's pressure-based shock sensor [66]

$$\eta_{c,x} = \frac{\langle Q \rangle_{c+\hat{e}_x} - 2 \langle Q \rangle_c + \langle Q \rangle_{c-\hat{e}_x}}{|\langle Q \rangle_{c+\hat{e}_x}| + 2|\langle Q \rangle_c| + |\langle Q \rangle_{c-\hat{e}_x}| + \epsilon}, \quad (35)$$

(similar expressions hold for $\eta_{c,y}$ and $\eta_{c,z}$) where $\epsilon$ is a small number to prevent negligible variations from triggering the switch. Eq. (35) is first evaluated with $\langle Q \rangle$ being the volume-averaged density, pressure, and magnetic pressure and the maximum over these three values is finally considered. Jameson's switch has the advantage of keeping the stencil compact (no additional ghost zones are needed), albeit three zones may not be sufficient to distinguish a discontinuity from a smooth extremum [39].

The second discontinuity sensor that we propose employs a derivatives ratio that attempts to identify high-frequency oscillations which typically arise near a discontinuous front. Since high wave-number Fourier modes will not be dumped around these critical points, we expect the contribution of higher order derivatives to be non-negligible. To this end, we first evaluate the normalized contribution of odd and even modes separately by evaluating the parameters $\eta_c^o$ and $\eta_c^e$

$$\eta_c^o = \frac{|\delta^{(3)} \langle Q \rangle_c|}{|\langle Q \rangle_{c,ref}| + |\delta^{(1)} \langle Q \rangle_c| + |\delta^{(3)} \langle Q \rangle_c| + \epsilon}, \quad \eta_c^e = \frac{|\delta^{(4)} \langle Q \rangle_c|}{|\langle Q \rangle_{c,ref}| + |\delta^{(2)} \langle Q \rangle_c| + |\delta^{(4)} \langle Q \rangle_c| + \epsilon}, \quad (36)$$

where $\delta^{(m)}Q_c$ (with $m = 1, .., 4$) are undivided 2nd-order accurate approximations to the $m$-th derivative. In the $x$-direction, for instance:

$$\begin{align}
\delta^{(1)} \langle Q \rangle_c &= \frac{1}{2} (\langle Q \rangle_{c+\hat{e}_x} - \langle Q \rangle_{c-\hat{e}_x}), \nonumber \\
\delta^{(2)} \langle Q \rangle_c &= \langle Q \rangle_{c+\hat{e}_x} - 2 \langle Q \rangle_c + \langle Q \rangle_{c-\hat{e}_x}, \nonumber \\
\delta^{(3)} \langle Q \rangle_c &= \frac{1}{2} (\langle Q \rangle_{i+2\hat{e}_x} - 2(\langle Q \rangle_{c+\hat{e}_x} - 2 \langle Q \rangle_{c-\hat{e}_x}) - \langle Q \rangle_{c-2\hat{e}_x}), \nonumber \\
\delta^{(4)} \langle Q \rangle_c &= \langle Q \rangle_{c+2\hat{e}_x} - 4 \langle Q \rangle_{c+\hat{e}_x} + 6 \langle Q \rangle_c - 4 \langle Q \rangle_{c-\hat{e}_x} + \langle Q \rangle_{c-2\hat{e}_x}. \tag{37}
\end{align}$$

In our experience, using primitive variables has shown to yield a more robust identification method. For this reason, we take $\langle Q \rangle_c = V(\langle U \rangle_c)$ at the beginning of the Runge-Kutta stage. The reference value $\langle Q \rangle_{c,ref}$ in Eq. (36) is equal to $\langle Q \rangle_c$ for density and pressure while we employ $\sqrt{\rho_c/\rho_c}$ (isothermal sound speed) for velocity components and $|B|_c$ for magnetic field components. Once $\eta_c^o$ and $\eta_c^e$ have been computed, we simply take the maximum of the two: $\eta_c = \max(\eta_c^o, \eta_c^e)$.

The derivative ratio detector needs a larger stencil thereby asking for one additional ghost zone. Along the same arguments, a more sophisticated approach that requires an even larger stencil, is shown by Bambozzi & Pires [67].

### 4. The 4th-order constrained transport

As we outlined in §2, the HD sub-system cannot be straightforwardly extended to the MHD system without introducing numerical strategies to control the rising of spurious magnetic monopoles. Failure to fulfill this requisite may lead to unphysical effects such as plasma acceleration parallel to the field, incorrect jump conditions, wrong propagation speed of discontinuities and odd-even decoupling [48,6,5]. In what follows we employ the upwind constrained transport method [45,3,5,46,47,18] in order to ensure that the divergence-free condition is maintained to machine accuracy.

In the FV-CT formalism the field's surface integrals $\hat{B}_{x,x_f}$, $\hat{B}_{y,y_f}$, $\hat{B}_{z,z_f}$ are the primary representation of the magnetic field evolved in time. Eq. (4) regulates the time variation of the magnetic field fluxes, implying that the formal correspondence
<figure><img src="image_2.png" /><figcaption>Fig. 1. Left: positioning of MHD variables in the UCT formalism. Staggered magnetic field components (blue) are face-centered, the electromotive force is edge-centered (red) while remaining HD quantities are located at the zone center (green). Right: top view of the intersection between four neighbor zones: N, S, E, and W indicate the four cardinal directions with respect to the zone edge (here represented by the intersection between four neighbor zones), $R_x(F_x)$ and $R_y(F_y)$ are 1-D reconstruction operators applied to each zone face, see Eq. (38). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

$$ 
\begin{aligned}
F_x^{[B_x]} &= 0, & F_y^{[B_x]} &= E_z, & F_z^{[B_x]} &= -E_y, \\
F_x^{[B_y]} &= -E_z, & F_y^{[B_y]} &= 0, & F_z^{[B_y]} &= E_x, \\
F_x^{[B_z]} &= E_y, & F_y^{[B_z]} &= -E_x, & F_z^{[B_z]} &= 0,
\end{aligned} 
\qquad (38) $$

holds when solving a Riemann problem. While this certainly holds for the fluxes computed at the face midpoint, Eq. (14) requires the EMF to become available at a zone edge to serve as a discrete version of Stokes theorem. In order to reconstruct the induction fluxes from the face center (where they are available as point values) to one of the four adjacent edges, we employ the UCT method of Mignone & Del Zanna [18]. This method provides a general formalism for systematic construction of EMF averaging procedures where dissipation terms depend on the left and right transverse magnetic field components.

Referring to the right panel of Fig. 1 and being consistent with the notation presented in Mignone & Del Zanna [18], we consider a top view of the intersection of four zones at a cell-edge. We denote, respectively, the left and right reconstructed states along the y-coordinate as south (S) and north (N), while on the x-direction the reconstructed states are labeled with west (W) or east (E) with respect to the intersection point. The edge-centered EMF with the desired upwind properties can be built up from

$$ E_{z,x_e} = - \left[ (a_x \bar{v}_x B_y)^W + (a_x \bar{v}_x B_y)^E \right] + \left[ (a_y \bar{v}_y B_x)^N + (a_y \bar{v}_y B_x)^S \right] + \left[ (d_x B_y)^E - (d_x B_y)^W \right] - \left[ (d_y B_x)^N - (d_y B_x)^S \right], \quad (39) $$

where $\bar{v}_{x,y}$ are the transverse velocities reconstructed from the interfaces, $a_{x,y}$ and $d_{x,y}$ are the flux and diffusion coefficients, and $B_{x,y}$ are the pointwise values of the magnetic field reconstructed from the corresponding interface. The other EMF components may be obtained by cyclic permutations. For the sake of briefness, we address the reader to find all the details of the different EMF averaging techniques that we employ in the present work in Mignone & Del Zanna [18]. Note that the employment of point value reconstruction methods (rather than 1D volume average-based ones) requires the electric field to be evaluated over 3 layers of ghost zones in the transverse direction rather than 4.

Finally, the line-averaged integrated electric field is computed as

$$ \bar{E}_{x,x_e} = E_{x,x_e} + \frac{\Delta^x E_{x,x_e}}{24}, \quad \bar{E}_{y,y_e} = E_{y,y_e} + \frac{\Delta^y E_{y,y_e}}{24}, \quad \bar{E}_{z,z_e} = E_{z,z_e} + \frac{\Delta^z E_{z,z_e}}{24}, \qquad (40) $$

and then employed to update the staggered magnetic field components using the discrete version of Stoke's theorem Eq. (14).

In the following we recap the steps of our 4th-order numerical scheme.

1. Start with the volume average of conserved quantities $\langle U \rangle_c$ and the face average of the staggered magnetic field components, $\hat{B}_{x,x_f}^x, \hat{B}_{y,y_f}^y, \hat{B}_{z,z_f}^z$.

2. Assign boundary conditions to these variables.

3. Obtain volume-averaged primitive variables by converting the volume averages $\langle V \rangle_c = V(\langle U \rangle_c)$. Notice that this array is 2nd-order accurate and it is only used in the limiting process, thus it does not replace the pointwise conversion later on.

4. Evaluate the $\theta_c$ parameter for order reduction in presence of discontinuities by means of either the Jameson's shock sensor (Eq. (35)) or by weighing up the contribution of higher-order derivatives (Eq. (36)) and flag troubled cells along with first neighbors.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Left: positioning of MHD variables in the UCT formalism. Staggered magnetic field components (blue) are face-centered, the electromotive force is edge-centered (red) while remaining HD quantities are located at the zone center (green). Right: top view of the intersection between four neighbor zones: N, S, E, and W indicate the four cardinal directions with respect to the zone edge (here represented by the intersection between four neighbor zones), $R_x(F_x)$ and $R_y(F_y)$ are 1-D reconstruction operators applied to each zone face, see Eq. (38). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

$$ 
\begin{aligned}
F_x^{[B_x]} &= 0, & F_y^{[B_x]} &= E_z, & F_z^{[B_x]} &= -E_y, \\
F_x^{[B_y]} &= -E_z, & F_y^{[B_y]} &= 0, & F_z^{[B_y]} &= E_x, \\
F_x^{[B_z]} &= E_y, & F_y^{[B_z]} &= -E_x, & F_z^{[B_z]} &= 0,
\end{aligned} 
\qquad (38) $$

holds when solving a Riemann problem. While this certainly holds for the fluxes computed at the face midpoint, Eq. (14) requires the EMF to become available at a zone edge to serve as a discrete version of Stokes theorem. In order to reconstruct the induction fluxes from the face center (where they are available as point values) to one of the four adjacent edges, we employ the UCT method of Mignone & Del Zanna [18]. This method provides a general formalism for systematic construction of EMF averaging procedures where dissipation terms depend on the left and right transverse magnetic field components.

Referring to the right panel of Fig. 1 and being consistent with the notation presented in Mignone & Del Zanna [18], we consider a top view of the intersection of four zones at a cell-edge. We denote, respectively, the left and right reconstructed states along the y-coordinate as south (S) and north (N), while on the x-direction the reconstructed states are labeled with west (W) or east (E) with respect to the intersection point. The edge-centered EMF with the desired upwind properties can be built up from

$$ E_{z,x_e} = - \left[ (a_x \bar{v}_x B_y)^W + (a_x \bar{v}_x B_y)^E \right] + \left[ (a_y \bar{v}_y B_x)^N + (a_y \bar{v}_y B_x)^S \right] + \left[ (d_x B_y)^E - (d_x B_y)^W \right] - \left[ (d_y B_x)^N - (d_y B_x)^S \right], \quad (39) $$

where $\bar{v}_{x,y}$ are the transverse velocities reconstructed from the interfaces, $a_{x,y}$ and $d_{x,y}$ are the flux and diffusion coefficients, and $B_{x,y}$ are the pointwise values of the magnetic field reconstructed from the corresponding interface. The other EMF components may be obtained by cyclic permutations. For the sake of briefness, we address the reader to find all the details of the different EMF averaging techniques that we employ in the present work in Mignone & Del Zanna [18]. Note that the employment of point value reconstruction methods (rather than 1D volume average-based ones) requires the electric field to be evaluated over 3 layers of ghost zones in the transverse direction rather than 4.

Finally, the line-averaged integrated electric field is computed as

$$ \bar{E}_{x,x_e} = E_{x,x_e} + \frac{\Delta^x E_{x,x_e}}{24}, \quad \bar{E}_{y,y_e} = E_{y,y_e} + \frac{\Delta^y E_{y,y_e}}{24}, \quad \bar{E}_{z,z_e} = E_{z,z_e} + \frac{\Delta^z E_{z,z_e}}{24}, \qquad (40) $$

and then employed to update the staggered magnetic field components using the discrete version of Stoke's theorem Eq. (14).

In the following we recap the steps of our 4th-order numerical scheme.

1. Start with the volume average of conserved quantities $\langle U \rangle_c$ and the face average of the staggered magnetic field components, $\hat{B}_{x,x_f}^x, \hat{B}_{y,y_f}^y, \hat{B}_{z,z_f}^z$.

2. Assign boundary conditions to these variables.

3. Obtain volume-averaged primitive variables by converting the volume averages $\langle V \rangle_c = V(\langle U \rangle_c)$. Notice that this array is 2nd-order accurate and it is only used in the limiting process, thus it does not replace the pointwise conversion later on.

4. Evaluate the $\theta_c$ parameter for order reduction in presence of discontinuities by means of either the Jameson's shock sensor (Eq. (35)) or by weighing up the contribution of higher-order derivatives (Eq. (36)) and flag troubled cells along with first neighbors.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Primitive variable profiles for the 3D rotated MHD shock tube problem at $t_f = t'_f \cos \alpha \cos \gamma$, where $t'_f = 0.2$ obtained with the WENOZ 4th-order scheme (red squares) on a uniform grid of 768 × 8 × 8 zones. All the variables are plotted along the non-rotated direction $x'$. On the left panels density (top) and pressure (bottom) are shown, while in the remaining panels, the three components of velocity (top) and magnetic field (bottom) are plotted. The black solid lines represent a reference 1D solution computed by employing a 2nd-order scheme with a resolution of $N_x = 16384$ grid points.</figcaption></figure>

5. Get pointwise value $U_c$ (Eq. (33)) at the cell center and $B_{x,f}$, $B_{y,f}$, $B_{z,f}$ (Eq. (21)) at the corresponding face centers. Also, during this step, interpolate the magnetic field from the face center to the cell center (Eq. (22)).

6. Convert pointwise conservative variables $U_c$ in primitive variables $V_c = \nabla(U_c)$.

7. Reconstruct L/R interface values (Eq. (23) or (24)) in primitive (or characteristic) variables.

8. Obtain upwind fluxes of hydro quantities (Eq. (29)) and average them over the corresponding faces (Eq. (30)). Concurrently, store the transverse velocities (e.g., $\bar{v}_y$ and $\bar{v}_z$ at an x-interface) as well as the coefficients $a^{L/R}$ and $d^{L/R}$ available with the chosen Riemann solver at a zone interface. These quantities are needed in the next step of the UCT scheme. For completeness we report the explicit form of these coefficients in Appendix B.

9. Reconstruct the transverse velocities obtained with the 1D Riemann solver and the normal components of the magnetic field from the face centers to the edges. This will yield the $N$, $S$ and $W$, $E$ values in Eq. (39). At an x-interface, for instance, one has to reconstruct $\bar{v}_{y,x_f}$ and $B_{x,x_f}$ along the y-direction, and then, again, $\bar{v}_{z,x_f}$ and $B_{x,x_f}$ along the z-direction.

10. Obtain the electric field at corners (Eq. (39)) by means of different averaging procedures (see Appendix B).

11. Line-average the electric field according to Eq. (40).

12. Add all contributions to the right-hand side and obtain the solution $\langle U \rangle_c$ at the next time level or RK stage (Eq. (31)).

## 5. Numerical benchmarks

We now test the capabilities of the novel method in a series of 2D and 3D numerical benchmarks. Each test is validated twice, by solving for the MHD and the RMHD case. As a reference, we will compare our results against the traditional 2nd- or 3rd-order schemes. As a default, Eq. (10) will be discretized by the 2nd-, 3rd-, or 4th-order schemes by means of a RK2, RK3, or RK4 time integrator. The WENOZ and MP5 spatial reconstructions will always use point values and, otherwise differently stated, each MHD and RMHD benchmark will employ, respectively, a HLLD solver coupled with the UCT-HLLD EMF average, and the GFORCE solver plus the UCT-GFORCE average that we suitably extended to relativistic frameworks.

Whenever the initial condition is a smooth analytic function our scheme integrates it by means of a Gaussian quadrature rule employing 4 Gaussian points.

We evaluate $L_1$ norm errors for a generic quantity $Q$ against a reference solution $Q_{\text{ref}}$ as $\epsilon_1(Q) = \sum_i |Q(x_i) - Q_{\text{ref}}(x_i)|/N$, where $x_i$ is a generic point of the domain and $N$ is the total number of grid zones. The convergence rate $O_{L_1}$ is thus calculated via a linear fit of $\log(\epsilon_1)$ against $\log(N_x)$, with $N_x$ being the linear resolution.

All computations have been carried out by means of the PLUTO astrophysical code [68,69], where the 4th-order method has been implemented.

### 5.1. 3D rotated shock tube

We begin this section by considering a 3D rotated version of a standard shock-tube problem. The initial condition is best expressed in the original (non-rotated) frame where $x'$ is the coordinate perpendicular to the plane of discontinuity. In the MHD case we have
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Primitive variable profiles for the 3D rotated MHD shock tube problem at $t_f = t'_f \cos \alpha \cos \gamma$, where $t'_f = 0.2$ obtained with the WENOZ 4th-order scheme (red squares) on a uniform grid of 768 × 8 × 8 zones. All the variables are plotted along the non-rotated direction $x'$. On the left panels density (top) and pressure (bottom) are shown, while in the remaining panels, the three components of velocity (top) and magnetic field (bottom) are plotted. The black solid lines represent a reference 1D solution computed by employing a 2nd-order scheme with a resolution of $N_x = 16384$ grid points.</figcaption></figure>

5. Get pointwise value $U_c$ (Eq. (33)) at the cell center and $B_{x,f}$, $B_{y,f}$, $B_{z,f}$ (Eq. (21)) at the corresponding face centers. Also, during this step, interpolate the magnetic field from the face center to the cell center (Eq. (22)).

6. Convert pointwise conservative variables $U_c$ in primitive variables $V_c = \nabla(U_c)$.

7. Reconstruct L/R interface values (Eq. (23) or (24)) in primitive (or characteristic) variables.

8. Obtain upwind fluxes of hydro quantities (Eq. (29)) and average them over the corresponding faces (Eq. (30)). Concurrently, store the transverse velocities (e.g., $\bar{v}_y$ and $\bar{v}_z$ at an x-interface) as well as the coefficients $a^{L/R}$ and $d^{L/R}$ available with the chosen Riemann solver at a zone interface. These quantities are needed in the next step of the UCT scheme. For completeness we report the explicit form of these coefficients in Appendix B.

9. Reconstruct the transverse velocities obtained with the 1D Riemann solver and the normal components of the magnetic field from the face centers to the edges. This will yield the $N$, $S$ and $W$, $E$ values in Eq. (39). At an x-interface, for instance, one has to reconstruct $\bar{v}_{y,x_f}$ and $B_{x,x_f}$ along the y-direction, and then, again, $\bar{v}_{z,x_f}$ and $B_{x,x_f}$ along the z-direction.

10. Obtain the electric field at corners (Eq. (39)) by means of different averaging procedures (see Appendix B).

11. Line-average the electric field according to Eq. (40).

12. Add all contributions to the right-hand side and obtain the solution $\langle U \rangle_c$ at the next time level or RK stage (Eq. (31)).

## 5. Numerical benchmarks

We now test the capabilities of the novel method in a series of 2D and 3D numerical benchmarks. Each test is validated twice, by solving for the MHD and the RMHD case. As a reference, we will compare our results against the traditional 2nd- or 3rd-order schemes. As a default, Eq. (10) will be discretized by the 2nd-, 3rd-, or 4th-order schemes by means of a RK2, RK3, or RK4 time integrator. The WENOZ and MP5 spatial reconstructions will always use point values and, otherwise differently stated, each MHD and RMHD benchmark will employ, respectively, a HLLD solver coupled with the UCT-HLLD EMF average, and the GFORCE solver plus the UCT-GFORCE average that we suitably extended to relativistic frameworks.

Whenever the initial condition is a smooth analytic function our scheme integrates it by means of a Gaussian quadrature rule employing 4 Gaussian points.

We evaluate $L_1$ norm errors for a generic quantity $Q$ against a reference solution $Q_{\text{ref}}$ as $\epsilon_1(Q) = \sum_i |Q(x_i) - Q_{\text{ref}}(x_i)|/N$, where $x_i$ is a generic point of the domain and $N$ is the total number of grid zones. The convergence rate $O_{L_1}$ is thus calculated via a linear fit of $\log(\epsilon_1)$ against $\log(N_x)$, with $N_x$ being the linear resolution.

All computations have been carried out by means of the PLUTO astrophysical code [68,69], where the 4th-order method has been implemented.

### 5.1. 3D rotated shock tube

We begin this section by considering a 3D rotated version of a standard shock-tube problem. The initial condition is best expressed in the original (non-rotated) frame where $x'$ is the coordinate perpendicular to the plane of discontinuity. In the MHD case we have
<figure><img src="image_4.png" /><figcaption>Fig. 3. Same of Fig. 2 but for the relativistic case. Here the numerical scheme adopted is using the MP5 algorithm and the “unrotated” final time is $t'_{f} = 0.55$.</figcaption></figure>

$$ egin{cases} V'_{L} = (1.08, 1.2, 0.01, 0.5, \frac{2}{\sqrt{4\pi}}, \frac{3.6}{\sqrt{4\pi}}, \frac{2}{\sqrt{4\pi}}, 0.95) & \text{for } x' < 0, \\ V'_{R} = (1, 0, 0, 0, \frac{2}{\sqrt{4\pi}}, \frac{4}{\sqrt{4\pi}}, \frac{2}{\sqrt{4\pi}}, 1) & \text{for } x' > 0, \end{cases} \qquad (41) $$

where $V' = (\rho, v'_{x}, v'_{y}, v'_{z}, B'_{x}, B'_{y}, B'_{z}, p)$ is an array of primitive quantities in the non-rotated frame. In the relativistic case we have, likewise,

$$ egin{cases} V'_{L} = (1.08, 0.4, 0.3, 0.2, 2, 0.3, 0.95) & \text{for } x' < 0, \\ V'_{R} = (1, -0.45, -0.2, 0.2, 2, -0.7, 0.5, 1) & \text{for } x' > 0. \end{cases} \qquad (42) $$

For each case, we rotate the initial condition by an angle $\alpha$ around the z-axis and then by an angle $\gamma$ around the y-axis. The coordinate transformation used for the 3D rotation is explained in detail in [15] (see §4.2.2 in that paper) and leads to a planar symmetry condition whereby, for any quantity $q(x)$, one must have $q(x + s) = q(x)$, where $s = (n_x \Delta x, n_y \Delta y, n_z \Delta z)$ corresponds to an integer shift of cells. We choose $\tan \alpha = 1/2$ and $\tan \gamma = \cos \alpha \tan \beta = -1/(2\sqrt{5})$ with $\tan \beta = -1/4$ corresponding to a translational invariance by an integer shift of cells, i.e., $(n_x, n_y, n_z) = (1, -2, 0)$ and $(n_x, n_y, n_z) = (1, 0, 4)$.

We solve the MHD / RMHD equations using our 4th-order method, adopting the WENOZ reconstruction in characteristic variables for the classical version and the MP5 reconstruction in primitive variables for the relativistic variant. The computational domain is the box $[-L/2, L/2] \times [0, L/96] \times [0, L/96]$ covered by a uniform grid of 768 × 8 × 8 zones. We evolve the solution until $t_f = t'_f \cos \alpha \cos \gamma = 4/\sqrt{21}$ where $t'_f = 0.2$ (for classical MHD), $t'_f = 0.55$ is the final time for the non-rotated (1D) versions of the problem, see, e.g., [70,71,15] for the classical MHD version and [72,58,73] (and references therein) for the relativistic variant. One dimensional plots of relevant fluid quantities are shown in the unrotated frame in Fig. (2) and (3) against a 1D high-resolution reference profile computed with a 2nd-order scheme with a resolution of $N_x = 16384$ grid points. The solution features a contact discontinuity in the middle ($x \approx 0$) separating two outermost fast magnetosonic shocks enclosing a pair of rotational waves followed by two slow waves. Profiles are correctly reproduced and oscillations in the normal component of magnetic field are of the same order as those of [15] and smaller than those in [71] who used a 2nd-order scheme. Our results demonstrate that our 4th-order method can correctly capture propagating discontinuous fronts although the order of the scheme degenerates to 1st-order.

## 5.2. 2D iso-density and iso-enthalpy vortex

We consider the MHD iso-density vortex advection problem (as introduced by Balsara [7] and later presented by Mignone [15] and Dumbser [74]) and also propose a new exact equilibrium solution of the ideal RMHD equations with constant enthalpy. For the MHD case, the initial condition consists of a magnetized vortex at constant density in force balance that propagates along the main diagonal of a computational domain of $[-5, 5] \times [-5, 5]$ under periodic boundary conditions. The equilibrium condition is best expressed in cylindrical coordinates $(R, \phi, z)$, from radial momentum balance

$$ \frac{d}{dR} \left( p + \frac{B_z^2}{2} \right) = \frac{\rho v_{\phi}^2}{R} - \frac{1}{2R^2} \frac{d}{dR} \left( R^2 B_{\phi}^2 \right). \qquad (43) $$

We assume constant density $\rho = 1$ and Eq. (43) is satisfied by adopting the following profiles:

$$ \mathbf{v} = \mathbf{v}_0 + (-y, x)ke^{q(1-R^2)}, \quad \mathbf{B} = (-y, x)\mu e^{q(1-R^2)} + B_z \hat{\mathbf{e}}_z, \qquad (44) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Same of Fig. 2 but for the relativistic case. Here the numerical scheme adopted is using the MP5 algorithm and the “unrotated” final time is $t'_{f} = 0.55$.</figcaption></figure>

$$ \left\{ \begin{aligned} V'_{L} &= \left( 1.08, 1.2, 0.01, 0.5, \frac{2}{\sqrt{4\pi}}, \frac{3.6}{\sqrt{4\pi}}, \frac{2}{\sqrt{4\pi}}, 0.95 \right) & \text{for } x' < 0, \\ V'_{R} &= \left( 1, 0, 0, 0, \frac{2}{\sqrt{4\pi}}, \frac{4}{\sqrt{4\pi}}, \frac{2}{\sqrt{4\pi}}, 1 \right) & \text{for } x' > 0, \end{aligned} \right. \qquad (41) $$

where $V' = (\rho, v'_{x}, v'_{y}, v'_{z}, B'_{x}, B'_{y}, B'_{z}, p)$ is an array of primitive quantities in the non-rotated frame. In the relativistic case we have, likewise,

$$ \left\{ \begin{aligned} V'_{L} &= (1.08, 0.4, 0.3, 0.2, 2, 0.3, 0.95) & \text{for } x' < 0, \\ V'_{R} &= (1, -0.45, -0.2, 0.2, 2, -0.7, 0.5, 1) & \text{for } x' > 0. \end{aligned} \right. \qquad (42) $$

For each case, we rotate the initial condition by an angle $\alpha$ around the z-axis and then by an angle $\gamma$ around the y-axis. The coordinate transformation used for the 3D rotation is explained in detail in [15] (see §4.2.2 in that paper) and leads to a planar symmetry condition whereby, for any quantity $q(x)$, one must have $q(x + s) = q(x)$, where $s = (n_x \Delta x, n_y \Delta y, n_z \Delta z)$ corresponds to an integer shift of cells. We choose $\tan \alpha = 1/2$ and $\tan \gamma = \cos \alpha \tan \beta = -1/(2\sqrt{5})$ with $\tan \beta = -1/4$ corresponding to a translational invariance by an integer shift of cells, i.e., $(n_x, n_y, n_z) = (1, -2, 0)$ and $(n_x, n_y, n_z) = (1, 0, 4)$.

We solve the MHD / RMHD equations using our 4th-order method, adopting the WENOZ reconstruction in characteristic variables for the classical version and the MP5 reconstruction in primitive variables for the relativistic variant. The computational domain is the box $[-L/2, L/2] \times [0, L/96] \times [0, L/96]$ covered by a uniform grid of 768 × 8 × 8 zones. We evolve the solution until $t_f = t'_f \cos \alpha \cos \gamma = 4/\sqrt{21}$ where $t'_f = 0.2$ (for classical MHD), $t'_f = 0.55$ is the final time for the non-rotated (1D) versions of the problem, see, e.g., [70,71,15] for the classical MHD version and [72,58,73] (and references therein) for the relativistic variant. One dimensional plots of relevant fluid quantities are shown in the unrotated frame in Fig. (2) and (3) against a 1D high-resolution reference profile computed with a 2nd-order scheme with a resolution of $N_x = 16384$ grid points. The solution features a contact discontinuity in the middle ($x \approx 0$) separating two outermost fast magnetosonic shocks enclosing a pair of rotational waves followed by two slow waves. Profiles are correctly reproduced and oscillations in the normal component of magnetic field are of the same order as those of [15] and smaller than those in [71] who used a 2nd-order scheme. Our results demonstrate that our 4th-order method can correctly capture propagating discontinuous fronts although the order of the scheme degenerates to 1st-order.

## 5.2. 2D iso-density and iso-enthalpy vortex

We consider the MHD iso-density vortex advection problem (as introduced by Balsara [7] and later presented by Mignone [15] and Dumbser [74]) and also propose a new exact equilibrium solution of the ideal RMHD equations with constant enthalpy. For the MHD case, the initial condition consists of a magnetized vortex at constant density in force balance that propagates along the main diagonal of a computational domain of $[-5, 5] \times [-5, 5]$ under periodic boundary conditions. The equilibrium condition is best expressed in cylindrical coordinates $(R, \phi, z)$, from radial momentum balance

$$ \frac{d}{dR} \left( p + \frac{B_z^2}{2} \right) = \frac{\rho v_{\phi}^2}{R} - \frac{1}{2R^2} \frac{d}{dR} \left( R^2 B_{\phi}^2 \right). \qquad (43) $$

We assume constant density $\rho = 1$ and Eq. (43) is satisfied by adopting the following profiles:

$$ \mathbf{v} = \mathbf{v}_0 + (-y, x)ke^{q(1-R^2)}, \quad \mathbf{B} = (-y, x)\mu e^{q(1-R^2)} + B_z \hat{\mathbf{e}}_z, \qquad (44) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Zoomed plot of the thermal pressure profile in the 2D MHD iso-density vortex (left panel) and density in the 2D RMHD iso-enthalpy vortex (right panel). In each panel the exact solution's profile (continuous black line) is compared with the profiles evolved by a 4th-order scheme adopting the WENOZ (red squares), the MP5 (blue circles), and the 2nd-order LINEAR scheme (green triangles). Filled symbols refer to data produced at a resolution of $N_x = 128$, while the empty ones are data at $N_x = 512$.</figcaption></figure>

<caption>Table 1</caption>

$L_1$ norm errors and corresponding convergence rates of the 2D MHD iso-density vortex (columns 2-5) and the 2D RMHD iso-enthalpy vortex (columns 6-9). Computations were carried out with the resolution of 32, 64, 128, 256 and 512 grid points along each direction.

<table><thead><tr><th rowspan="2">N<sub>x</sub></th><th colspan="2">MHD MP5</th><th colspan="2">MHD WENOZ</th><th colspan="2">RMHD MP5</th><th colspan="2">RMHD WENOZ</th></tr><tr><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th></tr></thead><tbody><tr><td>32</td><td>9.99E-04</td><td>-</td><td>9.01E-04</td><td>-</td><td>3.93E-04</td><td>-</td><td>3.11E-04</td><td>-</td></tr><tr><td>64</td><td>6.09E-05</td><td>4.04</td><td>6.00E-05</td><td>3.91</td><td>1.70E-05</td><td>4.53</td><td>1.83E-05</td><td>4.09</td></tr><tr><td>128</td><td>2.35E-06</td><td>4.70</td><td>2.37E-06</td><td>4.66</td><td>6.09E-07</td><td>4.80</td><td>6.41E-07</td><td>4.84</td></tr><tr><td>256</td><td>1.06E-07</td><td>4.47</td><td>1.06E-07</td><td>4.48</td><td>2.34E-08</td><td>4.71</td><td>2.36E-08</td><td>4.77</td></tr><tr><td>512</td><td>5.76E-09</td><td>4.20</td><td>5.76E-09</td><td>4.20</td><td>1.05E-09</td><td>4.47</td><td>1.05E-09</td><td>4.49</td></tr></tbody></table>

where $v_0$ is the translational velocity, $R = \sqrt{x^2 + y^2}$ is the cylindrical radius, $k = \mu = 1/2\pi$ and the parameter $q$ is set to 1 as in [74,15]. The pressure term on the left hand side of Eq. (43) is then readily obtained from

$$p_t = p + \frac{B_z^2}{2} = 1 + \frac{1}{4q} [\mu^2 (1 - 2qR^2) - fk^2] e^{2q(1-R^2)}, \quad (45)$$

where $f = \rho$ for the classical case. We equally distribute $p_t$ in Eq. (45) among the thermal and magnetic contributions, i.e. $p = \alpha p_t$ and $B_z^2/2 = (1 - \alpha)p_t$, where $\alpha \in [0, 1]$. In this setup, we use $\alpha = 1/2$ and $v_0 = (1, 1)$.

In the RMHD case (for which we also assume a purely radial dependence), the time-independent radial component of the momentum equation takes instead the form

$$\frac{dp}{dR} = \frac{w\gamma v_{\phi}^2}{R} + \frac{E_R}{R} \frac{d(RE_R)}{dR} - B_z \frac{dB_z}{dR} - \frac{B_{\phi}}{R} \frac{d(RB_{\phi})}{dR}, \quad (46)$$

where $w = 5$ is the constant enthalpy. We assume $B_z = E_R = 0$, so that Eq. (46) has solution given by Eq. (44) and (45) with $f = w$. Here we choose a static vortex ($v_0 = 0$). Density is recovered from the ideal equation of state (adiabatic index $\Gamma = 5/3$) as $\rho = w - 5p/2$, which is ensured to be strictly positive by our choice of the initial enthalpy.

Both the classical and relativistic cases are evolved until $t = 10$, which corresponds to one period for the MHD case. Table 1 reports the errors for $B_x$ measured in $L_1$ norm and the corresponding convergence rates both for the MHD and RMHD solutions by either using the WENOZ or the MP5 reconstructions.

Fig. 4 overlays the exact pressure profile (continuous black line) on the data produced with the 4th-order scheme adopting the WENOZ (red squares), the MP5 (blue circles) and the 2nd-order LINEAR scheme (green triangles). Filled symbols refer to data produced at a resolution of $N_x = 128$, while empty ones are produced at $N_x = 512$ and have been presented only for the LINEAR scheme, since high-order data do not display a sensible difference between the resolution of 128 and 512. The profiles of the 4th-order methods overlap almost perfectly on the exact solution at 128, while the 2nd-order data do not reach such accuracy neither at 128 nor at 512 grid points. Especially in the relativistic case, it is relevant to notice the significant amount of clipping introduced by the LINEAR scheme. On the other hand, the 4th-order method is capable of attaining significantly improved accuracy already at modest resolutions.

A close comparison with the results of the finite difference (FD) scheme of Mignone [15] reveals very similar errors for $N_x \lesssim 128$ even if the FD method, being 5th-order accurate, converges faster than the proposed FV method. We also point out that, owing to the employment of a 3rd-order Runge Kutta method, the time step of the FD scheme has to be re-scaled as $\Delta t \sim \Delta x^{5/3}$ in order to retain 5th-order accuracy (see Eq. [48] of that paper).
<figure><img src="image_6.png" /><figcaption>Fig. 4. Zoomed plot of the thermal pressure profile in the 2D MHD iso-density vortex (left panel) and density in the 2D RMHD iso-enthalpy vortex (right panel). In each panel the exact solution's profile (continuous black line) is compared with the profiles evolved by a 4th-order scheme adopting the WENOZ (red squares), the MP5 (blue circles), and the 2nd-order LINEAR scheme (green triangles). Filled symbols refer to data produced at a resolution of $N_x = 128$, while the empty ones are data at $N_x = 512$.</figcaption></figure>

<caption>Table 1</caption>

$L_1$ norm errors and corresponding convergence rates of the 2D MHD iso-density vortex (columns 2-5) and the 2D RMHD iso-enthalpy vortex (columns 6-9). Computations were carried out with the resolution of 32, 64, 128, 256 and 512 grid points along each direction.

<table><thead><tr><th rowspan="2">N<sub>x</sub></th><th colspan="2">MHD MP5</th><th colspan="2">MHD WENOZ</th><th colspan="2">RMHD MP5</th><th colspan="2">RMHD WENOZ</th></tr><tr><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>x</sub>)</th><th>O<sub>L<sub>1</sub></sub></th></tr></thead><tbody><tr><td>32</td><td>9.99E-04</td><td>-</td><td>9.01E-04</td><td>-</td><td>3.93E-04</td><td>-</td><td>3.11E-04</td><td>-</td></tr><tr><td>64</td><td>6.09E-05</td><td>4.04</td><td>6.00E-05</td><td>3.91</td><td>1.70E-05</td><td>4.53</td><td>1.83E-05</td><td>4.09</td></tr><tr><td>128</td><td>2.35E-06</td><td>4.70</td><td>2.37E-06</td><td>4.66</td><td>6.09E-07</td><td>4.80</td><td>6.41E-07</td><td>4.84</td></tr><tr><td>256</td><td>1.06E-07</td><td>4.47</td><td>1.06E-07</td><td>4.48</td><td>2.34E-08</td><td>4.71</td><td>2.36E-08</td><td>4.77</td></tr><tr><td>512</td><td>5.76E-09</td><td>4.20</td><td>5.76E-09</td><td>4.20</td><td>1.05E-09</td><td>4.47</td><td>1.05E-09</td><td>4.49</td></tr></tbody></table>

where $v_0$ is the translational velocity, $R = \sqrt{x^2 + y^2}$ is the cylindrical radius, $k = \mu = 1/2\pi$ and the parameter $q$ is set to 1 as in [74,15]. The pressure term on the left hand side of Eq. (43) is then readily obtained from

$$p_t = p + \frac{B_z^2}{2} = 1 + \frac{1}{4q} [\mu^2 (1 - 2qR^2) - fk^2] e^{2q(1-R^2)}, \quad (45)$$

where $f = \rho$ for the classical case. We equally distribute $p_t$ in Eq. (45) among the thermal and magnetic contributions, i.e. $p = \alpha p_t$ and $B_z^2/2 = (1 - \alpha)p_t$, where $\alpha \in [0, 1]$. In this setup, we use $\alpha = 1/2$ and $v_0 = (1, 1)$.

In the RMHD case (for which we also assume a purely radial dependence), the time-independent radial component of the momentum equation takes instead the form

$$\frac{dp}{dR} = \frac{w\gamma v_{\phi}^2}{R} + \frac{E_R}{R} \frac{d(RE_R)}{dR} - B_z \frac{dB_z}{dR} - \frac{B_{\phi}}{R} \frac{d(RB_{\phi})}{dR}, \quad (46)$$

where $w = 5$ is the constant enthalpy. We assume $B_z = E_R = 0$, so that Eq. (46) has solution given by Eq. (44) and (45) with $f = w$. Here we choose a static vortex ($v_0 = 0$). Density is recovered from the ideal equation of state (adiabatic index $\Gamma = 5/3$) as $\rho = w - 5p/2$, which is ensured to be strictly positive by our choice of the initial enthalpy.

Both the classical and relativistic cases are evolved until $t = 10$, which corresponds to one period for the MHD case. Table 1 reports the errors for $B_x$ measured in $L_1$ norm and the corresponding convergence rates both for the MHD and RMHD solutions by either using the WENOZ or the MP5 reconstructions.

Fig. 4 overlays the exact pressure profile (continuous black line) on the data produced with the 4th-order scheme adopting the WENOZ (red squares), the MP5 (blue circles) and the 2nd-order LINEAR scheme (green triangles). Filled symbols refer to data produced at a resolution of $N_x = 128$, while empty ones are produced at $N_x = 512$ and have been presented only for the LINEAR scheme, since high-order data do not display a sensible difference between the resolution of 128 and 512. The profiles of the 4th-order methods overlap almost perfectly on the exact solution at 128, while the 2nd-order data do not reach such accuracy neither at 128 nor at 512 grid points. Especially in the relativistic case, it is relevant to notice the significant amount of clipping introduced by the LINEAR scheme. On the other hand, the 4th-order method is capable of attaining significantly improved accuracy already at modest resolutions.

A close comparison with the results of the finite difference (FD) scheme of Mignone [15] reveals very similar errors for $N_x \lesssim 128$ even if the FD method, being 5th-order accurate, converges faster than the proposed FV method. We also point out that, owing to the employment of a 3rd-order Runge Kutta method, the time step of the FD scheme has to be re-scaled as $\Delta t \sim \Delta x^{5/3}$ in order to retain 5th-order accuracy (see Eq. [48] of that paper).
<figure><img src="image_6.png" /><figcaption>Fig. 5. L₁-norm errors (left panel) and numerical diffusion (right panel) for the 3D circularly polarized Alfvén wave test in the MHD regime (solid lines) and RMHD (dashed-dotted lines) as functions of the grid resolution Nₓ with 2<sup>nd</sup>- and 4<sup>th</sup>-order numerical schemes. The horizontal dotted line in the left panel is positioned at 10⁻⁴ in order to highlight the different resolution that achieves an accuracy of 10⁻⁴ with either a low or a high order scheme. The dashed lines give the ideal convergence slope (left panel) and a reference slope for the numerical diffusion (right panel).</figcaption></figure>

### 5.3. Circularly polarized Alfvén waves

We now consider the oblique propagation of circularly polarized Alfvén waves, which is an exact non-linear solution of the (R)MHD equations and was first presented by Tóth [6] in the context of classical MHD and later by Del Zanna [20] for the relativistic case.

In the MHD case we draw upon the setup of [15] by first constructing a 1D wave profile in the x-direction with angular frequency ω, wavenumber k and amplitude η as

$$ \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} = \begin{pmatrix} 0 \\ \eta \sin \phi \\ \eta \cos \phi \end{pmatrix}, \quad \begin{pmatrix} B_x \\ B_y \\ B_z \end{pmatrix} = \begin{pmatrix} c_a \sqrt{\rho} \\ \sqrt{\rho \eta \sin \phi} \\ \sqrt{\rho \eta \cos \phi} \end{pmatrix}, \qquad (47) $$

where φ = kx − ωt, cₐ = ω/k = 1 is the phase (Alfvén) velocity and η = 0.1 is the wave amplitude. Density and pressure are set to ρ = 1 and p = 2, respectively, while we employ the ideal EoS with Γ = 5/3. Such solution is then rotated by 45° around the y-axis in the 2D setup and by an additional 45° around the z-axis in the 3D case to have the wave propagating along the computational box's diagonal. The reader should refer to [15] for the details.

The RMHD case is set in a similar way, with a 1D wave profile defined by (see [20])

$$ \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} = \begin{pmatrix} 0 \\ -c_A \eta \cos \phi \\ -c_A \eta \sin \phi \end{pmatrix}, \quad \begin{pmatrix} B_x \\ B_y \\ B_z \end{pmatrix} = \begin{pmatrix} B_0 \\ B_0 \eta \cos \phi \\ B_0 \eta \sin \phi \end{pmatrix}, \qquad (48) $$

where we set B₀ = 1, η = 1, Γ = 4/3 and the Alfvén velocity is defined by

$$ c_A = \frac{B_0^2}{\rho h + B_0^2(1 + \eta^2)} \left[ \frac{1}{2} \left( 1 + \sqrt{1 - \left( \frac{2\eta B_0^2}{\rho h + B_0^2(1 + \eta^2)} \right)^2} \right) \right]^{-1} \qquad (49) $$

We then apply the same procedure as for the MHD case [15] to rotate the wave front.

The computational domain is the unit square (cube) in 2D (3D), while we follow the evolution of the wave for one period, after which its profile is supposed to return exactly on the initial values. We then estimate the error introduced by the numerical algorithms using the L₁ norm evaluated for the y-component of the magnetic field between the final and initial profiles.

Table 2 reports the errors and convergence rates evaluated for a series of 2D and 3D MHD and RMHD computations employing the MP5 and WENOZ reconstructions at increasing resolutions Nₓ = {8, 16, 32, 64, 128, 256, 512}. We use a CFL coefficient of 0.4 and 0.3 for the 2D and 3D case, respectively, and we impose periodic boundary conditions. The measured convergence rate O<sub>L₁</sub> well matches the expected value for all the configurations we considered (as shown in the first panel of Fig. 5). It is interesting to notice that, giving the same configuration, the L₁ norm errors are almost identical for the WENOZ or the MP5 schemes. This confirms that for such a smooth problem, the choice of the high-order spatial reconstruction algorithm does not change results significantly.

In the 3D MHD case, our results may be directly compared to those obtained with the 5th-order finite difference (FD) scheme of Mignone [15] (see Table 1 of that paper), indicating that our FV scheme yields smaller errors for grid sizes lower than 256³ zones. This difference lessens as the resolution grows (due to the additional order of accuracy of the FD scheme) so that similar errors are eventually found at the resolution of ~256³ zones, where ε₁(B_y) ~ 3.74 · 10⁻⁹ (for the FD scheme) while ε₁(B_y) ~ 2.91 · 10⁻⁹ (with the present method).

For problems without sharp gradients like the one considered here, the employment of a 4th-order (or higher) order scheme leads not only to an increase in accuracy, but also to a saving in computational time for fixed accuracy. It is therefore instructive to estimate
<figure><img src="image_7.png" /><figcaption>Fig. 5. L₁-norm errors (left panel) and numerical diffusion (right panel) for the 3D circularly polarized Alfvén wave test in the MHD regime (solid lines) and RMHD (dashed-dotted lines) as functions of the grid resolution Nₓ with 2<sup>nd</sup>- and 4<sup>th</sup>-order numerical schemes. The horizontal dotted line in the left panel is positioned at 10⁻⁴ in order to highlight the different resolution that achieves an accuracy of 10⁻⁴ with either a low or a high order scheme. The dashed lines give the ideal convergence slope (left panel) and a reference slope for the numerical diffusion (right panel).</figcaption></figure>

### 5.3. Circularly polarized Alfvén waves

We now consider the oblique propagation of circularly polarized Alfvén waves, which is an exact non-linear solution of the (R)MHD equations and was first presented by Tóth [6] in the context of classical MHD and later by Del Zanna [20] for the relativistic case.

In the MHD case we draw upon the setup of [15] by first constructing a 1D wave profile in the x-direction with angular frequency ω, wavenumber k and amplitude η as

$$ \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} = \begin{pmatrix} 0 \\ \eta \sin \phi \\ \eta \cos \phi \end{pmatrix}, \quad \begin{pmatrix} B_x \\ B_y \\ B_z \end{pmatrix} = \begin{pmatrix} c_a \sqrt{\rho} \\ \sqrt{\rho \eta \sin \phi} \\ \sqrt{\rho \eta \cos \phi} \end{pmatrix}, \qquad (47) $$

where φ = kx − ωt, cₐ = ω/k = 1 is the phase (Alfvén) velocity and η = 0.1 is the wave amplitude. Density and pressure are set to ρ = 1 and p = 2, respectively, while we employ the ideal EoS with Γ = 5/3. Such solution is then rotated by 45° around the y-axis in the 2D setup and by an additional 45° around the z-axis in the 3D case to have the wave propagating along the computational box's diagonal. The reader should refer to [15] for the details.

The RMHD case is set in a similar way, with a 1D wave profile defined by (see [20])

$$ \begin{pmatrix} v_x \\ v_y \\ v_z \end{pmatrix} = \begin{pmatrix} 0 \\ -c_A \eta \cos \phi \\ -c_A \eta \sin \phi \end{pmatrix}, \quad \begin{pmatrix} B_x \\ B_y \\ B_z \end{pmatrix} = \begin{pmatrix} B_0 \\ B_0 \eta \cos \phi \\ B_0 \eta \sin \phi \end{pmatrix}, \qquad (48) $$

where we set B₀ = 1, η = 1, Γ = 4/3 and the Alfvén velocity is defined by

$$ c_A = \frac{B_0^2}{\rho h + B_0^2(1 + \eta^2)} \left[ \frac{1}{2} \left( 1 + \sqrt{1 - \left( \frac{2\eta B_0^2}{\rho h + B_0^2(1 + \eta^2)} \right)^2} \right) \right]^{-1} \qquad (49) $$

We then apply the same procedure as for the MHD case [15] to rotate the wave front.

The computational domain is the unit square (cube) in 2D (3D), while we follow the evolution of the wave for one period, after which its profile is supposed to return exactly on the initial values. We then estimate the error introduced by the numerical algorithms using the L₁ norm evaluated for the y-component of the magnetic field between the final and initial profiles.

Table 2 reports the errors and convergence rates evaluated for a series of 2D and 3D MHD and RMHD computations employing the MP5 and WENOZ reconstructions at increasing resolutions Nₓ = {8, 16, 32, 64, 128, 256, 512}. We use a CFL coefficient of 0.4 and 0.3 for the 2D and 3D case, respectively, and we impose periodic boundary conditions. The measured convergence rate O<sub>L₁</sub> well matches the expected value for all the configurations we considered (as shown in the first panel of Fig. 5). It is interesting to notice that, giving the same configuration, the L₁ norm errors are almost identical for the WENOZ or the MP5 schemes. This confirms that for such a smooth problem, the choice of the high-order spatial reconstruction algorithm does not change results significantly.

In the 3D MHD case, our results may be directly compared to those obtained with the 5th-order finite difference (FD) scheme of Mignone [15] (see Table 1 of that paper), indicating that our FV scheme yields smaller errors for grid sizes lower than 256³ zones. This difference lessens as the resolution grows (due to the additional order of accuracy of the FD scheme) so that similar errors are eventually found at the resolution of ~256³ zones, where ε₁(B_y) ~ 3.74 · 10⁻⁹ (for the FD scheme) while ε₁(B_y) ~ 2.91 · 10⁻⁹ (with the present method).

For problems without sharp gradients like the one considered here, the employment of a 4th-order (or higher) order scheme leads not only to an increase in accuracy, but also to a saving in computational time for fixed accuracy. It is therefore instructive to estimate
<caption>Table 2</caption>

L₁ norm errors and corresponding convergence rates computed for the rotated two-dimensional (columns 2-5) and three-dimensional (columns 6-9) Alfvén wave propagation. Computations were carried out with the resolution of 8, 16, 32, 64, 128, 256, and 512 grid points along each direction.

<table><thead><tr><th rowspan="2">Scheme</th><th rowspan="2">N<sub>x</sub></th><th colspan="2">2D MP5</th><th colspan="2">2D WENOZ</th><th colspan="2">3D MP5</th><th colspan="2">3D WENOZ</th></tr><tr><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th></tr></thead><tbody><tr><td rowspan="6">MHD</td><td>8</td><td>2.76E-03</td><td>-</td><td>3.02E-03</td><td>-</td><td>3.61E-03</td><td>-</td><td>3.83E-03</td><td>-</td></tr><tr><td>16</td><td>1.47E-04</td><td>4.23</td><td>1.50E-04</td><td>4.33</td><td>2.05E-04</td><td>4.14</td><td>2.09E-04</td><td>4.19</td></tr><tr><td>32</td><td>8.39E-06</td><td>4.13</td><td>8.43E-06</td><td>4.16</td><td>1.24E-05</td><td>4.05</td><td>1.23E-05</td><td>4.08</td></tr><tr><td>64</td><td>5.07E-07</td><td>4.05</td><td>5.08E-07</td><td>4.05</td><td>7.57E-07</td><td>4.03</td><td>7.57E-07</td><td>4.03</td></tr><tr><td>128</td><td>3.13E-08</td><td>4.02</td><td>3.13E-08</td><td>4.02</td><td>4.69E-08</td><td>4.01</td><td>4.69E-08</td><td>4.01</td></tr><tr><td>256</td><td>1.94E-09</td><td>4.01</td><td>1.94E-09</td><td>4.01</td><td>2.92E-09</td><td>4.01</td><td>2.91E-09</td><td>4.00</td></tr><tr><td rowspan="6">RMHD</td><td>8</td><td>2.27E-02</td><td>-</td><td>2.60E-02</td><td>-</td><td>3.11E-02</td><td>-</td><td>3.48E-02</td><td>-</td></tr><tr><td>16</td><td>1.23E-03</td><td>4.20</td><td>1.29E-03</td><td>4.34</td><td>1.75E-03</td><td>4.15</td><td>1.79E-03</td><td>4.28</td></tr><tr><td>32</td><td>7.21E-05</td><td>4.10</td><td>7.25E-05</td><td>4.15</td><td>1.04E-04</td><td>4.07</td><td>1.04E-04</td><td>4.10</td></tr><tr><td>64</td><td>4.33E-06</td><td>4.06</td><td>4.33E-06</td><td>4.05</td><td>6.45E-06</td><td>4.01</td><td>6.45E-06</td><td>4.01</td></tr><tr><td>128</td><td>2.67E-07</td><td>4.02</td><td>2.67E-07</td><td>4.02</td><td>4.03E-07</td><td>4.00</td><td>4.03E-07</td><td>4.00</td></tr><tr><td>256</td><td>1.66E-08</td><td>4.01</td><td>1.66E-08</td><td>4.01</td><td>2.52E-08</td><td>4.00</td><td>2.52E-08</td><td>3.99</td></tr><tr><td></td><td>512</td><td>1.04E-09</td><td>4.00</td><td>1.04E-09</td><td>4.00</td><td>1.58E-09</td><td>4.00</td><td>1.58E-09</td><td>4.00</td></tr></tbody></table>

the net gain of using such approach. Let $\tau$ be the CPU time for updating a single cell with a chosen numerical scheme (e.g. 2nd or 4th-order). The total CPU time $T$ required to update $N^d$ cells ($d=1,2,3$ indicates the dimensionality of problem) scales as

$$T = \tau N^{d+1}, \qquad (50)$$

where the extra factor $N$ comes from the fact that the hyperbolic time step is inversely proportional to the grid resolution $N$. Now, let $\epsilon = CN^{-p}$ be the numerical error at the end of the computation, where $p$ is the order of the scheme and $C$ is an unknown constant depending, in general, on the chosen numerical method and the problem at hand. A 2nd- and 4th-order scheme will achieve the same accuracy $\epsilon_2 \sim \epsilon_4$ using, respectively, $N_2$ and $N_4$ grid points related by

$$N_4 \sim \left( \frac{C_4}{C_2} \right)^{1/4} \sqrt{N_2}. \qquad (51)$$

The previous expression demonstrates that a 4th-order scheme will require a reduced grid size of about $\sqrt{N_2}$. Furthermore, by expressing $N$ as a function of $T$ (from Eq. (50)) the previous relation can be adjusted in terms of CPU time:

$$\frac{T_4}{T_2} = \left( \frac{C_4}{C_2} \right)^{(d+1)/4} \frac{\tau_4}{\tau_2} \frac{1}{\sqrt{N_2^{d+1}}}. \qquad (52)$$

From our 3D tests we estimated $\tau_4/\tau_2 \sim \{4.3, 3.7\}$ and $C_4/C_2 \sim \{5.0, 8.7\}$ for the MHD and RMHD cases respectively, where we averaged the results produced with the MP5 and WENOZ spatial reconstructions. Plugging these parameters into Eq. (52) we obtain a net gain in computing time, for a given accuracy, that scales as $T_4/T_2 \sim K/N_2^2$ with $K \sim \{21, 32\}$. Our estimates show, for instance, that for a smooth MHD problem the accuracy obtained with a 2nd-order scheme on a grid with 512 × 512 × 512 points could be achieved by our high-order scheme with only $N_4 \sim 34$ per direction and a computing time reduced by a factor $\sim 1.2 \times 10^4$.

Another advantage of using more accurate schemes is the reduction of the numerical dissipation due to the round-off errors introduced by the grid discretization. We can estimate the diffusion intrinsic to the numerical scheme by looking at the damping of Alfvén waves (see [75]). The decay rate can be expressed as $\Gamma_A = \frac{1}{2}(\nu + \eta)k^2$, where $\nu$ is the fluid's kinematic viscosity and $\eta$ its resistivity. We then obtain an approximated value for $\Gamma_A$ by computing the quantity $\Delta \log(\delta B_z)/\Delta t$, where $\delta B_z = \sqrt{\iiint B_z^2 dx dy dz$ and the difference is evaluated after one period. From the right panel of Fig. 5 we can see how already at very small resolutions ($N \sim 8$) the high-order schemes produce an overall numerical dissipation which is more than one order of magnitude lower than the one delivered by the 2nd-order scheme. Moreover, the gap between the two schemes increases significantly at higher resolutions, as the decay rate scales as $\sim N^{-5}$ in the high-order setup, while we measured a slope of only $\sim -3$ in the low-order case.

## 5.4. Advection of a magnetic field loop

This problem consists of a weakly magnetized field loop advected by a uniform velocity field. Since the total pressure is dominated by the thermal contribution, the magnetic field is essentially transported as a passive scalar and the preservation of the initial circular shape tests the dissipative properties of the scheme [9, 71, 14, 15, 76].

Following Gardiner & Stone [9] and Mignone [15], the computational box is defined by $[-1, 1] \times [-0.5, 0.5]$ and discretized on $N_x \times N_x/2$ grid cells. Density and pressure are initially uniformly set to 1 in the MHD case, while $p=3$ for the RMHD setup (following Beckwith & Stone [76]). The velocity of the flow is given by
<caption>**Table 2**</caption>

L₁ norm errors and corresponding convergence rates computed for the rotated two-dimensional (columns 2-5) and three-dimensional (columns 6-9) Alfvén wave propagation. Computations were carried out with the resolution of 8, 16, 32, 64, 128, 256, and 512 grid points along each direction.

<table><thead><tr><th rowspan="2">Scheme</th><th rowspan="2">N<sub>x</sub></th><th colspan="2">2D MP5</th><th colspan="2">2D WENOZ</th><th colspan="2">3D MP5</th><th colspan="2">3D WENOZ</th></tr><tr><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th><th>ε<sub>1</sub>(B<sub>y</sub>)</th><th>O<sub>L<sub>1</sub></sub></th></tr></thead><tbody><tr><td rowspan="7">MHD</td><td>8</td><td>2.76E-03</td><td>-</td><td>3.02E-03</td><td>-</td><td>3.61E-03</td><td>-</td><td>3.83E-03</td><td>-</td></tr><tr><td>16</td><td>1.47E-04</td><td>4.23</td><td>1.50E-04</td><td>4.33</td><td>2.05E-04</td><td>4.14</td><td>2.09E-04</td><td>4.19</td></tr><tr><td>32</td><td>8.39E-06</td><td>4.13</td><td>8.43E-06</td><td>4.16</td><td>1.24E-05</td><td>4.05</td><td>1.23E-05</td><td>4.08</td></tr><tr><td>64</td><td>5.07E-07</td><td>4.05</td><td>5.08E-07</td><td>4.05</td><td>7.57E-07</td><td>4.03</td><td>7.57E-07</td><td>4.03</td></tr><tr><td>128</td><td>3.13E-08</td><td>4.02</td><td>3.13E-08</td><td>4.02</td><td>4.69E-08</td><td>4.01</td><td>4.69E-08</td><td>4.01</td></tr><tr><td>256</td><td>1.94E-09</td><td>4.01</td><td>1.94E-09</td><td>4.01</td><td>2.92E-09</td><td>4.01</td><td>2.91E-09</td><td>4.00</td></tr><tr><td>512</td><td>1.21E-10</td><td>4.00</td><td>1.21E-10</td><td>4.00</td><td>1.82E-10</td><td>4.00</td><td>1.82E-10</td><td>4.00</td></tr><tr><td rowspan="7">RMHD</td><td>8</td><td>2.27E-02</td><td>-</td><td>2.60E-02</td><td>-</td><td>3.11E-02</td><td>-</td><td>3.48E-02</td><td>-</td></tr><tr><td>16</td><td>1.23E-03</td><td>4.20</td><td>1.29E-03</td><td>4.34</td><td>1.75E-03</td><td>4.15</td><td>1.79E-03</td><td>4.28</td></tr><tr><td>32</td><td>7.21E-05</td><td>4.10</td><td>7.25E-05</td><td>4.15</td><td>1.04E-04</td><td>4.07</td><td>1.04E-04</td><td>4.10</td></tr><tr><td>64</td><td>4.33E-06</td><td>4.06</td><td>4.33E-06</td><td>4.05</td><td>6.45E-06</td><td>4.01</td><td>6.45E-06</td><td>4.01</td></tr><tr><td>128</td><td>2.67E-07</td><td>4.02</td><td>2.67E-07</td><td>4.02</td><td>4.03E-07</td><td>4.00</td><td>4.03E-07</td><td>4.00</td></tr><tr><td>256</td><td>1.66E-08</td><td>4.01</td><td>1.66E-08</td><td>4.01</td><td>2.52E-08</td><td>4.00</td><td>2.52E-08</td><td>3.99</td></tr><tr><td>512</td><td>1.04E-09</td><td>4.00</td><td>1.04E-09</td><td>4.00</td><td>1.58E-09</td><td>4.00</td><td>1.58E-09</td><td>4.00</td></tr></tbody></table>

the net gain of using such approach. Let $\tau$ be the CPU time for updating a single cell with a chosen numerical scheme (e.g. 2nd or 4th-order). The total CPU time $T$ required to update $N^d$ cells ($d=1,2,3$ indicates the dimensionality of problem) scales as

$$T = \tau N^{d+1}, \qquad (50)$$

where the extra factor $N$ comes from the fact that the hyperbolic time step is inversely proportional to the grid resolution $N$. Now, let $\epsilon = CN^{-p}$ be the numerical error at the end of the computation, where $p$ is the order of the scheme and $C$ is an unknown constant depending, in general, on the chosen numerical method and the problem at hand. A 2nd- and 4th-order scheme will achieve the same accuracy $\epsilon_2 \sim \epsilon_4$ using, respectively, $N_2$ and $N_4$ grid points related by

$$N_4 \sim \left( \frac{C_4}{C_2} \right)^{1/4} \sqrt{N_2}. \qquad (51)$$

The previous expression demonstrates that a 4th-order scheme will require a reduced grid size of about $\sqrt{N_2}$. Furthermore, by expressing $N$ as a function of $T$ (from Eq. (50)) the previous relation can be adjusted in terms of CPU time:

$$\frac{T_4}{T_2} = \left( \frac{C_4}{C_2} \right)^{(d+1)/4} \frac{\tau_4}{\tau_2} \frac{1}{\sqrt{N_2^{d+1}}}. \qquad (52)$$

From our 3D tests we estimated $\tau_4/\tau_2 \sim \{4.3, 3.7\}$ and $C_4/C_2 \sim \{5.0, 8.7\}$ for the MHD and RMHD cases respectively, where we averaged the results produced with the MP5 and WENOZ spatial reconstructions. Plugging these parameters into Eq. (52) we obtain a net gain in computing time, for a given accuracy, that scales as $T_4/T_2 \sim K/N_2^2$ with $K \sim \{21, 32\}$. Our estimates show, for instance, that for a smooth MHD problem the accuracy obtained with a 2nd-order scheme on a grid with 512 × 512 × 512 points could be achieved by our high-order scheme with only $N_4 \sim 34$ per direction and a computing time reduced by a factor $\sim 1.2 \times 10^4$.

Another advantage of using more accurate schemes is the reduction of the numerical dissipation due to the round-off errors introduced by the grid discretization. We can estimate the diffusion intrinsic to the numerical scheme by looking at the damping of Alfvén waves (see [75]). The decay rate can be expressed as $\Gamma_A = \frac{1}{2}(\nu + \eta)k^2$, where $\nu$ is the fluid's kinematic viscosity and $\eta$ its resistivity. We then obtain an approximated value for $\Gamma_A$ by computing the quantity $\Delta \log(\delta B_z)/\Delta t$, where $\delta B_z = \sqrt{\iiint B_z^2 dx dy dz$ and the difference is evaluated after one period. From the right panel of Fig. 5 we can see how already at very small resolutions ($N \sim 8$) the high-order schemes produce an overall numerical dissipation which is more than one order of magnitude lower than the one delivered by the 2nd-order scheme. Moreover, the gap between the two schemes increases significantly at higher resolutions, as the decay rate scales as $\sim N^{-5}$ in the high-order setup, while we measured a slope of only $\sim -3$ in the low-order case.

## 5.4. Advection of a magnetic field loop

This problem consists of a weakly magnetized field loop advected by a uniform velocity field. Since the total pressure is dominated by the thermal contribution, the magnetic field is essentially transported as a passive scalar and the preservation of the initial circular shape tests the dissipative properties of the scheme [9, 71, 14, 15, 76].

Following Gardiner & Stone [9] and Mignone [15], the computational box is defined by $[-1, 1] \times [-0.5, 0.5]$ and discretized on $N_x \times N_x/2$ grid cells. Density and pressure are initially uniformly set to 1 in the MHD case, while $p=3$ for the RMHD setup (following Beckwith & Stone [76]). The velocity of the flow is given by
<figure><img src="image_7.png" /><figcaption>Fig. 6. Normalized magnetic energy density for the 2D field loop after two periods at the resolution of 256 × 128 grid points in the classical (top) and relativistic (bottom) MHD. Simulations compare LINEAR (leftmost panel), WENO3 (center left panels), MP5 (center right panels) and WENOZ reconstructions (rightmost panels). Magnetic field lines are overlaid as black lines.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Plot of the normalized magnetic energy density decay of a 2D field loop advected over time at a resolution of 256 × 128 grid points. Each simulation has been performed with a 2nd-order LINEAR (green), a 3rd-order WENO3+RK3 (cyan), a 4th-order MP5 (blue), and a 4th-order WENOZ scheme (red). The top panels show the results obtained in the classical MHD regime, with the HLL (top left panel) and HLLD (top right panel) Riemann solvers. Similarly, the bottom panels show the results obtained in the relativistic MHD regime, with the HLL (bottom left panel) and GFORCE (bottom right panel) Riemann solvers. Notice that for the results obtained with the HLLD and GFORCE Riemann solvers, the time axes have been reversed.</figcaption></figure>

$$v = V_0(\cos \alpha, \sin \alpha), \qquad (53)$$

with $V_0 = \sqrt{5}$ in the classical case while $V_0 = \sqrt{5}/10$ for the relativistic case. The inclination angle is set to $\sin \alpha = 1/\sqrt{5}$, $\cos \alpha = 2/\sqrt{5}$. The magnetic field is defined by means of the vector potential as

$$A_z = \begin{cases} A_0(R-r) & \text{if } R_1 < r \le R, \\ 0 & \text{if elsewhere,} \end{cases} \qquad (54)$$

with $A_0 = 10^{-3}$, $R = 0.3$, $R_1 = 0.2R$, $r = \sqrt{x^2 + y^2}$. Simulations have been performed on a computational domain of 256 × 128 grid points with the LINEAR, WENO3, MP5 and WENOZ schemes up to $t_f = 2$ units for the MHD setup and $t_f = 20$ for the RMHD setup (both corresponding to two advection periods). The CFL was set to 0.4 and double periodic boundary conditions are imposed.

Fig. 6 shows the normalized magnetic energy density together with magnetic field lines with a resolution of 256 × 128 grid points for each configuration. As shown in both sequences of Fig. 6, smearing and diffusion are significantly reduced as we gradually increase the accuracy from 2nd to 4th-order. It is interesting to notice that the WENOZ algorithm performs the best in this problem while MP5 has some small amplitude overshoots (undershoots) at the loop tail (head), most likely caused by a small clipping of its limiting procedure.

Another extremely crucial aspect is the ability of the algorithm to conserve the energy of the system. Fig. 7 compares the normalized magnetic energy density over time for a set of 4 simulations employing a 2nd- (green), a 3rd- (cyan), and both the 4th-
<figure><img src="image_8.png" /><figcaption>Fig. 6. Normalized magnetic energy density for the 2D field loop after two periods at the resolution of 256 × 128 grid points in the classical (top) and relativistic (bottom) MHD. Simulations compare LINEAR (leftmost panel), WENO3 (center left panels), MP5 (center right panels) and WENOZ reconstructions (rightmost panels). Magnetic field lines are overlaid as black lines.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Plot of the normalized magnetic energy density decay of a 2D field loop advected over time at a resolution of 256 × 128 grid points. Each simulation has been performed with a 2nd-order LINEAR (green), a 3rd-order WENO3+RK3 (cyan), a 4th-order MP5 (blue), and a 4th-order WENOZ scheme (red). The top panels show the results obtained in the classical MHD regime, with the HLL (top left panel) and HLLD (top right panel) Riemann solvers. Similarly, the bottom panels show the results obtained in the relativistic MHD regime, with the HLL (bottom left panel) and GFORCE (bottom right panel) Riemann solvers. Notice that for the results obtained with the HLLD and GFORCE Riemann solvers, the time axes have been reversed.</figcaption></figure>

$$v = V_0(\cos \alpha, \sin \alpha), \qquad (53)$$

with $V_0 = \sqrt{5}$ in the classical case while $V_0 = \sqrt{5}/10$ for the relativistic case. The inclination angle is set to $\sin \alpha = 1/\sqrt{5}$, $\cos \alpha = 2/\sqrt{5}$. The magnetic field is defined by means of the vector potential as

$$A_z = \begin{cases} A_0(R-r) & \text{if } R_1 < r \le R, \\ 0 & \text{if elsewhere,} \end{cases} \qquad (54)$$

with $A_0 = 10^{-3}$, $R = 0.3$, $R_1 = 0.2R$, $r = \sqrt{x^2 + y^2}$. Simulations have been performed on a computational domain of 256 × 128 grid points with the LINEAR, WENO3, MP5 and WENOZ schemes up to $t_f = 2$ units for the MHD setup and $t_f = 20$ for the RMHD setup (both corresponding to two advection periods). The CFL was set to 0.4 and double periodic boundary conditions are imposed.

Fig. 6 shows the normalized magnetic energy density together with magnetic field lines with a resolution of 256 × 128 grid points for each configuration. As shown in both sequences of Fig. 6, smearing and diffusion are significantly reduced as we gradually increase the accuracy from 2nd to 4th-order. It is interesting to notice that the WENOZ algorithm performs the best in this problem while MP5 has some small amplitude overshoots (undershoots) at the loop tail (head), most likely caused by a small clipping of its limiting procedure.

Another extremely crucial aspect is the ability of the algorithm to conserve the energy of the system. Fig. 7 compares the normalized magnetic energy density over time for a set of 4 simulations employing a 2nd- (green), a 3rd- (cyan), and both the 4th-
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Colored maps for the density of the current sheet test problem for the classical (top panel) and relativistic (bottom panel) MHD equations. From left to right snapshots at, respectively, $t = 5, 7.5, 10$ are shown at a resolution of $600 \times 600$ grid points. Simulations are carried out with the HLLC Riemann solver in the non-relativistic case and with the GFORCE solver in the relativistic case.</figcaption></figure>

order accurate (blue and red) schemes, respectively, in the MHD (top panels) and RMHD (bottom panels) regimes. The top panels compare the results of the HLL + UCT-HLL solver combination with those of the HLLD + UCT-HLLD pair, where the time axis has been reversed to better compare the two sets of simulations. The RMHD case follows along but with the GFORCE + UCT-GFORCE solver, showing improved results, although the choice of the Riemann solver seems to affect only marginally the conservation of the magnetic energy, with the GFORCE achieving better performance. The magnetic energy decay is considerably faster for the 2^nd^-order scheme compared to the 3^rd^- and 4^th^-order ones, with the latter yielding the best performance.

Very similar results are also recovered with the FD scheme of Mignone [15] (see left panel in Fig. A.9 of that paper) when using either WENOZ or MP5, which dissipate just ~5% of total magnetic energy at the end of the computation.

## 5.5. 2D current sheet

This problem was first introduced by Hawley & Stone [77] and also tested by Núñez-de La Rosa & Munz both in a MHD [34] and RMHD [35] regime (see also references therein).

We use the computational box $[-1, 1] \times [-0.5, 0.5]$ discretized with $600 \times 600$ grid points. Density and pressure of an ideal, adiabatic gas with $\Gamma = 5/3$ are initially set to the constant values $\rho = 1$, $p = 0.5\beta$, with $\beta$ being the ratio between the thermal and magnetic pressure of the gas. The magnetic field initially lies in the $y$-direction and it is discontinuous with $B_y = 1$ for $|x| > 0.25$ and $B_y = -1$ otherwise. We set $v_x = \sin(2\pi y)$, $v_y = v_z = 0$, where $A$ is the amplitude of the perturbation. The parameters $A$ and $\beta$ are both set to the challenging values of $A = \beta = 0.1$ to test the robustness of the underlying scheme. Simulations have been run up to $t_f = 10$, with a CFL restriction of 0.4 and periodic boundary conditions.

The harmonic velocity perturbation produces non-linearly polarized Alfvén waves that grow into magnetosonic waves, triggering magnetic reconnection in the regions where the field flips sign. When $\beta < 1$, namely, in the regime of strongly magnetized fluid, magnetic reconnection forms overpressurized regions that develop magnetosonic waves orthogonal to the field. This leads to the conversion from magnetic to thermal energy and to the formation of steep gradients in the regions nearby the field. For this reason, this benchmark can be extremely severe in assessing the robustness of the underlying numerical scheme.

In Fig. 8 we display density snapshots at $t = 5, 7.5, 10$ for the MHD (top panels) and RMHD setups (bottom panels) obtained, respectively, with the HLLC + UCT-HLL, and with the GFORCE + UCT-GFORCE as solver plus EMF average combination. This highly non-linear problem tests our novel limiting strategy where the WENOZ scheme lowers to WENO3 whenever the order reduction procedure is activated, thus increasing the stability range of the scheme. A similar fallback strategy has been adopted by [34,35] who used a simpler Rusanov-Lax Friedrichs solver.

## 5.6. Orszag-Tang vortex

The Orszag-Tang MHD vortex [78] has traditionally been a standard numerical benchmark for inter-scheme comparison. The test consists of an initial vorticity distribution that spins the fluid counter-clockwise, leading eventually to the formation of a complex structure that includes shocks and turbulence. The initial density and pressure are set to the constant values $\rho_0 = 25/9$ and $p_0 = 5/3$, while the initial velocity and magnetic fields are $(v_x, v_y) = (-\sin y, \sin x)$ and $(B_x, B_y) = (-\sin y, \sin 2x)$, respectively, in the MHD setup, and $(v_x, v_y) = (-0.5 \sin y, 0.5 \sin x)$ in the RMHD setup (as presented in [79,80,76]). We ran computations up to $t_f = 2\pi$ with $N_x = N_y = \{64, 512\}$ grid zones for the MHD case, and $N_x = N_y = \{128, 1024\}$ for the RMHD case using periodic boundary conditions
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Colored maps for the density of the current sheet test problem for the classical (top panel) and relativistic (bottom panel) MHD equations. From left to right snapshots at, respectively, $t = 5, 7.5, 10$ are shown at a resolution of $600 \times 600$ grid points. Simulations are carried out with the HLLC Riemann solver in the non-relativistic case and with the GFORCE solver in the relativistic case.</figcaption></figure>

order accurate (blue and red) schemes, respectively, in the MHD (top panels) and RMHD (bottom panels) regimes. The top panels compare the results of the HLL + UCT-HLL solver combination with those of the HLLD + UCT-HLLD pair, where the time axis has been reversed to better compare the two sets of simulations. The RMHD case follows along but with the GFORCE + UCT-GFORCE solver, showing improved results, although the choice of the Riemann solver seems to affect only marginally the conservation of the magnetic energy, with the GFORCE achieving better performance. The magnetic energy decay is considerably faster for the 2^nd^-order scheme compared to the 3^rd^- and 4^th^-order ones, with the latter yielding the best performance.

Very similar results are also recovered with the FD scheme of Mignone [15] (see left panel in Fig. A.9 of that paper) when using either WENOZ or MP5, which dissipate just ~5% of total magnetic energy at the end of the computation.

## 5.5. 2D current sheet

This problem was first introduced by Hawley & Stone [77] and also tested by Núñez-de La Rosa & Munz both in a MHD [34] and RMHD [35] regime (see also references therein).

We use the computational box $[-1, 1] \times [-0.5, 0.5]$ discretized with $600 \times 600$ grid points. Density and pressure of an ideal, adiabatic gas with $\Gamma = 5/3$ are initially set to the constant values $\rho = 1$, $p = 0.5\beta$, with $\beta$ being the ratio between the thermal and magnetic pressure of the gas. The magnetic field initially lies in the $y$-direction and it is discontinuous with $B_y = 1$ for $|x| > 0.25$ and $B_y = -1$ otherwise. We set $v_x = \text{Asin}(2\pi y)$, $v_y = v_z = 0$, where $A$ is the amplitude of the perturbation. The parameters $A$ and $\beta$ are both set to the challenging values of $A = \beta = 0.1$ to test the robustness of the underlying scheme. Simulations have been run up to $t_f = 10$, with a CFL restriction of 0.4 and periodic boundary conditions.

The harmonic velocity perturbation produces non-linearly polarized Alfvén waves that grow into magnetosonic waves, triggering magnetic reconnection in the regions where the field flips sign. When $\beta < 1$, namely, in the regime of strongly magnetized fluid, magnetic reconnection forms overpressurized regions that develop magnetosonic waves orthogonal to the field. This leads to the conversion from magnetic to thermal energy and to the formation of steep gradients in the regions nearby the field. For this reason, this benchmark can be extremely severe in assessing the robustness of the underlying numerical scheme.

In Fig. 8 we display density snapshots at $t = 5, 7.5, 10$ for the MHD (top panels) and RMHD setups (bottom panels) obtained, respectively, with the HLLC + UCT-HLL, and with the GFORCE + UCT-GFORCE as solver plus EMF average combination. This highly non-linear problem tests our novel limiting strategy where the WENOZ scheme lowers to WENO3 whenever the order reduction procedure is activated, thus increasing the stability range of the scheme. A similar fallback strategy has been adopted by [34,35] who used a simpler Rusanov-Lax Friedrichs solver.

## 5.6. Orszag-Tang vortex

The Orszag-Tang MHD vortex [78] has traditionally been a standard numerical benchmark for inter-scheme comparison. The test consists of an initial vorticity distribution that spins the fluid counter-clockwise, leading eventually to the formation of a complex structure that includes shocks and turbulence. The initial density and pressure are set to the constant values $\rho_0 = 25/9$ and $p_0 = 5/3$, while the initial velocity and magnetic fields are $(v_x, v_y) = (-\sin y, \sin x)$ and $(B_x, B_y) = (-\sin y, \sin 2x)$, respectively, in the MHD setup, and $(v_x, v_y) = (-0.5 \sin y, 0.5 \sin x)$ in the RMHD setup (as presented in [79,80,76]). We ran computations up to $t_f = 2\pi$ with $N_x = N_y = \{64, 512\}$ grid zones for the MHD case, and $N_x = N_y = \{128, 1024\}$ for the RMHD case using periodic boundary conditions
<figure><img src="image_10.png" /><figcaption>Fig. 9. Thermal pressure of the non-relativistic Orszag-Tang vortex computed at $t = 2\pi$ with different algorithms and resolutions. Results are obtained with a 2^nd^-order LINEAR (left column), a 4^th^-order MP5 (middle column), and a 4^th^-order WENOZ scheme (right column). Each simulation has been computed with a Roe (1^st^ and 3^rd^ rows) and HLLD (2^nd^ and 4^th^ rows) Riemann solver with a grid resolution of $N_x = N_y = 64$ (1^st^ and 2^nd^ rows) and $N_x = N_y = 512$ (3^rd^ and 4^th^ rows).</figcaption></figure>

and a CFL number of 0.5. We repeated the computations by choosing two different Riemann solvers in the MHD case (Roe and HLLD, always with the UCT-HLLD average), and only the GFORCE Riemann solver coupled to the UCT-GFORCE EMF averaging algorithm in the relativistic case. Likewise, we also compared 3 different reconstruction schemes, namely, LINEAR, MP5, and WENOZ.

Multiple shock-vortex interactions regulate the dynamics of the system, leading to the formation of an inclined current sheet at the center of the domain. The magnetic energy is gradually dissipated and the current sheet twists, leading to the structures observed in the thermal pressure distribution shown in Fig. 9. For the MHD setup, even at modest resolutions (64 points) the 4^th^-order scheme with MP5 reconstruction and the HLLD Riemann solver allows to distinguish at the center of the computational domain the formation of a magnetic island (an O-point). It is well known (see [81,82]) that the tearing instability is expected to develop on the ideal (Alfvénic) timescales of an ideal MHD simulation only if the numerical dissipation is sufficiently low. At higher resolutions, all simulations produce high-pressure regions, albeit small-scale structures are best appreciated with the high-order runs.

Similarly to Mignone & Del Zanna [18], we can assess more quantitatively the improvements introduced by our 4^th^-order scheme by computing the parameter $\mu_p = \max(p)/\langle p \rangle$, i.e. the pressure maximum normalized by the average pressure at the end of a run. As this quantity generally increases with better resolution, we use it as a proxy to assess the decrease of numerical dissipation. For instance, the values of $\mu_p$ obtained with the MP5 reconstruction and HLLD solver against those obtained with the LINEAR one yield an average increase in effective resolution of ~2.23, which leads to an average speedup of ~2.58 (taking into account the extra computational cost of the 4^th^-order scheme estimated in §5.3). This is, however, only a rough estimate that cannot fully characterize the impact of our scheme on the typical numerical dissipation length-scale in a turbulent plasma, for which more dedicated tests and analysis are required.

For the RMHD case, the final stage of the simulation produces a central low-pressure cavity with a thin high-pressure filament, which exhibits finer structures for less dissipative schemes (Fig. 10). While at a resolution of 128 points there are no apparent deviations between the different setups, at higher resolution (1024 points) the numerical resistivity from the 4^th^-order scheme is low enough to lead to the production of multiple O-points in the high-pressure filament, although they remain still quite unresolved. Our results confirm the overall robustness of the scheme as well as its ability to capture small-scale structures that would otherwise be lost with a lower order scheme.

### 5.7. 3D blast wave

We finally consider the 3D magnetized blast wave problem in order to assess the robustness of the algorithms under strong magnetization conditions. Despite its simplicity, the blast wave problem is a particularly effective benchmark in testing the solver's ability at handling MHD wave degeneracies parallel and perpendicularly to the field orientation.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Thermal pressure of the non-relativistic Orszag-Tang vortex computed at $t = 2\pi$ with different algorithms and resolutions. Results are obtained with a 2^nd^-order LINEAR (left column), a 4^th^-order MP5 (middle column), and a 4^th^-order WENOZ scheme (right column). Each simulation has been computed with a Roe (1^st^ and 3^rd^ rows) and HLLD (2^nd^ and 4^th^ rows) Riemann solver with a grid resolution of $N_x = N_y = 64$ (1^st^ and 2^nd^ rows) and $N_x = N_y = 512$ (3^rd^ and 4^th^ rows).</figcaption></figure>

and a CFL number of 0.5. We repeated the computations by choosing two different Riemann solvers in the MHD case (Roe and HLLD, always with the UCT-HLLD average), and only the GFORCE Riemann solver coupled to the UCT-GFORCE EMF averaging algorithm in the relativistic case. Likewise, we also compared 3 different reconstruction schemes, namely, LINEAR, MP5, and WENOZ.

Multiple shock-vortex interactions regulate the dynamics of the system, leading to the formation of an inclined current sheet at the center of the domain. The magnetic energy is gradually dissipated and the current sheet twists, leading to the structures observed in the thermal pressure distribution shown in Fig. 9. For the MHD setup, even at modest resolutions (64 points) the 4^th^-order scheme with MP5 reconstruction and the HLLD Riemann solver allows to distinguish at the center of the computational domain the formation of a magnetic island (an O-point). It is well known (see [81,82]) that the tearing instability is expected to develop on the ideal (Alfvénic) timescales of an ideal MHD simulation only if the numerical dissipation is sufficiently low. At higher resolutions, all simulations produce high-pressure regions, albeit small-scale structures are best appreciated with the high-order runs.

Similarly to Mignone & Del Zanna [18], we can assess more quantitatively the improvements introduced by our 4^th^-order scheme by computing the parameter $\mu_p = \max(p)/\langle p \rangle$, i.e. the pressure maximum normalized by the average pressure at the end of a run. As this quantity generally increases with better resolution, we use it as a proxy to assess the decrease of numerical dissipation. For instance, the values of $\mu_p$ obtained with the MP5 reconstruction and HLLD solver against those obtained with the LINEAR one yield an average increase in effective resolution of ~2.23, which leads to an average speedup of ~2.58 (taking into account the extra computational cost of the 4^th^-order scheme estimated in §5.3). This is, however, only a rough estimate that cannot fully characterize the impact of our scheme on the typical numerical dissipation length-scale in a turbulent plasma, for which more dedicated tests and analysis are required.

For the RMHD case, the final stage of the simulation produces a central low-pressure cavity with a thin high-pressure filament, which exhibits finer structures for less dissipative schemes (Fig. 10). While at a resolution of 128 points there are no apparent deviations between the different setups, at higher resolution (1024 points) the numerical resistivity from the 4^th^-order scheme is low enough to lead to the production of multiple O-points in the high-pressure filament, although they remain still quite unresolved. Our results confirm the overall robustness of the scheme as well as its ability to capture small-scale structures that would otherwise be lost with a lower order scheme.

### 5.7. 3D blast wave

We finally consider the 3D magnetized blast wave problem in order to assess the robustness of the algorithms under strong magnetization conditions. Despite its simplicity, the blast wave problem is a particularly effective benchmark in testing the solver's ability at handling MHD wave degeneracies parallel and perpendicularly to the field orientation.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Thermal pressure of the relativistic Orszag-Tang vortex computed at $t = 2\pi$ with different algorithms and resolutions. From left to right, results are obtained with, respectively, a 2nd-order LINEAR (left column), a 4th-order MP5 (middle column), and a 4th-order WENOZ scheme (bottom column). The grid resolution adopted is $N_x = N_y = 128$ (top row) and $N_x = N_y = 1024$ (bottom row). The presence or absence of O-points emerging from the central current sheets (only for the high resolution cases) is highlighted by a zoom on the central region.</figcaption></figure>

<table><caption>**Table 3**
Parameters used in the 3D blast wave problem.</caption><thead><tr><th>Case</th><th>ρ<sub>in</sub></th><th>ρ<sub>out</sub></th><th>p<sub>in</sub></th><th>p<sub>out</sub></th><th>B<sub>0</sub></th><th>Γ</th><th>r<sub>0</sub></th><th>L</th><th>t<sub>stop</sub></th><th>Solver</th><th>En. Corr.</th></tr></thead><tbody><tr><td>C1 (MHD)</td><td>1</td><td>1</td><td>10</td><td>0.1</td><td>3</td><td>5/3</td><td>0.1</td><td>1</td><td>0.2</td><td>HLL</td><td>NO</td></tr><tr><td>C2 (MHD)</td><td>1</td><td>1</td><td>10<sup>3</sup></td><td>0.1</td><td>100/√4π</td><td>1.4</td><td>0.1</td><td>1</td><td>0.01</td><td>HLL</td><td>YES</td></tr><tr><td>R1 (RMHD)</td><td>10<sup>-2</sup></td><td>10<sup>-4</sup></td><td>1</td><td>5 × 10<sup>-3</sup></td><td>1</td><td>4/3</td><td>0.08</td><td>6</td><td>4</td><td>GFORCE</td><td>YES</td></tr></tbody></table>

The initial condition consists of a uniform medium with density and pressure values set, respectively, to $ρ_{out}$ and $p_{out}$ and threaded by an oblique constant magnetic field,

$$ \mathbf{B} = B_0 \left( \sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta \right) \qquad (55) $$

where $θ = π/2$ and $φ = π/4$ are used in what follows. Note also that oblique configurations make the computation considerably more challenging than in the grid-aligned cases (see, for instance the discussion in §4.4 of Mattia & Mignone [83]). The computational domain is a cube of length $L$ centered at the origin with zero-gradient boundary conditions. An over-pressurized spherical region where pressure and density can have different values, namely $ρ_{in}$ and $p_{in}$, is initially set for $r < r_0$.

We examine three different cases with parameters listed in Table 3. Cases C1 is taken from Felker & Stone [36], with the exception that we consider a magnetization three times as large here. Case C2 proposes a more severe configuration as initially presented by Balsara & Spicer [5] (see also [18]). Lastly, in case R1 (from [83]) we apply our 4th-order method to the solution of the relativistic MHD equations. Note that the energy correction step (see [5] and, e.g., see Sec. 6.4 of [18]) is applied to both cases C2 and R1, while case C1 has the largest magnetization our scheme can afford without it. This correction step is still performed at 2nd-order level and it was found to substantially improve the robustness of the code for larger magnetic field strengths. Computations are performed using 192 × 192 × 192 grid zones.

Results are shown in Fig. 11 (C1), 12 (C2) and 13 (R1) for the three selected cases. The explosion spawns an oval-shaped structure delimited by an outer fast forward shock and the presence of a magnetic field makes the propagation highly anisotropic by compressing the gas in the direction parallel to the field. Gas motion takes place along the field lines generating two oppositely moving plasma blobs. In the perpendicular direction, the outer fast shock becomes magnetically dominated with very weak compression. Discontinuous waves and smooth structures are well resolved by our scheme and the point-symmetry is preserved throughout the evolution. We point out that the fallback approach played an essential role for this test thus strengthening the benefits of a low/high order hybrid scheme while resolving complex flows featuring both smooth as well as discontinuous solutions.

## 5.8. 3D cloud-shock interaction

We propose a 3D version of the cloud-shock interaction test previously introduced by [72,84,85,29] in its MHD version, and by Mignone [69] for the RMHD case. For both configurations, the initial condition consists of a plane discontinuity propagating a fast shock and a rotational wave. Specifically, in the MHD case, the initial state $V = (ρ, v_x, v_y, v_z, B_x, B_y, B_z, p)$ is initialized as

$$ \begin{cases} V_L = (3.86859, 0, 0, 0, 0, 2.1826182, -2.1826182, 167.345) & \text{for } x < 0.6, \\ V_R = (1, -11.2536, 0, 0, 0, 0.56418958, 0.56418958, 1) & \text{for } x > 0.6, \end{cases} \qquad (56) $$
<figure><img src="image_12.png" /><figcaption>Fig. 10. Thermal pressure of the relativistic Orszag-Tang vortex computed at $t = 2\pi$ with different algorithms and resolutions. From left to right, results are obtained with, respectively, a 2nd-order LINEAR (left column), a 4th-order MP5 (middle column), and a 4th-order WENOZ scheme (bottom column). The grid resolution adopted is $N_x = N_y = 128$ (top row) and $N_x = N_y = 1024$ (bottom row). The presence or absence of O-points emerging from the central current sheets (only for the high resolution cases) is highlighted by a zoom on the central region.</figcaption></figure>

<table><caption>**Table 3**
Parameters used in the 3D blast wave problem.</caption><thead><tr><th>Case</th><th>ρ<sub>in</sub></th><th>ρ<sub>out</sub></th><th>p<sub>in</sub></th><th>p<sub>out</sub></th><th>B<sub>0</sub></th><th>Γ</th><th>r<sub>0</sub></th><th>L</th><th>t<sub>stop</sub></th><th>Solver</th><th>En. Corr.</th></tr></thead><tbody><tr><td>C1 (MHD)</td><td>1</td><td>1</td><td>10</td><td>0.1</td><td>3</td><td>5/3</td><td>0.1</td><td>1</td><td>0.2</td><td>HLL</td><td>NO</td></tr><tr><td>C2 (MHD)</td><td>1</td><td>1</td><td>10<sup>3</sup></td><td>0.1</td><td>100/√4π</td><td>1.4</td><td>0.1</td><td>1</td><td>0.01</td><td>HLL</td><td>YES</td></tr><tr><td>R1 (RMHD)</td><td>10<sup>-2</sup></td><td>10<sup>-4</sup></td><td>1</td><td>5 × 10<sup>-3</sup></td><td>1</td><td>4/3</td><td>0.08</td><td>6</td><td>4</td><td>GFORCE</td><td>YES</td></tr></tbody></table>

The initial condition consists of a uniform medium with density and pressure values set, respectively, to $ρ_{out}$ and $p_{out}$ and threaded by an oblique constant magnetic field,

$$ \mathbf{B} = B_0 \left( \sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta \right) \qquad (55) $$

where $θ = π/2$ and $φ = π/4$ are used in what follows. Note also that oblique configurations make the computation considerably more challenging than in the grid-aligned cases (see, for instance the discussion in §4.4 of Mattia & Mignone [83]). The computational domain is a cube of length $L$ centered at the origin with zero-gradient boundary conditions. An over-pressurized spherical region where pressure and density can have different values, namely $ρ_{in}$ and $p_{in}$, is initially set for $r < r_0$.

We examine three different cases with parameters listed in Table 3. Cases C1 is taken from Felker & Stone [36], with the exception that we consider a magnetization three times as large here. Case C2 proposes a more severe configuration as initially presented by Balsara & Spicer [5] (see also [18]). Lastly, in case R1 (from [83]) we apply our 4th-order method to the solution of the relativistic MHD equations. Note that the energy correction step (see [5] and, e.g., see Sec. 6.4 of [18]) is applied to both cases C2 and R1, while case C1 has the largest magnetization our scheme can afford without it. This correction step is still performed at 2nd-order level and it was found to substantially improve the robustness of the code for larger magnetic field strengths. Computations are performed using 192 × 192 × 192 grid zones.

Results are shown in Fig. 11 (C1), 12 (C2) and 13 (R1) for the three selected cases. The explosion spawns an oval-shaped structure delimited by an outer fast forward shock and the presence of a magnetic field makes the propagation highly anisotropic by compressing the gas in the direction parallel to the field. Gas motion takes place along the field lines generating two oppositely moving plasma blobs. In the perpendicular direction, the outer fast shock becomes magnetically dominated with very weak compression. Discontinuous waves and smooth structures are well resolved by our scheme and the point-symmetry is preserved throughout the evolution. We point out that the fallback approach played an essential role for this test thus strengthening the benefits of a low/high order hybrid scheme while resolving complex flows featuring both smooth as well as discontinuous solutions.

## 5.8. 3D cloud-shock interaction

We propose a 3D version of the cloud-shock interaction test previously introduced by [72,84,85,29] in its MHD version, and by Mignone [69] for the RMHD case. For both configurations, the initial condition consists of a plane discontinuity propagating a fast shock and a rotational wave. Specifically, in the MHD case, the initial state $V = (ρ, v_x, v_y, v_z, B_x, B_y, B_z, p)$ is initialized as

$$ \begin{cases} V_L = (3.86859, 0, 0, 0, 0, 2.1826182, -2.1826182, 167.345) & \text{for } x < 0.6, \\ V_R = (1, -11.2536, 0, 0, 0, 0.56418958, 0.56418958, 1) & \text{for } x > 0.6, \end{cases} \qquad (56) $$
<figure><img src="image_12.png" /><figcaption>Fig. 11. MHD 3D blast wave problem (case C1) at $t = 0.2$. The left and central panels show density (top left), pressure (top middle, in logarithmic scale), specific kinetic energy (bottom left) and magnetic energy (bottom right) at $z = 0$ for a simulation employing 192 × 192 × 192 grid points. In the rightmost panels we show the density (top) and magnetic energy density (bottom) on the major (red line) and minor (black line) diagonals.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Same as Fig. 11 but for the case C2 at $t = 0.01$.</figcaption></figure>

while in the RMHD test as

$$ \begin{cases} V_L = (39.5052, 0, 0, 0, 0, 1.9753, 129.72386) & \text{for } x < 0.6, \\ V_R = (1, -\sqrt{0.99}, 0, 0, 0, 0.5, 10^{-3}) & \text{for } x > 0.6. \end{cases} \tag{57} $$

At $x = (0.8, 0.5, 0.5)$ an overdense ($\rho = 10$) spherical cloud of radius 0.15 is initially in hydrostatic equilibrium with the ambient plasma under the adiabatic equation of state with $\Gamma = 5/3$ ($\Gamma = 4/3$ in the special relativistic case). For both setups the computational domain consists of the unit cube, with outflow (zero-gradient) boundary conditions employed everywhere, except at the lower $x$-boundary, where the initial condition is imposed.

During the evolution, the magnetic shock propagates to the right, and the high density cloud moves supersonically into the shock front. At the final evolution time ($t_f = 0.06$ MHD, $t_f = 0.6$ RMHD) the cloud has been completely shocked and it has developed its characteristic mushroom-like structure. Fig. 14 shows the initial and final density condition for a 256 × 256 × 256 grid points simulation. The 4th-order scheme properly recovers the expected structures. It is important to mention that, for this particular test problem, the employment of the order reduction procedure extended to the selective fluxes integration has been fundamental for the correct evolution of the simulations.
<figure><img src="image_13.png" /><figcaption>Fig. 11. MHD 3D blast wave problem (case C1) at $t = 0.2$. The left and central panels show density (top left), pressure (top middle, in logarithmic scale), specific kinetic energy (bottom left) and magnetic energy (bottom right) at $z = 0$ for a simulation employing 192 × 192 × 192 grid points. In the rightmost panels we show the density (top) and magnetic energy density (bottom) on the major (red line) and minor (black line) diagonals.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Same as Fig. 11 but for the case C2 at $t = 0.01$.</figcaption></figure>

while in the RMHD test as

$$ \begin{cases} V_L = (39.5052, 0, 0, 0, 0, 1.9753, 129.72386) & \text{for } x < 0.6, \\ V_R = (1, -\sqrt{0.99}, 0, 0, 0, 0.5, 10^{-3}) & \text{for } x > 0.6. \end{cases} \qquad (57) $$

At $x = (0.8, 0.5, 0.5)$ an overdense ($\rho = 10$) spherical cloud of radius 0.15 is initially in hydrostatic equilibrium with the ambient plasma under the adiabatic equation of state with $\Gamma = 5/3$ ($\Gamma = 4/3$ in the special relativistic case). For both setups the computational domain consists of the unit cube, with outflow (zero-gradient) boundary conditions employed everywhere, except at the lower $x$-boundary, where the initial condition is imposed.

During the evolution, the magnetic shock propagates to the right, and the high density cloud moves supersonically into the shock front. At the final evolution time ($t_f = 0.06$ MHD, $t_f = 0.6$ RMHD) the cloud has been completely shocked and it has developed its characteristic mushroom-like structure. Fig. 14 shows the initial and final density condition for a 256 × 256 × 256 grid points simulation. The 4th-order scheme properly recovers the expected structures. It is important to mention that, for this particular test problem, the employment of the order reduction procedure extended to the selective fluxes integration has been fundamental for the correct evolution of the simulations.
<figure><img src="image_14.png" /><figcaption>Fig. 13. Same as Fig. 11 but for the case R1 at $t = 4$. In the bottom left panel, the specific kinetic energy has been replaced by the relativistic counterpart.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Non-relativistic (top panels) and relativistic (bottom panels) cloud-shock interaction at initial ($t = 0$, left panels) and final time ($t = 0.06$ for the non-relativistic case, and $t = 0.6$ for the relativistic case, right panels). Shown is the density at $z = 0.5$ in logarithmic scale. Both simulations ran with 256 × 256 × 256 grid points. The colorscales span over a range of [1, 24] for the non-relativistic case and [1, 110] for the relativistic case.</figcaption></figure>

## 6. Conclusions

In this paper we have presented a genuinely 4^th^-order accurate finite volume method for the solution of the classical and special relativistic MHD equations employing the constrained transport algorithm. Our scheme has been implemented in the PLUTO code for astrophysical plasma dynamics [68] and is rooted over the method originally proposed by McCorquodale & Colella [28], adding a number of innovative aspects that yield an accurate, robust and efficient computational tool. In fact, our scheme maximizes the benefits of the FV formalism (e.g., robustness when dealing with shocks, discontinuities, and steep gradients) along with techniques inspired by the FD formalism, such as the introduction of pointwise to pointwise reconstruction operations that ease up the structure of the scheme. To this aim, we revisited standard spatial reconstruction schemes (such as WENOZ and MP5), traditionally based
<figure><img src="image_15.png" /><figcaption>Fig. 13. Same as Fig. 11 but for the case R1 at $t = 4$. In the bottom left panel, the specific kinetic energy has been replaced by the relativistic counterpart.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Non-relativistic (top panels) and relativistic (bottom panels) cloud-shock interaction at initial ($t = 0$, left panels) and final time ($t = 0.06$ for the non-relativistic case, and $t = 0.6$ for the relativistic case, right panels). Shown is the density at $z = 0.5$ in logarithmic scale. Both simulations ran with 256 × 256 × 256 grid points. The colorscales span over a range of [1, 24] for the non-relativistic case and [1, 110] for the relativistic case.</figcaption></figure>

## 6. Conclusions

In this paper we have presented a genuinely 4th-order accurate finite volume method for the solution of the classical and special relativistic MHD equations employing the constrained transport algorithm. Our scheme has been implemented in the PLUTO code for astrophysical plasma dynamics [68] and is rooted over the method originally proposed by McCorquodale & Colella [28], adding a number of innovative aspects that yield an accurate, robust and efficient computational tool. In fact, our scheme maximizes the benefits of the FV formalism (e.g., robustness when dealing with shocks, discontinuities, and steep gradients) along with techniques inspired by the FD formalism, such as the introduction of pointwise to pointwise reconstruction operations that ease up the structure of the scheme. To this aim, we revisited standard spatial reconstruction schemes (such as WENOZ and MP5), traditionally based
on the knowledge of 1D volume averages, so as to provide 4^th^-order accurate left and right states at zone interfaces by directly interpolating point values rather than volume averages. Pointwise reconstructions can be accomplished in primitive or characteristic variables and they lead to a compact stencil and a more cost-effective scheme. Moreover, our implementation is endowed with the generalized upwind constrained transport (UCT) of Mignone & Del Zanna [18], plus a new generalization of the UCT-GFORCE averaging scheme for special relativistic MHD that has been widely validated in the presented benchmarks.

The performance of the method has been assessed against selected 2D and 3D numerical test problems, demonstrating: i) genuine 4^th^-order accuracy in smooth problems, ii) reduced numerical dissipation iii) the code's ability to correctly capture propagating discontinuities and to resolve small-scale structures and turbulence. The robustness of the method is reinforced by a fallback approach that reverts the spatial order of the scheme to 2 in proximity of steep gradients. To identify such critical regions, we have proposed a new discontinuity detector - based on a derivatives ratio - which employs up to the fourth derivative to pinpoint high-frequency oscillations typically arising in proximity of a discontinuity. The order reduction procedure turned out to be, in our experience, a crucial element when undertaking severe configurations, allowing our high-order scheme to successfully carry out computations otherwise feasible only with 2^nd^-order schemes (see, e.g., configuration C2 §5.7).

We point out that, despite the augmented complexity of high-order algorithms, genuine 4^th^-order accuracy in smooth problems leads to a net saving of computational time. In §5.3 we demonstrated that, for a given accuracy, the high-order scheme will require, approximately $\sqrt{N_2}$ grid cells when compared to a 2^nd^-order scheme with a corresponding reduction of CPU time that scales as $\sim N_2^{-(d+1)/2}$, where d is the number of spatial dimensions. When applied to simple 3D computations, for instance, we found that the accuracy reached by a 2^nd^-order scheme with 512^3 grid points is equally obtained when employing our 4^th^-order scheme with only $\sim 34^3$ points, with a corresponding reduction of the total computational time by a few orders of magnitude.

Better accuracy leads also to reduced numerical dissipation and thus an enhanced potential at describing spatial scales that would otherwise be challenging to probe within a traditional 2^nd^-order framework. Our results indicate that the decay rate of numerical diffusion is $\sim N^{-5}$ for our novel scheme, and only $\sim N^{-3}$ for the traditional scheme. For all these reasons our 4^th^-order scheme proved to be a computational tool suited for the next generation of numerical simulations in computational fluid dynamics. In a companion paper we will extend this formalism to non-Cartesian geometries.

Compared to other high-order competitive methods, such as discontinuous Galerkin (DG), our FV scheme requires fewer equations to be evolved in time (e.g. 8 instead of 32 for a 4^th^-order method), does not reduce stability with increasing order and it is better suited for flows containing smooth and discontinuous features. Its parallel implementation requires, however, several more inter-processor communications (one per Runge-Kutta stage). A more thorough comparison, though, has to be considered in a separate paper.

### CRediT authorship contribution statement

**V. Berta:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing – original draft, Writing – review & editing. **A. Mignone:** Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Supervision, Validation, Writing – original draft, Writing – review & editing. **M. Bugli:** Data curation, Resources, Validation, Visualization, Writing – review & editing. **G. Mattia:** Visualization, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgements

We thank L. Del Zanna for useful discussions and comments. The PLUTO code is publicly available and the simulation data will be shared on reasonable request to the corresponding author.

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 Performance Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No 101093441 (SPACE).

This paper is 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).

We also acknowledge the OCCAM supercomputing facility available at the Competence Centre for Scientific Computing at the University of Torino.
on the knowledge of 1D volume averages, so as to provide 4^th^-order accurate left and right states at zone interfaces by directly interpolating point values rather than volume averages. Pointwise reconstructions can be accomplished in primitive or characteristic variables and they lead to a compact stencil and a more cost-effective scheme. Moreover, our implementation is endowed with the generalized upwind constrained transport (UCT) of Mignone & Del Zanna [18], plus a new generalization of the UCT-GFORCE averaging scheme for special relativistic MHD that has been widely validated in the presented benchmarks.

The performance of the method has been assessed against selected 2D and 3D numerical test problems, demonstrating: i) genuine 4^th^-order accuracy in smooth problems, ii) reduced numerical dissipation iii) the code's ability to correctly capture propagating discontinuities and to resolve small-scale structures and turbulence. The robustness of the method is reinforced by a fallback approach that reverts the spatial order of the scheme to 2 in proximity of steep gradients. To identify such critical regions, we have proposed a new discontinuity detector - based on a derivatives ratio - which employs up to the fourth derivative to pinpoint high-frequency oscillations typically arising in proximity of a discontinuity. The order reduction procedure turned out to be, in our experience, a crucial element when undertaking severe configurations, allowing our high-order scheme to successfully carry out computations otherwise feasible only with 2^nd^-order schemes (see, e.g., configuration C2 §5.7).

We point out that, despite the augmented complexity of high-order algorithms, genuine 4^th^-order accuracy in smooth problems leads to a net saving of computational time. In §5.3 we demonstrated that, for a given accuracy, the high-order scheme will require, approximately $\sqrt{N_2}$ grid cells when compared to a 2^nd^-order scheme with a corresponding reduction of CPU time that scales as $\sim N_2^{-(d+1)/2}$, where d is the number of spatial dimensions. When applied to simple 3D computations, for instance, we found that the accuracy reached by a 2^nd^-order scheme with 512^3 grid points is equally obtained when employing our 4^th^-order scheme with only $\sim 34^3$ points, with a corresponding reduction of the total computational time by a few orders of magnitude.

Better accuracy leads also to reduced numerical dissipation and thus an enhanced potential at describing spatial scales that would otherwise be challenging to probe within a traditional 2^nd^-order framework. Our results indicate that the decay rate of numerical diffusion is $\sim N^{-5}$ for our novel scheme, and only $\sim N^{-3}$ for the traditional scheme. For all these reasons our 4^th^-order scheme proved to be a computational tool suited for the next generation of numerical simulations in computational fluid dynamics. In a companion paper we will extend this formalism to non-Cartesian geometries.

Compared to other high-order competitive methods, such as discontinuous Galerkin (DG), our FV scheme requires fewer equations to be evolved in time (e.g. 8 instead of 32 for a 4^th^-order method), does not reduce stability with increasing order and it is better suited for flows containing smooth and discontinuous features. Its parallel implementation requires, however, several more inter-processor communications (one per Runge-Kutta stage). A more thorough comparison, though, has to be considered in a separate paper.

### CRediT authorship contribution statement

**V. Berta:** Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing – original draft, Writing – review & editing. **A. Mignone:** Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Supervision, Validation, Writing – original draft, Writing – review & editing. **M. Bugli:** Data curation, Resources, Validation, Visualization, Writing – review & editing. **G. Mattia:** Visualization, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgements

We thank L. Del Zanna for useful discussions and comments. The PLUTO code is publicly available and the simulation data will be shared on reasonable request to the corresponding author.

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 Performance Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece, Italy, Norway, and Spain under grant agreement No 101093441 (SPACE).

This paper is 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).

We also acknowledge the OCCAM supercomputing facility available at the Competence Centre for Scientific Computing at the University of Torino.
## Appendix A. Laplacian operators in Cartesian coordinates

We start by demonstrating the validity of Eq. (16) and deriving the expression of the 2nd-order accurate Laplacian operator for a given function $U(x)$ in a Cartesian mesh in one dimension. The possibility to straightforwardly generalize the treatment to multiple dimensions is guaranteed when considering Cartesian meshes due to the orthogonality of the standard basis. The (1D) cell-average is given by

$$ \langle U \rangle_c = \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} U(x) dx, \qquad (A.1) $$

where $h$ is the cell width. By expanding in Taylor series $U(x)$ near the cell-center $x_c$ up to 4th-order we obtain

$$ \langle U \rangle_c = \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx, \qquad (A.2) $$

and, since odd integrands vanish when integrated over symmetrical domains, from straightforward calculation we get

$$ \langle U \rangle_c = U(x_c) + \frac{h^2}{24} U''(x_c) + O(h^4). \qquad (A.3) $$

We now obtain a 2nd-order accurate expression of the 2nd-order derivative $U''$. By expanding in Taylor series up to 4th-order the nearby average values: $\langle U \rangle_{c-\hat{e}_x}$, $\langle U \rangle_{c+\hat{e}_x}$ we get

$$ \langle U \rangle_{c-\hat{e}_x} = \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx, \qquad (A.4) \\ \langle U \rangle_{c+\hat{e}_x} = \frac{1}{h} \int_{\frac{h}{2}}^{\frac{3h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx. $$

Since the interval is no longer symmetric, the odd powers no longer integrate to zero. The final result is

$$ \langle U \rangle_{c-\hat{e}_x} = U(x_c) - U'(x_c)h + \frac{13}{24} U''(x_c)h^2 - \frac{5}{24} U'''(x_c)h^3 + O(h^4), \qquad (A.5) $$

$$ \langle U \rangle_{c+\hat{e}_x} = U(x_c) + U'(x_c)h + \frac{13}{24} U''(x_c)h^2 + \frac{5}{24} U'''(x_c)h^3 + O(h^4). \qquad (A.6) $$

Adding Eq. (A.5) to Eq. (A.6) and subtracting twice Eq. (A.3)

$$ \langle U \rangle_{c+\hat{e}_x} + \langle U \rangle_{c-\hat{e}_x} - 2 \langle U \rangle_c = U''h^2 + O(h^4), \qquad (A.7) $$

we can express the 2nd-order accurate second derivative in terms of volume-averaged quantities as

$$ U'' = \frac{\langle U \rangle_{c+\hat{e}_x} - 2 \langle U \rangle_c + \langle U \rangle_{c-\hat{e}_x}}{h^2} + O(h^2), \qquad (A.8) $$

which defines the Laplacian operator in Eq. (16).

Conversely, Eq. (17) can be formally obtained applying a Simpson's quadrature rule by integrating over the cell width $h$

$$ \int_{-\frac{h}{2}}^{\frac{h}{2}} U(x) dx \approx \int_{-\frac{h}{2}}^{\frac{h}{2}} (\alpha x^2 + \beta x + \gamma) dx. \qquad (A.9) $$

The function $U(x)$ is approximated by a 2nd-order polynomial in the coefficients $\alpha, \beta, \gamma$. The polynomial in the right hand side of Eq. (A.9) is the best approximation of $U(x)$ when the coefficients $\alpha, \beta, \gamma$ are set by evaluating the function in the 3 nearby cell centers:

$$ \alpha = \frac{U_{c+\hat{e}_x} - 2U_c + U_{c-\hat{e}_x}}{2}, \quad \beta = \frac{U_{c+\hat{e}_x} - U_{c-\hat{e}_x}}{2}, \quad \gamma = U_c. \qquad (A.10) $$

By replacing these expressions in Eq. (A.9) and by direct calculation of the integral we obtain Eq. (17).
## Appendix A. Laplacian operators in Cartesian coordinates

We start by demonstrating the validity of Eq. (16) and deriving the expression of the 2nd-order accurate Laplacian operator for a given function $U(x)$ in a Cartesian mesh in one dimension. The possibility to straightforwardly generalize the treatment to multiple dimensions is guaranteed when considering Cartesian meshes due to the orthogonality of the standard basis. The (1D) cell-average is given by

$$ \langle U \rangle_c = \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} U(x) dx, \qquad (A.1) $$

where $h$ is the cell width. By expanding in Taylor series $U(x)$ near the cell-center $x_c$ up to 4th-order we obtain

$$ \langle U \rangle_c = \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx, \qquad (A.2) $$

and, since odd integrands vanish when integrated over symmetrical domains, from straightforward calculation we get

$$ \langle U \rangle_c = U(x_c) + \frac{h^2}{24} U''(x_c) + O(h^4). \qquad (A.3) $$

We now obtain a 2nd-order accurate expression of the 2nd-order derivative $U''$. By expanding in Taylor series up to 4th-order the nearby average values: $\langle U \rangle_{c-\hat{e}_x}$, $\langle U \rangle_{c+\hat{e}_x}$ we get

$$ \begin{aligned} \langle U \rangle_{c-\hat{e}_x} &= \frac{1}{h} \int_{-\frac{h}{2}}^{\frac{h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx, \\ \langle U \rangle_{c+\hat{e}_x} &= \frac{1}{h} \int_{\frac{h}{2}}^{-\frac{h}{2}} \left[ U(x_c) + U'(x_c)(x-x_c) + \frac{U''(x_c)}{2}(x-x_c)^2 + \frac{U'''(x_c)}{6}(x-x_c)^3 + O(x^4) \right] dx. \end{aligned} \qquad (A.4) $$

Since the interval is no longer symmetric, the odd powers no longer integrate to zero. The final result is

$$ \langle U \rangle_{c-\hat{e}_x} = U(x_c) - U'(x_c)h + \frac{13}{24} U''(x_c)h^2 - \frac{5}{24} U'''(x_c)h^3 + O(h^4), \qquad (A.5) $$

$$ \langle U \rangle_{c+\hat{e}_x} = U(x_c) + U'(x_c)h + \frac{13}{24} U''(x_c)h^2 + \frac{5}{24} U'''(x_c)h^3 + O(h^4). \qquad (A.6) $$

Adding Eq. (A.5) to Eq. (A.6) and subtracting twice Eq. (A.3)

$$ \langle U \rangle_{c+\hat{e}_x} + \langle U \rangle_{c-\hat{e}_x} - 2 \langle U \rangle_c = U''h^2 + O(h^4), \qquad (A.7) $$

we can express the 2nd-order accurate second derivative in terms of volume-averaged quantities as

$$ U'' = \frac{\langle U \rangle_{c+\hat{e}_x} - 2 \langle U \rangle_c + \langle U \rangle_{c-\hat{e}_x}}{h^2} + O(h^2), \qquad (A.8) $$

which defines the Laplacian operator in Eq. (16).

Conversely, Eq. (17) can be formally obtained applying a Simpson's quadrature rule by integrating over the cell width $h$

$$ \int_{-\frac{h}{2}}^{\frac{h}{2}} U(x) dx \approx \int_{-\frac{h}{2}}^{\frac{h}{2}} (\alpha x^2 + \beta x + \gamma) dx. \qquad (A.9) $$

The function $U(x)$ is approximated by a 2nd-order polynomial in the coefficients $\alpha, \beta, \gamma$. The polynomial in the right hand side of Eq. (A.9) is the best approximation of $U(x)$ when the coefficients $\alpha, \beta, \gamma$ are set by evaluating the function in the 3 nearby cell centers:

$$ \alpha = \frac{U_{c+\hat{e}_x} - 2U_c + U_{c-\hat{e}_x}}{2}, \quad \beta = \frac{U_{c+\hat{e}_x} - U_{c-\hat{e}_x}}{2}, \quad \gamma = U_c. \qquad (A.10) $$

By replacing these expressions in Eq. (A.9) and by direct calculation of the integral we obtain Eq. (17).
## Appendix B. Explicit expressions for the UCT scheme

Here, for the sake of completeness and to also amend some equations containing typos, we briefly recap the most relevant expressions needed in our UCT scheme, as originally presented by Mignone & Del Zanna [18]. The edge-centered EMF (Eq. (39)) requires the reconstructed values of the transverse velocities as well as the linear combination coefficients a and d (for the centered and diffusive flux contribution, respectively).

* **Transverse Velocities:** when using the HLL, HLLC or HLLD classical MHD Riemann solvers at a zone interface, the transverse velocities are obtained, e.g., at an x-face, as

$$ \bar{v}_{t,x_f} = \frac{\alpha_x^R v_t^L_{x_f} + \alpha_x^L v_t^R_{x_f}}{\alpha_x^R + \alpha_x^L}, \quad (t = y, z), \qquad (B.1) $$

where $v_{t,x_f}^{L/R}$ are the left / right transverse velocities reconstructed at an x-face,

$$ \alpha_x^R = \max(0, \lambda_{x_f}^R), \quad \alpha_x^L = -\min(0, \lambda_{x_f}^L) \qquad (B.2) $$

with $\lambda_{x_f}^{L/R}$ being the smallest and largest characteristic speeds (see Eq. 14 of [18]). In the case of the GFORCE Riemann solver instead, we use

$$ \bar{v}_t = \omega \left[ v_t^{\text{LW}} - \frac{\tau_x}{2} (v_t^R - v_t^L) v_x^{\text{LW}} \right] + (1 - \omega) \left( \frac{v_t^L + v_t^R}{2} \right) \qquad (B.3) $$

(note that Eq. 51 of [18] contains “y” instead of “t”), with $\omega = 1/(1+c)$ (c is the Courant factor, $c = 1$ yields the FORCE scheme), $v^{\text{LW}}$ is the velocity in the Lax-Wendroff state (see. Eq. 48 of [18]), and $\tau_x = 1/\max(|\lambda_{x_f}^L|, |\lambda_{x_f}^R|)$.

* **Combination Coefficients a and d.** For the HLL or HLLC Riemann solvers (away from degeneracies), the combination coefficients are evaluated as arithmetic averages across two adjacent faces (Eq. 34 and 35 of [18]) using the values available with the 1D Riemann solver:

$$ a_x^L = \frac{\alpha_x^R}{\alpha_x^R + \alpha_x^L}, \quad a_x^R = \frac{\alpha_x^L}{\alpha_x^R + \alpha_x^L}, \quad d_x^L = d_x^R = \frac{\alpha_x^R \alpha_x^L}{\alpha_x^R + \alpha_x^L}, \qquad (B.4) $$

where $a_x^{L/R}$ have been given in Eq. (B.2). In the limit $B_x \to 0$, the corresponding expressions for the HLLC solver are given by Eq. 38 of [18]. For the HLLD solver, instead, we employ

$$ a^L = \frac{1 + \nu^*}{2}, \quad a^R = \frac{1 - \nu^*}{2}, \quad d^s = \frac{1}{2}(\nu^s - \nu^*)\tilde{\chi}^s + \frac{1}{2}(|\lambda^{*s}| - \nu^*\lambda^{*s}), \quad (s = L, R) \qquad (B.5) $$

where $\tilde{\chi}^s = (\lambda^{*s} - \lambda^s)\chi^s$, while

$$ \nu^s = \frac{\lambda^{*s} + \lambda^s}{|\lambda^{*s}| + |\lambda^s|}, \quad \nu^* = \frac{\lambda^{*R} + \lambda^{*L}}{|\lambda^{*R}| + |\lambda^{*L}|}, \qquad (B.6) $$

and $\lambda^{*L/R}$ are the Alfvén velocities of the HLLD fan (see Eq. 40 of [18]).

## References

[1] S.A.E.G. Falle, S.S. Komissarov, P. Joarder, A multidimensional upwind scheme for magnetohydrodynamics, Mon. Not. R. Astron. Soc. 297 (1) (1998) 265-277, https://doi.org/10.1046/j.1365-8711.1998.01506.x.

[2] D.S. Balsara, Total variation diminishing scheme for adiabatic and isothermal magnetohydrodynamics, Astrophys. J. Suppl. Ser. 116 (1) (1998) 133-153, https://doi.org/10.1086/313093.

[3] W. Dai, P.R. Woodward, On the divergence-free condition and conservation laws in numerical simulations for supersonic magnetohydrodynamical flows, Astrophys. J. 494 (1) (1998) 317-335, https://doi.org/10.1086/305176.

[4] K.G. Powell, P.L. Roe, T.J. Linde, T.I. Gombosi, D.L. De Zeeuw, A solution-adaptive upwind scheme for ideal magnetohydrodynamics, J. Comput. Phys. 154 (2) (1999) 284-309, https://doi.org/10.1006/jcp.1999.6299.

[5] D.S. Balsara, D.S. Spicer, A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations, J. Comput. Phys. 149 (1999) 270-292, https://doi.org/10.1006/jcp.1998.6153.

[6] G. Tóth, The V·B=0 constraint in shock-capturing magnetohydrodynamics codes, J. Comput. Phys. 161 (2) (2000) 605-652, https://doi.org/10.1006/jcp.2000.6519.

[7] D.S. Balsara, Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction, Astrophys. J. Suppl. Ser. 151 (1) (2004) 149-184, https://doi.org/10.1086/381377, arXiv:astro-ph/0308249.

[8] M. Torrilhon, Locally divergence-preserving upwind finite volume schemes for magnetohydrodynamic equations, SIAM J. Sci. Comput. 26 (2005) 1166-1191, https://doi.org/10.1137/S1064827503426401.

[9] T.A. Gardiner, J.M. Stone, An unsplit Godunov method for ideal MHD via constrained transport, J. Comput. Phys. 205 (2005) 509-539, https://doi.org/10.1016/j.jcp.2004.11.016, arXiv:astro-ph/0501557.

[10] S. Fromang, P. Hennebelle, R. Teyssier, A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics, Astron. Astrophys. 457 (2) (2006) 371-384, https://doi.org/10.1051/0004-6361:20065371, arXiv:astro-ph/0607230.

[11] A. Mignone, G. Bodo, An HLLC Riemann solver for relativistic flows - II. Magnetohydrodynamics, Mon. Not. R. Astron. Soc. 368 (35) (2006) 1040-1054, https://doi.org/10.1111/j.1365-2966.2006.10162.x, arXiv:astro-ph/0601640.
## Appendix B. Explicit expressions for the UCT scheme

Here, for the sake of completeness and to also amend some equations containing typos, we briefly recap the most relevant expressions needed in our UCT scheme, as originally presented by Mignone & Del Zanna [18]. The edge-centered EMF (Eq. (39)) requires the reconstructed values of the transverse velocities as well as the linear combination coefficients a and d (for the centered and diffusive flux contribution, respectively).

* **Transverse Velocities:** when using the HLL, HLLC or HLLD classical MHD Riemann solvers at a zone interface, the transverse velocities are obtained, e.g., at an x-face, as

$$ \bar{v}_{t,x_f} = \frac{\alpha_x^R v_t^L_{x_f} + \alpha_x^L v_t^R_{x_f}}{\alpha_x^R + \alpha_x^L}, \quad (t = y, z), \qquad (B.1) $$

where $v_{t,x_f}^{L/R}$ are the left / right transverse velocities reconstructed at an x-face,

$$ \alpha_x^R = \max(0, \lambda_{x_f}^R), \quad \alpha_x^L = -\min(0, \lambda_{x_f}^L) \qquad (B.2) $$

with $\lambda_{x_f}^{L/R}$ being the smallest and largest characteristic speeds (see Eq. 14 of [18]). In the case of the GFORCE Riemann solver instead, we use

$$ \bar{v}_t = \omega \left[ v_t^{\text{LW}} - \frac{\tau_x}{2} (v_t^R - v_t^L) v_x^{\text{LW}} \right] + (1 - \omega) \left( \frac{v_t^L + v_t^R}{2} \right) \qquad (B.3) $$

(note that Eq. 51 of [18] contains “y” instead of “t”), with $\omega = 1/(1+c)$ (c is the Courant factor, $c = 1$ yields the FORCE scheme), $v^{\text{LW}}$ is the velocity in the Lax-Wendroff state (see. Eq. 48 of [18]), and $\tau_x = 1/\max(|\lambda_{x_f}^L|, |\lambda_{x_f}^R|)$.

* **Combination Coefficients a and d.** For the HLL or HLLC Riemann solvers (away from degeneracies), the combination coefficients are evaluated as arithmetic averages across two adjacent faces (Eq. 34 and 35 of [18]) using the values available with the 1D Riemann solver:

$$ a_x^L = \frac{\alpha_x^R}{\alpha_x^R + \alpha_x^L}, \quad a_x^R = \frac{\alpha_x^L}{\alpha_x^R + \alpha_x^L}, \quad d_x^L = d_x^R = \frac{\alpha_x^R \alpha_x^L}{\alpha_x^R + \alpha_x^L}, \qquad (B.4) $$

where $a_x^{L/R}$ have been given in Eq. (B.2). In the limit $B_x \to 0$, the corresponding expressions for the HLLC solver are given by Eq. 38 of [18]. For the HLLD solver, instead, we employ

$$ a^L = \frac{1 + \nu^*}{2}, \quad a^R = \frac{1 - \nu^*}{2}, \quad d^s = \frac{1}{2}(\nu^s - \nu^*)\tilde{\chi}^s + \frac{1}{2}(|\lambda^{*s}| - \nu^*\lambda^{*s}), \quad (s = L, R) \qquad (B.5) $$

where $\tilde{\chi}^s = (\lambda^{*s} - \lambda^s)\chi^s$, while

$$ \nu^s = \frac{\lambda^{*s} + \lambda^s}{|\lambda^{*s}| + |\lambda^s|}, \quad \nu^* = \frac{\lambda^{*R} + \lambda^{*L}}{|\lambda^{*R}| + |\lambda^{*L}|}, \qquad (B.6) $$

and $\lambda^{*L/R}$ are the Alfvén velocities of the HLLD fan (see Eq. 40 of [18]).

## References

[1] S.A.E.G. Falle, S.S. Komissarov, P. Joarder, A multidimensional upwind scheme for magnetohydrodynamics, Mon. Not. R. Astron. Soc. 297 (1) (1998) 265-277, https://doi.org/10.1046/j.1365-8711.1998.01506.x.

[2] D.S. Balsara, Total variation diminishing scheme for adiabatic and isothermal magnetohydrodynamics, Astrophys. J. Suppl. Ser. 116 (1) (1998) 133-153, https://doi.org/10.1086/313093.

[3] W. Dai, P.R. Woodward, On the divergence-free condition and conservation laws in numerical simulations for supersonic magnetohydrodynamical flows, Astrophys. J. 494 (1) (1998) 317-335, https://doi.org/10.1086/305176.

[4] K.G. Powell, P.L. Roe, T.J. Linde, T.I. Gombosi, D.L. De Zeeuw, A solution-adaptive upwind scheme for ideal magnetohydrodynamics, J. Comput. Phys. 154 (2) (1999) 284-309, https://doi.org/10.1006/jcp.1999.6299.

[5] D.S. Balsara, D.S. Spicer, A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations, J. Comput. Phys. 149 (1999) 270-292, https://doi.org/10.1006/jcp.1998.6153.

[6] G. Tóth, The V·B=0 constraint in shock-capturing magnetohydrodynamics codes, J. Comput. Phys. 161 (2) (2000) 605-652, https://doi.org/10.1006/jcp.2000.6519.

[7] D.S. Balsara, Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction, Astrophys. J. Suppl. Ser. 151 (1) (2004) 149-184, https://doi.org/10.1086/381377, arXiv:astro-ph/0308249.

[8] M. Torrilhon, Locally divergence-preserving upwind finite volume schemes for magnetohydrodynamic equations, SIAM J. Sci. Comput. 26 (2005) 1166-1191, https://doi.org/10.1137/S1064827503426401.

[9] T.A. Gardiner, J.M. Stone, An unsplit Godunov method for ideal MHD via constrained transport, J. Comput. Phys. 205 (2005) 509-539, https://doi.org/10.1016/j.jcp.2004.11.016, arXiv:astro-ph/0501557.

[10] S. Fromang, P. Hennebelle, R. Teyssier, A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics, Astron. Astrophys. 457 (2) (2006) 371-384, https://doi.org/10.1051/0004-6361:20065371, arXiv:astro-ph/0607230.

[11] A. Mignone, G. Bodo, An HLLC Riemann solver for relativistic flows - II. Magnetohydrodynamics, Mon. Not. R. Astron. Soc. 368 (35) (2006) 1040-1054, https://doi.org/10.1111/j.1365-2966.2006.10162.x, arXiv:astro-ph/0601640.