A well-balanced semi-implicit IMEX finite volume scheme for ideal Magnetohydrodynamics at all Mach numbers

We propose a second-order accurate semi-implicit and well-balanced finite volume scheme for the equations of ideal magnetohydrodynamics (MHD) including gravitational source terms. The scheme treats all terms associated with the acoustic pressure implicitly while keeping the remaining terms part of the explicit sub-system. This semi-implicit approach makes the method particularly well suited for problems in the low Mach regime. We combine the semi-implicit scheme with the deviation well-balancing technique and prove that it maintains equilibrium solutions for the magnetohydrostatic case up to rounding errors. In order to preserve the divergence-free property of the magnetic field enforced by the solenoidal constraint, we incorporate a constrained transport method in the semi-implicit framework. Second order of accuracy is achieved by means of a standard spatial reconstruction technique with total variation diminishing (TVD) property, and by an asymptotic preserving (AP) time stepping algorithm built upon the implicit-explicit (IMEX) Runge-Kutta time integrators. Numerical tests in the low Mach regime and near magnetohydrostatic equilibria support the low Mach and well-balanced properties of the numerical method.


Introduction
The equations of magnetohydrodynamics (MHD) describe electrically conducting fluids in the presence of a magnetic field.They can be used to model a wide range of physical phenomena, including astrophysical systems such as stars and galaxies, plasma physics, and fusion reactors.Their simplest form is given by the ideal MHD equations, in which the influence of fluid viscosity is neglected.This model consists of a system of nonlinear hyperbolic partial differential equations (PDEs) that involve the conservation of mass, momentum, and total energy, along with Faraday's law for the magnetic field.In certain applications, such as flows in the deep convective layers of stars, the acoustic Mach number, which is the ratio of the fluid speed to the speed of sound, can become very small.Small Mach numbers pose a challenge for explicit finite volume methods, which are typically used to simulate compressible flows, since in this case the flow approaches the incompressible regime and the acoustic waves propagate much faster than the fluid motion.This gives rise to the following problems for standard finite volume methods: the numerical dissipation increases significantly and the time step condition becomes more restrictive.Basically, there are two strategies in the literature to tackle these challenges.
1.One approach is to rescale the dissipation term in the numerical flux to make it independent of the Mach number and combine this with a fully implicit time integration to circumvent the restrictive time step condition [35,26].The drawback of this strategy is that solving large nonlinear systems of equations is required due to the fully implicit method, which can be computationally expensive.2. Another approach is to split the PDE system into a convective part that is solved explicitly and an acoustic part that is solved implicitly [16,8,14].This makes the CFL condition only coupled to the convective sub-system and therefore independent of the Mach number.At the same time, the implicit part becomes smaller and easier to invert, making it less computationally expensive to solve.Furthermore, no numerical dissipation is embedded in the implicit solver.
For our scheme, we pursue the latter approach and treat all terms associated with the acoustic pressure implicitly.All remaining terms are part of the explicit sub-system.For this purpose, a time integration inspired by the class of implicit-explicit (IMEX) schemes is used [5,8].Concerning the development of semi-implicit numerical schemes applied to compressible magnetized plasma flows governed by the ideal MHD equations, only very little work has been done so far.Indeed, semiimplicit schemes for MHD flows have been firstly designed in [27,1,22], however they are concerned with the incompressible case.In [31] the anelastic case of the MHD equations is discussed, while in [16] a semi-implicit finite difference scheme is propsed, inspired by the flux splitting approach originally forwarded for compressible gas dynamics [17].An implicit treatment of the magnetic field has been recently devised in [18], and in [14] a semi-implicit IMEX numerical method for the ideal MHD has been proposed with a finite difference spatial discretization.Our work is strictly related to what presented in [14], however there are some important differences: i) we use a finite volume discretization for the convective terms; ii) no numerical dissipation is added in the implicit part even in the case of shock waves; iii) we also include gravitational source terms.Regarding the discretisation in the case of low Mach numbers, this has no effect and the source terms can easily be assigned to the explicit part.More challenging is its impact on magnetohydrostatic solutions.These are stationary solutions, where the velocity of the fluid is zero and the magnetic field and pressure are in balance with the gravitational force.If one wants to resolve small perturbations of this equilibrium on a coarse grid, the method should be able to maintain the equilibrium exactly.Usually, this is prevented by different discretizations of the flux and the source term.Numerical schemes that are still able to maintain this balance are called well-balanced.There are different strategies in the literature to construct well-balanced schemes, e.g.hydrostatic reconstructions [23,24,34], path conservative methods [28,12] or special approximate Riemann solvers based on relaxation systems [15,32,4].In this work, we combine the semi-implicit scheme with the method of deviation well-balancing [3].Instead of evolving the actual solution, the deviation from an a priori known equilibrium solution is evolved.Through special discretization, the numerical fluxes and source terms cancel out in the case of equilibrium, so that the discrete solution is preserved up to rounding errors.A similar approach has been recently forwarded in [20] in the context of general relativity.
In simulations of the MHD equations, it is important to consider the solenoidal constraint, which states that the magnetic field is divergence-free.Numerically, we achieve this structure preserving property by incorporating a constrained transport (CT) method [21] into the semi-implicit scheme.
The CT method corrects the magnetic field directly after the explicit step, as the induction equations are entirely assigned to the explicit sub-system of the splitting.By correcting the magnetic field, a specific discrete definition of divergence is kept within machine precision, thereby increasing the stability of the overall scheme.
The paper is structured as follows.In Sect.2, we provide an overview of the MHD equations and their associated properties and propose a splitting into convective and pressure sub-systems.On the basis of this splitting, we develop in Sect. 3 a second order accurate semi-implicit and well-balanced numerical scheme.The numerical results of this scheme for different test problems are presented in Sect.4, and finally, Sect. 5 concludes the article providing an outlook to future developments.

