# A structure-preserving semi-implicit IMEX finite volume scheme for ideal magnetohydrodynamics at all Mach and Alfvén numbers

Walter Boscheri, Andrea Thomann

► To cite this version:

Walter Boscheri, Andrea Thomann. A structure-preserving semi-implicit IMEX finite volume scheme for ideal magnetohydrodynamics at all Mach and Alfvén numbers. Journal of Scientific Computing, 2024, 100 (3), pp.67. 10.1007/s10915-024-02606-1. hal-04495185v2

HAL Id: hal-04495185

https://hal.science/hal-04495185v2

Submitted on 28 Jun 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
# A structure-preserving semi-implicit IMEX finite volume scheme for ideal magnetohydrodynamics at all Mach and Alfvén numbers

Walter Boscheri · Andrea Thomann

Received: date / Accepted: date

**Abstract** We present a divergence-free semi-implicit finite volume scheme for the simulation of the ideal magnetohydrodynamics (MHD) equations which is stable for large time steps controlled by the local transport speed at all Mach and Alfvén numbers. An operator splitting technique allows to treat the convective terms explicitly while the hydrodynamic pressure and the magnetic field contributions are integrated implicitly, yielding two decoupled linear implicit systems. The linearity of the implicit part is achieved by means of a semi-implicit time linearization. This structure is favorable as second-order accuracy in time can be achieved relying on the class of semi-implicit IMplicit-EXplicit Runge-Kutta (IMEX-RK) methods. In space, implicit cell-centered finite difference operators are designed to discretely preserve the divergence-free property of the magnetic field on three-dimensional Cartesian meshes. The new scheme is also particularly well suited for low Mach number flows and for the incompressible limit of the MHD equations, since no explicit numerical dissipation is added to the implicit contribution and the time step is scale independent. Likewise, highly magnetized flows can benefit from the implicit treatment of the magnetic fluxes, hence improving the computational efficiency of the novel method. The convective terms undergo a shock-capturing second order finite volume discretization to guarantee the effectiveness of the proposed method even for high Mach number flows. The new scheme is benchmarked against a series of test cases for the ideal MHD equations addressing different acoustic and Alfvén Mach number regimes where the performance and the stability of the new scheme is assessed.

**Keywords** semi-implicit · divergence-free · low Mach and low Alfvén number flows · magnetohydrodynamics · finite volume scheme

**Mathematics Subject Classification (2020)** MSC 65 · MSC 68

## 1 Introduction

Compressible magnetized plasma flows are mathematically modeled by the equations of magnetohydrodynamics (MHD), that constitute a first order hyperbolic system of nonlinear partial differential equations (PDE). Potential practical applications concern astrophysics to describe plasma flows in the magnetosphere of neutron stars, or inertial and magnetic confinement fusion, in order to calibrate or design fusion reactors for civil energy production.

W. Boscheri

Laboratoire de Mathématiques UMR 5127 CNRS, Université Savoie Mont Blanc, 73376 Le Bourget du Lac, France
E-mail: walter.boscheri@univ-smb.fr

A. Thomann

Université de Strasbourg, CNRS, Inria, IRMA, F-67000 Strasbourg, France
E-mail: andrea.thomann@inria.fr
The dynamics of plasma flows can involve different time scales, that correspond to the characteristics of the governing PDE system. The multi-scale nature of the MHD equations can be analyzed by considering the Mach number $M_c$, which is the ratio between the fluid velocity and the sound speed, and the Alfvén number $M_b$, that represents the ratio between the fluid velocity and the speed of the Alfvén wave related to the magnetic field. Different values of $M_c$ and $M_b$ yield different fluid regimes, and, consequently, the nature of the governing equations can change passing for instance from hyperbolic to elliptic. A classical example is given by low Mach number flows, which are dominated by acoustic waves, hence recovering an incompressible fluid regime in the stiff relaxation limit for $M_c \to 0$. On the other hand, the material speed is of the same order than the acoustic wave speeds in high Mach number flows, thus leading to compressible fluids with shock waves. Another source of stiffness arises for highly magnetized flows, i.e. for $M_b \to 0$, that are typically encountered in Tokamak scenarios for magnetic confinement fusion applications.

The numerical resolution of the MHD equations is therefore very challenging if the numerical method is required to deal with the entire spectrum of the different time scales. Explicit Godunov-type finite volume schemes [40, 42, 4] exhibit excellent shock-capturing properties that can perfectly handle high Mach number flows. However, if different regimes coexist in the flow, two main problems are generated by the different time scales that can be assumed by the model: i) the computational time step obeys a stability condition that is constrained by the magnitude of the fastest characteristic speeds, hence leading to very small time steps in the stiff relaxation limits of the model which make the simulation extremely slow, and, ultimately, even impossible to be accomplished; ii) the numerical viscosity does not have the correct scaling but is of order $O(1/M_c)$ which is especially severe for low Mach number flows and therefore the accuracy of the solution is compromised [29, 41].

To tackle these challenges, a first strategy consists in rescaling the dissipation term in the numerical flux so that the numerical viscosity becomes independent of the Mach number [22, 8, 9, 11]. A fully implicit time integration has then to be adopted to circumvent the restrictive time step condition [47, 61], which leads to a highly nonlinear system to be solved with a large number of unknowns, or local time stepping techniques can be foreseen [1, 19]. Another possibility is instead given by the implicit treatment of only one part of the system to be solved while keeping the remaining part explicit [26, 32, 20, 25, 31, 37, 10, 49]. The terms responsible to a possible stiffness of the governing model undergo an implicit discretization in which no numerical diffusion is needed, since stability follows from implicitness, so that the above mentioned problems of explicit solvers can be simultaneously solved. This technique makes often use of the class of IMplicit-EXplicit (IMEX) time integrators forwarded e.g. in [2, 13, 16, 53, 14, 52], that have been developed with the specific purpose of dealing with multiple time scales, typically a slow and a fast time scale. Furthermore, these IMEX schemes are proven to be Asymptotic Preserving (AP) and, some even Asymptotically Accurate (AA) meaning that they recover a consistent discretization of the limit model also with high order accuracy.

Besides the multi-scale character, the MHD equations embed another constraint on the divergence of the magnetic field, which must remain zero at all times because no magnetic monopoles exist in nature. This is an involution that has to be carefully taken into account by the numerical scheme, since in principle it is not satisfied at the discrete level. In the literature, there are three main approaches to satisfy the solenoidal property of magnetic field in a numerical scheme: i) the use of staggered grids is typically adopted to cancel the divergence errors, which requires the knowledge of the electric field [62, 7]; ii) divergence cleaning methods [51, 28] add one hyperbolic equation to the MHD system in order to transport away the divergence errors, which require the determination of the cleaning speed as a new parameter; iii) constrained transport schemes on collocated grids that also make use of the electric field [38, 43, 44].

The aim of this paper is to design a numerical scheme for the ideal MHD equations that can handle fast time scales which might be induced by both the stiff limits $M_c \to 0$ and $M_b \to 0$, while retaining shock-capturing properties in the case of high Mach number flows. To face the multi-scale nature of the model, we rely on the implicit-explicit approach, whereas the divergence-free involution is preserved by means of mimetic finite difference operators [45, 50, 48] on collocated meshes. For the low Mach number regime, we follow the schemes forwarded in [20] for compressible viscous flows. This approach has been
extended to magnetohydrodynamics in [10,25], thus devising asymptotic preserving methods for the hydrodynamic sub-system of the MHD model up to third order space-time accuracy. On staggered meshes, semi-implicit finite volume schemes have been developed in [31], and hybrid finite volume/finite element methods are presented in [37]. These schemes are second order accurate in space and first order accurate in time. For the low Alfvén regime, the only semi-implicit scheme currently available in the literature is presented in [36] with first order in time accuracy, where staggered grids are used, and the implicit sub-system results to be non-linearly coupled, hence requiring to resort to a fixed point algorithm to linearize the problem combined with a nested Newton algorithm which is very challenging in practice.

In this work, we propose a novel semi-implicit finite volume scheme that is second order accurate in both space and time. The semi-implicit IMEX strategy of [12] is adopted to devise an overall second-order all Mach and Alfvén number scheme that only involves linear implicit systems for the fast scales. Consequently, neither Newton solvers nor Picard iteration techniques will be needed. A standard explicit TVD finite volume method is employed for the nonlinear convective terms. To respect the involution on the magnetic field, simple collocated grids are used onto which we design a structure-preserving div-curl operator that is retrieved as the second order version of the div-curl discontinuous Galerkin operators recently forwarded in [17]. To this end, the MHD model is reformulated in terms of the magnetic potential which generates the magnetic field by its curl.

The paper is structured as follows. In Section 2, we briefly recall the model of ideal MHD equations and its non-dimensional formulation in terms of Mach and Alfvén number. Then, the numerical scheme is introduced in Section 3 where we first focus on the operator splitting in a one-dimensional framework and subsequently describe the multi-dimensional time semi-discrete scheme based on the magnetic potential instead of evolving the magnetic field. In the multi-dimensional case, we especially focus on the cell-centered operators that are designed to preserve exactly the div-curl involution constraint for the magnetic potential. This finally yields a divergence-free fully discrete multi-dimensional finite volume scheme. In Section 4, a series of test problems in up to three space dimensions is presented to numerically verify the performance of the numerical schemes at various Mach and Alfvén numbers. A section providing a summary and an outlook to future developments concludes the paper.

## 2 Governing equations

We consider the ideal MHD equations on a three-dimensional domain $\Omega \subset \mathbb{R}^3$ where $\mathbf{x} = (x, y, z) \in \Omega$ is the vector of spatial coordinates and $t \in \mathbb{R}_0^+$ denotes the time coordinate. They write

$$ \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{v} \\ \rho E \\ \mathbf{B} \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{v} \\ \mathbf{v} \otimes \mathbf{v} + \left( p + \frac{\|\mathbf{B}\|^2}{8\pi} \right) \mathbf{I} - \frac{1}{4\pi} \mathbf{B} \otimes \mathbf{B} \\ \left( \rho E + p + \frac{\|\mathbf{B}\|^2}{8\pi} \right) \mathbf{v} - \frac{1}{4\pi} \mathbf{B}(\mathbf{v} \cdot \mathbf{B}) \\ \mathbf{B} \otimes \mathbf{v} - \mathbf{v} \otimes \mathbf{B} \end{pmatrix} = \mathbf{0}. \quad (1) $$

Therein, the density is denoted by $\rho$, the pressure by $p$, the velocity field by $\mathbf{v} = (u, v, w)$ and the magnetic field by $\mathbf{B} = (B_x, B_y, B_z)$ which obeys

$$ \nabla \cdot \mathbf{B} = 0, \quad (2) $$

hence imposing that the magnetic field lines are always closed. The identity matrix is denoted by $\mathbf{I}$. The total energy density is composed of the internal energy $\rho e$, the kinetic energy $\rho k$ and the magnetic energy density $m$ as follows:

$$ \rho E = \rho e + \rho k + m \quad \text{with} \quad k = \frac{1}{2} \|\mathbf{v}\|^2, \quad m = \frac{\|\mathbf{B}\|^2}{8\pi}. \quad (3) $$
The system is closed by the ideal gas equation of state (EOS)

$$ \rho e = \frac{p}{\gamma - 1}, \qquad (4) $$

with $\gamma$ representing the ratio of specific heats. Thus, the pressure is a linear function of the total energy density $\rho E$ and can be computed as $p = (\gamma - 1)(\rho E - \rho k - m)$. Further, the sound speed, the so-called total and directional Alfvén speeds are defined respectively by

$$ c = \sqrt{\frac{\gamma p}{\rho}}, \quad b = \sqrt{\frac{\|\mathbf{B}\|^2}{4\pi\rho}}, \quad b_n = \sqrt{\frac{(\mathbf{B} \cdot \mathbf{n})^2}{4\pi\rho}}, \qquad (5) $$

where $\mathbf{n}$ denotes a unit vector pointing normal to the flow direction. The ideal MHD system (1) exhibits 8 characteristic speeds $\lambda_{\{1...8\}}$ in normal direction $\mathbf{n}$ given by

$$ \lambda_{1,8} = \mathbf{v} \cdot \mathbf{n} \mp c_f, \quad \lambda_{2,7} = \mathbf{v} \cdot \mathbf{n} \mp b_n, \quad \lambda_{3,6} = \mathbf{v} \cdot \mathbf{n} \mp c_s, \quad \lambda_{4,5} = \mathbf{v} \cdot \mathbf{n}, \qquad (6) $$

where $c_{s,f}$ denote the slow and fast magneto-sonic wave speeds which read

$$ c_{s,f} = \sqrt{\frac{1}{2} \left( c^2 + b^2 \mp \sqrt{(c^2 + b^2)^2 - 4c^2 b_n^2} \right)}. \qquad (7) $$

For a detailed discussion of the ideal MHD equations (1) and its structure see e.g. [55].

## 2.1 Scaling of the ideal MHD equations

The MHD system (1) is characterized by multiple time scales, which depend on the physical phenomena under consideration. A sub-system of the ideal MHD equations is indeed given by the Euler equations of gas dynamics, in which a stiffness might be originated in the momentum equations by the pressure waves. In this case, the sound speed is much faster than the local fluid velocity, thus weakly compressible flows are retrieved. The scaling parameter is the Mach number, defined as $M_c = \|\mathbf{v}\|/c$, and these flows are typically affected by a low Mach regime. Another interesting case is given by large magnetic fields, where the Alfvén speeds $b$ are significantly faster than the local velocity or even faster than the sound speed $c$. Analogously to the Mach number, we can define a parameter $M_b = \|\mathbf{v}\|/b$ associated to the Alfvén speed $b$, which is referred to as Alfvén Mach number.

To analyze the influence of those parameters on the dynamics of the model, we consider the non-dimensional formulation of the ideal MHD equations (1). To that end, we employ a decomposition of a dimensional variable $q$ in the associated non-dimensional variable $\tilde{q}$ and a reference value $q_{\text{ref}}$ that contains the units of measurements as well. Thus we can write $q = \tilde{q}q_{\text{ref}}$. The convective scale is given by $v_{\text{ref}} = x_{\text{ref}}/t_{\text{ref}}$, where $x_{\text{ref}}$ and $t_{\text{ref}}$ denote the reference length and time respectively. Further, we have a reference density $\rho_{\text{ref}}$, sound speed $c_{\text{ref}}$, pressure $p_{\text{ref}} = \rho_{\text{ref}}c_{\text{ref}}^2$, specific internal energy $e_{\text{ref}} = c_{\text{ref}}^2$ and a reference Alfvén speed $b_{\text{ref}}$ associated to the magnetic field. Inserting this decomposition into the governing equations, the rescaled non-dimensional MHD system writes

