# A multi-dimensional, robust, and cell-centered finite-volume scheme for the ideal MHD equations

Pascal Tremblin, Rémi Bourgeois, Solène Bulteau, Samuel Kokh, Thomas Padioleau, Maxime Delorme, Antoine Strugarek, Matthias González, Allan Sacha Brun

► To cite this version:

Pascal Tremblin, Rémi Bourgeois, Solène Bulteau, Samuel Kokh, Thomas Padioleau, et al. A multi-dimensional, robust, and cell-centered finite-volume scheme for the ideal MHD equations. Journal of Computational Physics, 2024, 519, 10.1016/j.jcp.2024.113455. insu-04843503

HAL Id: insu-04843503

https://insu.hal.science/insu-04843503v1

Submitted on 17 Dec 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_2.png" /></figure>

# A multi-dimensional, robust, and cell-centered finite-volume scheme for the ideal MHD equations

Pascal Tremblin<sup>a,*</sup>, Rémi Bourgeois<sup>a</sup>, Solène Bulteau<sup>a</sup>, Samuel Kokh<sup>b</sup>, Thomas Padioleau<sup>a</sup>, Maxime Delorme<sup>c</sup>, Antoine Strugarek<sup>c</sup>, Matthias González<sup>d</sup>, Allan Sacha Brun<sup>c</sup>

## ARTICLE INFO

**Keywords:**

Finite volume
Magnetohydrodynamics
Relaxation
Splitting
Entropy satisfying

## ABSTRACT

We present a new multi-dimensional, robust, and cell-centered finite-volume scheme for the ideal MHD equations. This scheme relies on relaxation and splitting techniques and can be easily used at high order. A fully conservative version is not entropy satisfying but is observed experimentally to be more robust than standard constrained transport schemes at low plasma beta. At very low plasma beta and high Alfvén number, we have designed an entropy-satisfying version that is not conservative for the magnetic field but preserves admissible states and we switch locally a-priori between the two versions depending on the regime of plasma beta and Alfvén number. This strategy is robust in a wide range of standard MHD test cases, all performed at second order with a classic MUSCL-Hancock scheme.

## 1. Introduction

Developing a robust multi-dimensional numerical scheme for the ideal MHD equations remains a challenge that is of great importance for astrophysics and plasma physics applications. A MHD flow is characterized by an exact zero-divergence magnetic field, and by using terms that are proportional to the divergence of the magnetic field, the MHD equations can be put in a fully conservative form, with density, momentum, energy and magnetic field conservation. However this form introduces a source term proportional to the divergence of the magnetic field on the entropy evolution equation, leading to an unstable scheme for multi-dimensional test cases, because of discretization errors on this source term.

A solution to this problem is to remove the divergence errors so that the source term in the entropy evolution equation is as small as possible. Such a solution encompass the divergence-cleaning method (see [1–4]) and the constrained transport method (see [5–8]). These methods greatly improve the stability of MHD numerical schemes and have been used in numerous applications in astrophysics and plasma physics. However, they are not entropy satisfying and may fail with negative energies especially in the low plasma beta regime. This problem is mitigated by using a threshold value for the internal energy, effectively breaking the energy conservation of the numerical scheme (ref divB cleaning with threshold). Another solution is to design an entropy satisfying numerical

<footnote><sup>a</sup> Université Paris-Saclay, UVSQ, CNRS, CEA, Maison de la Simulation, 91191, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Université Paris-Saclay, CEA, Service de Génie Logiciel pour la Simulation, 91191, Gif-sur-Yvette, France</footnote>

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

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

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: pascal.tremblin@cea.fr (P. Tremblin).</footnote>
scheme for any value of the divergence of the magnetic field. This solution has been explored using relaxation methods in [9–11]. Originally, it has been shown that a multi-dimensional solver with the introduction of non-conservative Powell source terms in the momentum, induction and energy evolution equations allows to obtain a symmetric form of the MHD equations [12,13], but [11] has demonstrated that it is also possible to obtain a symmetric form with a source term only on the induction equation, therefore preserving energy and momentum conservation with an entropy satisfying numerical scheme.

In recent years, significant advancements have been made in splitting strategies for designing numerical solvers for Euler equations. In the works by [14–16], the approximation algorithm is divided into two steps: an acoustic step and a transport step. For one-dimensional cases, these methods resemble the explicit Lagrange-Projection approach [17–19]. However, this new splitting technique avoids the use of a moving Lagrangian mesh and is readily adaptable to multi-dimensional problems. On the other hand, significant progresses have also been made on approximate Riemann solvers based on relaxation strategies [20–25].

In this paper, we build on the proposition of a relaxation approximation for the MHD system [10,11] by taking advantage of splitting techniques introduced in [14] to design a fully-conservative multi-dimensional MHD solver in regions of high plasma beta / low Alfvén number, and an entropy satisfying version with an entropy correction in regions of low plasma beta / high Alfvén number. The resulting solver, therefore, allies a robust entropy-satisfying and a fully-conservative scheme depending on the regime of the flow. In Sect. 2, we introduce the different systems of equations for MHD (conservative and non-conservative) and the corresponding entropy evolution equation. In Sect. 3, we present the splitting of the equations in a magneto-acoustic and transport subsystems. Sect. 4 and 5 introduce the numerical methods used to solve the evolution of these subsystems and Sect. 6 provides the global fully-conservative scheme for the MHD system of equations. Sect. 7 is devoted to the entropy analysis of the numerical method, showing that the fully-conservative solver is not entropy-satisfying and introduces an entropy correction on the induction equation in order to obtain an entropy-satisfying method at the price of loosing the magnetic field conservation. In Sect. 8, we provide numerical tests in 1D and 2D at second order by leveraging the advantages of the fully-conservative and entropy-satisfying solvers depending on the regime of the flow. We provide our conclusions and a discussion in Sect. 9.

## 2. MHD equations

The ideal MHD equations are given by the evolution equations of the fluid density $ρ$, momentum $ρu$, energy $ρ(e + u^2/2)$, and the Faraday's law of induction describing the evolution of the magnetic field $B$

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + \mathbf{j} \times \mathbf{B}, $$

$$ \partial_t (\rho(e + \mathbf{u}^2/2)) + \nabla \cdot (\rho(e + \mathbf{u}^2/2)\mathbf{u}) = -\nabla \cdot (\rho\mathbf{u}) + (\mathbf{j} \times \mathbf{B}) \cdot \mathbf{u}, $$

$$ \partial_t \mathbf{B} + \nabla \times \mathbf{E} = 0. \qquad (1) $$

The term $\mathbf{j} \times \mathbf{B}$ is the Lorentz force. This system of equations is closed with the ideal Ohm's law $\mathbf{E} = -\mathbf{u} \times \mathbf{B}$, the low frequency Maxwell equation $\mathbf{j} = \nabla \times \mathbf{B}$ assuming a system of units in which the vacuum permeability is one, and an equation of state connecting the pressure $p$ to the density $ρ$ (or specific volume $\tau = 1/ρ$) and internal energy $e$. The equation of state also defines the specific physical entropy $s(\tau, e)$ assuming that $-s$ is a convex function of $(\tau, e)$, and satisfies

$$ de + pd\tau = Tds. \qquad (2) $$

This equivalently means that the internal energy is convex with respect to specific volume and entropy, hence the sound speed $c_s$ defined by

$$ c_s^2 = \left( \frac{\partial p}{\partial \rho} \right)_s \qquad (3) $$

is positive and ensures the hyperbolicity of the system. Assuming smooth solutions of (1), one can show that they satisfy the following equation of conservation for the entropy

$$ \partial_t (\rho s) + \nabla \cdot (\rho s \mathbf{u}) = 0. \qquad (4) $$

For the non-conservative form of the MHD equations, this holds for any value of the divergence of the magnetic field $\nabla \cdot \mathbf{B}$. Assuming that the divergence of the magnetic field is zero at an initial time $\nabla \cdot \mathbf{B} = 0$, it remains zero at all time following the divergence of the induction equation,

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

The free divergence constraint is therefore a consequence of the induction equation and not a dynamical constraint.

Equivalently, by adding terms proportional to $\nabla \cdot \mathbf{B}$ in the momentum and energy equations (see Appendix A), one can obtain a conservative form for the MHD equations

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + \sigma - \mathbf{B} \otimes \mathbf{B}) = 0, $$
$$ \begin{aligned} \partial_t(\rho E) + \nabla \cdot (\rho E \mathbf{u} + \sigma \cdot \mathbf{u} - (\mathbf{B} \cdot \mathbf{u})\mathbf{B}) &= 0, \\ \partial_t \mathbf{B} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) &= 0, \end{aligned} \tag{6} $$

with $\sigma = (p + \mathbf{B}^2/2)\mathbf{I}$ and $E = e + \mathbf{u}^2/2 + \mathbf{B}^2/(2\rho)$. Assuming smooth solutions of (6), one can show that they satisfy the following equation for the evolution of the entropy by subtracting the evolution of the kinetic and magnetic energy from the evolution of the total energy

$$ \partial_t(\rho s) + \nabla \cdot (\rho s \mathbf{u}) = -\frac{\mathbf{u} \cdot \mathbf{B}}{T} \nabla \cdot \mathbf{B}, \tag{7} $$

which is compatible with entropy conservation only when $\nabla \cdot \mathbf{B} = 0$ in contrast to the non-conservative form presented above [26]. This shows that the entropy balance is closely related to the free divergence constraint for the conservative MHD equations.

In the case of discontinuities such as shocks and in order to ensure dissipation, the second law of thermodynamics must be enforced and implies the entropy inequality

$$ \partial_t(\rho s) + \nabla \cdot (\rho s \mathbf{u}) \geq 0, \tag{8} $$

After discretization, truncation errors on the $\nabla \cdot \mathbf{B}$ source term in Eq. (7) therefore leads to some issues in order to obtain an entropy satisfying numerical scheme ensuring a discrete version of Eq. (8).

In the next sections (3, 4, 5 and 6), we introduce a new fully-conservative solver relying on a splitting between a magneto-acoustic and a transport subsystem. This solver is entropy satisfying and is not compatible with Eq. (8). We then introduce in Sect. 7 an entropy correction following [11] ensuring that the modified scheme is compatible with Eq. (8) while breaking the magnetic field conservation but maintaining the momentum and energy conservation.

### 3. Magneto-acoustic/transport splitting

Similarly to [14], we propose the following splitting of the conservative MHD equations into a magneto-acoustic sub-system

$$ \begin{aligned} \partial_t \rho + \rho \nabla \cdot \mathbf{u} &= 0, \\ \partial_t (\rho \mathbf{u}) + \rho \mathbf{u} \nabla \cdot \mathbf{u} + \nabla \cdot (\sigma \mathbf{B} \otimes \mathbf{B}) &= 0, \\ \partial_t (\rho E) + \rho E \nabla \cdot \mathbf{u} + \nabla \cdot (\sigma \mathbf{u} - (\mathbf{B} \cdot \mathbf{u})\mathbf{B}) &= 0, \\ \partial_t \mathbf{B} + \mathbf{B} \nabla \cdot \mathbf{u} - \nabla \cdot (\mathbf{B} \otimes \mathbf{u}) &= 0, \end{aligned} \tag{9} $$