Governing equations
Let us consider a three-dimensional computational domain Ω( , , ) ∈ R 3 and let the time coordinate be denoted with ∈ R + 0 .The ideal MHD equations with a gravitational source term constitute a system of balance laws of the form In this system, denotes the density, v = ( , , ) the velocity field, the pressure and B = ( , , ) the magnetic field.The total energy is defined as the sum of internal ( ), kinetic ( ) and magnetic energy ( ), i.e.
The internal energy is computed from the ideal gas equation of state (EOS) with the ratio of specific heats given by .The vector g = ( , , ) in the source term on the right hand side of (1) denotes the time-independent gravitational acceleration.Using these defintions we can write the MHD equations in the equivalent form with the state vector q, the flux in -direction F(q) and the source S(q) that explicitly write The fluxes G(q) and H(q) can be expressed in similar forms.

Eigenstructure of the MHD system
In order to analyze the eigenstructure of the MHD system, we consider its one-dimensional homogeneous version.The resulting system is hyperbolic since the eigenvalues ={1,...,8} of the associated Jacobian matrix A = F/ q with = are The Alfvén wave speed ( ) and the slow ( ) and fast ( ) magnetosonic wave speeds herein are defined by where the adiabatic sound speed = / is computed from the ideal gas equation of state.Furthermore, we use the abbreviation 2 = B 2 /(4 ).

Low Mach number regime
The hydrodynamics behavior of the fluid can be studied by considering the Mach number M = |v|/ .In the low Mach number limit, the sound speed becomes very high compared to the fluid velocity, hence the terms related to the pressure are dominant.Consequently, larger values of the fast and slow magnetosonic wave speeds are retrieved, and fully explicit numerical methods suffer from both an excessive amount of numerical viscosity, which is proportional to the eigenvalues, and a drastic reduction of the admissible time step Δ to ensure stability under a classical CFL condition of the type with (Δ , Δ , Δ ) denoting the characteristic mesh spacings in the three spatial directions and the CFL ≤ 1 being the CFL number.Therefore, we propose to discretize the pressure gradient in the momentum equation and the enthalpy term in the energy equation implicitly, while keeping an explicit discretization for the nonlinear convective fluxes and the terms related to the magnetic field.To that aim, let the fluxes in -direction be split into a convective-type flux F (q) and a pressure-type flux F (q), that is We obtain two sub-systems with the following eigenvalues.
-Convective sub-system: -Pressure sub-system: The fluxes G and H can be split in a similar way with analogous sub-systems and eigenvalues.In the above splitting we have added the source term to the convective sub-system, since this has no effect on the eigenvalues and thus it does not pose a numerical problem in the low Mach limit.It is clear that, by taking the pressure sub-system implicitly, the maximum admissible time step of the scheme becomes hence making the scheme particularly well suited for low Mach number flows (M ≪ 1) where the pressure terms are dominant.On the other hand, for strongly convected flows with shocks, the convective eigenvalues lead the computation of the time step granting stability.