$$ \frac{\partial}{\partial t} \begin{pmatrix} \tilde{\rho} \\ \tilde{\rho} \tilde{\mathbf{v}} \\ \tilde{\rho} \tilde{E} \\ \tilde{\mathbf{B}} \end{pmatrix} + \nabla \cdot \begin{pmatrix} \tilde{\rho} \tilde{\mathbf{v}} \\ \tilde{\rho} \tilde{\mathbf{v}} \otimes \tilde{\mathbf{v}} + \left( \frac{\tilde{p}}{M_c^2} + \frac{1}{M_b^2} \frac{\|\tilde{\mathbf{B}}\|^2}{2} \right) \mathbb{I} - \frac{1}{M_b^2} \tilde{\mathbf{B}} \otimes \tilde{\mathbf{B}} \\ \left( \tilde{\rho} \tilde{E} + \frac{\tilde{p}}{M_c^2} + \frac{1}{M_b^2} \frac{\|\tilde{\mathbf{B}}\|^2}{2} \right) \tilde{\mathbf{v}} - \frac{1}{M_b^2} \tilde{\mathbf{B}} (\tilde{\mathbf{v}} \cdot \tilde{\mathbf{B}}) \\ \tilde{\mathbf{B}} \otimes \tilde{\mathbf{v}} - \tilde{\mathbf{v}} \otimes \tilde{\mathbf{B}} \end{pmatrix} = 0, \qquad (8) $$
where the total energy is given by

$$ \tilde{E} = \frac{1}{M_c^2} \tilde{e} + \frac{1}{M_b^2} \frac{\tilde{m}}{\tilde{\rho}} + \tilde{k}. \qquad (9) $$

Introducing the plasma-beta parameter $\beta = \tilde{p}/\tilde{m}$, which is given by ratio between (reference) pressure and (reference) magnetic pressure, the non-dimensional formulation (8) can be rewritten as

$$ \frac{\partial}{\partial t} \begin{pmatrix} \tilde{\rho} \\ \tilde{\rho} \tilde{\mathbf{v}} \\ \tilde{\rho} \tilde{E} \\ \tilde{\mathbf{B}} \end{pmatrix} + \nabla \cdot \begin{pmatrix} \tilde{\rho} \tilde{\mathbf{v}} \otimes \tilde{\mathbf{v}} + \left( \frac{\tilde{p}}{M_c^2} + \frac{1}{M_c^2 \beta} \frac{\|\tilde{\mathbf{B}}\|^2}{2} \right) \mathbf{I} - \frac{1}{M_c^2 \beta} \tilde{\mathbf{B}} \otimes \tilde{\mathbf{B}} \\ \tilde{\rho} \tilde{\mathbf{v}} \otimes \tilde{\mathbf{v}} + \left( \tilde{\rho} \tilde{E} + \tilde{p} + \frac{1}{\beta} \frac{\|\tilde{\mathbf{B}}\|^2}{2} \right) \tilde{\mathbf{v}} - \frac{1}{\beta} \tilde{\mathbf{B}}(\tilde{\mathbf{v}} \cdot \tilde{\mathbf{B}}) \\ \tilde{\mathbf{B}} \otimes \tilde{\mathbf{v}} - \tilde{\mathbf{v}} \otimes \tilde{\mathbf{B}} \end{pmatrix} = 0, \qquad (10) $$

where the total energy is given by

$$ \tilde{E} = \tilde{e} + \frac{1}{\beta} \frac{\tilde{m}}{\tilde{\rho}} + M_c^2 \tilde{k}, \qquad (11) $$

which is closer to the formulation known from the Euler equations with $M_c$ in front of the kinetic energy contribution [20].

## 2.2 Novel 3-split form of the ideal MHD equations

Let us consider only one space dimension, so that the MHD system (1) reduces to

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}(\mathbf{q})}{\partial x} = 0, \qquad (12) $$

with the definitions

$$ \mathbf{q} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho E \\ B_x \\ B_y \\ B_z \end{pmatrix}, \quad \mathbf{f}(\mathbf{q}) = \begin{pmatrix} \rho u \\ \rho u^2 + p + \frac{\|\mathbf{B}\|^2}{8\pi} - \frac{1}{4\pi} B_x^2 \\ \rho uv - \frac{1}{4\pi} B_x B_y \\ \rho uw - \frac{1}{4\pi} B_x B_z \\ \left(\rho E + p + \frac{\|\mathbf{B}\|^2}{8\pi}\right) u - \frac{1}{4\pi} B_x (\mathbf{v} \cdot \mathbf{B}) \\ 0 \\ u B_y - v B_x \\ u B_z - w B_x \end{pmatrix}. \qquad (13) $$

To handle multiple time scales, we resort to a flux splitting technique that permits to design different time discretizations according to the stiffness of each resulting sub-system. Following [36], the fluxes of the MHD system (1) are split into three contributions: i) convective fluxes $\mathbf{f}^c$ related to $u$; ii) pressure
fluxes $f^p$ related to $p$; iii) magnetic fluxes $f^b$ related to $\mathbf{B}$. Here, we propose the following definitions

$$ \mathbf{f}^c(\mathbf{q}) = \begin{pmatrix} \rho u \\ \rho u^2 \\ \rho uv \\ \rho uw \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \quad \mathbf{f}^p(\mathbf{q}) = \begin{pmatrix} 0 \\ p \\ 0 \\ 0 \\ (\rho E + p) u \\ 0 \\ 0 \\ 0 \end{pmatrix}, \quad \mathbf{f}^b(\mathbf{q}) = \begin{pmatrix} 0 \\ \frac{\|\mathbf{B}\|^2}{8\pi} - \frac{1}{4\pi}B_x^2 \\ -\frac{1}{4\pi}B_x B_y \\ -\frac{1}{4\pi}B_x B_z \\ \frac{\|\mathbf{B}\|^2}{8\pi}u - \frac{1}{4\pi}B_x(\mathbf{v} \cdot \mathbf{B}) \\ 0 \\ uB_y - vB_x \\ uB_z - wB_x \end{pmatrix}. \quad (14) $$

Note that the convective flux $f^c$ is convex in the first two components with respect to $(\rho, \rho u)$, and the remaining non-zero components are passive transport with $u$. We further remark that the above splitting is different from the formulation proposed in [59] for the Euler equations which was applied in [31,36] on the MHD system. The difference here is that in (14) the hydrodynamic component of the energy flux is not split into a convective and a pressure part, but is entirely assigned to the pressure fluxes $f^p$. In particular, the magnetic energy contribution $\rho m$ will also play a role in the fluxes $f^p$. The novel 3-split form of the one-dimensional system (12) is then given in compact notation by

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}^c(\mathbf{q})}{\partial x} + \frac{\partial \mathbf{f}^p(\mathbf{q})}{\partial x} + \frac{\partial \mathbf{f}^b(\mathbf{q})}{\partial x} = \mathbf{0}. \quad (15) $$

In the following we detail for each sub-system the associated set of eigenvalues.

1. For the convective sub-system

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}^c(\mathbf{q})}{\partial x} = \mathbf{0}, \quad (16a) $$

we find the following real eigenvalues

$$ \lambda_{1,2,3,4}^c = 0, \quad \lambda_{5,6,7,8}^c = u. \quad (16b) $$

Thus sub-system (16a) is weakly hyperbolic since one eigenvector associated to $\lambda^c = 0$ is missing. The characteristic speeds depend only on the fluid flow.

2. For the pressure sub-system

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}^p(\mathbf{q})}{\partial x} = \mathbf{0}, \quad (17a) $$

we find the following eigenvalues

$$ \lambda_{1,2,3,4,5,6}^p = 0, \quad \lambda_{7,8}^p = \frac{1}{2} \left( u \mp \sqrt{u^2 + 4(c^2 - (\gamma - 1)(m + k + u^2))} \right), \quad (17b) $$

with a full set of corresponding eigenvectors. Compared to the eigenvalues of the pressure sub-system derived in [31,36], we notice the additional contribution of the specific magnetic and kinetic energy. We remark that the eigenvalues $\lambda_{7,8}^p$ might become complex in the case of highly magnetized fluids, thus the pressure sub-system may result to be elliptic.
3. For the magnetic sub-system

$$ \frac{\partial \mathbf{q}}{\partial t} + \frac{\partial \mathbf{f}^b(\mathbf{q})}{\partial x} = \mathbf{0}, \qquad (18a) $$

we find the following real eigenvalues

$$ \lambda_{1,2,3,4}^b = 0, \quad \lambda_{5,6}^b = \frac{1}{2} \left( u \mp \sqrt{u^2 + 4 \left( \frac{B_x}{\sqrt{4\pi\rho}} \right)^2} \right), \quad \lambda_{7,8}^b = \frac{1}{2} \left( u \mp \sqrt{u^2 + 4 \left( \frac{\|\mathbf{B}\|^2}{\sqrt{4\pi\rho}} \right)^2} \right). \qquad (18b) $$

Associated with $\lambda^b = 0$ there are two eigenvectors missing thus the magnetic sub-system is weakly hyperbolic.

Looking at the eigenvalues of the three sub-systems, it is clear that for low Mach number flows, the acoustic sound speed $c$ in (17b) becomes dominant compared to the material speed $u$ of the flow, hence needing an implicit discretization of the terms related to the pressure. This is also confirmed by analyzing the non-dimensional system (8) for $M_c \to 0$. Likewise, if the magnitude of the magnetic field $\|\mathbf{B}\|$ increases, i.e. when $\beta \to 0$ in (10), the MHD model is dominated by the magnetic sub-system (18a) with eigenvalues $\lambda^b$. As such, an implicit treatment of the terms related to the magnetic field would allow for a larger stability region of the resulting numerical method. The numerical scheme which will be designed in this work starts from these considerations.

## 3 Numerical scheme

Our objective is to construct a numerical scheme for the ideal MHD equations (1) which is stable independently of the parameters $M_c$ and $\beta$ which are associated to the scaling of the pressure (and internal energy) and the strength of the magnetic field, respectively. To achieve this goal, we rely on the novel 3-split formulation (15), dividing the fluxes into explicitly and implicitly treated terms. In particular, the pressure and magnetic sub-systems with $f^P$ and $f^b$ are discretized implicitly, while we retain an explicit scheme for the convective sub-system with $f^c$.

This approach simultaneously yields two main advantages: i) no degradation of accuracy occurs which would be caused for explicit schemes by excessive numerical viscosity for low Mach number flows as it scales with the inverse of the Mach number, see e.g. [29, 46]; ii) improvement of the computational efficiency of the scheme in the low Mach number or low plasma-beta regime, since the stability of the scheme is guaranteed by a time step imposed only by the velocity of the local fluid flow $v$.

However, the resulting implicit system is in general nonlinear due to the non-linearity of the implicitly treated flux terms. Moreover, this implies a huge computational overhead since nonlinear solvers would be required to solve large coupled systems. To overcome this problem, the non-linearity can be treated by adopting a fixed point iteration algorithm, that has been successfully applied to low Mach number flows for the Euler and MHD systems [32, 31]. If the magnetic sub-system is discretized implicitly as done in [36], the nonlinear coupling is even stronger, thus a fixed point method combined with a nested Newton solver [24] has to be employed, which makes the resulting method quite demanding in terms of implementation and computational efficiency. Furthermore, there is no guarantee of convergence since the contraction constant which satisfies the Banach-Cacciopoli theorem is actually unknown.

Therefore, our goal is to give a numerical scheme which only includes linear systems to be solved reducing the computational costs without modifying the model equations (1). Furthermore, a 3-split form as given in (15) and proposed in [36] is not suitable for the class of IMEX Runge-Kutta methods [53, 3] that allows only the splitting according to two time scales, namely a slow and a fast scale. Hence dealing with a 2-split form of the governing equations is favorable and we will detail how to resort to a semi-implicit IMEX discretization of the fluxes, following the time integrators proposed in [12] yielding an overall second-order time discretization.

Besides the treatment in time, the MHD equations require to be compliant with the divergence-free property of the magnetic field (2). A widespread technique consists in the use of staggered grids [62],
since it has been shown that an electric field defined on a mesh that is dual with respect to the mesh of the magnetic field allows the divergence-free condition to be preserved up to machine accuracy. Nevertheless, the staggering of the meshes can become cumbersome, especially in higher dimensions. Consequently, we prefer to adopt a collocated approach utilizing a reformulation of the governing equations in terms of the magnetic potential **A** [43, 44]. Since the magnetic field **B** is the curl of **A**, a mimetic finite difference div-curl operator will then be designed following the general approach introduced in [17].

Summarizing, we will first derive the time semi-discrete numerical scheme in one space dimension where the focus lies on the semi-implicit formulation of the flux splitting. Since the divergence-free constraint reduces to $\partial_x(B_x) = 0$, we solve for the magnetic field directly. Secondly, the time semi-discrete numerical scheme in multiple space dimensions is considered, where the magnetic field **B** is replaced by the curl of the magnetic potential, i.e. $\mathbf{B} = \nabla \times \mathbf{A}$. Thirdly, the extension to second order time accuracy is carried out relying on semi-implicit IMEX time marching schemes. Then, the details of the spatial discretization are provided, including the structure-preserving div-curl operator needed to respect the solenoidal property of the magnetic field at the discrete level upon which the section is concluded with the description of the fully discrete numerical scheme.

*Space and time computational domain.* The three-dimensional computational domain is defined in terms of the spatial intervals as $\Omega = [x_{\min}; x_{\max}] \times [y_{\min}; y_{\max}] \times [z_{\min}; z_{\max}]$, and it is discretized using a Cartesian grid composed of $N_x \times N_y \times N_z$ cells. The characteristic mesh sizes are then given by

$$ \Delta x = \frac{x_{\max} - x_{\min}}{N_x}, \quad \Delta y = \frac{y_{\max} - y_{\min}}{N_y}, \quad \Delta z = \frac{z_{\max} - z_{\min}}{N_z}. \qquad (19) $$

The time interval $T = [0; t_f]$ is discretized by time steps $\Delta t = t^{n+1} - t^n$ which are always subject to a CFL-type stability condition that is only based on the characteristic speeds of the convective sub-system (16a). Thus, with $|\lambda^c|^n$ denoting the maximum eigenvalue of the convective sub-system (16b) at the time level $t^n$, the CFL condition is given by