and a transport sub-system

$$ \begin{aligned} \partial_t \rho + \mathbf{u} \cdot \nabla \rho &= 0, \\ \partial_t (\rho \mathbf{u}) + \mathbf{u} \cdot \nabla (\rho \mathbf{u}) &= 0, \\ \partial_t (\rho E) + \mathbf{u} \cdot \nabla (\rho E) &= 0, \\ \partial_t \mathbf{B} + \mathbf{u} \cdot \nabla (\mathbf{B}) &= 0. \end{aligned} \tag{10} $$

We emphasize that all the components of the magnetic field are transported at velocity $\mathbf{u}$ in the transport sub-system. We then propose to approximate the solution of Eq. (6) by approximating the solutions of the two sub-systems (9) and (10), i.e. for a discrete state $\mathbf{U}_i^n = (\rho, \rho \mathbf{u}, \rho E, \mathbf{B}_i^n)$ in a cell $\Omega_i$ at time $t^n$, the update to $\mathbf{U}_i^{n+1}$ is first an update from $\mathbf{U}_i^n$ to $\mathbf{U}_i^{n+1-}$ by approximating the solution of (9), then an update from $\mathbf{U}_i^{n+1-}$ to $\mathbf{U}_i^{n+1}$ by approximating the solution of (10). We present in Sect. 4 and in Sect. 5 the discretization and the entropy analysis for each sub-system respectively.

### 4. Relaxation approximation of the magneto-acoustic sub-system

The relaxation approximation of the magneto-acoustic sub-system and the associated entropy analysis in Sect. 7 heavily relies on earlier works by [10,11]. We highlight two main differences in our approach: we keep in the analysis gradients of the magnetic field perpendicular to the interface that appears in the multi-dimensional case and we propose a different choice of relaxation parameters in the 5-wave solver to ensure the strict hyperbolicity of the relaxed system.

The multi-dimensional scheme will be obtained by taking advantage of the rotational invariance of the magneto-acoustic sub-system, following the lines of [18]. We, therefore, rewrite sub-system (9) in 1D, and simplify it by using the density evolution equation

$$ \begin{aligned} \rho \partial_t \tau - \partial_x u &= 0, \\ \rho \partial_t \mathbf{u} + \partial_x (\sigma \mathbf{e}_x - \mathbf{B}_x \mathbf{B}) &= 0, \\ \rho \partial_t E + \partial_x (\sigma u_x - (\mathbf{B} \cdot \mathbf{u}) \mathbf{B}_x) &= 0, \\ \rho \partial_t (\tau \mathbf{B}) - \partial_x (\mathbf{B}_x \mathbf{u}) &= 0, \end{aligned} \tag{11} $$
with $e_x$, the unit vector normal to the interface, $B_x$, $B_y$, and $B_z$ the components of the magnetic field and $u_x$, $u_y$, and $u_z$ the components of the velocity field. The eigenvalues of this sub-system are given by

$$-u, 0, \pm c_{ms}, \pm c_{ma}, \pm c_{mf} \qquad (12)$$

with $c_{ma}$, the magnetic Alfvén speed, $c_{ms}$, the slow magnetosonic speed, $c_{mf}$, the fast magnetosonic speed defined by

$$c_{ma} = \frac{|B_x|}{\sqrt{\rho}}, \\ c_{ms}^2 = \frac{1}{2} \left( c_s^2 + \frac{\mathbf{B}^2}{\rho} - \sqrt{\left(c_s^2 + \frac{\mathbf{B}^2}{\rho}\right)^2 - 4c_s^2 c_{ma}} \right), \\ c_{mf}^2 = \frac{1}{2} \left( c_s^2 + \frac{\mathbf{B}^2}{\rho} + \sqrt{\left(c_s^2 + \frac{\mathbf{B}^2}{\rho}\right)^2 - 4c_s^2 c_{ma}} \right). \qquad (13)$$

We then introduce a relaxation procedure [10,14] with the relaxation pressures $\pi_u$ playing the role of the fluxes in the impulsion equation and the relaxation variable $r$ playing the role of the density in front of the time derivatives

$$r \partial_t \tau - \partial_x u = 0, \\ r \partial_t \mathbf{u} + \partial_x \pi_u = 0, \\ r \partial_t E + \partial_x (\pi_u \cdot \mathbf{u}) = 0, \\ r \partial_t (\tau \mathbf{B}) - \partial_x (B_x \mathbf{u}) = 0, \qquad (14)$$

with the following equations for the relaxation variables

$$\partial_t r = \frac{\rho - r}{\epsilon}, \\ r \partial_t \pi_u + (c_b^2 + b_y^2 + b_z^2) \partial_x u - c_a b_y \partial_x v - c_a b_z \partial_x w + d_x \partial_x B_x = \frac{\sigma - B_x^2 - \pi_u}{\epsilon}, \\ r \partial_t \pi_v - c_a b_y \partial_x u + c_a^2 \partial_x v + d_y \partial_x B_x = \frac{-B_x B_y - \pi_v}{\epsilon}, \\ r \partial_t \pi_w - c_a b_z \partial_x u + c_a^2 \partial_x w + d_z \partial_x B_x = \frac{-B_x B_z - \pi_w}{\epsilon}. \qquad (15)$$

The parameters $c_a$, $c_b$, $b_y$, $b_z$ play the role of approximations of $\sqrt{\rho}|B_x|$, $\rho c_s$, sign($B_x$) $\sqrt{\rho}B_y$, sign($B_x$) $\sqrt{\rho}B_z$, respectively, as in [10]. The extra parameters $d_x$, $d_y$, $d_z$ are linked to the possibility of a non-constant $B_x$ in the magneto-acoustic sub-system and play the role of approximations of $2B_x u/\tau + \mathbf{u} \cdot \mathbf{B}(\partial_c p - 1/\tau)$, $(B_x v + B_y u)/\tau$, and $(B_x w + B_z u)/\tau$, respectively. If these extra parameters are fixed to zero, the relaxation equations for $\pi_u$ is the Lagrangian form of the relaxation equations used in [10]. By replacing all these parameters exactly by the quantities they approximate, Eq. (15) reduces to the evolution equation of $\sigma - B_x^2$, $-B_x B_y$, and $-B_x B_z$ in the limit $\epsilon \to \infty$. In order to obtain the same Riemann invariants as [10], we fix $d_x$, $d_y$, and $d_z$ to zero and the other constants are evolved with

$$\partial_t c_a = \partial_t c_b = \partial_t b_y = \partial_t b_z = 0. \qquad (16)$$

In the limit $\epsilon \to 0$, the relaxation equations in Eq. (15) ensures that $r \to \rho$, $\pi_u \to \sigma - B_x^2$, $\pi_v \to -B_x B_y$, and $\pi_w \to -B_x B_z$. In this limit, Eq. (14) is then equivalent to Eq. (11). A classical approach to achieve the limit $\epsilon \to 0$ numerically is to first enforce the equilibrium relations $r = \rho$ and $\pi_u = \sigma e_x - B_x \mathbf{B}$ at time $t^n$ and then solve (14) and (15) without the relaxation source terms. Using $L \equiv r/\rho$, the full system without the relaxation source term is

$$\begin{aligned} & \partial_t L - \partial_x u = 0, \\ & \partial_t (\rho L \mathbf{u}) + \partial_x \pi_u = 0, \\ & \partial_t (\rho L E) + \partial_x (\pi_u \cdot \mathbf{u}) = 0, \\ & \partial_t (L \mathbf{B}) - \partial_x (B_x \mathbf{u}) = 0, \\ & \partial_t (\rho L) = 0, \\ & \partial_t (\rho L \pi_u) + (c_b^2 + b_y^2 + b_z^2) \partial_x u - c_a b_y \partial_x v - c_a b_z \partial_x w = 0, \\ & \partial_t (\rho L \pi_v) - c_a b_y \partial_x u + c_a^2 \partial_x v = 0, \\ & \partial_t (\rho L \pi_w) - c_a b_z \partial_x u + c_a^2 \partial_x w = 0. \end{aligned} \qquad (17)$$
After some tedious algebra, one can compute the eigenvalues of this system of 16 equations (including $\partial_t c_a = \partial_t c_b = \partial_t b_y = \partial_t b_z = 0$),

$$-u/L, 0, \pm c_{rs}/(\rho L), \pm c_{ra}/(\rho L), \pm c_{rf}/(\rho L) \quad (18)$$

with

$$c_{ra} = c_a, \\ c_{rs}^2 = \frac{1}{2} \left( c_b^2 + c_a^2 + b_y^2 + b_z^2 - \sqrt{(c_b^2 + c_a^2 + b_y^2 + b_z^2)^2 - 4c_a^2 c_b^2} \right), \\ c_{rf}^2 = \frac{1}{2} \left( c_b^2 + c_a^2 + b_y^2 + b_z^2 + \sqrt{(c_b^2 + c_a^2 + b_y^2 + b_z^2)^2 - 4c_a^2 c_b^2} \right). \quad (19)$$

The central wave at zero velocity has multiplicity 9. All the waves are linearly degenerate. Similarly to [10], $c_{rs} \le c_a \le c_{rf}$, $c_{rs} \le c_b \le c_{rf}$ and the eigenvalues of (17) match the eigenvalues of (11) for $c_a = \sqrt{\rho}|B_x|$, $c_b = \rho c_s$, $b_y = \text{sign}(B_x)\sqrt{\rho}B_y$, $b_z = \text{sign}(B_x)\sqrt{\rho}B_z$. Similarly to [10], a Chapman-Enskog analysis can be performed on the relaxation equations which leads to the following stability conditions

$$\frac{1}{\rho} - \frac{B_x^2}{c_a^2} \ge 0, \\ c_b^2 - \rho^2 c_s^2 \ge 0, \\ (c_b^2 - \rho^2 c_s^2) \left( \frac{1}{\rho} - \frac{B_x^2}{c_a^2} \right) \ge \left( B_y - \frac{B_x b_y}{c_a} \right)^2 + \left( B_z - \frac{B_x b_z}{c_a} \right)^2, \quad (20)$$

in order to ensure positive eigenvalues of the entropy diffusion matrix.

**The 3+1 and 5+1 wave solver** The solution of the Riemann problem associated to (17) contains 7+1 waves in the general case, 7 waves that are identical to a Lagrangian version of the 1D relaxation solver presented in [10] to which we add a wave at $-u/L$ associated to $B_x$. Similarly to [11] we can design an approximate Riemann solver with 5+1 waves by choosing $b_y = b_z = 0$, or with 3+1 waves by choosing in addition $c_a = c_b = c$. The 5+1 wave solver is a good compromise between accuracy and computational cost and we will use this approximation for now on.

We now look for strong Riemann invariants for the different waves by finding quantities transported at the corresponding wave speed [18]. $B_x$ is a strong Riemann invariant associated to the wave at $-u/L$. Note that $B_x$ is not constant but advected at velocity $-u/L$. $B_x$ has to be understood as evaluated locally, upwind relative to the wave $-u/L$. $c_a$ and $c_b$ are strong Riemann invariants for the central wave with