Divergence constraint
The induction equations, which describe the evolution of the magnetic field in (1), can be rewritten in the equivalent form where E = −v × B denotes the electric field.Applying the divergence operator to (12) yields This implies that the divergence of the magnetic field is always equal to zero if it is zero for the initial data.Therefore, the MHD equations in multiple spatial dimensions naturally contain the condition which is a physically meaningful condition, as it expresses that there should be no magnetic monopoles.
At the discrete level, on the other hand, the divergence of the numerically computed magnetic field can increase in time.This can lead to severe stability problems [11].Therefore, particular care must be taken in the discretisation to ensure that the divergence in the numerical solution remains sufficiently small [2,7].In this work, we rely on a constrained transport method that keeps the divergence at the level of machine accuracy [21].Sect.3.5 explains how this method works and how it is incorporated into the overall method.

Magnetohydrostatic solutions
When considering the MHD equations with gravitational source term, magnethydrostatic solutions can play an important role.These solutions fulfil the following two conditions As a consequence, for magnetohydrostatic solutions all fluxes except the momentum flux drop out.
The momentum flux, however, fulfils an equilibrium with the source term, so that magnetohydrostatic solutions are stationary.Different ways of discretising the hyperbolic flux and the graviatational source term prevent in general the exact conservation of discrete stationary magnetohydrostatic solutions of the form (15).In order to ensure that the numerical method described in this work is nevertheless exact (up to machine precision) on this type of stationary solutions, we resort to the idea of deviation well-balancing [3].The key idea of this approach is to evolve the deviation from an a priori known equilibrium solution q instead of the actual solution q.Therefore, we do not directly discretize equation ( 3), but the following equation for the deviation Δq = q − q: In the equilibrium case (q = q), the spatial discretizations of ( 16) cancel so that the equation reduces to Δq = 0 and the solution remains stationary.

Numerical scheme
For reasons of clarity and comprehensibility we restrict ourselves in the following to one spatial dimension before returning to three spatial dimensions in Section 3.4.The computational domain Ω = [ ; ] is discretized using a total number of equidistant cells of volume Δ = ( − )/ .The cell centers are indicated with and the cell interfaces are referred to with +1/2 .The time coordinate is bounded in the interval ∈ [0; ], and the final time is reached performing a sequence of time steps Δ = +1 − that are computed according to the one-dimensional version of the CFL stability condition (11).