$$ \Delta t \le \text{CFL} \frac{(\Delta x \Delta y \Delta z)^{1/3}}{\max_{\Omega} |\lambda^c|^n}. \qquad (20) $$

### 3.1 Time semi-discrete scheme in one space dimension

We consider the one-dimensional MHD system (12)-(13). The divergence free condition on the magnetic field in one dimension reduces to $\partial_x(B_x) = 0$, thus the first component of the magnetic field needs to be constant, that is $B_x = \text{const}$.

Applying a forward Euler method on the interval $[t^n; t^{n+1}]$ to the convective sub-system (16a), we obtain the intermediate states $\mathbf{q}^*$ for density and momentum as

$$ \mathbf{q}^* = \mathbf{q}^n - \Delta t \frac{\partial \mathbf{f}^c(\mathbf{q}^n)}{\partial x}. \qquad (21) $$

The remaining flux terms will be considered implicitly, leading to a coupling of the magnetic field with the momentum and energy equation. The two implicit sub-systems will be solved sequentially.

1. The coupling of the magnetic field with the momentum equation provides the magnetic field at the new time level $\mathbf{B}^{n+1}$.

2. The coupling of the energy with the momentum equation yields the total energy density at the new time level $(\rho E)^{n+1}$.
The time discretization of the implicit magnetic subsystem reads

$$ (\rho u)^{n+1} = \rho u^* - \Delta t \frac{\partial}{\partial x} \left( p^n + \frac{\mathbf{B}^n \cdot \mathbf{B}^{n+1}}{8\pi} - \frac{1}{4\pi} B_x^2 \right), \quad (22a) $$

$$ (\rho v)^{n+1} = \rho v^* - \Delta t \frac{\partial}{\partial x} \left( -\frac{1}{4\pi} B_x B_y^{n+1} \right), \quad (22b) $$

$$ (\rho w)^{n+1} = \rho w^* - \Delta t \frac{\partial}{\partial x} \left( -\frac{1}{4\pi} B_x B_z^{n+1} \right), \quad (22c) $$

$$ B_y^{n+1} = B_y^n - \Delta t \frac{\partial}{\partial x} \left( \frac{(\rho u)^{n+1}}{\rho^{n+1}} B_y^n - \frac{(\rho v)^{n+1}}{\rho^{n+1}} B_x^n \right), \quad (22d) $$

$$ B_z^{n+1} = B_z^n - \Delta t \frac{\partial}{\partial x} \left( \frac{(\rho u)^{n+1}}{\rho^{n+1}} B_z^n - \frac{(\rho w)^{n+1}}{\rho^{n+1}} B_x^n \right). \quad (22e) $$

Note that the density at the new time $\rho^{n+1}$ is known after the explicit step (21), and $B_x$ is constant. Furthermore, the pressure $p^n$ in the momentum equation (22a) accounts for more physical information in the computation of the new magnetic field. In addition, this contribution is treated explicitly inspired by the approach forwarded in [21] for an implicit treatment of the viscous terms in the compressible Navier-Stokes equations. We remark that we could even not consider the pressure contribution at this stage, since here the momentum equation is needed only to provide a flux for the magnetic sub-system, and we are not solving for the momentum at the new time level. Formal substitution of the evolution equations (22a)-(22c) into the equation for the magnetic field components (22d)-(22e) yields the following sub-system for the unknowns $B_y^{n+1}$ and $B_z^{n+1}$:

$$ B_y^{n+1} = B_y^* + \Delta t^2 \frac{\partial}{\partial x} \left( \frac{B_y^n}{\rho^{n+1}} \frac{\partial}{\partial x} \left( \frac{B_y^n B_y^{n+1} + B_z^n B_z^{n+1}}{8\pi} \right) + \frac{B_x}{\rho^{n+1}} \frac{\partial}{\partial x} \left( \frac{B_x B_y^{n+1}}{4\pi} \right) \right), \quad (23a) $$

$$ B_z^{n+1} = B_z^* + \Delta t^2 \frac{\partial}{\partial x} \left( \frac{B_z^n}{\rho^{n+1}} \frac{\partial}{\partial x} \left( \frac{B_y^n B_y^{n+1} + B_z^n B_z^{n+1}}{8\pi} \right) + \frac{B_x}{\rho^{n+1}} \frac{\partial}{\partial x} \left( \frac{B_x B_z^{n+1}}{4\pi} \right) \right), \quad (23b) $$

where

$$ B_y^* = B_y^n - \Delta t \frac{\partial}{\partial x} \left( \frac{B_y^n}{\rho^{n+1}} \left( \rho u^* - \Delta t \frac{\partial}{\partial x} \left( p^n - \frac{B_x^2}{8\pi} \right) \right) + \frac{B_x}{\rho^{n+1}} \rho v^* \right), \quad (24) $$

$$ B_z^* = B_z^n - \Delta t \frac{\partial}{\partial x} \left( \frac{B_z^n}{\rho^{n+1}} \left( \rho u^* - \Delta t \frac{\partial}{\partial x} \left( p^n - \frac{B_x^2}{8\pi} \right) \right) + \frac{B_x}{\rho^{n+1}} \rho w^* \right). \quad (25) $$

Thanks to the semi-implicit discretization of the term $||\mathbf{B}||^2 \approx \mathbf{B}^n \cdot \mathbf{B}^{n+1}$ in the $x$-momentum equation (22a), the above system is linear.

Finally, we can write the implicit system for the energy. Since the magnetic field at the new time level $\mathbf{B}^{n+1}$ is known from the solution of (23), the approaches introduced in implicit-explicit schemes for the Euler equations can be readily applied. Following [3], the energy sub-system reads

$$ (\rho u)^{n+1} = (\rho u)^* - \Delta t \frac{\partial}{\partial x} \left( p^{n+1} + \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} - \frac{1}{4\pi} B_x^2 \right) \quad (26a) $$

$$ (\rho E)^{n+1} = (\rho E)^n - \Delta t \frac{\partial}{\partial x} \left( (\rho E^n + p^n) \frac{(\rho u)^{n+1}}{\rho^{n+1}} + \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} u^n - \frac{1}{4\pi} B_x (\mathbf{v}^n \cdot \mathbf{B}^{n+1}) \right). \quad (26b) $$

Let us notice that a semi-implicit scheme has again been applied to the split fluxes (14), in which the chosen discretization for the terms $(\rho E^n + p^n) u^{n+1}$ in $f^p$ and $\frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} u^n$ in $f^b$ allows to retrieve a linear sub-system. Indeed, defining

$$ p^{n+1} = (\rho E^{n+1} - \rho k^n - m^{n+1})(\gamma - 1), \quad (27) $$
and inserting the x-momentum equation (26a) into (26b), the following elliptic equation on the total energy is obtained

$$ (\rho E)^{n+1} = (\rho E)^{\star} + (\gamma - 1)\Delta t^2 \frac{\partial}{\partial x} \left( \frac{\rho E^n + p^n}{\rho^{n+1}} \frac{\partial}{\partial x} \left( (\rho E)^{n+1} \right) \right), \quad (28) $$

where

$$ \rho E^{\star} = \rho E^n - \Delta t \frac{\partial}{\partial x} \left( \frac{\rho E^n + p^n}{\rho^{n+1}} (\rho u)^{\star\star} + \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} u^n - \frac{1}{4\pi} B_x (\mathbf{v}^n \cdot \mathbf{B}^{n+1}) \right), \quad (29) $$

$$ (\rho u)^{\star\star} = (\rho u)^{\star} - \Delta t \frac{\partial}{\partial x} \left( (2 - \gamma) \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} - (\gamma - 1) \rho^n k^n + \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} - \frac{1}{4\pi} B_x^2 \right). \quad (30) $$

We recall that the density $\rho^{n+1}$ at the new time level is already known from the explicit step (21) as the mass flux is contained in the explicitly treated convective sub-system given in (14). The implicit energy equation (28) is thus a linear problem for the scalar unknown $(\rho E)^{n+1}$. This is because of: i) the choice of the ideal gas EOS (4); ii) the semi-implicit time discretization in the pressure and magnetic fluxes $\mathbf{f}^p$ and $\mathbf{f}^b$ of (14).

Once the total energy at the new time level is obtained, the x-component of the momentum can be updated according to (26a) and the y- and z-components according to (22b) and (22c), respectively.

**Remark 1 (On the stability of the time semi-discrete scheme)** In the semi-implicit discretizations designed for both the magnetic and the energy sub-systems, given by (23) and (26), respectively, we make use of a time linearization of the nonlinear terms such that we only have to solve linear systems. Consequently, the proposed schemes cannot be expected to be unconditionally stable, as it would be granted by solving the corresponding original nonlinear sub-systems. However, the time step of our schemes is always limited by the convective sub-system, where we have a CFL stability condition oriented towards the fastest material wave as given in (20). Therefore, we do not use arbitrarily large time steps in our numerical scheme, but they are always limited by the fastest material wave speed. For high Mach number flows, the flow speed is indeed responsible for the most severe restriction, while for low Mach number flows only weak shocks might arise. In those regimes, our scheme is stable as is also numerically verified in the numerical results presented in Section 4. Nevertheless, a thorough study and analysis of the stability condition of the proposed semi-implicit linearization would be necessary.

## 3.2 Time semi-discrete scheme in multiple space dimensions

In multiple space dimension, the preservation of the divergence-free property of the magnetic field at the discrete level is extremely important to ensure physically admissible solutions. According to [23], violating the $\nabla \cdot \mathbf{B} = 0$ constraint leads to non-physical phenomena, e.g spurious numerical monopoles may be generated or magnetic field lines may break into open lines. To guarantee that condition (2) is obeyed, instead of solving for $\mathbf{B}$ directly, as done in the one-dimensional case, we use the magnetic vector potential $\mathbf{A}$ which generates the magnetic field by its curl. Thus, since $\nabla \cdot \mathbf{B} = \nabla \cdot (\nabla \times \mathbf{A}) = 0$, the objective of the numerical scheme is to preserve the div-curl property exactly up to machine precision. The details will be discussed later in Section 3.4. Here, we provide the semi-discrete scheme in time for the ideal MHD system in multiple space dimensions equipped with $\mathbf{A}$ instead of $\mathbf{B}$.

The system of equations is non-conservative for the time evolution of $\mathbf{A}$, and it is given by

$$ \frac{\partial}{\partial t} \begin{pmatrix} \rho \\ \rho \mathbf{v} \\ \rho E \\ \mathbf{A} \end{pmatrix} + \nabla \cdot \begin{pmatrix} \rho \mathbf{v} \\ \rho \mathbf{v} \otimes \mathbf{v} + \left( p + \frac{\|\nabla \times \mathbf{A}\|^2}{8\pi} \right) \mathbf{I} - \frac{1}{4\pi} (\nabla \times \mathbf{A}) \otimes (\nabla \times \mathbf{A}) \\ \left( \rho E + p + \frac{\|\nabla \times \mathbf{A}\|^2}{8\pi} \right) \mathbf{v} - \frac{1}{4\pi} (\nabla \times \mathbf{A})(\mathbf{v} \cdot (\nabla \times \mathbf{A})) \\ \mathbf{0} \end{pmatrix} + \begin{pmatrix} 0 \\ \mathbf{0} \\ 0 \\ (\nabla \times \mathbf{A}) \times \mathbf{v} \end{pmatrix} = \mathbf{0}, \quad (31) $$
with the identity

$$ B - \nabla \times A = 0. \qquad (32) $$

The structure of the time semi-discrete scheme is the same as in the one-dimensional case, but the magnetic sub-system is now written in terms of **A**. Consequently, we focus only on the magnetic sub-system where $(\rho\mathbf{v})^*$ and $\rho^{n+1}$ are known from the convective step (21). Then, assuming that $\mathbf{B}^n = \nabla \times \mathbf{A}^n$ is fulfilled at the current time step, the magnetic sub-system writes

$$ (\rho \mathbf{v})^{n+1} = (\rho \mathbf{v})^* - \Delta t \nabla \cdot \left( \left( p^n + \frac{\mathbf{B}^n \cdot (\nabla \times \mathbf{A}^{n+1})}{8\pi} \right) \mathbb{I} - \frac{\mathbf{B}^n \otimes (\nabla \times \mathbf{A}^{n+1})}{4\pi} \right), \qquad (33a) $$

$$ \mathbf{A}^{n+1} = \mathbf{A}^n - \Delta t \mathbf{B}^n \times \frac{(\rho \mathbf{v})^{n+1}}{\rho^{n+1}}. \qquad (33b) $$

Analogously to the one-dimensional case, plugging the momentum equations (33a) into (33b) yields

$$ \mathbf{A}^{n+1} = \mathbf{A}^* + \Delta t \mathbf{B}^n \times \frac{\Delta t}{\rho^{n+1}} \nabla \cdot \left( \left( \frac{\mathbf{B}^n \cdot (\nabla \times \mathbf{A}^{n+1})}{8\pi} \right) \mathbb{I} - \frac{\mathbf{B}^n \otimes (\nabla \times \mathbf{A}^{n+1})}{4\pi} \right) \qquad (34) $$

where

$$ \mathbf{A}^* = \mathbf{A}^n - \Delta t \mathbf{B}^n \times \frac{1}{\rho^{n+1}} ((\rho \mathbf{v})^* - \Delta t \nabla p^n). \qquad (35) $$

We remark that a semi-implicit discretization has been employed for the terms $\|\nabla \times \mathbf{A}\|^2 \approx (\nabla \times \mathbf{A}^n) \cdot (\nabla \times \mathbf{A}^{n+1}) = \mathbf{B}^n \cdot (\nabla \times \mathbf{A}^{n+1})$ and $(\nabla \times \mathbf{A}) \otimes (\nabla \times \mathbf{A}) \approx (\nabla \times \mathbf{A}^n) \otimes (\nabla \times \mathbf{A}^{n+1}) = \mathbf{B}^n \otimes (\nabla \times \mathbf{A}^{n+1})$. The magnetic potential at the new time level permits to easily evaluate the new magnetic field as $\mathbf{B}^{n+1} = \nabla \times \mathbf{A}^{n+1}$.

Analogously to the one-dimensional case, we obtain the implicit energy update (28) in a multi-dimensional setting, which remains a linear system with the scalar unknown $(\rho E)^{n+1}$. After solving the energy sub-system, the momentum equations can be updated according to