$$\frac{1}{\rho} + \frac{\pi_u}{c_b^2} + \frac{B_y}{\rho} + \frac{B_x}{c_a^2} \pi_v, \frac{B_z}{\rho} + \frac{B_x}{c_a^2} \pi_w, e + \frac{\mathbf{B}^2}{2\rho} - \frac{\pi_u^2}{2c_b^2} - \frac{\pi_v^2 + \pi_w^2}{2c_a^2}. \quad (21)$$

Similarly to [11], there are six strong Riemann invariants for the left and right waves $\pi_u + c_u u$ and $\pi_u - c_u u$, respectively, in which we have defined $c_u = (c_b, c_a, c_a)$. Strong Riemann invariants for a given wave are weak Riemann invariants for the other waves. They are, therefore, weak Riemann invariants for the central wave, hence, $u$ and $\pi_u$ take the same value on the left and right of this wave that we shall define as $u^*$ and $\pi_u^*$ respectively. By using the weak Riemann invariants, we get

$$\mathbf{u}^* = \frac{c_{u,l}\mathbf{u}_l + c_{u,r}\mathbf{u}_r + \pi_{u,l} - \pi_{u,r}}{c_{u,l} + c_{u,r}}, \\ \pi_{u}^{*} = \frac{c_{u,r}\pi_{u,l} + c_{u,l}\pi_{u,r} + c_{u,l}c_{u,r}(\mathbf{u}_l - \mathbf{u}_r)}{c_{u,l} + c_{u,r}}. \quad (22)$$

Then one has

$$B_x(x, t) = \begin{cases} B_{x,l} & \text{if } x/t < -u/L \\ B_{x,r} & \text{if } x/t > -u/L, \end{cases} \quad (23)$$

hence, at the interface, we define $B_x^{-u^*} = B_x(0, t)$ with

$$B_x^{-u^*} = \begin{cases} B_{x,l} & \text{if } u^* < 0 \\ B_{x,r} & \text{if } u^* > 0. \end{cases} \quad (24)$$

The other intermediate states, e.g. $\tau_{l,r}^*$ and $e_{l,r}^*$ can be obtained by using (21), but are not needed for deriving the update of the numerical scheme. The discrete numerical scheme for the magneto-acoustic sub-system is then given by

$$L_i^{n+1-} = 1 + \frac{\Delta t}{\Delta x}(u_{i+1/2}^* - u_{i-1/2}^*),$$
$$
\begin{align}
\rho_i^{n+1} L_i^{n+1} &= \rho_i^n, \nonumber \\
\rho_i^{n+1} \mathbf{u}_i^{n+1} L_i^{n+1} &= \rho_i^n \mathbf{u}_i^n - \frac{\Delta t}{\Delta x} (\pi_{\mathbf{u},i+1/2}^* - \pi_{\mathbf{u},i-1/2}^*), \nonumber \\
\rho_i^{n+1} E_i^{n+1} L_i^{n+1} &= \rho_i^n E_i^n - \frac{\Delta t}{\Delta x} (\pi_{\mathbf{u},i+1/2}^* \cdot \mathbf{u}_{i+1/2}^* - \pi_{\mathbf{u},i-1/2}^* \cdot \mathbf{u}_{i-1/2}^*), \nonumber \\
\mathbf{B}_i^{n+1} L_i^{n+1} &= \mathbf{B}_i^n + \frac{\Delta t}{\Delta x} (B_{x,i+1/2}^{-u^*} \mathbf{u}_{i+1/2}^* - B_{x,i-1/2}^{-u^*} \mathbf{u}_{i-1/2}^*), \tag{25}
\end{align}
$$

with the CFL condition for this scheme

$$
\max_{i \in \mathbb{Z}} \left( \frac{c_{rf,i}}{\rho_i} \right) \Delta t \le \frac{\Delta x}{2}. \qquad (26)
$$

5. Transport sub-system

The transport sub-system is a quasi-hyperbolic system that only involves the transport of conservative variables with the velocity
u. We choose to approximate the solution of the 1D version of (10) thanks to a standard upwind Finite-Volume approximation for
U = (ρ, ρu, ρE, B) by discretizing

$$
\frac{\partial \mathbf{U}}{\partial t} + u \frac{\partial \mathbf{U}}{\partial x} = \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial (u\mathbf{U})}{\partial x} - \mathbf{U} \frac{\partial u}{\partial x} = 0, \qquad (27)
$$

with

$$
\mathbf{U}_i^{n+1} = \mathbf{U}_i^{n+1-} - \frac{\Delta t}{\Delta x} (u_{i+1/2}^* \mathbf{U}_{i+1/2} - u_{i-1/2}^* \mathbf{U}_{i-1/2}) + \frac{\Delta t}{\Delta x} \mathbf{U}_i^{n+1-} (u_{i+1/2}^* - u_{i-1/2}^*), \quad (28)
$$

with two possible choices of discretization for the interface states $\mathbf{U}_{i-1/2}$ and $\mathbf{U}_{i+1/2}$. The first choice

$$
\mathbf{U}_{i+1/2} = \begin{cases} \mathbf{U}_i^{n+1-} & \text{if } u_{i+1/2}^* \ge 0, \\ \mathbf{U}_i^{n+1-} & \text{if } u_{i+1/2}^* \le 0, \end{cases} \qquad (29)
$$

leads to a magneto-acoustic+transport scheme of stencil 2 similar to [14]. The second choice

$$
\mathbf{U}_{i+1/2} = \begin{cases} \mathbf{U}_i^n & \text{if } u_{i+1/2}^* \ge 0, \\ \mathbf{U}_{i+1}^n & \text{if } u_{i+1/2}^* \le 0, \end{cases} \qquad (30)
$$

leads to a magneto-acoustic+transport scheme of stencil 1 similar to [16]. We will refer to these choices of discretization as “stencl 1”
and “stencl 2” in the rest of the paper. In both cases and using the notation $u^{\pm} = \frac{u \pm |u|}{2}$, the CFL condition of the transport sub-system
is given by

$$
\max_{i \in \mathbb{Z}} ((u_{i-1/2}^*)^+ - (u_{i+1/2}^*)^-) \Delta t \le \Delta x. \qquad (31)
$$

The transport can also be written in the form

$$
\mathbf{U}_i^{n+1} = \mathbf{U}_i^{n+1-} L_i^{n+1-} - \frac{\Delta t}{\Delta x} (u_{i+1/2}^* \mathbf{U}_{i+1/2} - u_{i-1/2}^* \mathbf{U}_{i-1/2}). \quad (32)
$$

6. Magneto-acoustic+transport scheme

The global scheme is given by

$$
\begin{align}
\rho_i^{n+1} &= \rho_i^n - \frac{\Delta t}{\Delta x} (\rho_{i+1/2} u_{i+1/2}^* - \rho_{i-1/2} u_{i-1/2}^*), \nonumber \\
(\rho \mathbf{u})_i^{n+1} &= (\rho \mathbf{u})_i^n - \frac{\Delta t}{\Delta x} ((\rho \mathbf{u})_{i+1/2} u_{i+1/2}^* + \pi_{\mathbf{u},i+1/2}^* \nonumber \\
&\quad - (\rho \mathbf{u})_{i-1/2} u_{i-1/2}^* - \pi_{\mathbf{u},i-1/2}^*), \nonumber \\
(\rho E)_i^{n+1} &= (\rho E)_i^n - \frac{\Delta t}{\Delta x} ((\rho E)_{i+1/2} u_{i+1/2}^* + \pi_{\mathbf{u},i+1/2}^* \cdot \mathbf{u}_{i+1/2}^* \nonumber \\
&\quad - (\rho E)_{i-1/2} u_{i-1/2}^* - \pi_{\mathbf{u},i-1/2}^* \cdot \mathbf{u}_{i-1/2}^*), \nonumber \\
\mathbf{B}_i^{n+1} &= \mathbf{B}_i^n - \frac{\Delta t}{\Delta x} (\mathbf{B}_{i+1/2} u_{i+1/2}^* - \mathbf{B}_{x,i+1/2}^{-u^*} \mathbf{u}_{i+1/2}^* \nonumber \\
&\quad - \mathbf{B}_{i-1/2} u_{i-1/2}^* + \mathbf{B}_{x,i-1/2}^{-u^*} \mathbf{u}_{i-1/2}^*). \tag{33}
\end{align}
$$

The global scheme of stencil 2 is stable under the most restrictive CFL condition between the magneto-acoustic and transport sub-systems. The scheme of stencil 1 is stable under a CFL condition involving the sum of the speeds of the magneto-acoustic and transport subsystem as demonstrated in [16] and in Sect. 7.
## 7. Entropy analysis

In this section, we first introduce under which conditions the 1D relaxation solver is entropy-satisfying. For a non-constant $B_x$ in a multi-dimensional setup, it is clear that the fully-conservative solver is not entropy-satisfying: on the $-u/L$ wave, $B_x$ is the only quantity that jumps, hence, induces a jump in internal energy because of the last Riemann invariant in (21). Similarly to [11], an entropy satisfying solver will require the introduction of an entropic correction on the induction equation to get a symmetric version of the MHD equations. We will present the multi-dimensional entropy-satisfying solver at the end of the section.

The choice of the relaxation parameter $c = c_a = c_b$ for the 3+1 wave approximate Riemann solver and $c_a, c_b$ for the 5+1 wave solver is made to ensure that the solver is entropy satisfying for a constant $B_x$ in 1D. If for all intermediate states $U_{i,r}^*$, one has $\tau_{i,r}^* > 0$ and

$$ (\rho^2 c_s^2)_{*,l,r} \le c_b^2, \\ \tau_{l,r}^* - \frac{B_x^2}{c_a^2} \ge 0, \\ (B_{y,l,r}^2 + B_{z,l,r}^2) \le (c_b^2 - (\rho^2 c_s^2)_{*,l,r}) \left( \tau_{l,r}^* - \frac{B_x^2}{c_a^2} \right), \qquad (34) $$

with $(\rho^2 c_s^2)_{*,l,r} \equiv \sup_{\rho \in (\rho_*, \rho_l, \rho_r)} (\rho^2 c_s^2(\rho, s_{l,r}))$, there exists a numerical flux function $q_{i+1/2}^n = q(U_i^n, U_{i+1}^n)$, consistent with zero (see [14]) such that

$$ \rho_i^{n+1} s(U_i^{n+1}) - \rho_i^n s(U_i^n) + \frac{\Delta t}{\Delta x} (q_{i+1/2}^n + (\rho s)_{i+1/2} u_{i+1/2}^* - q_{i-1/2}^n - (\rho s)_{i-1/2} u_{i-1/2}^*) \ge 0. \qquad (35) $$

Following [11], optimal choices of $c_a$ and $c_b$ for smooth solutions are given by

$$ c_a^2 = \rho(B_x^2 + |B_x|\sqrt{B_y^2 + B_z^2}) \\ c_b^2 = \rho^2 c_s^2 + \rho(B_y^2 + B_z^2 + |B_x|\sqrt{B_y^2 + B_z^2}) \qquad (36) $$