First order semi-discrete scheme in time
The method constructed in this paper should be able to preserve magnetohydrostatic solutions of the form (15). Therefore, we assume to know an equilibrium solution q a priori and discretize equation (16).The splitting ( 8) is therefore applied to equation (16).To simplify the notation, we write ¯ = Δ + ˜ .For the convective sub-system, which is treated explicitly, we then get: Written out for the single components, this results in In the above equations, we have included terms containing ˜ , ˜ or ˜ .In the following, however, we will omit these terms, since the equilibrium velocities are always zero.The above definitions are then employed to obtain a first order semi-implicit time discretization [5,8,10] of the MHD equations ( 3), which writes Here, we use the known enthalpy ℎ at time level in the energy equation to avoid nonlinear terms in the implicit part, differently form the schemes proposed in [16,18].Additionally, we have made use of the fact that in the implicit fluxes we only work with point values at the cell centers, i.e. we use an implicit finite difference discretization, thus we can write Therefore, the implicit fluxes are simplified by According to the definition given in the energy equation ( 2), we split the deviation of the total energy at the new time level as follows The deviation of the kinetic energy therein is set to where a semi-implicit strategy is adopted for the term ( ) +1 .Indeed, the split of the momentum contribution into an explicit and an implicit part is again done in order to avoid nonlinear implicit terms.The density +1 and the deviation in the magnetic energy Δ +1 are known because both continuity and induction equations are fully explicit.With the help of the ideal gas law, the deviation in the internal energy can be expressed in terms of the pressure In order to derive a preliminary discretization of the total energy equation we transform (19b) into and insert the result into the energy equation (19e).This leads to an elliptic equation for the pressure with the known right-hand-side given by The pressure equation ( 25) constitutes a linear system for the scalar unknown +1 that is solved using the iterative GMRES solver [30] up to a prescribed tolerance (we typically set tol = 10 −12 ).Differently from [16,18], this approach does not need any fixed point method thanks to the semiimplicit splitting of the enthalpy flux and the kinetic energy in the energy equation.Once the new pressure is known, the deviation in the momentum (Δ ) +1 is updated with (19b), and then the deviation in the total energy is updated using (19e).Notice that the scheme is written in flux form, therefore it is locally and globally conservative.

Discrete spatial operators
The convective sub-system ( 17) is discretized by a Godunov-type finite volume method which writes In this notation, the hat symbol indicates a numerical flux or source term.For the numerical flux, we decide to use a simple and robust Rusanov-type flux of the form where the numerical dissipation = max |( ) +1/2 |, |( ) +1/2 | only considers the convective eigenvalues and is therefore independent of the Mach number.The numerical flux is computed in the states q , +1/2 = Δq , where the superscripts , denote the left and right extrapolated data at the interface.It is essential for the well-balancing that only the deviation Δq is reconstructed, while the equilibrium solution is evaluated at the cell interface.The third term in the right hand side of ( 27) simply computes the physical fluxes based on the equilibrium solution at the cell interface, that is F (q , +1/2 ).For the source term, we substitute the volume-averaged quantity with its cell centered value, which is accurate up to second order: The implicit terms appearing in the pressure sub-system (10) are approximated by means of finite difference operators with no numerical dissipation [8]: 3.3 Second order extension in time and space

Second order in time
The scheme achieves second order of accuracy in time by using the implicit-explicit Runge Kutta method LSDIRK(2,2,2) [29].The triplet (2,2,2) refers to the number of stages of the implicit part, the number of stages of the explicit part, and the order of the IMEX scheme, respectively.Let us denote the spatial operator in the MHD equations by H (q ( ), q ( )), where the first argument is discretized explicitly, while the second argument is discretized implicitly.According to the splitting (8), H is defined by At the beginning of each time step we initialize q = q = q and then compute for each stage its stage fluxes k in the following way: Notice that the only implicit evaluation lies in (33c), which indeed correspond to the solution of the elliptic equation on the pressure (25).Using these stage fluxes, the updated solution at the new time level is computed by The coefficients are given by a double Butcher tableau of the form with matrices ( ˆ , ) ∈ R × and vectors ( ˆ , , ˆ , ) ∈ R .In this notation, the hat symbol indicates the coefficients for the explicit scheme.The coefficients for LSDIRK2(2,2,2) are listed in Table 1.

Second order in space
The extension to second order in space in the explicit part is based on a piecewise linear reconstruction of the conservative variables.The values at the cell interface, which serve as initial data for the Riemann problems, are computed by evaluating the function at the cell boundaries −1/2 and +1/2 .The slopes are computed by applying a minmod limiter to the left and right slope, i.e.
We compute a separate slope for each of the conservative variables in q.At this point, it should be recalled that the reconstruction is only applied to the deviation Δq in (29), and not to the equilibrium solution q.With this procedure, the spatial discretization of the explicit part becomes second order.In the implicit part, no further changes are needed, since only point values are used and the operators in (31) are already second order.