$$ (\rho \mathbf{v})^{n+1} = (\rho \mathbf{v})^* - \Delta t \nabla \cdot \left( \left( p^{n+1} + \frac{\|\mathbf{B}^{n+1}\|^2}{8\pi} \right) \mathbb{I} - \frac{\mathbf{B}^{n+1} \otimes \mathbf{B}^{n+1}}{4\pi} \right), \qquad (36) $$

where the new magnetic field $\mathbf{B}^{n+1}$ is used, and the new pressure $p^{n+1}$ is computed by (27).

We remark at this point, that the scheme is asymptotic preserving (AP) for $M_b = 1$ and $M_c$ tending to zero for well-prepared initial data towards the incompressible Euler equations. For a proof see [49], where the AP property is proven in the context of an IMEX finite volume scheme for single-temperature two-fluid all speed flows. The numerical scheme for the Euler sub-system coincides with the scheme presented in this work, thus the proof can be applied directly on the time semi-discrete scheme presented in this section for $M_b = O(1)$.

### 3.3 Second order time discretization

The semi-discrete scheme presented in the previous section is first order accurate in time. To achieve second order of accuracy, we rely on the class of semi-implicit IMPLICIT-EXPLICIT (IMEX) time integrators proposed in [12]. These Runge-Kutta methods perfectly fit to our semi-discrete scheme since we have applied a time linearization to the 3-split form of the fluxes (14). Such linearization is typically referred to as semi-implicit technique.

Without loss of generality, let us consider the one-dimensional system (12)-(13) and the associated time semi-discrete scheme presented in Section 3.1. According to [12], the governing equations can be written in the form of an autonomous system

$$ \frac{\partial \mathbf{q}}{\partial t} = \mathcal{H}(\mathbf{q}_E(t), \mathbf{q}_I(t)), \qquad (37) $$
with initial condition $\mathbf{q}_0 = \mathbf{q}(t=0)$. The function $\mathcal{H}$ represents the spatial approximation of the flux terms (13) that will be detailed later, and it is a function of two arguments, namely the explicit and the implicit state vector $\mathbf{q}_E$ and $\mathbf{q}_I$, respectively. To integrate the system (37) in time, Implicit-Explicit Runge-Kutta schemes are used, see e.g. [53] for Stiffly Accurate (SA) IMEX-RK methods widely used for the construction of multi-scale schemes. An IMEX-RK scheme is described by a double Butcher tableau

$$ \frac{\tilde{c} | \tilde{A} |}{b^{\top}} = \frac{c | A |}{b^{\top}}, \qquad (38) $$

with the matrices $(\tilde{A}, A) \in \mathbb{R}^{s \times s}$ and the vectors $(\tilde{c}, c, \tilde{b}, b) \in \mathbb{R}^s$. The tilde symbol refers to the explicit scheme and matrix $\tilde{A} = (\tilde{a}_{ij})$ is a lower triangular matrix with zero elements on the diagonal, while $A = (a_{ij})$ is a triangular matrix which accounts for the implicit scheme, thus having non-zero elements on the diagonal. Here, we adopt the LSDIRK2 scheme from [3], which is represented by the following tableaux

$$ \frac{\begin{array}{c|cc} 0 & 0 & 0 \\ \tilde{c} & \tilde{c} & 0 \\ \hline 1 - \alpha & \alpha & 0 \end{array}}{1 - \alpha \alpha} = \frac{\begin{array}{c|cc} \alpha & \alpha & 0 \\ 1 & 1 - \alpha & \alpha \\ \hline 1 - \alpha & \alpha & 0 \end{array}}{1 - \alpha \alpha}, \qquad (39) $$

with $\alpha = 1 - 1/\sqrt{2}$ and $\tilde{c} = 1/(2\alpha)$. Since $\tilde{b} = b$ and $b$ coincides with the last stage in $A$ in (39), the LSDIRK2 scheme is SA which is a crucial property for assuring asymptotic consistency and accuracy of the scheme, see also [15]. To obtain a semi-implicit IMEX-RK method, let us first set $\mathbf{q}_E^n = \mathbf{q}^n$ and $\mathbf{q}_I^n = \mathbf{q}^n$, then the stage fluxes for $i=1, \dots, s$ are calculated as

$$ \mathbf{q}_E^i = \mathbf{q}_E^n + \Delta t \sum_{j=1}^{i-1} \tilde{a}_{ij} k_j, \quad 2 \le i \le s, \qquad (40a) $$

$$ \tilde{\mathbf{q}}_I^i = \mathbf{q}_E^n + \Delta t \sum_{j=1}^{i-1} a_{ij} k_j, \quad 2 \le i \le s, \qquad (40b) $$

$$ k_i = \mathcal{H}(\mathbf{q}_E^i, \tilde{\mathbf{q}}_I^i + \Delta t a_{ii} k_i), \quad 1 \le i \le s. \qquad (40c) $$

Note that the evaluation of the stage fluxes (40c) requires the solution of an implicit system, that corresponds to the implicit magnetic and energy sub-systems given by (23) and (28), respectively.

Finally, the numerical solution is updated with $\mathbf{q}^{n+1} = \mathbf{q}^s$, i.e. the solution coincides with the last Runge-Kutta stage because of the Stiffly Accurate property of the chosen IMEX scheme (39).

## 3.4 Discrete spatial operators

A triple index $ijk$ referred to each space direction allows a cell to be uniquely identified, whereas half indices are adopted to label cell interfaces in $x-$, $y-$ and $z$-direction, namely $(i + \frac{1}{2}jk)$, $(ij + \frac{1}{2}k)$ and $(ijk + \frac{1}{2})$, respectively. The cell barycenter has coordinates $(x_i, y_j, z_k)$, and $|\omega^{ijk}| = \Delta x \Delta y \Delta z$ is the cell volume. Collocated grids are adopted for the spatial discretization of the numerical solution, thus the discrete representation of a generic scalar quantity $q(\mathbf{x})$ at a fixed time is expressed in each cell in the finite volume sense as the cell average

$$ q_{ijk} = \frac{1}{|\omega_{ijk}|} \int_{\omega_{ijk}} q(\mathbf{x}) \, d\mathbf{x}, \qquad (41) $$

Likewise, we denote by $\mathbf{q}(\mathbf{x}) = (q_x, q_y, q_z)$ a generic vector, whose components undergo the same cell centered discretization (41), hence implying $\mathbf{q}_{ijk} = (q_{x,ijk}, q_{y,ijk}, q_{z,ijk})$.

We introduce a set of discrete spatial operators that will be employed in the numerical scheme.
Numerical flux operator $\mathbf{F}(f(q))$. To compute numerical fluxes, we rely on the numerical flux operator $\mathbf{F}(f(q))$ which is defined in the x-direction by

$$ \mathbf{F}_x(f(q)) = \frac{\mathcal{F}_{i+\frac{1}{2}jk}(f(q)) - \mathcal{F}_{i-\frac{1}{2}jk}(f(q))}{\Delta x}, \qquad (42) $$

where $\mathcal{F}(f(q))$ is a Rusanov-type numerical flux function which reads

$$ \mathcal{F}_{i+\frac{1}{2}jk}(q) = \frac{1}{2} \left( f(q_{i+\frac{1}{2}jk}^{+}) + f(q_{i+\frac{1}{2}jk}^{-}) \right) - \frac{1}{2} \alpha_{i+\frac{1}{2}jk} \left( q_{i+\frac{1}{2}jk}^{+} - q_{i+\frac{1}{2}jk}^{-} \right), \qquad (43) $$

with $f(q)$ representing the physical convex flux related to variable $q$, and $q_{i+\frac{1}{2}jk}^{\pm}$ denoting the boundary extrapolated values of the quantity $q$ at the cell interface $(i+\frac{1}{2}jk)$. The numerical dissipation is computed as

$$ \alpha_{i+\frac{1}{2}jk} = \max(|u_{i+1jk}|, |u_{ijk}|), \qquad (44) $$

thus it is only proportional to the material speed $u$ in x-direction. The above operator can be easily extended to the other spatial directions, namely

$$ \mathbf{F}(f(q)) = \mathbf{F}_x(f(q)) + \mathbf{F}_y(f(q)) + \mathbf{F}_z(f(q)). \qquad (45) $$

A first order accurate numerical flux is obtained by choosing $q_{i+\frac{1}{2}jk}^{+} = q_{i+1jk}$ and $q_{i+\frac{1}{2}jk}^{-} = q_{ijk}$ in (43). To attain second order of accuracy, a TVD reconstruction is carried out for the quantity $q(\mathbf{x})$ in each spatial direction, hence obtaining a piecewise polynomial representation $r(\mathbf{x})$ which, in x-direction for cell $\omega_{ijk}$, writes

$$ r_{ijk}(\mathbf{x}) = c_0 q_{ijk} + c_1 (\mathbf{x} - \mathbf{x}_i), \qquad (46) $$

with the expansion coefficients computed using a minmod limiter [58] as

$$ c_0 = 1, \qquad c_1 = \minmod \left( \frac{q_{i+1jk} - q_{ijk}}{\Delta x}, \frac{q_{ijk} - q_{i-1jk}}{\Delta x} \right). \qquad (47) $$

A second order numerical flux is then obtained by setting $q_{i+\frac{1}{2}jk}^{+} = r_{i+1jk}(x_{i+\frac{1}{2}})$ and $q_{i+\frac{1}{2}jk}^{-} = r_{ijk}(x_{i-\frac{1}{2}})$ in (43).

Central flux operator $\mathbf{B}(f(q))$. The operator $\mathbf{B}(f(q))$ is the same as $\mathbf{F}(f(q))$ but zero numerical dissipation is considered, i.e. $\alpha_{i+\frac{1}{2}jk} = 0$ in (43), hence obtaining a second order central flux operator.

Differential operator $\mathbf{H}(h, q)$. This operator applies to two generic quantities $h(\mathbf{x})$ and $q(\mathbf{x})$. In x-direction, it computes a second derivative with second order of accuracy of the form $\partial_x(h \partial_x q)$ by

$$ \mathbf{H}_x(h, q) = \frac{1}{\Delta x} \left( h_{i+\frac{1}{2}jk} \frac{q_{i+1jk} - q_{ijk}}{\Delta x} - h_{i-\frac{1}{2}jk} \frac{q_{ijk} - q_{i-1jk}}{\Delta x} \right), \qquad (48) $$

with the interface quantities $h_{i\pm 1/2jk}$ given by the arithmetic average

$$ h_{i+\frac{1}{2}jk} = \frac{1}{2}(h_{i+1jk} + h_{ijk}), \qquad h_{i-\frac{1}{2}jk} = \frac{1}{2}(h_{ijk} + h_{i-1jk}). \qquad (49) $$

The multidimensional extension writes $\mathbf{H}(h, q) = \mathbf{H}_x(h, q) + \mathbf{H}_y(h, q) + \mathbf{H}_z(h, q)$.

Gradient operator $\mathbf{G}(q)$. The gradient operator is computed at the aid of a second order central finite difference discretization:

$$ \mathbf{G}(q) = \begin{pmatrix} \mathbf{G}_x(q) \\ \mathbf{G}_y(q) \\ \mathbf{G}_z(q) \end{pmatrix} = \begin{pmatrix} \frac{q_{i+1jk} - q_{i-1jk}}{2\Delta x} \\ \frac{q_{ij+1k} - q_{ij-1k}}{2\Delta y} \\ \frac{q_{ijk+1} - q_{ijk-1}}{2\Delta z} \end{pmatrix}. \qquad (50) $$
Curl operator $\mathbf{C}(\mathbf{q})$. The curl operator, directly applied to the vector $\mathbf{q}$, follows from the computation of the gradient operator. Indeed, it simply writes

$$ \mathbf{C}(\mathbf{q}) = \begin{pmatrix} \mathbf{C}_x(\mathbf{q}) \\ \mathbf{C}_y(\mathbf{q}) \\ \mathbf{C}_z(\mathbf{q}) \end{pmatrix} = \begin{pmatrix} \mathbf{G}_x(q_z) - \mathbf{G}_z(q_y) \\ \mathbf{G}_z(q_x) - \mathbf{G}_x(q_z) \\ \mathbf{G}_x(q_y) - \mathbf{G}_y(q_x) \end{pmatrix}. \qquad (51) $$

Divergence operator $\mathbf{D}(\mathbf{q})$. The divergence operator also makes use of the gradient operator (50) as follows

$$ \mathbf{D}(\mathbf{q}) = \mathbf{G}_x(q_x) + \mathbf{G}_y(q_y) + \mathbf{G}_z(q_z). \qquad (52) $$

It is easy to check that the chosen divergence and curl operators are structure preserving, meaning that they satisfy the algebraic identity $\nabla \cdot (\nabla \times \mathbf{q}) = 0$ at the discrete level, namely

$$ \mathbf{D}(\mathbf{C}(\mathbf{q})) = \mathbf{G}_x(\mathbf{C}_x(\mathbf{q})) + \mathbf{G}_y(\mathbf{C}_y(\mathbf{q})) + \mathbf{G}_z(\mathbf{C}_z(\mathbf{q})) = 0. \qquad (53) $$

This relation is crucial for ensuring the solenoidal property of the magnetic field at the discrete level. A proof of the div-curl relation (53) can be found in Appendix A of [17].

Laplace operator $\mathbf{L}(\kappa, \mathbf{q})$. This operator approximates the Laplacian of a scalar function with a second order central finite difference discretization for a space-dependent diffusion coefficient $\kappa$ as

$$ \mathbf{L}(\nu, \mathbf{q}) = \begin{pmatrix} \frac{1}{\Delta x^2} \left( \kappa_{i+\frac{1}{2}jk} (q_{i+1jk} - q_{ijk}) - \kappa_{i-\frac{1}{2}jk} (q_{ijk} - q_{i-1jk}) \right) \\ \frac{1}{\Delta y^2} \left( \kappa_{ij+\frac{1}{2}k} (q_{ij+1k} - q_{ijk}) - \kappa_{ij-\frac{1}{2}k} (q_{ijk} - q_{ij-1k}) \right) \\ \frac{1}{\Delta z^2} \left( \kappa_{ijk+\frac{1}{2}} (q_{ijk+1} - q_{ijk}) - \kappa_{ijk-\frac{1}{2}} (q_{ijk} - q_{ijk-1}) \right) \end{pmatrix}. \qquad (54) $$

### 3.5 Fully discrete scheme in multiple space dimensions

A single time step of the new semi-implicit method from time $t^n$ to $t^{n+1}$ in a fully discrete set-up is made of the following steps.

1. Solve the explicit convective sub-system based on the flux $\mathbf{f}^c$ for density and momentum as