for the 5+1 wave solver and $c = \rho c_{mf}$ for the 3+1 wave solver. Optimal choices for discontinuous solutions are given in [11], however, in all the tests performed in Sect. 8 the smooth version has been sufficient to ensure stability and is therefore preferred for its low computational cost. As noted in [11], the diffusion of the 5+1 solver is zero when $B_x = 0$ or $B_y^2 + B_z^2 = 0$ which means that the solver is exact in these conditions. We, however, point out that this is exactly where the MHD system is not strictly hyperbolic with $c_{ma} = 0$ for $B_x = 0$ and $c_{ma} = c_{ms}$ for $B_y^2 + B_z^2 = 0$. Therefore, in practice, we employ a more diffusing approximation for the choices of $c_a$ and $c_b$ by using the following inequality $|B_x|\sqrt{B_y^2 + B_z^2} \le (B_x^2 + B_y^2 + B_z^2)/2$:

$$ c_a^2 = \rho(B_x^2 + (B_x^2 + B_y^2 + B_z^2)/2) \\ c_b^2 = \rho^2 c_s^2 + \rho(B_y^2 + B_z^2 + (B_x^2 + B_y^2 + B_z^2)/2) \qquad (37) $$

to ensure the use of a stable strictly hyperbolic approximation even when $B_x$ or $B_y^2 + B_z^2$ vanishes. It also helps with the isotropy of the numerical diffusion whenever there is a large difference between the normal and transverse magnetic intensity, avoiding the generation of spurious patterns. We decompose the proof of the entropy analysis of the global scheme into an entropy analysis of each sub-system, magneto-acoustic and transport, respectively.

### 7.1. Entropy analysis of the magneto-acoustic sub-system in 1D

**Proposition 1.** Let $s_{l,r} = s(\tau_{l,r}, e_{l,r})$. If the inequality

$$ e_{l,r}^* \ge e(\tau_{l,r}^*, s_{l,r}) \qquad (38) $$

is verified, there exists a numerical flux function $q_{i+1/2}^n = q(U_i^n, U_{i+1}^n)$, consistent with zero such that

$$ L_i^{n+1} - \rho_i^{n+1} s(\tau_i^{n+1}, e_i^{n+1}) - \rho_i^n s(\tau_i^n, e_i^n) + \frac{\Delta t}{\Delta x} (q_{i+1/2}^n - q_{i-1/2}^n) \ge 0 \qquad (39) $$

**Proof.** According to Eq. (2), at fixed $\tau$, $e(\tau, s)$ is an increasing function of $s$, hence $e(\tau_{l,r}^*, s_{l,r}^*) \ge e(\tau_{l,r}^*, s_{l,r})$ implies $s_{l,r}^* \ge s_{l,r}$. This inequality then implies that for any $c > 0$

$$ 0 \ge -c(s_l^* - s_l) + c(s_r - s_r^*) \qquad (40) $$
which is consistent with the integral form of the entropy inequality $\partial_t(s(\tau, e)) \ge 0$. As in [14], this implies the existence of $q_{i+1/2}^n = q(U_i^n, U_{i+1}^n)$ such that

$$s(\tau_i^{n+1-}, e_i^{n+1-}) - s(\tau_i^n, e_i^n) + \tau_i^n \frac{\Delta t}{\Delta x} (q_{i+1/2}^n - q_{i-1/2}^n) \ge 0 \quad (41)$$

The inequality (39) follows from $L_i^{n+1-} \rho_i^{n+1-} = \rho_i^n$. $\square$

**Proposition 2.** The 5+1 wave approximate Riemann solver associated to the relaxation (17) of the magneto-acoustic sub-system is positive and satisfies all discrete entropy inequalities whenever for all intermediate states $U_{l,r}^*$, $\tau_{l,r}^*$ are positive and the inequalities (34) are verified.

**Proof.** According to (21), the 5+1 wave relaxation Riemann problem has the same Riemann invariants as [11] apart from the addition of $B_x$ as a strong Riemann invariant of the $-u/L$ wave. $B_x$ has therefore to be understood as evaluated locally according to (23). By introducing the decomposition into elementary dissipation terms similarly as in [27], using the Riemann invariants (21) and defining $\sigma(\mathbf{U}) = p(\tau, s = s_{l,r}) + \mathbf{B}^2/2$, one can show that

$$e(\tau_{l,r}^*, s_{l,r}) - e_{l,r}^* = D_0(\mathbf{U}_{l,r}^*, \mathbf{U}_{l,r}) - \frac{1}{2} \left| \frac{\sigma(\mathbf{U}_{l,r}^*)\mathbf{n} - \mathbf{B}_x\mathbf{B}^* - \tau_u^*}{c_u} \right|^2, \quad (42)$$

with $D_0$ the dissipation associated to the central wave given by

$$D_0(\mathbf{U}_{l,r}^*, \mathbf{U}_{l,r}) = e(\tau_{l,r}^*, s_{l,r}) - e(\tau_{l,r}, s_{l,r}) + p(\tau_{l,r}^*, s_{l,r}) (\tau_{l,r}^* - \tau_{l,r}) \\ + \frac{1}{2c_b^2} (\sigma(\mathbf{U}_{l,r}^*) - \sigma(\mathbf{U}_{l,r}))^2 \\ - (\tau_{l,r} - B_x^2/c_a^2) \frac{1}{2} |\mathbf{B}^* - \mathbf{B}_{l,r}|^2. \quad (43)$$

The proof of Proposition 2 then follows directly from the entropy analysis of [10] who showed that under (34) and by using (42), the inequality (38) is verified. $\square$

The final part of the analysis requires to give the conditions under which the relaxation approximation is positive for the intermediate states of the specific volume $\tau_{l,r}^* > 0$. These conditions for the relaxation parameters are provided in proposition 3.3 of [11], however we do not explicitly specify them here because we will use a less restrictive choice with Eq. (37) which seems sufficient in practice in all the numerical tests performed in Sect. 8.

## 7.2. Entropy analysis of the transport sub-system in 1D

By using $u^\pm = \frac{u\pm|u|}{2}$, the transport step of the global scheme of stencil 2 can be written in the form

$$\mathbf{U}_i^{n+1} = \frac{\Delta t}{\Delta x} u_{i-1/2}^{*,+} \mathbf{U}_{i-1}^{n+1-} - \frac{\Delta t}{\Delta x} u_{i+1/2}^{*,-} \mathbf{U}_{i+1}^{n+1-} + \left(1 - \frac{\Delta t}{\Delta x} (u_{i-1/2}^{*,+} - u_{i+1/2}^{*,-})\right) \mathbf{U}_i^{n+1-}, \quad (44)$$

hence $\mathbf{U}_i^{n+1}$ is a convex combination of $\mathbf{U}_{i-1}^{n+1-}$, $\mathbf{U}_i^{n+1-}$ and $\mathbf{U}_{i+1}^{n+1-}$ as their pre-factors are positive and sum to 1. By convexity of the function $\mathbf{U} \to -\rho s(\mathbf{U})$

$$\rho_i^{n+1} s(\mathbf{U}_i^{n+1}) \ge \rho_i^{n+1-} L_i^{n+1-} s(\mathbf{U}_i^{n+1-}) - \frac{\Delta t}{\Delta x} ((\rho s)_{i+1/2} u_{i+1/2}^* - (\rho s)_{i-1/2} u_{i-1/2}^*). \quad (45)$$

By combining, the inequalities (39) and (45) we obtain the inequality (35). Following [16], the global scheme of stencil 1 can be written in the form

$$\mathbf{U}_i^{n+1} = \alpha_i \mathbf{U}_i^A + (1 - \alpha_i) \mathbf{U}_i^T \quad (46)$$

for any $\alpha_i \in ]0, 1[$ and

$$\begin{aligned} \mathbf{U}_i^A &= \mathbf{U}_i^n + \frac{1}{\alpha_i \Delta x} (\mathbf{U}_i^{n+1-} L_i^{n+1-} - \mathbf{U}_i^n), \\ \mathbf{U}_i^T &= \mathbf{U}_i^n - \frac{1}{1 - \alpha_i \Delta x} \frac{\Delta t}{\Delta x} (u_{i+1/2}^* \mathbf{U}_{i+1/2} - u_{i-1/2}^* \mathbf{U}_{i-1/2}), \end{aligned} \quad (47)$$

with $\mathbf{U}_i^A$ corresponding to a magneto-acoustic update with $\Delta t^A = \frac{1}{\alpha_i} \Delta t$ and $\mathbf{U}_i^T$ corresponding to a conservative transport update also with $\Delta t^T = \frac{1}{1-\alpha_i} \Delta t$. Following [16], $\mathbf{U}_i^T / \rho_i^{n+1}$ can be written as a convex combination of $\mathbf{U}_i^n / \rho_i^n$. Thus, we can also obtain (35) by using the convexity of (46) under the CFL conditions:
$$ \max_{i \in \mathcal{Z}} ((u_{i-1/2}^*)^+ - (u_{i+1/2}^*)^-) \frac{1}{\alpha_i} \Delta t \le \Delta x. \quad (48) $$

$$ \max_{i \in \mathcal{Z}} \left( \frac{c_{rf,i}}{\rho_i} \right) \frac{1}{1-\alpha_i} \Delta t \le \frac{\Delta x}{2}. \quad (49) $$

As the local choice of $\alpha_i$ is free, we can pick it so that both conditions coincide, giving the following condition for the stencil 1 scheme:

$$ \max_{i \in \mathcal{Z}} ((u_{i-1/2}^*)^+ - (u_{i+1/2}^*)^-) + 2 \frac{c_{rf,i}}{\rho_i} \Delta t \le \Delta x. \quad (50) $$

### 7.3. Entropic correction for multi-dimensional MHD

Similarly to [11], we introduce an entropic correction on the induction equation proportional to $\nabla \cdot \mathbf{B}$,

$$ \partial_t \mathbf{B} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) + \mathbf{u} \nabla \cdot \mathbf{B} = 0. \quad (51) $$

The rest of the MHD system is not changed and, of course, (51) is equivalent to the standard form when $\nabla \cdot \mathbf{B} = 0$. Smooth solutions follow the entropy evolution

$$ \partial_t(\rho s) + \nabla \cdot (\rho s \mathbf{u}) = 0. \quad (52) $$

We recall the discretization of the source term as [11] which results in two different values of $B_x^{-u^*}$ at an interface, $B_{x,i+1/2,i}^{-u^*} = B_{x,i}^n$ and $B_{x,i+1/2,r}^{-u^*} = B_{x,i+1}^n$, hence giving a non-conservative discretization of the induction equation with

$$ \mathbf{B}_i^{n+1} = \mathbf{B}_i^n - \frac{\Delta t}{\Delta x} (\mathbf{B}_{i+1/2} u_{i+1/2}^* - \mathbf{B}_{x,i}^n u_{i+1/2}^* - \mathbf{B}_{i-1/2} u_{i-1/2}^* + \mathbf{B}_{x,i}^n u_{i-1/2}^*). \quad (53) $$