Multi-dimensional extension
In three spatial dimensions we discretize the domain Ω = [ , ] × [ , ] × [ , ] by a Cartesian grid with × × cells, which have the uniform cell size Δ × Δ × Δ .The discretization of the explicit part is extended to three spatial dimensions by using an unsplit finite volume method according to [33], which writes The numerical fluxes F , Ĝ and Ĥ have the form of the Rusanov flux and are constructed as in (28).For the source term Ŝ we again use the cell centered value.The updated density Δ +1 , , and magnetic field ΔB +1 , , are equal to their explicit update, since in the splitting the complete flux of these components is explicit.The update of the momentum components, on the other hand, contains implicit parts and reads in the fully discrete and three-dimensional form as Implicit terms also appear in the update of the total energy: The pressure +1 , which is needed for the updates (39) and (40), can be determined by solving the following elliptic equation: , , , , , , , , − ˜ , , with the right-hand side , , , , , , (42)

Constrained transport method
When simulating the ideal MHD equations in multiple space dimensions, the divergence constraint ( 14) must also be taken into account in order to maintain stability and accuracy.Furthermore, the solenoidal property of the magnetic field represents a physical constraint.The scheme described up to this point does not do this in the sense that the divergence can increase significantly during the simulation.Therefore, an additional correction of the magnetic field is required after each time step.In order to retrieve at the discrete level a solenoidal magnetic field, we rely on a second order accurate CT method [21].As it is typical in staggered CT methods, this method is based on the integration of the induction equations ( 12) over the cell boundaries using Stokes theorem.This results in an update formula for the magnetic field at the face centers of the cell based on the electric field at the cell corners.In our case, we update the deviation at the face centers: The deviation of the magnetic field at -and -faces is updated by similar formulae.The discretizations of the electric field at the cell corners can be taken from [21] and transferred to the deviations without any problems, since the physical fluxes evaluated in the equilibrium solution are omitted due to the zero velocity.The cell averages of the magnetic field is then computed via the arithmetic mean from the values at the faces, which still maintains second order of accuracy: In the semi-implicit scheme, this whole procedure of constained transport is carried out directly after the explicit update.Thus, for solving the pressure equation ( 41) we already use the corrected and divergence-free magnetic field.Overall, by incorporating the CT method, the semi-implicit scheme keeps the following discrete definition of the divergence within the order of machine accuracy.This is checked and validated at each time step of the computation.

Well-balanced property
The following theorem states that the presented fully discrete scheme preserves magnetohydrostatic equilibria of the form (15) up to rounding errors.
Theorem 1 Let us assume that the numerical solution q is equal to the discrete equilibrium solution q, i.e.
Proof From the assumption (46) it follows that Δq = 0 at every grid point.Thus, the input of the Rusanov fluxes in the convective sub-system is reduced to the equilibrium solution at the interface, i.e. q +1/2, , = q +1/2, , = q +1/2, , , q , , +1/2 = q , , +1/2 = q , , +1/2 .(47c) It follows by the consistency of the Rusanov flux that the numerical fluxes are equal to the physical fluxes and these truncate away with the physical fluxes in (38).Likewise, the source terms are truncated away.Therefore, the contributions from the explicit part are zero: (48) Consequently, also the deviation in the magnetic energy Δ +1 , , is zero.The constrained transport step after the explicit update does not change anything, because the electric field at the corners in the update (43) is equal to zero due to the zero velocity in the equilibrium.Under these conditions, the elliptic equation for the pressure (41) simplifies to +1 , , , , + which admits the unique solution This means that the right-hand side of the momentum update (39) becomes zero and consequently the right-hand side of the energy update (40) also becomes zero.Thus the updated deviations for all components of the state vector remain zero, which means that the solution q does not change during one time step, i.e. q +1 , , = q , , = q , , ∀( , , ) ∈ {1, ..., } × {1, ..., } × {1, ..., }.
3.7 Summary of the scheme Since various methods are combined in our scheme, we provide an overview of the individual steps of the scheme at this point.Let us assume that we know the time-independent equilibrium solution q at every point in the computational domain Ω.Then a single time-step in the well-balanced semi-implicit method from time = to time = +1 consists of the following sub-steps: 1. Start with the values at the current time level: q , , , B , +1/2, , , B , , +1/2, , B , , , +1/2 .