$$ \rho_{ijk}^{\star} = \rho_{ijk}^{n} - \Delta t (\mathbf{F}_{x}((\rho u)^{n}) + \mathbf{F}_{y}((\rho v)^{n}) + \mathbf{F}_{z}((\rho w)^{n})), \qquad (55a) $$

$$ (\rho u)_{ijk}^{\star} = (\rho u)_{ijk}^{n} - \Delta t (\mathbf{F}_{x}((\rho u u)^{n}) + \mathbf{F}_{y}((\rho u v)^{n}) + \mathbf{F}_{z}((\rho u w)^{n})), \qquad (55b) $$

$$ (\rho v)_{ijk}^{\star} = (\rho v)_{ijk}^{n} - \Delta t (\mathbf{F}_{x}((\rho v u)^{n}) + \mathbf{F}_{y}((\rho v v)^{n}) + \mathbf{F}_{z}((\rho v w)^{n})), \qquad (55c) $$

$$ (\rho w)_{ijk}^{\star} = (\rho w)_{ijk}^{n} - \Delta t (\mathbf{F}_{x}((\rho w u)^{n}) + \mathbf{F}_{y}((\rho w v)^{n}) + \mathbf{F}_{z}((\rho w w)^{n})). \qquad (55d) $$

Let us notice that $\rho_{ijk}^{\star} = \rho_{ijk}^{n+1}$, thus the density at the next time level is already available. Furthermore, the explicit contribution of momentum is given by the vector

$$ (\rho \mathbf{v})_{ijk}^{\star} = ((\rho u)_{ijk}^{\star}, (\rho v)_{ijk}^{\star}, (\rho w)_{ijk}^{\star})^{\top}. $$
2. Solve the implicit magnetic sub-system (34) and obtain the vector potential $\mathbf{A}_{ijk}^{n+1}$ from

$$ 
\begin{aligned}
\mathbf{A}_{ijk}^{n+1} - \Delta t \mathbf{B}_{ijk}^n &\times \frac{\Delta t}{\rho_{ijk}^{n+1}} \mathbf{D}\left(\left(\frac{\mathbf{B}_{ijk}^n \cdot \mathbf{C}(\mathbf{A}_{ijk}^{n+1})}{8\pi}\right) \mathbf{I} - \frac{\mathbf{B}_{ijk}^n \otimes \mathbf{C}(\mathbf{A}_{ijk}^{n+1})}{4\pi}\right) \\
&+ \mathbf{L}(\lambda^b \Delta x, A_{x,ijk}^{n+1}) + \mathbf{L}(\lambda^b \Delta y, A_{y,ijk}^{n+1}) + \mathbf{L}(\lambda^b \Delta z, A_{z,ijk}^{n+1}) \\
&= \mathbf{A}_{ijk}^*,
\end{aligned}
\quad (56) $$

with the right hand side

$$ \mathbf{A}_{ijk}^* = \mathbf{A}_{ijk}^n - \Delta t \mathbf{B}_{ijk}^n \times \frac{1}{\rho_{ijk}^{n+1}} \left( (\rho \mathbf{v})_{ijk}^* - \Delta t \mathbf{G}(p_{ijk}^n) \right). \quad (57) $$

Therein, following [44], we have added numerical dissipation by means of the Laplace operator (54) applied to all components of $\mathbf{A}$ to stabilize the numerical scheme especially in high Mach and Alfvén regimes. The diffusion coefficients are proportional to the eigenvalues of the magnetic sub-system (18b), e.g. $\kappa_x = |\lambda^b|\Delta x$. The resulting Laplace operators are discretized implicitly, thus the extra viscosity does not affect the stability of the overall scheme based on the material CFL condition (20). To compute the solution of this linear system (56),(57), we use the GMRES method [57] up to a prescribed tolerance (typically set to $10^{-14}$).

3. Compute the magnetic field and the magnetic energy density at the next time level as

$$ \mathbf{B}_{ijk}^{n+1} = \mathbf{C}(\mathbf{A}_{ijk}^{n+1}), \quad m_{ijk}^{n+1} = \frac{\|\mathbf{B}_{ijk}^{n+1}\|^2}{8\pi}. \quad (58) $$

Thanks to the structure preserving property (53), the magnetic field $\mathbf{B}_{ijk}^{n+1}$ is divergence-free up to machine accuracy.

4. After inserting the definition for $p^{n+1}$ via (27) in the momentum equation (36), one can summarize all known terms, including the implicit magnetic fluxes $\mathbf{f}^b$, in the intermediate quantities $(\rho\mathbf{v})^{**}$ and $(\rho E)^{**}$. Therein, all magnetic contributions are evaluated at the new time level $\mathbf{B}^{n+1}$. Using notation (58) for the magnetic energy density, they read

$$ (\rho u)_{ijk}^{**} = (\rho u)_{ijk}^* - \Delta t \left( \mathbf{B}_x \left( -(\gamma - 1) \rho^n k^n + (2 - \gamma) m^{n+1} - \frac{B_x^{n+1} B_y^{n+1}}{4\pi} \right) + \mathbf{B}_y \left( -\frac{B_x^{n+1} B_y^{n+1}}{4\pi} \right) + \mathbf{B}_z \left( -\frac{B_x^{n+1} B_z^{n+1}}{4\pi} \right) \right), \quad (59a) $$

$$ (\rho v)_{ijk}^{**} = (\rho v)_{ijk}^* - \Delta t \left( \mathbf{B}_x \left( -\frac{B_y^{n+1} B_x^{n+1}}{4\pi} \right) + \mathbf{B}_y \left( -(\gamma - 1) \rho^n k^n + (2 - \gamma) m^{n+1} - \frac{B_y^{n+1} B_z^{n+1}}{4\pi} \right) + \mathbf{B}_z \left( -\frac{B_y^{n+1} B_z^{n+1}}{4\pi} \right) \right), \quad (59b) $$

$$ (\rho w)_{ijk}^{**} = (\rho w)_{ijk}^* - \Delta t \left( \mathbf{B}_x \left( -\frac{B_z^{n+1} B_x^{n+1}}{4\pi} \right) + \mathbf{B}_y \left( -\frac{B_z^{n+1} B_y^{n+1}}{4\pi} \right) + \mathbf{B}_z \left( -(\gamma - 1) \rho^n k^n + (2 - \gamma) m^{n+1} - \frac{B_z^{n+1} B_z^{n+1}}{4\pi} \right) \right), \quad (59c) $$

$$ (\rho E)_{ijk}^* = (\rho E)_{ijk}^n - \Delta t \left( \mathbf{B}_x \left( u^n m^{n+1} - \frac{B_x^{n+1}}{4\pi} \mathbf{v}^n \cdot \mathbf{B}^{n+1} \right) + \mathbf{B}_y \left( v^n m^{n+1} - \frac{B_y^{n+1}}{4\pi} \mathbf{v}^n \cdot \mathbf{B}^{n+1} \right) + \mathbf{B}_z \left( w^n m^{n+1} - \frac{B_z^{n+1}}{4\pi} \mathbf{v}^n \cdot \mathbf{B}^{n+1} \right) \right). \quad (59d) $$
5. Solve the elliptic equation for the total energy $(\rho E)_{ijk}^{n+1}$ according to

$$ (\rho E)_{ijk}^{n+1} - (\gamma - 1) \Delta t^2 (\mathbf{H}_x(h^n, (\rho E)_{ijk}^{n+1}) + \mathbf{H}_y(h^n, (\rho E)_{ijk}^{n+1}) + \mathbf{H}_z(h^n, (\rho E)_{ijk}^{n+1})) = b_{ijk}^{**}, \quad (60) $$

with the right hand side

$$ b_{ijk}^{**} = (\rho E)_{ijk}^{*} - \Delta t (\mathbf{B}_x(h^n(\rho u)^{**}) + \mathbf{B}_y(h^n(\rho v)^{**}) + \mathbf{B}_z(h^n(\rho w)^{**})). \quad (61) $$

Even in this case, the linear system is solved with the GMRES method up to a residual of the order of $10^{-14}$.

6. Update the momentum at the next time level $(\rho \mathbf{v})_{ijk}^{n+1}$ according to

$$ (\rho u)_{ijk}^{n+1} = (\rho u)_{ijk}^{**} - (\gamma - 1) \Delta t \mathbf{B}_x ((\rho E)^{n+1}), \quad (62a) $$

$$ (\rho v)_{ijk}^{n+1} = (\rho v)_{ijk}^{**} - (\gamma - 1) \Delta t \mathbf{B}_y ((\rho E)^{n+1}), \quad (62b) $$

$$ (\rho w)_{ijk}^{n+1} = (\rho w)_{ijk}^{**} - (\gamma - 1) \Delta t \mathbf{B}_z ((\rho E)^{n+1}), \quad (62c) $$

which is equivalent to (36) after the introduction of the new pressure (27).

Summarizing, to stabilize the numerical scheme, we have introduced two sources of artificial viscosity. The first one is the dissipative part of the explicit Rusanov flux proportional to the local fluid velocity $\mathbf{v}$, and the second one via an implicit Laplacian operator in the magnetic sub-system proportional to the local eigenvalue $\lambda^b$ of the magnetic-subsystem. In particular, the magnetic dissipation contributes to the stability of the numerical scheme especially in high Mach and Alfvén Mach number regimes. Since it is treated implicitly, the material time step restriction does not need to be restricted to guarantee the stability of the numerical scheme.

## 4 Numerical results

In the sequel, we present a large suite of test cases for ideal MHD aiming at demonstrating the accuracy, robustness and effectiveness of our new scheme in different regimes of both the Mach and the Alfvén number. If not otherwise specified, we set $\gamma = 3/5$ and the time step $\Delta t$ is computed according to the stability condition (20) with CFL = 0.9, which only depends on the material waves. The novel scheme is labeled with SIFV-EB, that stands for Semi-Implicit Finite Volume method with implicit discretization for the energy (E) and the magnetic field (B). If not stated otherwise, it is always used the second order version of the SIFV-EB scheme in space and time.

*Boundary conditions for the magnetic potential A.* In all the following test problems we always perform a linear extrapolation to determine the boundary conditions for the magnetic potential **A**. Without loss of generality, let us consider the component $A_{x,i}$ for a cell attached to the left boundary $x_{i-1/2} = x_{\min}$ of the computational domain, i.e. $i = 1$. The sought boundary state $A_{x,b}$ in $x$-direction is formally located at $x_b = x_{\min} - \Delta x/2$ and is computed as

$$ A_{x,b} = A_{x,i} - \alpha \Delta x \quad \text{with} \quad \alpha = \frac{A_{x,i+1} - A_{x,i}}{\Delta x}. \quad (63) $$

This procedure applies for all the components of **A** and all the spatial directions $\mathbf{x}$. The Laplace operator $\mathbf{L}$ in (56), which accounts for artificial viscosity in the implicit update of the magnetic potential, can thus be applied without any modification.
## 4.1 Accuracy of the numerical scheme

We test the accuracy of the numerical scheme on a slight modification of the analytical vortex solution forwarded in [4] on the two-dimensional computational domain $\Omega = [-5; 5] \times [-5; 5]$. The vortex is constructed as a perturbation $\delta q$ on top of an unperturbed magnetohydrodynamics background flow given by the constant states $q_0 = (\rho_0, v_0, p_0, 0)$. The perturbations of the velocity and magnetic field are given by

$$ (\delta u, \delta v) = \frac{\tilde{v}}{2\pi} \exp\left(\frac{1-r^2}{2}\right) \cdot (-r \sin(\theta), r \cos(\theta)), \quad (64) $$

$$ (\delta B_x, \delta B_y) = \frac{\tilde{B}}{2\pi} \exp\left(\frac{1-r^2}{2}\right) \cdot (-r \sin(\theta), r \cos(\theta)), \quad (65) $$

with the constants $\tilde{v} = \sqrt{2\pi}$ and $\tilde{B} = \sqrt{4\pi}$. The perturbation of the magnetic field is generated by the magnetic potential perturbation in z-direction $\delta A_z = \frac{\tilde{B}}{2\pi} \exp(1-r^2)$. The pressure is set such that it balances the centrifugal force generated by the circular motion against the centripetal force stemming from the tension in the magnetic field lines. In radial coordinates, with the generic radial position $r = \sqrt{x^2 + y^2}$, the equation for the pressure perturbation writes

$$ \frac{\partial \delta p}{\partial r} = \left( \frac{\rho v_{\theta}^2}{r} - \frac{\delta B_{\theta}^2}{4\pi r} - \frac{\partial (\delta B_{\theta})^2}{8\pi r} \right), \quad (66) $$

with the definitions

$$ \delta v_{\theta} = \frac{\tilde{v}}{2\pi} r \exp\left(\frac{1-r^2}{2}\right), \quad \delta B_{\theta} = \frac{\tilde{B}}{2\pi} r \exp\left(\frac{1-r^2}{2}\right). \quad (67) $$

Integration of (66) yields

$$ \delta p = \frac{1}{2} e^{1-r^2} \left( \frac{1}{8\pi} \frac{\tilde{B}^2}{(2\pi)^2} (1-r^2) - \frac{1}{2} \rho \left( \frac{\tilde{v}}{2\pi} \right)^2 \right). \quad (68) $$

We set the background velocity field to $v_0 = (1, 1, 0)$, thus the vortex is moving diagonally through the computational domain, and the background pressure is $p_0 = 1$. Since the novelty of the scheme is the stability under a material time step which is independent of the sound and Alfvén speeds, we study the experimental order of convergence (EOC) under a variation of the Alfvén speeds. To this end, let us note that simply adding a constant background magnetic field does not yield an analytic solution of the ideal MHD equations in this setup. Thus, since the Alfvén speeds depend on the inverse of the density, the Alfvén Mach number can be decreased by lowering the background density $\rho_0$. We underline that decreasing the density will increase the sound speed as well, and the kinetic influence on the pressure perturbation (68) will decrease. However, the contribution of the magnetic field remains the same for all the cases considered. We remark that, due to a zero background magnetic field, it is not possible to compute the Alfvén Mach number $M_b$ over the whole computational domain. Therefore, the maximal acoustic Mach number $M_c$ and Alfvén speed $b$ are reported in Table 1 for different background densities $\rho_0 = 10^{-k}, k = 0, \dots, 5$.

The simulations are run on a sequence of successively refined meshes until the final time $t_f = 1$. In Table 2, the $L^2$ error in $\rho, u, p, B_x$ and $A_z$ are summarized. We observe that, for all regimes, the expected EOC of two is achieved. Furthermore, the scheme performs well for small values of the density. We recall, however, that the scheme is not naturally positivity preserving.