With this non-conservative source term, the evolution equation of $B_x$ is simply $\partial_t B_x = 0$ and the system becomes symmetric with an additional wave centered at 0 instead of the $-u/L$ wave [12]. The strong Riemann invariant $B_x$ jumps at 0, similarly to the other Riemann invariant (21). As in [11], the 3+1 and 5+1 approximate Riemann solvers with the non-conservative source term are entropy satisfying with the same proof presented above, $B_x$ simply needs to be understood as evaluated locally with a jump on the central wave.

We emphasize that the normal component of the magnetic field for $B_x^{-u^*}$ in (53) is always the value at cell center $B_{x,i}^n$, both at first and second order. As noted by [28], the source term vanishes for smooth solutions at second order if one uses the reconstructed values at interfaces. The proposed discretization in (53) avoids this problem and can be employed for both 1st and 2nd order.

## 8. Numerical results

In order to take advantage of the fully conservative and entropy-satisfying solvers, we use an hybrid strategy in our simulations by switching between both depending on the regime of the flow. On cells where the plasma beta number $\beta = p/\frac{B^2}{2}$ is inferior to a tunable threshold $\beta_{min}$ or where the local Alfvén number $Al = \sqrt{\rho |u|/|B|}$ is superior to another tunable threshold $Al_{max}$, we use locally the entropy-satisfying solver instead of the fully conservative solver. In all our experiments, we set $\beta_{min} = 10^{-3}$ and $Al_{max} = 10$. The entropy correction is only activated in the specifically designed low-plasma-beta blast problem (see Sect. 8.2.4) and the field loop advection test case (see Sect. 8.2.6).

All the simulations performed in this section are using the stencil 1 solver and a MUSCL-Hancock scheme [29], with second order accuracy in space with states reconstructions at interfaces and second order accuracy in time with a predictor-corrector step at half time-step. We perform the extrapolation on the primitive variables $(\rho, p, \mathbf{u}, \mathbf{B})$ and use the classical minmod limiter in order to ensure the admissibility of the Riemann states. The time-step is computed with $\Delta t = \text{CFL} \times \Delta x/(c_{mf} + |u|)$ with $c_{mf}$ the speed of the fast magneto-acoustic waves. With a standard MUSCL-Hancock scheme the CFL condition with CFL < 0.5 ensures a positive numerical scheme. In practice, this CFL condition is often observed to be too restrictive and we use in all simulations a fixed CFL number of 0.8. We also use an ideal gas equation of state. All numerical experiments were conducted using the one step 5 + 1 waves solver with $c_a$ and $c_b$ given by Eq. (37) to avoid the loss of hyperbolicity of the relaxation whenever $B_x$ or $B_y$ and $B_z$ vanish. We also use the 3 + 1 waves solver for the 2D rotated shock tube test to provide a comparison with the 5 + 1 waves solver.

### 8.1. 1D tests cases

In this section, we reproduce several 1D Riemann problems that were used in [11]. The values of the left and right states, the final time, length of the domain and adiabatic indexes are given in table (54). The simulations were all performed with $\Delta x = 10^{-2}$. The reference solutions were all generated with the 5 + 1 waves solver using $\Delta x = 5 \times 10^{-4}$.
<figure><img src="image_3.png" /><figcaption>Fig. 1. $\rho$ and $B_y$ for the Dai and Woodward shock tube at $t = 0.2, 5 + 1$ waves solver against a reference solution.</figcaption></figure>

<table><thead><tr><th>Test case name, (γ, t<sub>end</sub>, L)</th><th>ρ</th><th>(u, v, w)</th><th>p</th><th>(B<sub>x</sub>, B<sub>y</sub>, B<sub>z</sub>)</th></tr></thead><tbody><tr><td>Dai & Woodward, (5/3, 0.2, 1.1)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.08</td><td>(1.2, 0.01, 0.5)</td><td>0.95</td><td>(4/√4π, 3.6/√4π, 2/√4π)</td></tr><tr><td>R state</td><td>1.0</td><td>(0.0, 0.0, 0.0)</td><td>1.0</td><td>(4/√4π, 4/√4π, 2/√4π)</td></tr><tr><td>Brio & Wu I, (2.0, 0.2, 1.0)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.0</td><td>(0.0, 0.0, 0.0)</td><td>1.0</td><td>(0.65, 1.0, 0.0)</td></tr><tr><td>R state</td><td>0.125</td><td>(0.0, 0.0, 0.0)</td><td>0.1</td><td>(0.65, -1.0, 0.0)</td></tr><tr><td>Brio & Wu II, (2.0, 0.012, 1.4)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.0</td><td>(0.0, 0.0, 0.0)</td><td>1000.0</td><td>(0.0, 1.0, 0.0)</td></tr><tr><td>R state</td><td>0.125</td><td>(0.0, 0.0, 0.0)</td><td>0.1</td><td>(0.0, -1.0, 0.0)</td></tr><tr><td>Slow rarefaction, (5/3, 0.2, 1.0)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.0</td><td>(0.0, 0.0, 0.0)</td><td>2.0</td><td>(1.0, 0.0, 0.0)</td></tr><tr><td>R state</td><td>0.2</td><td>(1.186, 2.967, 0.0)</td><td>0.1368</td><td>(1.0, 1.6405, 0.0)</td></tr><tr><td>Expansion I, (5/3, 0.15, 1.4)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.0</td><td>(-3.1, 0.0, 0.0)</td><td>0.45</td><td>(0.0, 0.5, 0.0)</td></tr><tr><td>R state</td><td>1.0</td><td>(3.1, 0.0, 0.0)</td><td>0.45</td><td>(0.0, 0.5, 0.0)</td></tr><tr><td>Expansion II, (5/3, 0.15, 1.4)</td><td></td><td></td><td></td><td></td></tr><tr><td>L state</td><td>1.0</td><td>(-3.1, 0.0, 0.0)</td><td>0.45</td><td>(1.0, 0.5, 0.0)</td></tr><tr><td>R state</td><td>1.0</td><td>(-3.1, 0.0, 0.0)</td><td>0.45</td><td>(1.0, 0.5, 0.0)</td></tr></tbody></table>

### 8.1.1. Dai-Woodward shock tube

This shock tube configuration was introduced in [30]. During the computation, the solution displays the full eigen-structure of the MHD system as it generates shocks and discontinuities on all fields. We observe in Fig. 1 that our method captures the density and transverse magnetic field robustly, without spurious oscillations. We observe the effect of numerical diffusion smoothing the various waves. A density undershoot is observed at $x \approx 0.7$ and is due to the choice of CFL number 0.8, higher than what the 0.5 allowed by the stability analysis of MUSCL methods. These results are very similar to what is obtained in [11].

### 8.1.2. Brio-Wu shock tube, configuration I

The Brio-Wu shock tube was first introduced in [31]. The solution of this shock tube is composed of shocks, rarefactions, contact discontinuities and a compound wave, in this case a discontinuity attached to a slow rarefaction. In Fig. 2, we can see that our solver captures all features of the solution of this Riemann problem. The effect of diffusion is mainly observed on the $x \approx 0.6$ shock and the density peak around $x \approx 0.45$ as it is a very fine feature. At the same location, the low-resolution result does present a smoothed bump. These results are very similar to what is obtained in [11]. Note that as in [11], the slow shock position does not seem to be well captured at low resolution.

### 8.1.3. Brio-Wu shock tube, configuration II

The second Riemann problem from [31] also involves a complex wave structure but with a high magneto-acoustic Mach number. In Fig. 3, we observe that our solver captures all features of the shock tube, similarly to the results of [11]. The effect of diffusion is
<figure><img src="image_4.png" /><figcaption>Fig. 2. $\rho$ and $B_y$ for the Brio and Wu -I- shock tube at $t = 0.2$, 5 + 1 waves solver against a reference solution.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. $\rho$ and $B_y$ for the Brio and Wu -II- shock tube at $t = 0.012$, 5 + 1 waves solver against a reference solution.</figcaption></figure>

mainly observed at $x \approx 1.05$ where a discontinuity and an undershoot are observed on the high resolution plot. This corresponds to the smoothed dip observed in the low-resolution solution.

### 8.1.4. Slow rarefaction tube

This test has been first proposed in [32]. It involves a sonic point, where the slow magneto-acoustic speed equals the fluid velocity. This feature is problematic for linearized method like the Roe solver, but our scheme is stable as we can see in Fig. 4, just like the resolution shown in [11]. The $x \approx 0.75$ dip and $x \approx 0.85$ bump present on the high-resolution line are smoothed but still present on the low-resolution solution.

### 8.1.5. Expansion problem, configuration I

This test is taken from [33]. It consists of two out-going rarefaction separating a low density region that is difficult to tackle in a stable manner. Our solver is able to simulate this region as we can see in Fig. 5. The effect of numerical diffusion on the sharpness of the $x = 0.5$ density and magnetic field dip is visually enhanced by the use of the log scale. Similar results are found in [11].

### 8.1.6. Expansion problem, configuration II

This test is a modification of 8.1.5 suggested by [11] where we simply set $B_x = 1.0$ instead of 0. Taking $B_x$ nonzero causes the thermal pressure to be low in the central region which can be hard to tackle robustly. Nevertheless, we can see in Fig. 6 that our method is stable and provides results that are very similar to the ones presented in [11].
<figure><img src="image_6.png" /><figcaption>Fig. 4. $\rho$ and $B_y$ for the slow rarefaction tube at $t = 0.2, 5 + 1$ waves solver against a reference solution.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. $\rho$ and $B_y$ for the expansion-I-tube at $t = 0.15, 5 + 1$ waves solver against a reference solution. logscale on the y-axis.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. $\rho$ and $B_y$ for the expansion-II-tube at $t = 0.15, 5 + 1$ waves solver against a reference solution. logscale on the y-axis.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7. Density map of the Orszag-Tang vortex at t = 0.5 s. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

## 8.2. 2D tests cases

All 2D test cases are using $\Delta x = \Delta y = \frac{1}{256}$. We also have tested all the resolutions between 64 and 2048 without any issue to report. In all 2D setups, the quantity r always refers to the distance from the center of the domain.

### 8.2.1. Orszag-Tang vortex

The Orszag-Tang vortex test case was first introduced by [34] and has become a standard multi-dimensional benchmark case for ideal MHD. The dynamic of this vortex involves the formation of shocks as well as interactions between them which are challenging to simulate robustly. For instance, 1D solvers like HLLD straightforwardly extended to 2D fail at this task. We recall that this problem takes place in the $[0 : 1]^2$ periodic domain with initial data:

$$ \rho(x, y) = 25/36\pi, $$

$$ \rho(x, y) = 5/12\pi, $$

$$ \mathbf{u}(x, y) = \begin{pmatrix} -\sin 2\pi y \\ \sin 2\pi x \end{pmatrix}, $$

$$ \mathbf{B}(x, y) = 1/\sqrt{4\pi} \begin{pmatrix} -\sin 2\pi y \\ \sin 4\pi x \end{pmatrix}, $$

$$ \gamma = 5/3. $$