Numerical results
In this section, the second order semi-implicit well-balanced method is applied to various numerical test problems for the Euler and ideal MHD equations.The time step of the method is calculated in each case based on the CFL condition (11) with CFL = 0.9.Partially, the results are compared with a non-well-balanced method.This is the same semi-implicit scheme, except for the fact that the equilibrium solution q is set to zero.All tests are performed on an uniform Cartesian grid.

Shock tube under gravitational field
The first test case is the standard Sod shock tube for the one-dimensional Euler equations, to which a gravitational source term with constant acceleration = −1 is added [13].The initial conditions in the domain Ω = [0, 1] are given by ( , v, ) = (1, 0, 1), < 0.5, (0.125, 0, 0.1), > 0.5, with solid wall boundary conditions.The magnetic field B is set to zero and the ratio of specific heats is = 1.4.The solution at final time = 0.2 is computed by the second order semi-implicit scheme on 100 cells.Since we have no flow around an equilibrium in this test, we set the equilibrium solution q to zero.The numerical solution is compared to a reference solution, which is computed by a fully explicit second order finite volume method on 20000 cells.The results in Fig. 1 show good agreement with the reference solution and are also consistent with solutions in the literature [13].This test demonstrates the capability of the semi-implicit scheme to deal with shocks, thus flows which are not in the low Mach regime.

General Euler Steady State in 1D
The following test investigates the method's behaviour in the presence of a hydrostatic equilibrium [15].On the domain Ω = [0, 1] we define the gravitational acceleration in -direction by = 2 cos(2 ) and use the initial values ( , v, ) = (3 + 2 sin(2 ), 0, 3 + 3 sin(2 ) − 0.5 cos(4 )) . (53) The magnetic field B is set to zero and = 1.4.Periodic boundary conditions are imposed on all sides.It can be easily checked that the values in (53) satisfy the equilibrium (15) and thus should be preserved until any final time .Here, we set = 1.In a first step, we perform this test without defining an equilibrium solution q.Tab. 2 lists the 1 error as well as the experimental order of convergence ( ) for different numbers of grid points .Clearly, in this case, the equilibrium is not preserved exactly, but only with a second order of accuracy.In this form, the method is therefore not well-balanced.In a second step we use the initial values as equilibrium solution q.The results in Tab. 3 show that the equilibrium is now maintained exactly.No iterations are needed to solve the linear system for the pressure and the residual is of the size of the machine accuracy. .The maximum number of the iterations and the corresponding residuals needed in the linear solver for the pressure system are also reported.

Isothermal atmosphere for Euler in 2D
After performing one-dimensional tests, the well-balanced property is now verified for the twodimensional scheme.For this purpose we set as initial values an isothermal equilibrium [13] of the form ( , v, ) = 1.21 −1.21( + ) , 0, −1.21( + ) .
The magnetic field B is set to zero and the gravitational acceleration is defined by g = (−1, −1, 0).The ratio of specific heats is = 1.4.Clearly, this initial data satisfies the equilibrium (15).We run this test until a final time = on a domain Ω = [0, 1] 2 with different numbers of cells ( , ).As boundary conditions we enforce the exact solution.The results in Tab. 4 show that indeed no iterations are needed to solve the pressure linear system (41) because the right hand side is perfectly balanced by the pressure terms.As a consequence, the residual has the order of the machine accuracy and the solution is preserved exactly throughout the simulation.).The maximum number of the iterations and the corresponding residuals needed in the linear solver for the pressure system are also reported.