To demonstrate the advantage of having a scale independent CFL stability condition, we plot in Figure 1 the ratio between a fully explicit time step restricted by the fastest characteristic speeds against the time step given by the CFL condition of our scheme (20). The SIFV-EB scheme allows, especially for low acoustic and Alfvén Mach numbers, to use a time step which is several orders of magnitude larger than the one of an explicit scheme while maintaining second order accuracy in space and time.
<table><caption>**Table 1** MHD traveling vortex. Maximum acoustic Mach number $M_c$ and Alfvén speed $b$ at the initial time.</caption><thead><tr><th rowspan="2"></th><th colspan="6">Background density ρ₀ = 10⁻ᵏ</th></tr><tr><th>k = 0</th><th>k = 1</th><th>k = 2</th><th>k = 3</th><th>k = 4</th><th>k = 5</th></tr></thead><tbody><tr><th>M_c</th><td>1.606E+01</td><td>4.489E-01</td><td>1.529E-01</td><td>4.832E-02</td><td>1.529E-02</td><td>4.832E-03</td></tr><tr><th>b</th><td>1.549E-01</td><td>4.900E-01</td><td>1.549E+00</td><td>4.900E+00</td><td>1.549E+01</td><td>4.900E+01</td></tr></tbody></table>

<table><caption>**Table 2** MHD traveling vortex. $L^2$ error and EOC for the second order SIFV-EB scheme with varying background density $\rho_0$ at final time $t_f = 1$.</caption><thead><tr><th>ρ₀</th><th>Nₓ = Nᵧ</th><th>ρ</th><th>u</th><th>p</th><th>Bₓ</th><th>A_z</th></tr></thead><tbody><tr><td rowspan="4">10<sup>0</sup></td><td>32</td><td>7.13e-02</td><td>4.40e-02</td><td>5.86e-02</td><td>8.26e-02</td><td>4.94e-02</td></tr><tr><td>64</td><td>1.63e-02</td><td>2.13</td><td>9.16e-03</td><td>2.26</td><td>1.65e-02</td><td>1.83</td><td>2.20e-02</td><td>1.91</td><td>1.25e-02</td><td>1.99</td></tr><tr><td>128</td><td>3.79e-03</td><td>2.10</td><td>2.18e-03</td><td>2.07</td><td>4.19e-03</td><td>1.98</td><td>5.66e-03</td><td>1.96</td><td>3.14e-03</td><td>1.99</td></tr><tr><td>256</td><td>9.26e-04</td><td>2.03</td><td>6.79e-04</td><td>1.69</td><td>1.05e-03</td><td>1.99</td><td>1.54e-03</td><td>1.88</td><td>8.12e-04</td><td>1.95</td></tr><tr><td rowspan="4">10<sup>-1</sup></td><td>32</td><td>2.33e-03</td><td>3.50e-02</td><td>6.84e-03</td><td>8.21e-02</td><td>4.84e-02</td></tr><tr><td>64</td><td>6.39e-04</td><td>1.87</td><td>8.55e-03</td><td>2.04</td><td>1.91e-03</td><td>1.84</td><td>2.17e-02</td><td>1.92</td><td>1.23e-02</td><td>1.98</td></tr><tr><td>128</td><td>1.62e-04</td><td>1.98</td><td>2.14e-03</td><td>2.00</td><td>4.90e-04</td><td>1.96</td><td>5.53e-03</td><td>1.97</td><td>3.09e-03</td><td>1.99</td></tr><tr><td>256</td><td>4.08e-05</td><td>2.00</td><td>6.71e-04</td><td>1.67</td><td>1.24e-04</td><td>1.99</td><td>1.51e-03</td><td>1.87</td><td>7.99e-04</td><td>1.95</td></tr><tr><td rowspan="4">10<sup>-2</sup></td><td>32</td><td>4.12e-04</td><td>9.14e-02</td><td>3.90e-03</td><td>9.05e-02</td><td>5.02e-02</td></tr><tr><td>64</td><td>7.64e-05</td><td>2.43</td><td>2.45e-02</td><td>1.90</td><td>9.82e-04</td><td>1.99</td><td>2.32e-02</td><td>1.96</td><td>1.24e-02</td><td>2.02</td></tr><tr><td>128</td><td>1.63e-05</td><td>2.23</td><td>6.27e-03</td><td>1.97</td><td>2.69e-04</td><td>1.87</td><td>5.85e-03</td><td>1.99</td><td>3.08e-03</td><td>2.00</td></tr><tr><td>256</td><td>3.92e-06</td><td>2.06</td><td>1.63e-03</td><td>1.95</td><td>6.83e-05</td><td>1.98</td><td>1.58e-03</td><td>1.89</td><td>7.95e-04</td><td>1.95</td></tr><tr><td rowspan="4">10<sup>-3</sup></td><td>32</td><td>2.80e-04</td><td>2.07e-01</td><td>3.45e-03</td><td>9.44e-02</td><td>4.80e-02</td></tr><tr><td>64</td><td>4.50e-05</td><td>2.64</td><td>5.88e-02</td><td>1.81</td><td>9.63e-04</td><td>1.84</td><td>2.06e-02</td><td>2.20</td><td>8.76e-03</td><td>2.46</td></tr><tr><td>128</td><td>5.95e-06</td><td>2.92</td><td>1.72e-02</td><td>1.77</td><td>1.98e-04</td><td>2.28</td><td>4.45e-03</td><td>2.21</td><td>1.79e-03</td><td>2.29</td></tr><tr><td>256</td><td>8.10e-07</td><td>2.88</td><td>4.69e-03</td><td>1.88</td><td>4.74e-05</td><td>2.06</td><td>1.19e-03</td><td>1.90</td><td>4.61e-04</td><td>1.96</td></tr><tr><td rowspan="4">10<sup>-4</sup></td><td>32</td><td>1.55e-04</td><td>1.67e-00</td><td>1.28e-02</td><td>2.60e-01</td><td>2.15e-01</td></tr><tr><td>64</td><td>4.98e-05</td><td>1.64</td><td>2.25e-01</td><td>2.89</td><td>4.49e-03</td><td>1.51</td><td>8.01e-02</td><td>1.70</td><td>4.92e-02</td><td>2.12</td></tr><tr><td>128</td><td>5.87e-06</td><td>3.09</td><td>3.44e-02</td><td>2.71</td><td>8.21e-04</td><td>2.45</td><td>1.55e-02</td><td>2.37</td><td>7.77e-03</td><td>2.66</td></tr><tr><td>256</td><td>7.22e-07</td><td>3.02</td><td>9.02e-03</td><td>1.93</td><td>1.30e-04</td><td>2.66</td><td>2.45e-03</td><td>2.66</td><td>1.14e-03</td><td>2.76</td></tr><tr><td rowspan="4">10<sup>-5</sup></td><td>32</td><td>1.85e-05</td><td>8.14e-00</td><td>1.32e-02</td><td>2.85e-01</td><td>2.35e-01</td></tr><tr><td>64</td><td>9.85e-06</td><td>0.91</td><td>8.62e-01</td><td>3.24</td><td>6.75e-03</td><td>0.96</td><td>1.24e-01</td><td>1.20</td><td>8.47e-02</td><td>1.47</td></tr><tr><td>128</td><td>1.41e-06</td><td>2.81</td><td>8.42e-02</td><td>3.36</td><td>1.64e-03</td><td>2.04</td><td>3.04e-02</td><td>2.03</td><td>1.70e-02</td><td>2.32</td></tr><tr><td>256</td><td>1.61e-07</td><td>3.12</td><td>2.05e-02</td><td>2.04</td><td>2.71e-04</td><td>2.60</td><td>4.85e-03</td><td>2.65</td><td>2.54e-03</td><td>2.74</td></tr></tbody></table>

## 4.2 1D Riemann Problems

Next, we consider a series of Riemann Problems (RP) which are standard one-dimensional test cases for ideal MHD solvers. The setup is detailed in Table 3 and it is taken from [30]. They have already been used to validate for instance the Osher-type scheme [33] or the HLLEM Riemann solver presented in [30]. We compare the results obtained with our new SIFV-EB scheme against the exact solution, that has been computed with an exact Riemann solver that has kindly been provided by Falle and Komissarov [34,35] and which is labeled in Figures 2 and 3 as reference solution. For an alternative exact Riemann solver of the MHD equations, see also [60]. These test cases verify the robustness of the new SIFV-EB scheme in the case of moderate and high Mach number flows. The computational domain is given by $\Omega = [-0.5; 0.5]$ and the grid is made of 2000 cells along the x-direction with Dirichlet boundaries.

For each RP, depicted in Figures 2 and 3. The SIFV-EB scheme is able to capture the correct shock position and all waves accurately for RP2-RP7. In particular, the overshoot in density and pressure
<figure><img src="image_2.png" /><figcaption>**Fig. 1** MHD traveling vortex. Time evolution of the ratio between the material time step of the SIFV-EB scheme and the magneto-sonic time step of a fully explicit finite volume scheme for different values of the background density $\rho_0$.</figcaption></figure>

<table><caption>**Table 3** Initialization of Riemann problems. Initial states left (L) and right (R) are reported as well as the final time of the simulation $t_f$. In all cases $\gamma = 5/3$ and the position of the initial discontinuity is $x_d = 0$, except for RP2 and RP3 where it is $x_d = -0.1$.</caption><thead><tr><th>Case</th><th></th><th>ρ</th><th>u</th><th>v</th><th>w</th><th>p</th><th>B<sub>x</sub></th><th>B<sub>y</sub></th><th>B<sub>z</sub></th><th>t<sub>f</sub></th></tr></thead><tbody><tr><td rowspan="2">RP1</td><td>L:</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0.75√4π</td><td>+√4π</td><td>0.0</td><td rowspan="2">0.10</td></tr><tr><td>R:</td><td>0.125</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.1</td><td>0.75√4π</td><td>-√4π</td><td>0.0</td></tr><tr><td rowspan="2">RP2</td><td>L:</td><td>1.08</td><td>1.2</td><td>0.01</td><td>0.5</td><td>0.95</td><td>2.0</td><td>3.6</td><td>2.0</td><td rowspan="2">0.2</td></tr><tr><td>R:</td><td>0.9891</td><td>-0.0131</td><td>0.0269</td><td>0.010037</td><td>0.97159</td><td>2.0</td><td>4.0244</td><td>2.0026</td></tr><tr><td rowspan="2">RP3</td><td>L:</td><td>1.7</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.7</td><td>3.899398</td><td>3.544908</td><td>0.0</td><td rowspan="2">0.15</td></tr><tr><td>R:</td><td>0.2</td><td>0.0</td><td>0.0</td><td>-1.496891</td><td>0.2</td><td>3.899398</td><td>2.785898</td><td>2.192064</td></tr><tr><td rowspan="2">RP4</td><td>L:</td><td>1.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>1.0</td><td>1.3√4π</td><td>+√4π</td><td>0.0</td><td rowspan="2">0.16</td></tr><tr><td>R:</td><td>0.4</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.4</td><td>1.3√4π</td><td>-√4π</td><td>0.0</td></tr><tr><td rowspan="2">RP5</td><td>L:</td><td>0.15</td><td>21.55</td><td>1.0</td><td>1.0</td><td>0.28</td><td>0.05</td><td>-2.0</td><td>-1.0</td><td rowspan="2">0.04</td></tr><tr><td>R:</td><td>0.10</td><td>-26.45</td><td>0.0</td><td>0.0</td><td>0.10</td><td>0.05</td><td>+2.0</td><td>+1.0</td></tr><tr><td rowspan="2">RP6</td><td>L:</td><td>1.0</td><td>36.87</td><td>-0.115</td><td>-0.0386</td><td>1.0</td><td>4.0</td><td>4.0</td><td>1.0</td><td rowspan="2">0.03</td></tr><tr><td>R:</td><td>1.0</td><td>-36.87</td><td>0.0</td><td>0.0</td><td>1.0</td><td>4.0</td><td>4.0</td><td>1.0</td></tr><tr><td rowspan="2">RP7</td><td>L:</td><td>1/μ<sub>0</sub></td><td>-1.0</td><td>+1.0</td><td>-1.0</td><td>1.0</td><td>1.0</td><td>-1.0</td><td>1.0</td><td rowspan="2">0.25</td></tr><tr><td>R:</td><td>1/μ<sub>0</sub></td><td>-1.0</td><td>-1.0</td><td>-1.0</td><td>1.0</td><td>1.0</td><td>+1.0</td><td>1.0</td></tr></tbody></table>
in RP1 is an artifact that also appears in approximate solutions produced with the HLLEM scheme from [30], the Osher-type scheme [33] or the Rusanov scheme. Moreover, RP7 contains an isolated steady Alfvén wave which is well resolved. However, small perturbations of density and pressure appear around the jump position. To further investigate the nature of these errors, we compare the solution to an explicit second order MUSCL-FV scheme using a TVD limiter in Figure 4 with the same mesh resolution. We can observe that these artifacts also arise with the same order of magnitude using the explicit solver. This is comparable to the errors that appear in the density and pressure for RP1 which have been reproduced with different approximate Riemann solvers, see e.g. [30,33].

Overall, the new SIFV-EB scheme is clearly more diffusive than an explicit approximate Riemann solver. This is due to the semi-implicit nature and the material time stepping, based on the fluid velocity $u$, that allows much larger time steps. This is illustrated in Figure 5 where a comparison between the material time stepping $\Delta t$ of the SIFV-EB scheme with respect to a fully explicit time step $\Delta t_{\text{MHD}}$, enforced by the fast magneto-sonic wave speeds $\lambda_{1,8}$, is given. Even though most RP are placed in a compressible regime, our scheme allows for larger time steps that is reflected by slightly more diffusive fast traveling waves, while the material waves are captured accurately, as expected. Let us note that the first time step is always given by a magneto-sonic explicit time step $\Delta t = \Delta t_{\text{MHD}}$, since tests RP1, RP3 and RP4 are initialized with a zero velocity $u$ which would lead to an infinitely large first time step and thus to inaccurate numerical results. Since this yields $\Delta t/\Delta t_{\text{MHD}} = 1$ for the first time step, we have plotted the ratio starting from the second time step in Figure 5. Moreover, the time step ratio remains almost constant for those Riemann problems which satisfy a maximum principle for the velocity $u$.

## 4.3 Blast wave