We show the density map at $t = 0.5$ in Fig. 7. We observe that the shocks and discontinuities are well captured without spurious numerical artifacts. We also notice the usual “eye-shape” high frequency feature at the center of the domain, demonstrating the accuracy of our solver. Note that this test does not show any low $\beta$ zone. Thus, the solver is fully conservative with respect to $\mathbf{B}$ as the entropic correction is never activated.

### 8.2.2. Rotated shock tube

The rotated shock tube problem has been proposed in [7]. It consists of a 1D shock tube rotated by an angle $\theta$ in order to obtain a 2D shock propagation that is not aligned with the grid. The test takes place in the $[0 : 1]^2$ square with Neumann boundary conditions. The setup is given by:

$$ \theta = \arctan(-2), $$

$$ \mathbf{R}(\theta) = \begin{pmatrix} \sin \theta & \cos \theta \\ \cos \theta & -\sin \theta \end{pmatrix}, $$

$$ \mathbf{u}_0 = \begin{pmatrix} 0 \\ 10 \end{pmatrix}, $$

$$ \mathbf{B}_0 = \frac{5}{\sqrt{4\pi}} \begin{pmatrix} 1 \\ 1 \end{pmatrix}, $$
<figure><img src="image_10.png" /><figcaption>Fig. 8. Parallel component of the magnetic field along the rotated shock propagation at t=0.03.</figcaption></figure>

$$ (x_\theta, y_\theta) = (\tan \theta(x - 0.5), y - 0.5), $$

$$ \rho(x, y) = 1, $$

$$ \mathbf{B}(x, y) = \mathbf{R}(\theta)\mathbf{B}_0, $$

$$ \mathbf{u}(x, y) = \begin{cases} \mathbf{R}(\theta)\mathbf{u}_0 & \text{for } x_\theta < y_\theta, \\ -\mathbf{R}(\theta)\mathbf{u}_0 & \text{elsewhere.} \end{cases} $$

$$ p(x, y) = \begin{cases} 20 & \text{for } x_\theta < y_\theta, \\ 1 & \text{elsewhere.} \end{cases}. $$

Note that the magnetic field is initialized as a constant on the whole domain, hence the condition $\nabla \cdot \mathbf{B} = 0$ is verified at the beginning of the computation. Our solver is able to robustly and accurately simulate this rotated shock propagation. A quantity of interest in this problem is the component of the magnetic field that is parallel to the shock propagation. Without discretization error, this quantity should remain constant similarly to $B_x$ in a purely 1D setup. In Fig. 8, we show the component of the magnetic field that is parallel to the shock propagation, with both 3+1 and 5+1 solvers. Both schemes produces discretization errors at the location of discontinuities, the errors with the 5+1 waves solver are larger than the errors with the 3+1 waves solver. These errors can be compared with [7] for constrained transport schemes and we point out that the 3+1 and 5+1 waves solvers produce less oscillations around the discontinuities.

### 8.2.3. MHD blast - standard configuration

The Blast test case was introduced in [35]. The setup takes place in the periodic $[0 : 1]^2$ square. A circular region of radius $r_c = 0.1$ is initialized with a greater pression than the rest of the domain. As the computation starts, the blast expands outwards in an elliptical shape due to the presence of a magnetic field. We recall the exact setup:

$$ p(x, y) = \begin{cases} 10 & \text{for } r < r_c, \\ 0.1 & \text{for } r \ge r_c, \end{cases} $$

$$ \mathbf{B}(x, y) = \begin{pmatrix} \sqrt{2\pi} \\ \sqrt{2\pi} \end{pmatrix}, $$

$$ \gamma = 5/3, $$

$$ \rho(x, y) = 1, $$

$$ \mathbf{u}(x, y) = 0. $$

Our numerical method is able to simulate the expansion of this blast wave accurately and is stable as demonstrated in Fig. 9 where we show the density map at $t = 0.2$. We can see that the expanding wave is well captured. Note that this test does not show any low $\beta$ zone. Thus, the solver is fully conservative with respect to $\mathbf{B}$.

### 8.2.4. MHD blast - low $\beta$ configuration

This test case is inspired from [36]. It consists of the same setup as section 8.2.3 with a lower $\beta \approx 10^{-6}$:

$$ p(x, y) = \begin{cases} 1000 & \text{for } r < r_c, \\ 0.1 & \text{for } r \ge r_c, \end{cases} $$
<figure><img src="image_11.png" /><figcaption>Fig. 9. Density map of the MHD Blast at t = 0.2 s.</figcaption></figure>

$$ \mathbf{B}(x, y) = \begin{pmatrix} 250/\sqrt{2} \\ 250/\sqrt{2} \end{pmatrix}, \\ \gamma = 1.4, $$

$$ \rho(x, y) = 1, $$

$$ \mathbf{u}(x, y) = 0. $$

The dynamic of the low $\beta$ blast wave is the same as in 8.2.3 but is harder to tackle as the simulation reaches the limit of the admissibility domain ($e \approx 0$) and develops strong $\mathbf{B}$ gradients. Note that the 5 + 1 wave solver and the constrained transport method [37] fail to produce an admissible result as the computation presents negative internal energies (directly after few iterations). We point out that the 5 + 1 solver seems, however, more robust than the constrained transport method on such problems: for lower values of the magnetic field $25/\sqrt{2}$, the relaxation solver is stable while the constrained transport method fails after few iterations. It is possible to still get an admissible result by artificially forcing the internal energy to stay above a small threshold (hence losing energy conservation), a solution used here with the constrained transport method, or by using an entropic correction on the induction equation (hence losing the magnetic field conservation), a solution used here with the 5 + 1 waves relaxation solver. In Figs. 10, we show the density map of this test case at $t = 0.02$ with our method and the energy-fixed constrained transport solver from the Heracles code [38]. Both methods are able to capture the low $\beta$ Blast propagation, however, we point out that the 5 + 1 waves solver is less diffusing as it reaches higher values for the magnetic field up (+18%).

### 8.2.5. MHD rotor

The MHD Rotor test case was first introduced in [6]. The setup consists of launching a rapidly spinning cylinder in a light ambient fluid. This rotation sends strong torsional Alfvén waves in the surrounding fluid. We initialize the solution in the $[0 : 1]^2$ periodic square as follows:

$$ \rho(x, y) = 1.0, $$

$$ \rho(x, y) = \begin{cases} 10 & \text{for } r < r_0, \\ 1 + 9f & \text{for } r \ge r_1 \text{ & } r \le r_0, \\ 1 & \text{elsewhere} \end{cases} $$

$$ \mathbf{u}(x, y) = \begin{cases} \frac{u_0}{r_0} (0.5 - y, x - 0.5) & \text{for } r < r_0, \\ \frac{f u_0}{r_0} (0.5 - y, x - 0.5) & \text{for } r \ge r_1 \text{ & } r \le r_0, \\ (0, 0) & \text{elsewhere} \end{cases} $$

$$ \mathbf{B}(x, y) = \begin{pmatrix} 5/\sqrt{4\pi} \\ 0 \end{pmatrix}, \\ \gamma = 1.4, $$

$$ (r_0, r_1) = (0.1, 0.115), $$
<figure><img src="image_12.png" /><figcaption>Fig. 10. Density map of the low β MHD blast at t = 0.2 s with our solver and the Heracles code's constrained transport method [38], [37].</figcaption></figure>

$$f = (r_1 - r)/(r_1 - r_0),$$

$$u_0 = 2.$$

We show the result of our simulation in Fig. 11. We observe that the central shear ring as well as the torsional waves is well captured by our solver. Note that this simulation does not require the use of entropic correction. Thus, the solver is fully conservative with respect to **B**.

### 8.2.6. Field loop advection

This test was introduced in [39] and involves advecting a field loop (a cylindrical current distribution) diagonally across the grid. One can choose any arbitrary angle. For the 2D results presented here, the problem domain is defined as $-1 < x < 1$ and $-0.5 < y < 0.5$. The flow has an inclination with $V_x = 2$ and $V_y = 1$. Both the density and pressure are set to 1.0, with the gas constant given by $\gamma = 5/3$. Periodic boundary conditions are applied across the domain. The magnetic field is initialized using an arbitrary vector potential. We set $A_z = \max([A_0(r_0 - r)], 0)$. This results in $(B_x, B_y)(r) = \frac{A_0}{r}(-x, y)$ if $r < r_0$, and $(0, 0)$ otherwise. We chose $A_0 = 0.001$ and set the radius for the loop as $r_0 = 0.3$. After a duration of $t = 2.0s$, the field loop is expected to have been advected and returned to its initial state. The quality of the solution can be assessed by comparing it to the initial solution shown in Fig. 12. The magnetic intensity, defined as $I = \sqrt{B_x^2 + B_y^2}$, obtained with our 5 + 1 waves solver, is illustrated in Fig. 13. One can observe that the entropic correction helps with preserving the shape of the cylinder and suppresses the spurious patterns observed with the conservative method. The source terms are activated here as the Alfvén number is above $Al_{max} = 10$ in this test.
<figure><img src="image_13.png" /><figcaption>Fig. 11. Density map of the MHD Rotor at $t = 0.15$ s.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Magnetic intensity of the field loop advection at time $t = 0$.</figcaption></figure>

## 9. Conclusion and discussion

In this paper, we have developed a new multi-dimensional, robust, and cell-centered finite-volume solver for ideal MHD. The solver is based on splitting and relaxation techniques, and can easily be extended to higher orders because of its reduced stencil. A symmetric version of the solver has been developed by introducing an entropic correction on the induction equation, in order to obtain an entropy-satisfying (but non-conservative for the magnetic field) scheme robust in low plasma beta regions and accurate in high Alfvén number regions. Another solution could be to use a floor value for the internal energy as classically done with constrained transport or divergence cleaning schemes that are not entropy satisfying. We, however, point out that the fully conservative relaxation solver is observed to be more robust than constrained transport schemes on low plasma beta test cases.

This cell-centered scheme could be coupled to a divergence cleaning or constrained transport method. We, however, highlight that all the tests we have performed do not seem to require a specific treatment of the divergence of the magnetic field, and a divergence consistent with zero with errors proportional to $\Delta x$ and $\Delta t$ at the power of the order of the spatial and temporal reconstructions seem sufficient. It is a common belief that the stability of MHD numerical schemes is closely tied to errors in magnetic field divergence. However, our research, as presented in this paper, suggests that this may not always be the case. To illustrate, we have successfully designed an entropy-satisfying MHD solver using the symmetric form of MHD equations without specifically addressing divergence issues. Furthermore, we have found that constrained transport schemes, while maintaining zero divergence at machine precision, do not necessarily satisfy entropy conditions and can fail to maintain positive internal energy in areas of low plasma beta.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Magnetic intensity of the field loop advection at time $t = 2.0$. Top: With the entropic correction. Bottom: Without the entropic correction.</figcaption></figure>