Perturbation of an isothermal atmosphere for Euler in 2D
The aim of well-balanced methods is not only to maintain hydrostatic equilibria exactly, but in particular to resolve small perturbations of these equilibria.In this sense we now add a perturbation in the pressure to the equilibrium (54), i.e.
The parameter regulates the size of the perturbation.We perform the test on a 64 × 64-grid with the non-well-balanced method as well as with the well-balanced method.For the well-balanced scheme, the unperturbed initial data (54) is used for q.First, we run the test with a rather large perturbation ( = 10 −1 ).The top row in Fig. 2 shows the perturbation = ( ) − (0) at time = 0.15.For this perturbation, both methods provide a good resolution and perform qualitatively equally well.In a second simulation, we choose a smaller perturbation ( = 10 −10 ).In this case, the non-well-balanced method is no longer able to resolve the perturbation, while the well-balanced method resolves the small perturbation as good as the large one before.These results are consistent with the results from the literature [13,25,4].
This test is also carried out with a large ( = 10 −1 ) and a small ( = 10 −10 ) perturbation.Fig. 3 shows the results for the well-balanced and the non-well-balanced method at final time = 0.15.The behavior is similar to the one for the Euler equilibrium in Sect.4.4.The non-well-balanced the numerical scheme adds a certain amount of numerical dissipation to the approximate solution in comparison to the initial data, which dampens the local Mach number.The vortex, however, is well resolved in each case and there is no qualitative difference between the solutions for different M .This indicates that the numerical dissipation is independent of the Mach number thanks to the implicit discretisation of the acoustic terms in the flux splitting.This finding is underlined by the time evolution of the total kinetic energy.Since the domain is a closed setup due to the periodic boundary conditions, the total kinetic energy should be conserved over time in the low Mach limit.Therefore, the loss of kinetic energy is a good measure for the amount of numerical dissipation.Fig. 5 shows the time evolution of the kinetic energy for different M and grid resolutions = = .It turns out that the loss of kinetic energy only depends on the grid resolution and not on the Mach number.The results are consistent with those in [32].

Conclusions
In this work we have presented a well-balanced semi-implicit scheme for the ideal MHD equations with gravitational source terms.The method is based on splitting the equations into the convective part, which is discretized explicitly in time, and the acoustic part, which is discretized implicitly.An elliptic equation for the pressure is obtained, which is solved in the implicit part without adding any numerical dissipation.This makes the method particularly suitable for problems in the low Mach regime.Numerical results show indeed that the numerical dissipation is independent of the Mach number.At the same time, the explicit finite volume method for the nonlinear convective terms ensures stability in the presence of shocks.A constrained transport method is embedded in the numerical scheme to ensure the solenoidal property of the magnetic field at the discrete level.In order to maintain magnetohydrostatic equilibria exactly and to handle small perturbations of these equilibria even on coarse grids, the method utilizes the deviation well-balancing technique.The resulting scheme solves the equations for the deviation of the solution from an a priori known equilibrium solution.We have proven that, in the equilibrium, the pressure equation has only one solution, which is the equilibrium pressure itself.As a result, the solution is preserved exactly in the case of a magnetohydrostatic equilibrium.Thus, the scheme is well-balanced.Numerical tests demonstrate that arbitrary equilibria are preserved exactly and even very small perturbations are resolved accurately.Future developments will concern the extension of the method to unstructured meshes along the lines of [9,6], in order to address more challenging and physically relevant problems in plasma flows as tokamak simulations [19].Conflicts of interest.The authors declare that they have no conflict of interest.

Fig. 1
Fig. 1 Shock tube under gravitational field at time = 0.2.Comparison against the reference solution for density (top left), horizontal velocity (top right), pressure (bottom left) and internal energy (bottom right).

Fig. 5
Fig. 5 Time evolution of the total kinetic energy for different maximum Mach numbers M and grid resolutions = = .

Table 2
General Steady State for the non well-balanced scheme. 1 error at time = 1 and the experimental order of convergence for the density, horizontal velocity and pressure using different numbers of cells .

Table 3
General Steady State for the well-balanced scheme.The errors are measured in 1 norm at time = 1 for the density, horizontal velocity and pressure using different numbers of cells

Table 4
Isothermal atmosphere.The errors are measured in 1 norm at time = 1 for the density, velocity and pressure using different numbers of cells ( ,