The MHD blast wave problem introduced in [7] represents a challenging test since it contains a strong magnetic field and a jump in the pressure of four orders of magnitude. The initial density is constant $\rho = 1$, the initial velocity field is set to zero $(u, v, w) = (0, 0, 0)$ and the magnetic field is $\mathbf{B} = (70, 0, 0)$. The computational domain $\Omega = [-0.5; 0.5] \times [-0.5; 0.5]$ is paved with $N_x \times N_y = 1024 \times 1024$ cells, and Dirichlet boundary conditions are imposed everywhere. Inside the circular region $r < 0.1$ with $r = \sqrt{x^2 + y^2}$, the pressure is set to $p = 1000$ whereas outside it is assigned with $p = 0.1$. Thus initially, the Alfvén waves are faster than the sound speed inside the circular domain, that are given by $b = 70$ and $c = 37.4166$, respectively. Here, we set $\gamma = 1.4$ and the final time is chosen to be $t_f = 0.01$. The results at final time for density, pressure, magnitude of the velocity computed as $v^2$ as well as the magnetic energy density $m$ are displayed in Figure 6 for the second-order SIFV-EB scheme. We observe that there are no oscillations arising in the simulation and all waves are captured accurately. The results compare qualitatively well with the literature [18,6,36].

## 4.4 Ideal MHD rotor problem

The next test case concerns a benchmark for the ideal MHD equations, the so-called rotor problem proposed in [7], which is characterized by a big density difference inside and outside of the rotor. The computational domain is given by the square $\Omega = [-0.5; 0.5] \times [-0.5; 0.5]$ with Dirichlet boundaries, and it is discretized using $N_x \times N_y = 1024 \times 1024$ control volumes. The background pressure and magnetic field are constant and set to $p = 1$ and $\mathbf{B} = (2.5, 0, 0)^T$, respectively. The corresponding magnetic vector potential in $z$-direction is $A_z = 2.5y$. The density and velocity field are defined differently inside and outside of the rotor, i.e.

$$ \rho = \begin{cases} 10 & \text{if } 0 \le r \le 0.1, \\ 1 & \text{else,} \end{cases} \qquad (u, v, w) = \begin{cases} (0, 0, 10) \times (x, y, z) & \text{if } 0 \le r \le 0.1, \\ (0, 0, 0) & \text{else.} \end{cases} \qquad (69) $$

The ratio of specific heats is set to $\gamma = 1.4$ and the simulation is run until the final time $t_f = 0.25$. In Figure 7, the results for the second-order SIFV-EB scheme are plotted for the density, pressure, acoustic
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Riemann problems. RP1 at time $t_f = 0.1$ (top), RP2 at time $t_f = 0.2$ (middle) and RP3 at time $t_f = 0.15$ (bottom). Left: density $\rho$. Center: pressure $p$. Right: magnetic field component $B_y$.</figcaption></figure>

Mach number $M_c$ and the magnetic energy density $m$. The scheme, though spanning a Mach number regime from weakly compressible to compressible with $M_c > 1$, yields results that compare well with simulations reported in literature see e.g. [31,18].

Furthermore, in Figure 8 we also monitor the error in the divergence of the magnetic field, which is measure in the $L_\infty$ norm over the whole computational domain at all time instants. For both the blast wave and the rotor problem we achieve a discrete divergence-free condition up to machine accuracy. Alongside, the new semi-implicit schemes allows the time step to be still larger than the one of a fully explicit scheme, even if strong shocks occur in these test cases. Looking at the right panel of Figure 8, we clearly see that the SIFV-EB method can lead to a time marching scheme with a time step that is up to four times larger than a classical explicit solver.
<figure><img src="image_4.png" /><figcaption>Fig. 3 Riemann problems. RP4 at time $t_f = 0.16$, RP5 at time $t_f = 0.04$, RP6 at time $t_f = 0.03$ and RP7 at time $t_f = 0.2$ (from top to bottom). Left: density $\rho$. Center: pressure $p$. Right: magnetic field component $B_y$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4 Riemann problem RP7 at time $t_f = 0.2$. Comparison between the second order accurate SIFV-EB scheme and an explicit MUSCL-TVD scheme. Left: density $ρ$. Right: pressure $p$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Riemann problems. Time evolution of the ratio between the material time step of the SIFV-EB scheme and the magneto-sonic time step of a fully explicit finite volume scheme for all the Riemann problems RP1–RP7.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6 Blast problem at time $t_f = 0.01$. Results for the logarithm (base 10) of density and pressure, magnitude of the velocity and magnetic pressure (from top left to bottom right).</figcaption></figure>

## 4.5 Orszag-Tang vortex

We consider the well-known Orszag-Tang vortex for the two-dimensional ideal MHD equations [56, 27, 54]. The initial condition is smooth and reads

$$ (\rho, u, v, w, p, B_x, B_y, B_z) = \left( \frac{25}{36\pi}, -\sin(2\pi y), \sin(2\pi x), 0, \frac{5}{12\pi}, -\frac{\sin(2\pi y)}{\sqrt{4\pi}}, \frac{\sin(4\pi y)}{\sqrt{4\pi}}, 0 \right). \quad (70) $$

The magnetic field **B** is generated by the z-component of the magnetic potential

$$ A_z = \frac{1}{8\pi^{3/2}} (\cos(2\pi x) + 2 \cos(\pi y)). \quad (71) $$

Though initially smooth, the dynamics develop shocks along the diagonal direction in combination with a vortex located at the center of the computational domain, which is defined by $\Omega = [0; 1] \times [0; 1]$ with periodic boundaries. We use a mesh made of $N_x \times N_y = 1024 \times 1024$ cells. Figure 9 depicts snapshots of the pressure computed with our second-order SIFV-EB scheme at times $t = \{0.5, 1, 2, 4\}$. The features
<figure><img src="image_8.png" /><figcaption>Fig. 7 Rotor problem at time $t_f = 0.25$. Results for density, pressure, Mach number and magnetic energy density (from top left to bottom right).</figcaption></figure>

of the flows and the development of the shock are captured well and the results are qualitatively in very good agreement with those reported in the literature [31, 6, 5].

## 4.6 Field loop advection

We solve the magnetic field loop advection problem originally proposed in [39], in the low Mach number version of [31] with $\rho = 1$, an initial velocity field $(u, v, w) = (2, 1, 0)$ and the pressure $p = 10^5$. The computational domain is given by $\Omega = [-1; 1] \times [-0.5; 0.5]$ with periodic boundaries, and we use a mesh composed of $N_x \times N_y = 512 \times 256$ cells. The magnetic field is prescribed by the magnetic vector potential in z-direction given by

$$A_z = \begin{cases} A_0(R-r) & \text{if } r \le R, \\ 0 & \text{else,} \end{cases} \qquad (72)$$

where $A_0 = 10^{-3}$, $R = 0.3$ and $r = \sqrt{x^2 + y^2}$. This test case is very challenging because of a singular point in the magnetic field B at the center of the computational domain. This setting places the flow in a low
<figure><img src="image_9.png" /><figcaption>Fig. 8 Blast problem and rotor problem. left: time evolution of the infinity norm of the divergence of the magnetic field $|\nabla \cdot \mathbf{B}|_{\infty}$ for 1000 time steps. Right: Time evolution of the ratio between the material time step of the SIFV-EB scheme and the magneto-sonic time step of a fully explicit finite volume scheme</figcaption></figure>

acoustic Mach number regime $M_c \approx 6 \cdot 10^{-3}$ and high Alfvén Mach number regime $M_b \approx 7.9 \cdot 10^3$. To make the test more challenging, we perform four different simulations with a progressively increasing value of $A_0$, reaching a configuration with an Alfvén Mach number $M_b \approx 7.9$ corresponding to $A_0 = 1$. The final time is chosen to be $t_f = 1$, and the results are depicted in Figure 10. The shape of the magnitude of the magnetic field remains qualitatively unchanged despite the difference in the Alfvén Mach number, hence demonstrating that the new SIFV-EB scheme is independent of the magnetic scales. Indeed, the time step for all the simulations is the same and is only determined by the convective velocity $(u, v) = (2, 1)$ which transports the loop diagonally through the domain.

In Figure 11 we show the errors related to the divergence of the magnetic field for all four simulations, demonstrating that the discrete divergence-free condition is preserved up to machine accuracy. The infinity norm of the divergence of $\mathbf{B}$ is considered over the whole domain and at each time level.

## 4.7 Magnetized Kelvin Helmholtz instability

We consider the setup for the Kelvin-Helmholtz instability recently proposed in [47] for ideal MHD. The computational domain is the rectangular box $\Omega = [0; 2] \times [-0.5; 0.5]$ with periodic boundaries. It is paved with a computational grid made of $N_x \times N_y = 512 \times 256$ control volumes. A uniform density $\rho = \gamma$ with $\gamma = 1.4$ and a uniform pressure $p = 1$ is initially assigned to the fluid. The velocity and the magnetic field are prescribed as a function of the parameter $M_x$, which represents the maximum acoustic Mach number of the horizontal flow. The velocity field is given by

$$u = M_x (1 - 2\eta(x)), \quad v = 0.1 M_x \sin(2\pi x), \qquad (73)$$

with the perturbation $\eta(x)$

$$\eta(x) = \begin{cases} 0.5 (1 + \sin(16\pi(y + 0.25))) & \text{if } y \ge -\frac{9}{32} \text{ and } y < -\frac{7}{32}, \\ 1 & \text{if } y \ge -\frac{7}{32} \text{ and } y < \frac{7}{32}, \\ 0.5 (1 - \sin(16\pi(y - 0.25))) & \text{if } y \ge \frac{7}{32} \text{ and } y < \frac{9}{32}, \\ 0 & \text{otherwise.} \end{cases} \qquad (74)$$
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Orszag-Tang vortex at time $t = 0.5$, $t = 2$, $t = 3$ and $t = 4$ (from top left to bottom right). Results are given for the pressure $p$ with 31 equidistant contour lines in the interval $[0.5, 5.0]$.</figcaption></figure>

The magnetic field at the initial time is uniform and horizontal, i.e. $B_x = 0.1 M_x$, and the associated $z$-component of the magnetic potential is given by $A_0 = B_x y$. This configuration at time $t = 0$ leads to a unit adiabatic sound speed and to a minimum Alfvén mach number $M_{b,min} = 11.82$. The simulation is run with different values of $M_x$ until the final time $t_f = t_{max}/6$ with $t_{max} = 4.8/M_x$. The flow develops a smooth and resolved interface across the shear flow, hence yielding to convergent results in the early stages of the dynamics. Figure 12 plots the acoustic Mach number at the end of the simulation. It is rescaled with the prescribed value of $M_x$. Note that all the results are almost identical since excessive numerical dissipation is avoided thanks to the central implicit discretization of the pressure and the magnetic terms. Fully explicit schemes would not preserve the vortical structures as the value of $M_x$ gets smaller, as noticed in [47]. Furthermore, our SIFV-EB scheme can run these test cases with a time step that is up to 10,000 times larger compared to a fully explicit method, as shown in Figure 13.
<figure><img src="image_11.png" /><figcaption>Fig. 10 Field loop advection problem at time $t_f = 1$. Magnitude of the magnetic field with $A_0 = \{10^{-3}, 10^{-2}, 10^{-1}, 10^0\}$ (from top to bottom).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Field loop advection problem with $A_0 = \{10^{-3}, 10^{-2}, 10^{-1}, 10^0\}$. Evolution of the infinity norm of the divergence of the magnetic field $|\nabla \cdot \mathbf{B}|_\infty$ for 1000 time steps of each simulation.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12** Kelvin-Helmholtz instability test at time $t/t_{max} = 1/6$ for $M_x = \{10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}\}$. Distribution of the acoustic Mach number rescaled by the corresponding value of $M_x (M_c/M_x)$ with 31 equidistant contour lines in the interval $[0, 1.6]$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13 Kelvin-Helmholtz instability test. Time evolution of the ratio between the material time step of the SIFV-EB scheme and the magneto-sonic time step of a fully explicit finite volume scheme for different values of $M_x$.</figcaption></figure>

## 4.8 3D cloud-shock interaction problem

As last test case we consider the 3D cloud-shock interaction problem. The computational domain is given by $\Omega = [0; 1]^3$ with Dirichlet boundary conditions, and we employ a computational mesh with a total number of $N_x = N_y = N_z = 180^3$ cells, corresponding to approximately $16.5 \cdot 10^6$ degrees of freedom. The initial condition is taken from [44] where the initial shock front is located at $x_0 = 0.0$ with

$$ (\rho, u, v, w, p, B_x, B_y, B_z) = \begin{cases} (3.86859, 11.2536, 0, 0, 167.345, 0, 2.1826182, -2.1826182) & \text{if } x < x_0 \\ (1, 0, 0, 0, 1, 0, 0.56418958, 0.56418958) & \text{if } x \ge x_0 \end{cases} . \quad (75) $$

The associated magnetic vector potential writes

$$ \mathbf{A} = \begin{cases} ( 2.1826182y, 0, -2.1826182(x - 0.05) )^T & \text{if } x < x_0 \\ ( -0.56418958y, 0, 0.56418956(x - 0.05) )^T & \text{if } x \ge x_0 \end{cases} . \quad (76) $$

This shock front travels towards a spherical cloud of density $\rho = 10$ situated at $\mathbf{x} = (0.25, 0.5, 0.5)$ with radius $r = 0.15$. The final time is set to $t_f = 0.06$. In Figure 14, we compare the first against the second order SIFV-EB scheme. The density profiles along the planes cutting through the center of the initial cloud position are displayed. The first snapshot is taken shortly after the shock front hits the density cloud, the second one after the shock front passes through the cloud, and the third one shortly before the shock front leaves the computational domain. One can observe that the second-order scheme greatly improves the resolution of the flow structures, which are qualitatively in very good agreement with the results presented in [44]. This test case demonstrates the capability of the SIFV-EB scheme to handle
genuinely 3D magnetized flows at reasonable compute times due to the linear implicit structure of our
novel SIFV-EB scheme.

5 Conclusions