Additionally, there is a prevalent view that errors in magnetic divergence significantly impact the physical accuracy of simulations, potentially leading to artificial magnetic monopoles. We offer several arguments to challenge this perspective. Even in constrained transport schemes, certain terms involving divergence in the conservative forms of the Lorentz force and the energy evolution equation do not achieve zero at machine precision, despite a zero divergence. These residual terms in the entropy evolution equation are indeed the reason why constrained transport schemes are not entropy satisfying. Moreover, it can be demonstrated that constrained transport schemes are not immune to divergence errors. For example, the rotated shock tube case detailed in [7] shows that at the continuous level, a zero divergence equates to a constant magnetic field parallel to the shock tube. However, constrained transport schemes do not maintain this constant magnetic field at machine precision, thus resulting in “divergence errors” that are significant for the physics at play.

In conclusion, while ensuring zero magnetic divergence at machine precision in simulations is physically relevant, this is only feasible when aligning the grid to a specific magnetic field configuration. This issue is akin to preserving angular momentum in a rotating structure, achievable at machine precision only in a polar grid. Consequently, for simulations with highly dynamic magnetic fields, maintaining zero divergence at machine precision on a Cartesian grid may not be as critical with a solver that is entropy-satisfying.

The MHD relaxation solver presented in this paper is a direct extension of 1D relaxation solvers already used for the Euler equations and can be implemented in a one-step flux-update algorithm, that can easily be extended to higher orders and to non-ideal MHD. Because of its simplicity, this solver should also have improved performances compared to other multi-dimensional MHD solvers (constrained transport and divergence cleaning) and offers interesting possibilities for large-scale physical applications on the next generation of exascale supercomputers.

CRediT authorship contribution statement

**Pascal Tremblin:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Methodology, Conceptualization.
**Rémi Bourgeois:** Writing – original draft, Visualization, Validation, Software, Methodology.
**Solène Bulteau:**
Methodology, Investigation. **Samuel Kokh**: Writing – review & editing, Methodology, Conceptualization. **Thomas Padioleau**: Software. **Maxime Delorme**: Investigation. **Antoine Strugarek**: Writing – review & editing, Investigation. **Matthias González**: Writing – review & editing. **Allan Sacha Brun**: Writing – review & editing.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Pascal Tremblin reports financial support was provided by European Research Council. If there are other authors, they 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

No data was used for the research described in the article.

## Acknowledgements

P. Tremblin and S. Bulteau would like to acknowledge and thank the ERC for funding this work under the Horizon 2020 program project ATMO (ID: 757858).

## Appendix A. Deriving the conservative MHD equations

### A.1. Useful vector identities

#### A.1.1. Lorentz's force in conservative form

The first identity we derive is

$$j \times \mathbf{B} = -(\nabla \cdot \mathbf{B})\mathbf{B} - \nabla \cdot \left(\frac{\mathbf{B}^2}{2}\mathbf{I} - \mathbf{B} \otimes \mathbf{B}\right). \qquad (A.1)$$

We only verify this equality for the x component as the relationship for the two other components are checked by rotational invariance. We have $j \times \mathbf{B} = (\nabla \times \mathbf{B}) \times \mathbf{B}$. Expanding the first component, we get $[(\nabla \times \mathbf{B}) \times \mathbf{B}]_x = B_z(\partial_z B_x - \partial_x B_z) - B_y(\partial_x B_y - \partial_y B_x)$. Moreover, $[\nabla \cdot (\frac{\mathbf{B}^2}{2}\mathbf{I})]_x = B_x \partial_x B_x + B_y \partial_y B_y + B_z \partial_z B_z$. Lastly, $[\nabla \cdot (\mathbf{B} \otimes \mathbf{B})]_x = (\nabla \cdot \mathbf{B})B_x + B_x \partial_x B_x + B_y \partial_y B_x + B_z \partial_z B_x$. Collecting the right hand side terms, we get $-(\nabla \cdot \mathbf{B})B_x - B_x \partial_x B_x - B_y \partial_y B_y - B_z \partial_z B_z + (\nabla \cdot \mathbf{B})B_x + B_x \partial_x B_x + B_y \partial_y B_x + B_z \partial_z B_x$ where both the terms proportional to the divergence of $\mathbf{B}$ and $B_x \partial_x B_x$ cancel out and provide the desired result.

#### A.1.2. Fully developed Lorentz force

Using $\nabla \cdot (\mathbf{B} \otimes \mathbf{B}) = \mathbf{B}(\nabla \cdot \mathbf{B}) + (\mathbf{B} \cdot \nabla)\mathbf{B}$, we get:

$$j \times \mathbf{B} = (\mathbf{B} \cdot \nabla)\mathbf{B} - \nabla \left(\frac{\mathbf{B}^2}{2}\right) \qquad (A.2)$$

#### A.1.3. Curl of a cross product

$$\nabla \times (\mathbf{u} \times \mathbf{B}) = \nabla \cdot (\mathbf{B} \otimes \mathbf{u} - \mathbf{u} \otimes \mathbf{B}). \qquad (A.3)$$

$$\nabla \times (\mathbf{u} \times \mathbf{B}) = \mathbf{u}(\nabla \cdot \mathbf{B}) - \mathbf{B}(\nabla \cdot \mathbf{u}) + (\mathbf{B} \cdot \nabla)\mathbf{u} - (\mathbf{u} \cdot \nabla)\mathbf{B} \qquad (A.4)$$

#### A.1.4. Transport of a squared quantity

$$((\mathbf{u} \cdot \nabla)\mathbf{A}) \cdot \mathbf{A} = (\mathbf{u} \cdot \nabla)\frac{\mathbf{A}^2}{2} = \nabla \left(\frac{\mathbf{A}^2}{2}\right) \cdot \mathbf{u} \qquad (A.5)$$

### A.2. Full system

Our goal is now to go from the non conservative MHD system:

$$\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0,$$

$$\partial_t(\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + j \times \mathbf{B},$$

$$\partial_t(\rho e) + \nabla \cdot (\rho e \mathbf{u}) = -p \nabla \cdot \mathbf{u},$$

$$\partial_t \mathbf{B} - \nabla \times (\mathbf{u} \times \mathbf{B}) = 0, \qquad (A.6)$$
to the conservative MHD system.

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + \sigma \mathbf{B} \otimes \mathbf{B}) = 0, $$

$$ \partial_t (\rho E) + \nabla \cdot (\rho E \mathbf{u} + \sigma \mathbf{u} - (\mathbf{B} \cdot \mathbf{u}) \mathbf{B}) = 0, $$

$$ \partial_t \mathbf{B} + \nabla \cdot (\mathbf{u} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{u}) = 0. \qquad (A.7) $$

Where $e_{mag} = \frac{B^2}{2\rho}$ and $\sigma = p + \frac{B^2}{2}$. Obtaining the conservative momentum equation is straightforward using (A.1), substituting for $j \times \mathbf{B}$ and assuming $\nabla \cdot \mathbf{B} = 0$. Obtaining the conservative induction equation is also straightforward using (A.3) (note that using the $\nabla \cdot \mathbf{B} = 0$ hypothesis is not necessary to obtain the induction equation). This leaves us with deriving the total energy equation.

### A.2.1. Kinetic energy evolution equation

From the non conservative momentum equation, we can deduce the evolution equation of the velocity $\partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{u} + \nabla p / \rho = j \times \mathbf{B} / \rho$. Dotting this equation against $\mathbf{u}$, we get $\partial_t (\frac{u^2}{2}) + ((\mathbf{u} \cdot \nabla)\mathbf{u}) \cdot \mathbf{u} + \nabla p \cdot \mathbf{u} / \rho = (j \times \mathbf{B}) \cdot \mathbf{u} / \rho$. Using (A.5), we have that $((\mathbf{u} \cdot \nabla)\mathbf{u}) \cdot \mathbf{u} = (\mathbf{u} \cdot \nabla) (\frac{u^2}{2})$. Substituting this transport term and multiplying by $\rho$, we get $\rho \partial_t (\frac{u^2}{2}) + \rho(\mathbf{u} \cdot \nabla) (\frac{u^2}{2}) + \nabla p \cdot \mathbf{u} = (j \times \mathbf{B}) \cdot \mathbf{u}$. Adding $\frac{u^2}{2} (\partial_t \rho + \nabla \cdot (\rho \mathbf{u})) = 0$, we get: $\partial_t (\rho \frac{u^2}{2}) + \rho(\mathbf{u} \cdot \nabla) (\frac{u^2}{2}) + \frac{u^2}{2} \nabla \cdot (\rho \mathbf{u}) + \nabla p \cdot \mathbf{u} = (j \times \mathbf{B}) \cdot \mathbf{u}$. Since $\rho(\mathbf{u} \cdot \nabla) (\frac{u^2}{2}) + \frac{u^2}{2} \nabla \cdot (\rho \mathbf{u}) = \nabla \cdot (\frac{\rho u^2 u}{2})$, noting $e_{kin} = \frac{u^2}{2}$, we get:

$$ \partial_t (\rho e_{kin}) + \nabla \cdot (\rho e_{kin} \mathbf{u}) + \nabla p \cdot \mathbf{u} = (j \times \mathbf{B}) \cdot \mathbf{u}. \qquad (A.8) $$

Summing this with the internal energy evolution equation, we get: $\partial_t (\rho(e + e_{kin})) + \nabla \cdot (\rho(e + e_{kin})\mathbf{u} + p\mathbf{u}) = (j \times \mathbf{B}) \cdot \mathbf{u}$. Replacing the right hand side using (A.2), we get:

$$ \partial_t (\rho(e + e_{kin})) + \nabla \cdot (\rho(e + e_{kin})\mathbf{u} + p\mathbf{u}) = ((\mathbf{B} \cdot \nabla)\mathbf{B}) \cdot \mathbf{u} - \nabla \left(\frac{\mathbf{B}^2}{2}\right) \cdot \mathbf{u}. \qquad (A.9) $$

### A.2.2. Magnetic energy evolution equation

Using the identity (A.4), we get $\partial_t \mathbf{B} - \mathbf{u}(\nabla \cdot \mathbf{B}) + \mathbf{B}(\nabla \cdot \mathbf{u}) - (\mathbf{B} \cdot \nabla)\mathbf{u} + (\mathbf{u} \cdot \nabla)\mathbf{B} = 0$. Dotting against $\mathbf{B}$, we get

$$ \partial_t (\rho e_{mag}) - (\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}) + (\nabla \cdot \mathbf{u}) \mathbf{B}^2 - ((\mathbf{B} \cdot \nabla)\mathbf{u}) \cdot \mathbf{B} + ((\mathbf{u} \cdot \nabla)\mathbf{B}) \cdot \mathbf{B} = 0. \qquad (A.10) $$

### A.2.3. Total energy evolution equation