In this work, we have introduced a new divergence-free semi-implicit finite volume method for the simulation of ideal magnetized flows in multiple space dimensions. To tackle the time multi-scale nature of the governing equations, a flux splitting technique is adopted, which allows to separate the convective terms, the hydrodynamic pressure contribution and its magnetic counterpart. More specifically, the nonlinear convective contribution is discretized explicitly with a classical finite volume scheme, while we use a finite difference semi-implicit strategy to compute the hydrodynamic pressure and the magnetic fluxes, that ultimately yields two linear systems to be solved at each time step of the algorithm. We underline that these two systems are uncoupled, thus firstly we solve for the unknown magnetic field, and secondly we determine the unknown pressure contribution. This approach makes the scheme very efficient, especially in the low acoustic or Alfvén regimes, since the time step is restricted by a milder CFL stability condition that is only based on the material waves. Furthermore, no numerical dissipation is added in the discretization of the implicit pressure fluxes, hence the scheme does not pollute the solution in the low acoustic Mach number regime, contrarily to standard explicit finite volume methods. To satisfy the solenoidal property of the magnetic field, the equations are written in terms of the magnetic vector potential in a multi-dimensional setting, and a set of finite difference operators are specifically designed to mimic the algebraic property of the zero div-curl at the discrete level, thus preserving the structure of the underlying model. The entire spatial discretization is cell-centered, and second-order of accuracy is achieved by means of a TVD reconstruction procedure in the finite volume scheme. The class of semi-implicit IMEX Runge-Kutta schemes is employed to reach second-order accuracy also in time. The new schemes have been thoroughly tested in terms of accuracy and robustness against a wide set of test cases, addressing different acoustic and Alfvén Mach number regimes, demonstrating the capability of our novel methods to span a wide range of physical settings in ideal MHD.

In the future, we plan to carry out a theoretical and numerical study of the limit models in the
low Alfvén Mach number, which is still rather unexplored. The extension of the presented approach
to unstructured meshes is also foreseen, with particular care in the development of divergence-free
operators that can guarantee the solenoidal property of the magnetic field. Finally, future research
lines will consider the application to viscous and resistive MHD equations in the context of Tokamak
simulations.

Acknowledgments

WB received financial support by Fondazione Cariplo and Fondazione CDP (Italy) under the project No. 2022-1895 and by the Italian Ministry of University and Research (MUR) with the PRIN Project 2022 No. 2022N9BM3N. WB is member of the GNCS-INdAM (*Istituto Nazionale di Alta Matematica*) group.

Declarations

*Funding.* WB received financial support by Fondazione Cariplo and Fondazione CDP (Italy) under the project No. 2022-1895 and by the Italian Ministry of University (MUR) with the PRIN Project 2022 No. 2022N9BM3N.

*Conflicts of interest.* The authors declare that they have no conflict of interest.
<figure><img src="image_15.png" /><figcaption>**Fig. 14** 3D cloud-shock interaction problem at time $t = 0.02$ (top), $t = 0.04$ (middle) and $t = 0.06$ (bottom). Distribution of the logarithm of density obtained with first (left) and second (right) order SIFV-EB schemes. Slices along the planes $x = 0.25$, $y = 0.5$ and $z = 0.5$.</figcaption></figure>
*Availability of data and material.* Data and material are available upon reasonable request addressed to the corresponding author.

*Code availability.* The code is written in Fortran programming language and is available upon reasonable request addressed to the corresponding author.

*Ethics approval.* Not applicable.

*Consent to participate.* Not applicable.

*Consent for publication.* Not applicable.

## Data availability

The datasets generated during the current study are available from the corresponding author upon reasonable request.

## References

1. C. Altmann, T. Belat, M. Gutnic, P. Helluy, H. Mathis, E. Sonnendruecker, W. Angulo, and J-M. Herard. A local time-stepping discontinuous Galerkin algorithm for the MHD system. In *ESAIM: proceedings*, volume 28, pages 33–54, 2009.

2. U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. *Appl. Numer. Math.*, 25:151–167, 1982.

3. S. Avgerinos, F. Bernard, A. Iollo, and G. Russo. Linearly implicit all Mach number shock capturing schemes for the Euler equations. *J. Comput. Phys.*, 393:278–312, 2019.

4. D. Balsara. Second-Order Accurate Schemes for Magnetohydrodynamics with Divergence-Free Reconstruction. *Astrophys. J. Suppl. Ser.*, 151:149–184, 2004.

5. D. S. Balsara. Multidimensional HLLE Riemann solver: application to Euler and magnetohydrodynamic flows. *J. Comput. Phys.*, 229(6):1970–1993, 2010.

6. D. S. Balsara and M. Dumbser. Divergence-free mhd on unstructured meshes using high order finite volume schemes based on multidimensional riemann solvers. *J. Comput. Phys.*, 299:687–715, 2015.

7. D. S. Balsara and 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):270–292, 1999.

8. W. Barsukow. Truly multi-dimensional all-speed schemes for the Euler equations on Cartesian grids. *J. Comput. Phys.*, 435:110216, 2021.

9. W. Barsukow. All-speed numerical methods for the Euler equations via a sequential explicit time integration. *J. Sci. Comput.*, 95(2):53, 2023.

10. C. Birke, W. Boscheri, and C. Klingenberg. A well-balanced semi-implicit IMEX finite volume scheme for ideal Magnetohydrodynamics at all Mach numbers. *J. Sci. Comput.*, 98(2):34, 2024.

11. C. Birke, C. Chalons, and C. Klingenberg. A low mach two-speed relaxation scheme for the compressible euler equations with gravity. *Commun. Math. Sci.*, 21, 2023.

12. S. Boscarino, F. Filbet, and G. Russo. High Order Semi-implicit Schemes for Time Dependent Partial Differential Equations. *J. Sci. Comput.*, 68:975–1001, 2016.

13. S. Boscarino and L. Pareschi. On the asymptotic properties of IMEX Runge-Kutta schemes for hyperbolic balance laws. *J. Comput. Appl. Math.*, 316:60–73, 2017.

14. S. Boscarino, L. Pareschi, and G. Russo. A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation. *SIAM Journal of Numerical Analysis*, 55(4):2085–2109, 2017.

15. S. Boscarino, J.-M. Qiu, G. Russo, and T. Xiong. A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system. *Journal Of Computational Physics*, 392:594–618, 2019.

16. S. Boscarino and G. Russo. On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolic systems with relaxation. *SIAM J. Sci. Comput.*, 31:1926–1945, 2009.

17. W. Boscheri, G. Dimarco, and L. Pareschi. Locally Structure-Preserving div-curl operators for high order Discontinuous Galerkin schemes. *J. Comput. Phys.*, 486:112130, 2023.

18. W. Boscheri, M. Dumbser, and D. S. Balsara. High-order ADER-WENO ALE schemes on unstructured triangular meshes—application of several node solvers to hydrodynamics and magnetohydrodynamics. *Int. J. Numer. Methods Fluids*, 76(10):737–778, 2014.
19. W. Boscheri, M. Dumbser, and O. Zanotti. High order cell-centered lagrangian-type finite volume schemes with time-accurate local time stepping on unstructured triangular meshes. *J. Comput. Phys.*, 291:120–150, 2015.

20. W. Boscheri and L. Pareschi. High order pressure-based semi-implicit IMEX schemes for the 3D Navier-Stokes equations at all Mach numbers. *Journal Of Computational Physics*, 434:110206, 2021.

21. W. Boscheri and M. Tavelli. High order semi-implicit schemes for viscous compressible flows in 3d. *Appl. Math. Comput.*, 434:127457, 2022.

22. F. Bouchut, C. Klingenberg, and K. Waagan. A multiwave approximate Riemann solver for ideal MHD based on relaxation. I: theoretical framework. *Numer. Math.*, 108:7–42, 2007.

23. J.U. Brackbill and D.C. Barnes. The effect of nonzero $\nabla \cdot \mathbf{B}$ on the numerical solution of the magnetohydrodynamic equations. *J. Comput. Phys.*, 35(3):426–430, 1980.

24. L. Brugnano and V. Casulli. Iterative solution of piecewise linear systems. *SIAM J. Sci. Comput.*, 30:463–472, 2007.

25. W. Chen, K. Wu, and T. Xiong. High order asymptotic preserving finite difference WENO schemes with constrained transport for MHD equations in all sonic Mach numbers. *Journal Of Computational Physics*, 488:112240, 2023.

26. F. Cordier, P. Degond, and A. Kumbaro. An Asymptotic-Preserving all-speed scheme for the Euler and Navier-Stokes equations. *J. Comput. Phys.*, 231:5685–5704, 2012.

27. R. B. Dahlburg and J. M. Picone. Evolution of the Orszag–Tang vortex system in a compressible medium. I. Initial average subsonic flow. *Phys. Fluids B*, 1(11):2153–2171, 1989.

28. A. Dedner, F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, and M. Wesenberg. Hyperbolic divergence cleaning for the MHD equations. *J. Comput. Phys.*, 175(2):645–673, 2002.

29. S. Dellacherie. Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number. *J. Comput. Phys.*, 229:978–1016, 2010.

30. M. Dumbser and D. S. Balsara. A new efficient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems. *J. Comput. Phys.*, 304:275–319, 2016.

31. M. Dumbser, D.S. Balsara, M. Tavelli, and F. Fambri. A divergence-free semi-implicit finite volume scheme for ideal, viscous, and resistive magnetohydrodynamics. *Int. J. Numer. Methods Fluids*, 89:16–42, 2019.

32. M. Dumbser and V. Casulli. A conservative, weakly nonlinear semi-implicit finite volume scheme for the compressible Navier-Stokes equations with general equation of state. *Appl. Math. Comput.*, 272:479–497, 2016.

33. M. Dumbser and E. F. Toro. On universal Osher-type schemes for general nonlinear hyperbolic conservation laws. *Commun. Comput. Phys.*, 10(3):635–671, 2011.

34. S. A. E. G. Falle. Rarefaction Shocks, Shock Errors, and Low Order of Accuracy in ZEUS. *ApJ*, 577(2):L123–L126, October 2002.

35. S. A. E. G. Falle and S. S. Komissarov. On the inadmissibility of non-evolutionary shocks. *J. Plasma Phys.*, 65(1):29–58, 2001.

36. F. Fambri. A novel structure preserving semi-implicit finite volume method for viscous and resistive magnetohydrodynamics. *Int. J. Numer. Methods Fluids*, 93:3447–3489, 2021.

37. F. Fambri, E. Zampa, S. Busto, L. Río-Martín, F. Hindenlang, E. Sonnendrücker, and M. Dumbser. A well-balanced and exactly divergence-free staggered semi-implicit hybrid finite volume / finite element scheme for the incompressible MHD equations. *J. Comput. Phys.*, 493:112493, 2023.

38. M. Fey and M. Torrilhon. A constrained transport upwind scheme for divergence-free advection. In *Hyperbolic Problems: Theory, Numerics, Applications: Proceedings of the Ninth International Conference on Hyperbolic Problems held in CalTech, Pasadena, March 25–29, 2002*, pages 529–538, 2003.

39. T.A. Gardiner and J.M. Stone. An unsplit Godunov method for ideal MHD via constrained transport. *J. Comput. Phys.*, 205(2):509–539, 2005.

40. S. Godunov. Finite difference methods for the computation of discontinuous solutions of the equations of fluid dynamics. *Mat. Sb.*, 47:271–306, 1959.

41. H. Guillard and A. Murrone. On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes. *Comput. & Fluids*, 33:655–675, 2004.

42. A. Harten, P. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Rev.*, 25:35–61, 1983.

43. C. Helzel, J. A. Rossmanith, and B. Taetz. An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations. *J. Comput. Phys.*, 230(10):3803–3829, 2011.

44. C. Helzel, J. A. Rossmanith, and B. Taetz. A high-order unstaggered constrained-transport method for the three-dimensional ideal magnetohydrodynamic equations based on the method of lines. *SIAM J. Sci. Comput.*, 35(2):A623–A651, 2013.

45. J.M. Hyman and M. Shashkov. Natural discretizations for the divergence, gradient, and curl on logically rectangular grids. *Computers and Mathematics with Applications*, 33:81–104, 1997.

46. R. Klein. Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow. *J. Comp. Phys.*, 121:213–237, 1995.

47. G. Leidi, C. Birke, R. Andrassy, J. Higl, P.V.F. Edelmann, G. Wiest, C. Klingenberg, and F.K.Röpke. A finite-volume scheme for modeling compressible magnetohydrodynamic flows at low Mach numbers in stellar interiors. *A&A*, 668:A143, 2022.

48. K. Lipnikov, G. Manzini, and M. Shashkov. Mimetic finite difference method. *J. Comput. Phys.*, 257:1163–1227, 2014.
49. M. Lukáčová-Medvid'ová, I. Peshkov, and A. Thomann. An implicit-explicit solver for a two-fluid single-temperature model. *J. Comput. Phys.*, 498:112696, 2024.

50. G. Margolin, M. Shashkov, and P.K. Smolarkiewicz. A discrete operator calculus for finite difference approximations. *Comput. Methods Appl. Mech. Eng.*, 187:365-383, 2000.

51. C-D. Munz, P. Omnes, R. Schneider, E. Sonnendrücker, and U. Voss. Divergence correction techniques for Maxwell solvers based on a hyperbolic model. *J. Comput. Phys.*, 161(2):484-511, 2000.

52. S. Osher and F. Solomon. A partially implicit method for large stiff systems of ODE's with only few equations introducing small time-constants. *SIAM Journal of Numerical Analysis*, 13:645-663, 1976.

53. L. Pareschi and G. Russo. Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation. *J. Sci. Comput.*, 25:129-155, 2005.

54. J. M. Picone and R. B. Dahlburg. Evolution of the Orszag-Tang vortex system in a compressible medium. II. Supersonic flow. *Phys. Fluids B*, 3(1):29-44, 1991.

55. P. L. Roe and D. S. Balsara. Notes on the eigensystem of magnetohydrodynamics. *SIAM J. Appl. Math.*, 56(1):57-67, 1996.

56. C-M. Tang S. A. Orszag. Small-scale structure of two-dimensional magnetohydrodynamic turbulence. *J. Fluid Mech.*, 90(1):129-143, 1979.

57. Y. Saad and M. Schultz. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. *SIAM J. Sci. Stat. Comput.*, 7:856-869, 1986.

58. E.F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics: a Practical Introduction*. Springer, 2009.

59. E.F. Toro and M.E. Vázquez-Cendón. Flux splitting schemes for the Euler equations. *Comput. & Fluids*, 70:1-12, 2012.

60. M. Torrilhon. Non-uniform convergence of finite volume schemes for Riemann problems of ideal magnetohydrodynamics. *J. Comput. Phys.*, 192(1):73-94, 2003.

61. M. Viallet, I. Baraffe, and R. Walder. Towards a new generation of multi-dimensional stellar evolution models: development of an implicit hydrodynamic code. *A&A*, 531:A86, 2011.

62. K. Yee. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. *IEEE Transactions on antennas and propagation*, 14(3):302-307, 1966.