Summing (A.9) and (A.10), we get $\partial_t (\rho E) + \nabla \cdot (\rho(e + e_{kin})\mathbf{u} + p\mathbf{u}) = ((\mathbf{B} \cdot \nabla)\mathbf{B}) \cdot \mathbf{u} - \nabla \left(\frac{\mathbf{B}^2}{2}\right) \cdot \mathbf{u} + (\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}) - (\nabla \cdot \mathbf{u}) \mathbf{B}^2 + ((\mathbf{B} \cdot \nabla)\mathbf{u}) \cdot \mathbf{B} - ((\mathbf{u} \cdot \nabla)\mathbf{B}) \cdot \mathbf{B}$. Using (A.5), we have $-((\mathbf{u} \cdot \nabla)\mathbf{B}) \cdot \mathbf{B} - \nabla \left(\frac{\mathbf{B}^2}{2}\right) \cdot \mathbf{u} = \nabla(\mathbf{B}^2) \cdot \mathbf{u}$. Moreover, since $\nabla(\mathbf{B}^2) \cdot \mathbf{u} + (\nabla \cdot \mathbf{u}) \mathbf{B}^2 = \nabla \cdot (\mathbf{B}^2 \mathbf{u}) = \nabla \cdot (\rho e_{mag} \mathbf{u} + \mathbf{B}^2 / 2 \mathbf{u})$, we can show that:

$$ \partial_t (\rho E) + \nabla \cdot (\rho E \mathbf{u} + \sigma \mathbf{u}) = ((\mathbf{B} \cdot \nabla)\mathbf{B}) \cdot \mathbf{u} + (\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}) + ((\mathbf{B} \cdot \nabla)\mathbf{u}) \cdot \mathbf{B}. \qquad (A.11) $$

As $((\mathbf{B} \cdot \nabla)\mathbf{B}) \cdot \mathbf{u} + ((\mathbf{B} \cdot \nabla)\mathbf{u}) \cdot \mathbf{B} = (\mathbf{B} \cdot \nabla)(\mathbf{u} \cdot \mathbf{B}) = \nabla(\mathbf{B} \cdot \mathbf{u}) \cdot \mathbf{B}$ and $\nabla(\mathbf{B} \cdot \mathbf{u}) \cdot \mathbf{B} + (\nabla \cdot \mathbf{B})(\mathbf{u} \cdot \mathbf{B}) = \nabla \cdot ((\mathbf{B} \cdot \nabla)\mathbf{B})$, we get the desired result. Note that it is not required to assume $\nabla \cdot \mathbf{B} = 0$ to obtain the conservative total energy equation.

## A.3. Entropy inequality

### A.3.1. Entropy inequality of the non conservative MHD system

We start with the classical result of the entropy inequality of the MHD system (A.6), starting from the evolution equation of the internal energy. We note $D_t = \partial_t + u \nabla$. We have $D_t e = -p(\nabla \cdot \mathbf{u})\tau$ where $\tau = 1/\rho$. From the density evolution equation, we have that $D_t \tau = \tau(\nabla \cdot \mathbf{u})$. Therefore, $D_t e + p D_t \tau = 0$. Using the first principle of thermodynamics $de + pd\tau = Tds$, we get

$$ D_t s = 0. \qquad (A.12) $$

### A.3.2. Entropy inequality of the conservative MHD system

To go from the non conservative system to the conservative system, we only had to cancel one term in the momentum equation, using the $\nabla \cdot \mathbf{B} = 0$ hypothesis. This means that if we are discretizing the conservative momentum equation and that the numerical value of the divergence is not zero, we are in fact discretizing $\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + j \times \mathbf{B} + (\nabla \cdot \mathbf{B})\mathbf{B}$. We want to derive the corresponding internal energy equation. We dot the momentum equation against $\mathbf{u}$ and subtract it to the conservative total energy equation. Doing this, we get $\partial_t (\rho e) + \nabla \cdot (\rho e \mathbf{u}) = -p \nabla \cdot \mathbf{u} - (\nabla \cdot \mathbf{B})(\mathbf{B} \cdot \mathbf{u})$. Performing the same steps as above, we get $D_t e + p D_t \tau = -\tau(\nabla \cdot \mathbf{B})(\mathbf{B} \cdot \mathbf{u})$ thus:

$$ D_t s = -\frac{\tau}{T} (\nabla \cdot \mathbf{B})(\mathbf{B} \cdot \mathbf{u}) \qquad (A.13) $$
References

[1] J.U. Brackbill, D.C. Barnes, The effect of nonzero V · B on the numerical solution of the magnetohydrodynamic equations, J. Comput. Phys. 35 (3) (1980) 426-430, https://doi.org/10.1016/0021-9991(80)90079-0.

[2] D. Ryu, F. Miniati, T.W. Jones, A. Frank, A divergence-free upwind code for multidimensional magnetohydrodynamic flows, Astrophys. J. 509 (1) (1998) 244–255, https://doi.org/10.1086/306481, arXiv:astro-ph/9807228.

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

[4] A. Dedner, F. Kemm, D. Kröner, C. Munz, T. Schnitzer, M. Wessenberg, Hyperbolic divergence cleaning for the MHD equations, J. Comput. Phys. 175 (2002) 645–673.

[5] C.R. Evans, J.F. Hawley, Simulation of magnetohydrodynamic flows: a constrained transport model, Astrophys. J. 332 (1988) 659, https://doi.org/10.1086/166684.

[6] 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 (2) (1999) 270–292, https://doi.org/10.1006/jcp.1998.6153.

[7] 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.

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

[9] G. Gallice, Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates, Numer. Math. 94 (4) (2003) 673–713.

[10] F. Bouchut, C. Klingenberg, K. Waagan, A multiwave approximate Riemann solver for ideal MHD based on relaxation: I: theoretical framework, Numer. Math. 108 (1) (2007) 7–42.

[11] F. Bouchut, C. Klingenberg, K. Waagan, A multiwave approximate Riemann solver for ideal MHD based on relaxation II: numerical implementation with 3 and 5 waves, Numer. Math. 115 (4) (2010) 647–679.

[12] S. Godunov, Symmetric form of the magnetohydrodynamic equation, Tech. rep., Computer Center, Novosibirsk, USSR, 1972.

[13] S. Busto, M. Dumbser, A new thermodynamically compatible finite volume scheme for magnetohydrodynamics, SIAM J. Numer. Anal. 61 (1) (2023) 343–364, https://doi.org/10.1137/22M147815X.

[14] C. Chalons, M. Girardin, S. Kokh, An all-regime Lagrange-projection-like scheme for the gas dynamics equations on unstructured meshes, Commun. Comput. Phys. 20 (1) (2016) 188–233, https://doi.org/10.4208/cicp.260614.061115a.

[15] T. Padioleau, P. Tremblin, E. Audit, P. Kestener, S. Kokh, A high-performance and portable all-Mach regime flow solver code with well-balanced gravity. Application to compressible convection, Astrophys. J. 875 (2) (2019) 128, https://doi.org/10.3847/1538-4357/ab0f2c.

[16] R. Bourgeois, P. Tremblin, S. Kokh, T. Padioleau, Recasting an operator splitting solver into a standard finite volume flux-based algorithm. The case of a Lagrange-projection-type method for gas dynamics, J. Comput. Phys. 496 (2024) 112594, https://doi.org/10.1016/j.jcp.2023.112594.

[17] C. Hirt, A. Amsden, J. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (3) (1974) 227–253, https://doi.org/
10.1016/0021-9991(74)90051-5.

[18] E. Godlewski, P. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Applied Mathematical Sciences, vol. 118, Springer, 1996, https://books.google.fr/books?id = 9BwMIDMmTmcC.

[19] B. Després, Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques, Mathématiques et Applications, Springer Berlin Heidelberg, 2010, https://
books.google.cn/books?id = yVbnQbtq1JcC.

[20] S. Jin, Z. Xin, The relaxation schemes for systems of conservation laws in arbitrary space dimensions, Commun. Pure Appl. Math. 48 (1995) 235–276, https://
apisemanticscholar.org/CorpusID:13245844.

[21] I. Suliciu, On the thermodynamics of rate-type fluids and phase transitions. I. Rate-type fluids, Int. J. Eng. Sci. 36 (1998) 921–947, https://apisemanticscholar.
org/CorpusID:121171038.

[22] F. Coquel, E. Godlewski, B. Perthame, A. In, P. Rascle, Some new Godunov and relaxation methods for two-phase flow problems, https://apisemanticscholar.
org/CorpusID:115535063, 2001.

[23] F. Bouchut, Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws: And Well-Balanced Schemes for Sources, Springer Science & Business Media, 2004.

[24] C. Chalons, J.-F. Coulombel, Relaxation approximation of the Euler equations, J. Math. Anal. Appl. 348 (2) (2008) 872–893, https://doi.org/10.1016/j.jmaa.
2008.07.034.

[25] F. Coquel, Q.L. Nguyen, M. Postel, Q.H. Tran, Entropy-satisfying relaxation method with large time-steps for Euler ibvps, Math. Comput. 79 (2010) 1493–1533.

[26] B. Després, A new Lagrangian formulation of ideal magnetohydrodynamics, J. Hyperbolic Differ. Equ. 08 (01) (2011) 21–35, https://doi.org/10.1142/
S0219891611002329.

[27] F. Bouchut, Entropy satisfying flux vector splittings and kinetic BGK models, Numer. Math. 94 (4) (2002) 623–672.

[28] C. Klingenberg, K. Waagan, Relaxation solvers for ideal mhd equations - a review, Acta Math. Sci. 30 (2010) 621–632.

[29] B. van Leer, On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe, SIAM J. Sci. Stat. Comput. 5 (1984) 1–20, https://
doi.org/10.1137/0905001.

[30] W. Dai, P. Woodward, An approximate Riemann solver for ideal magnetohydrodynamics, J. Comput. Phys. 111 (2) (1994) 354–372, https://doi.org/10.1006/
jcp.1994.1069.

[31] M. Brio, C. Wu, An upwind differencing scheme for the equations of ideal magnetohydrodynamics, J. Comput. Phys. 75 (2) (1988) 400–422, https://doi.org/10.
1016/0021-9991(88)90120-9.

[32] 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.

[33] T. Miyoshi, K. Kusano, A multi-state hll approximate Riemann solver for ideal magnetohydrodynamics, J. Comput. Phys. 208 (2005) 315–344, https://doi.org/
10.1016/j.jcp.2005.02.017.

[34] S.A. Orszag, C.-M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, J. Fluid Mech. 90 (1) (1979) 129–143, https://doi.org/10.
1017/S002211207900210X.

[35] J.M. Stone, T. Gardiner, A simple unsplit Godunov method for multidimensional mhd, New Astron. 14 (2) (2009) 139–148, https://doi.org/10.1016/j.newast.
2008.06.003.

[36] D.S. Balsara, Self-adjusting, positivity preserving high order schemes for hydrodynamics and magnetohydrodynamics, J. Comput. Phys. 231 (22) (2012) 7504–7517, https://doi.org/10.1016/j.jcp.2012.01.032.

[37] J. Vides, E. Audit, H. Guillard, B. Nkonga, Divergence-free mhd simulations with the heracles code, ESAIM Proc. 43 (2013) 180–194, https://doi.org/10.1051/
proc/201343012.

[38] M. González, E. Audit, P. Huynh, HERACLES: a three-dimensional radiation hydrodynamics code, Astron. Astrophys. 464 (2) (2007) 429–435, https://doi.org/
10.1051/0004-6361:20065486.

[39] G. Tóth, D. Odstrčil, Comparison of some flux corrected transport and total variation diminishing numerical schemes for hydrodynamic and magnetohydrodynamic problems, J. Comput. Phys. 128 (1) (1996) 82–100.