High order asymptotic preserving finite difference WENO schemes with constrained transport for MHD equations in all sonic Mach numbers

In this paper, a high-order semi-implicit (SI) asymptotic preserving (AP) and divergence-free finite difference weighted essentially nonoscillatory (WENO) scheme is proposed for magnetohydrodynamic (MHD) equations. We consider the sonic Mach number $\varepsilon$ ranging from $0$ to $\mathcal{O}(1)$. High-order accuracy in time is obtained by SI implicit-explicit Runge-Kutta (IMEX-RK) time discretization. High-order accuracy in space is achieved by finite difference WENO schemes with characteristic-wise reconstructions. A constrained transport method is applied to maintain a discrete divergence-free condition. We formally prove that the scheme is AP. Asymptotic accuracy (AA) in the incompressible MHD limit is obtained if the implicit part of the SI IMEX-RK scheme is stiffly accurate. Numerical experiments are provided to validate the AP, AA, and divergence-free properties of our proposed approach. Besides, the scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime, while it is also a good incompressible solver with uniform large-time step conditions in the low sonic Mach limit.

For MHD equations, divergence free of the magnetic field is a very important property. If the initial magnetic field is divergence-free, this condition will be maintained for all later times. Numerically, if this divergence-free condition is violated seriously, an unphysical force will be created to parallel the magnetic field [16], which can produce some numerical instabilities or nonphysical features in simulation. The constrained transport (CT) methodology, which was originally introduced by Evans and Hawley [26], is an effective way to maintain the (discrete) divergence-free condition up to machine precision. Many improvements and extensions [1, 4, 15-18, 20, 21, 29] are followed.
Recently, the MHD system in the low sonic Mach limit has attracted a lot of interests [19,25,33,34,42,45,48]. Jiang, Ju, and Li [33] have shown that a weak solution of compressible MHD equations will converge to a strong solution of their corresponding incompressible MHD equations in this low sonic Mach limit. They have also investigated the low Mach limit of full MHD equations with a heat conductivity [34]. Cui, Ou, and Ren [19] have established a uniform convergence from full compressible MHD equations to isentropic incompressible MHD equations with well-prepared initial conditions in the three-dimensional case. Numerically it is also very attractive to develop schemes for MHD flows at any speed. For the above mentioned explicit shock capturing schemes, when applied to the MHD system in the low sonic Mach regime, the time step is subject to a very strict CFL condition which is proportional to the sonic Mach number, making it very undesirable for all speed flows. An implicit time discretization can release the small time step restrictions from an explicit one, but usually it results in a highly nonlinear system, which is not easy to be solved efficiently with an iterative method and it is hard to guarantee a fast convergence. Furthermore, numerical viscosities for such schemes are inversely proportional to the sonic Mach number, which introduce excessive numerical dissipations. Instead, in recent years some researchers are in pursuit of semi-implicit (SI) schemes with low dissipations. In [25], Dumbser et. al. developed a pressure-based SI finite volume all Mach number flow solver for compressible MHD equations. Minoshima and Miyoshi [48] proposed a multistate low-dissipation advection upstream splitting method for MHD with both high and low Mach numbers, by using a Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver. Leidi et al. [42] presented an SI finite volume solver with a 5-wave HLLD and a well-balanced method to efficiently simulate MHD flows at low Mach numbers with a gravitational source. For the MHD system with all sonic Mach numbers, it is important to design schemes with asymptotic stability and consistency as the sonic Mach number goes to 0 in the incompressible limit, namely, asymptotic preserving (AP). AP schemes have been widely used in a variety of areas, such as hydrodynamic or diffusive limits of kinetic models, relaxation limits of hyperbolic models, and low Mach number limits of compressible fluid models. We refer readers to the papers [23,24,27,[35][36][37][38] for more information. To achieve high order AP schemes in time, SI implicit-explicit Runge-Kutta (IMEX-RK) methods are widely used, which apply an explicit time discretization for nonstiff terms, while an implicit time discretization for stiff terms [5-11, 13, 39, 40, 47, 50]. With a suitable choice of explicit and implicit discretizations, both high efficiency and uniform stability independent of the stiff parameter can be obtained. As related, in [57] Tavelli and Dumbser developed an SI space-time discontinuous Galerkin method for all Mach-number flow, which involves staggered meshes and a Picard iteration for solving a large nonlinear system. In [8,10], Boscarino, Qiu, Russo, and Xiong presented high-order SI IMEX weighted essentially nonoscillatory (WENO) schemes for all-Mach isentropic Euler equations and all-Mach full Euler system, respectively. In [11], Boscheri and Pareschi developed high order pressure-based semi-implicit IMEX schemes for the three dimensional Navier-Stokes equations at all Mach numbers. In [30], Huang, Xing, and Xiong designed a high-order well-balanced SI AP scheme for shallow water equations in all Fraude numbers. We refer to many other related works from the references therein.
In this paper, we would like to develop a high-order AP finite difference WENO scheme with CT for MHD equations in all sonic Mach numbers. To our best knowledge, there are very few AP schemes for MHD systems in literature. The new contributions and innovations of this work are outlined as follows: • In the low sonic Mach number regime, SI IMEX-RK methods are adopted for time discretization to obtain a uniform time step independent of the sonic Mach number. Besides, to avoid a nonlinearity from the equation of state (EOS), an SI approach is used, leading to a linearized elliptic equation for the pressure. With carefully designed explicit and implicit discretizations, only a linear system needs to be solved. The resulting scheme is more efficient than explicit shock capturing schemes, especially in the low sonic Mach regime.
• Numerical viscosities are carefully designed to avoid excessive numerical dissipations, and a CT method is applied to maintain a discrete divergence-free property. The scheme can well capture discontinuities such as shocks in an essentially non-oscillatory fashion in the compressible regime. Meanwhile, the scheme is also a good incompressible solver as ε → 0. Due to less dissipation, our scheme can be shown to perform better than explicit schemes for low sonic Mach problems. • Formal AP and asymptotically accurate (AA) properties in the stiff limit as the sonic Mach number ε → 0 are proved, by assuming the implicit part of an SI IMEX-RK scheme is stiffly accurate (SA). The rest of the paper is as follows. In Section 2, we will briefly review the ideal MHD system, corresponding to its low sonic Mach limit based on asymptotic expansions, and a CT methodology for the divergence-free condition. A first-order SI scheme in time will first be introduced in Section 3, and then a high-order SI scheme with IMEX-RK in time and finite difference WENO in space is followed. Formal proofs of AP and AA properties are given in Section 4. Numerical experiments are performed in Section 5. A brief conclusion will be drawn in the last section.

Equations of compressible ideal MHD
The ideal MHD equations in a conservative form can be written as where ρ is the fluid mass density, u = (u, v, w) is the velocity, B = (B x , B y , B z ) is the magnetic field, E is the total energy density, p is the gas pressure, and · is the Euclidean vector norm. The total energy is given by where γ is the specific heat ratio. The system is subject to a divergence-free condition of the magnetic field If initially, the divergence-free condition holds, from (2.1), it also holds for all later times [18]. For this reason, (2.3) is usually not regarded as a constraint (like the ∇·u = 0 constraint for the incompressible Navier-Stokes equations), but rather an involution [16,18,29].
2.1. Eigenstructure of the MHD system. For (2.1), the eigenvalues of the Jacobian for the flux function along a normal direction n are given as: where λ 1,8 , λ 2,7 , λ 3,6 , λ 4 , and λ 5 are the left/right fast magnetosonic waves, left/right Alfvén waves, left/right slow magnetosonic waves, entropy wave, and divergence wave, correspondingly. The sound speed a, Alfvén speed c a , fast magnetosonic speed c f , and slow magnetosonic speed c s are defined as A fixed ordering for the MHD eigenvalues is Spectral decomposition for this system has already been described in [12,14,32,51,52,55], and we will follow the one in [32] for all simulations in the current paper.
2.2. Non-dimensionalized MHD. To derive a low sonic Mach limit for the ideal MHD system, we start with rewriting (2.1), (2.2), and (2.3) into a dimensionless form. We choose some reference or characteristic values, such as a length x 0 , a time t 0 , a fluid mass density ρ 0 , a velocity u 0 , a magnetic field B 0 , and a gas pressure p 0 , with u 0 = x 0 /t 0 . For the sake of simplicity, we consider a regime where the characteristic Alfvén Mach number u 0 √ ρ 0 /B 0 = 1. Then we define the following dimensionless Inserting these new variables into (2.1), we obtain the following non-dimensionalised MHD equations [42]: where we have dropped the hats of the dimensionless variables for ease of presentation. The EOS, divergence-free condition and the characteristic sonic Mach number ε are now given by In this dimensionless MHD system (2.6) with (2.7), the expressions of those eigenvalues (2.4) remain unchanged, except the sound speed a, Alfvén speed c a , fast magnetosonic speed c f , and slow magnetosonic speed c s from (2.5) now become Following the eigenstructure of the dimensionless MHD system discussed above, it is natural to find that if the background flow velocity u 0 is slow, then the sonic Mach number ε = u 0 √ ρ 0 / √ p 0 1, and the fast magnetosonic speed c f becomes very large. We can see that large magnetosonic speed c f will make the fast magnetosonic waves λ 1,8 to be much faster as compared to other waves. In this situation, for explicit shock-capturing schemes, the time step subjects to ∆t = CFL ∆x max where ∆t is the time step size, ∆x is the mesh size, and CFL is the time stability CFL number. On one hand, highly inaccurate solutions will appear due to the excessive numerical viscosity (scales as ε −1 ) in standard upwind schemes. On the other hand, this restriction results in an increasingly large computational cost for low sonic Mach fluid flows. A straight approach is to apply an implicit time discretization to avoid such severe time step constraints. However, a fully implicit scheme will result in a complicated nonlinear system which is usually very hard to be solved efficiently. Besides, it may not be able to guarantee a correct asymptotic limit. Thus, it is important to design efficient schemes with both asymptotic stability and consistency in the low sonic Mach limit, namely, the AP property.
2.3. Low sonic Mach limit. Next, we recall a formal derivation of the incompressible MHD system from the non-dimensionalized compressible MHD equations (2.6) with (2.7). We consider an asymptotic expansion ansatz for the following variables [46]: Inserting (2.8) into (2.6) and (2.7), equating to zero for different orders of ε, corresponding to the first three leading orders, we have where we have assumed a similar expansion of E as p in (2.8), and take E 0 = p 0 /(γ − 1). Denoting the material derivative d/dt = ∂/∂t + u 0 · ∇, from the energy equation in (2.9), E 0 = p 0 /(γ − 1), with ∇p 0 = 0, we have Taking into account periodic or no-slip u · n = 0 boundary conditions on the spatial domain Ω, and integrating (2.10) over Ω, we get Ω ∇ · u 0 dx = ∂Ω u 0 · nds = 0, where n is the unit outward normal vector along ∂Ω. This implies p 0 is constant not only in space but also in time, that is p 0 = Const. A direct result ∇ · u 0 = 0 can be obtained through (2.10) with p 0 = Const. Inserting ∇ · u 0 = 0 into the mass-continuity equation in (2.9), we get dρ 0 /dt = 0 which means that the density is constant along the characteristic trajectories. For the order O(ε) in (2.6), we have: Due to E 1 = p 1 /(γ − 1) from the EOS in (2.7), similarly as above, we get p 1 = Const and ∇ · u 1 = 0. In this case, we include εu 1 into u 0 , and εp 1 into p 0 , the asymptotic expansions of u and p become (2.11) u = u 0 + ε 2 u 2 + · · · , and (2.12) p = p 0 + ε 2 p 2 + · · · .
Using such expansions for u and p, as ε → 0, the leading order of the system (2.6) yields Numerically, in order to have a correct asymptotic limit from (2.6) to (2.13), it is important to impose well-prepared initial conditions, which are consistent with the expansions given in (2.8), with u and p from (2.11) and (2.12). According to [22,41], the well-prepared initial conditions for (2.6) are given as where ρ 0 (x) is a strictly positive function, p 0 = Const, ∇ · u 0 = 0, and ∇ · B 0 = 0.
2.4. CT for ∇ · B = 0. Many CT methods have been proposed to satisfy the discrete divergence-free condition ∇ · B = 0 for the ideal MHD equations [15][16][17][18]. Here we consider the one described in [18], in which the magnetic field is written as the curl of a magnetic vector potential A: , the equation of the magnetic field B in the MHD system can be rewritten as Inserting (2.15) into (2.16), we further get which implies the existence of a scalar function ψ such that Stable solutions can be obtained by introducing a Weyl gauge, i.e. ψ ≡ 0 [29], and the evolution equation for the vector potential becomes It is easy to check that ∇ · B = ∇ · (∇ × A) = 0, i.e. the divergence-free condition will be satisfied for all times. In summary, we get the following MHD system: Corresponding to (2.17), with well-prepared initial conditions (2.14), the low sonic Mach limit incompressible MHD system becomes

Numerical schemes
In this section, we will construct and analyze a high-order finite difference scheme with AP and divergence-free properties for the non-dimensionalized MHD system (2.17). The features of our scheme include the following: we design an SI IMEX-RK time discretization so that the scheme is stable with a time step constraint independent of the sonic Mach number ε; the AP property is preserved in the zero sonic Mach number limit; the scheme can be implemented efficiently without solving a nonlinear system. We will adopt the high order finite difference WENO scheme in [8,10,30], and the CT method in [18]. The final scheme can well capture discontinuities such as shocks in the compressible regime, while keeping a discrete divergence-free condition in the low sonic Mach regime, with a uniform time step condition independent of ε. We start with a first-order SI scheme. The high-order IMEX-RK scheme in time and a finite difference WENO reconstruction in space will be followed. After that, an overall scheme will be summarized.

3.2.
High order SI IMEX-RK scheme. In this subsection, we will follow a similar procedure as in [8,10] to extend the first-order SI IMEX scheme (3.1) to high order. For simplicity, we only consider (2.6) with CT methods, namely, (2.17). First, (2.17) can be written as an autonomous system: where U = (ρ, q, A, E) T and H : R n × R n × R m → R n is a sufficiently regular mapping. We use subscript "E" to express an explicit treatment for the first argument U and the third one B, and subscript "I" to express an implicit treatment for the second U . For the system (2.17), the function To obtain high-order accuracy in time, we can apply a multi-stage IMEX-RK time discretization with a double Butcher tableau,cÃ is an s × s matrix for an explicit scheme, withã ij = 0 for j ≥ i, and A = (a ij ) is an s × s matrix for an implicit scheme. For the implicit part of the methods, we use a diagonally implicit scheme, i.e. a ij = 0, for j > i, to guarantee simplicity and efficiency in solving the algebraic equations corresponding to the implicit part of the discretization. The vectorsc = ( Later, for the consideration of AP property, we consider implicit schemes with SA property, i.e., the implicit part of the Butcher table satisfies the condition b T = e T s A, with e s = (0, · · · 0, 1) T . We will see that the SA property guarantees that the numerical solution is identical to the last internal stage value of the scheme. Now, we may update the solutions as follows. Starting from U • Update the known values for the implicit part U and then solve The first order scheme (3.1) is the same as applying the following Butcher table to (2.17): namely U E = U n and U I = U n+1 .
3.3. High order finite difference WENO scheme. Below we will describe our spatial discretization strategies which incorporate a WENO mechanism in order to capture shocks in the compressible regime and produce a high order incompressible solver for the flow in the zero Mach limit. In particular, we propose to apply a fifth order characteristic-wise WENO procedure to F 1 E of the ideal MHD system so that in the compressible regime numerical oscillations could be well controlled, and to apply a component-wise WENO procedure for the flux functions of F 2 E and F SI . As for G E , which is a Hamiltonian, a WENO procedure described in [31] will be applied.
• Component-wise WENO ∇ W for ∇ · F 2 E and ∇ · F SI : Let v = (ρ, q, B, E) T and f = f (v) be the unknown variables and flux of the ideal MHD system (2.1). For simplicity, we consider the 1D case as an example, and we take a uniform mesh with N + 1 grid points: a = x 0 < x 1 < · · · < x N = b. A conservative finite difference WENO scheme for system (2.1) can be written in the following flux-difference form: are the numerical fluxes, which can be reconstructed by the following procedures: (1) Compute the physical flux at each grid point: (2) Perform a Lax-Friedrichs flux vector splitting for each component of the physical variables to get (3) Perform a finite difference WENO reconstruction to obtain upwind and downwind numerical fluxes: Here Φ WENO5 denotes a 5th-order WENO reconstruction [31], while other orders can also be used. Then we setf • Characteristic-wise WENO ∇ CW for ∇ · F 1 E : The Jacobian matrix for the flux function ∂f /∂v, has a spectral decomposition of the form ∂f /∂v = RDL [18], where D is a diagonal matrix of real eigenvalues, R is the matrix of right eigenvectors, and L = R −1 is the matrix of left eigenvectors. Then, for a characteristic-wise WENO reconstruction ∇ CW , we have the following procedures: (1) Compute the physical flux at each grid point: (2) Compute the average state v i+ 1 2 by using an arithmetic mean: (3) Compute the right and left eigenvectors of the flux Jacobian matrix ∂f /∂v by taking ).
(4) Project the solution and physical flux into the characteristic space: for all j in the numerical stencil associated with x = x i+ 1 2 . In the case of a 5th-order finite difference WENO scheme: j = i − 2, i − 1, i, i + 1, i + 2, i + 3.
where α has the same expression in the component-wise WENO (3.11). (6) Perform a WENO reconstruction Φ WENO5 on each component of g ± j to obtain their corresponding numerical fluxes, and then set (7) Project the numerical fluxĝ i+ 1 2 back to the physical space: However, when the sonic Mach number ε is small, the fast eigenvalue c f will make such a spectral characteristic decomposition to be unstable. Since a characteristic decomposition is mainly required in the compressible regime where ε = O(1), while it does not make any significant differences for small Mach numbers; in this case we directly use the eigenvectors in the expression of (2.1) by taking ε = 1, which works well from our numerical results.
• A Hamiltonian WENO reconstruction WENO-HJ (·) HJ for G E : We take a 1D Hamilton-Jacobi equation as an example where HT is the Hamiltonian. As described in [31], a semi-discrete form for (3.13) is where HT is a numerical Hamiltonian which approximates (3.13), and v − x,i , v + x,i are upwind and downwind approximations of v x at x = x i . A Lax-Friedrichs Hamiltonian from [31] is given as follows (3.14) HT Such a Hamiltonian WENO reconstruction helps us to control unphysical oscillations both in v x and v. For G E in 2D, the Hamiltonian is and the numerical Hamiltonian (3.14) can be applied along each direction for u(t, x, y) ∂Az ∂x and v(t, x, y) ∂A 3 ∂y , respectively. Here A z is the third component of A.

3.4.
Flowchart of a high order SI IMEX-RK finite difference WENO scheme. With the above space and time discretizations, we now summarize our high order scheme as follows: • Step 1. Starting from U n at time t n , for an intermediate stage (2) We can obtain ρ /2). Notice that in the process of substitution to obtain (3.16) or (3.18), the gradient and the divergence operators are kept to be continuous, obtaining the second order operator ∇ · (H (i) ∇p (i) I,2 ). The second-order spatial derivative is then discretized by a compact discretization as proposed in [10].
and successively update E Step 4. Finally, we update the numerical solution U n+1 = U (s) I , and B n+1 = ∇ × A n+1 . 4. Asymptotic preserving (AP) and Asymptotically Accurate (AA) properties 4.1. AP property. In this section, we will prove the AP property for the first order SI IMEX scheme (3.1), and the AA property for the high order SI IMEX-RK scheme (3.7)-(3.10). We focus on time discretizations while keeping the space continuous when we discuss the AP or AA property. First, we have the following theorem: Theorem 4.1. The first order SI IMEX scheme (3.1) with space continuous is AP, in the sense that, with well-prepared initial conditions (2.14), at the leading order asymptotic expansions, the scheme (3.1) is a consistent approximation of the incompressible MHD equations (2.18).

AA property.
Similarly to the AP property, here we focus on the AA analysis on time discretizations, while keeping space continuous. Then we have the following theorem: We consider a high order SI IMEX-RK scheme (3.7)-(3.10) of order r applied to system (2.17) on a bounded domain Ω ⊂ R with periodic or compact support boundary conditions, and assume that the implicit part of the SI IMEX-RK scheme is SA and the initial conditions (ρ 0 (x), ρ 0 (x)u 0 (x), p 0 , A 0 (x), B 0 (x)) T are well-prepared, namely in the form of (4.1). If we denote (ρ 1 (x; ε), ρ 1 (x; ε)u 1 (x; ε), p 1 (x; ε), A 1 (x; ε), B 1 (x; ε)) T as the numerical solution after one-time step, then, with p * being a constant, we have: ) T be the exact solution of the incompressible MHD equations (2.18) with the same initial data, one has the following one-step error estimate i.e. the scheme is AA.
Proof. We consider the first step from t 0 = 0 to t 1 = ∆t for the SI IMEX-RK scheme (3.7)-(3.10) of order r applied to system (2.17) with well-prepared initial data (4.1): where ρ 0 inc := ρ inc (x, 0), u 0 inc := u inc (x, 0), A 0 inc := A inc (x, 0), B 0 inc := B inc (x, 0). By well-prepared assumptions we have: p inc (x, 0) := p * is a constant independent of time and space, ∇ · u 0 inc = 0 and ∇ · B 0 (x) = 0. Now we consider a formal ε-expansion of the quantities: For example, the expansions of the density and pressure are To prove the theorem, we use mathematical induction.
• AA property for internal stages i = 1, · · · , s. When i = 1, it leads to the same AP analysis for the scheme (3.1) with ∆t replaced by a 11 ∆t. To prove the result for i > 1, we make use of the induction hypothesis, assuming the property holds for j ≤ i − 1, and prove that it still holds for j = i. For j = 1, · · · , i − 1, we have Now we insert the expansions (4.7) into the explicit step in (3.7), and from the energy equation we get with E 0 inc = p * /(γ − 1) and for j = 1, · · · , i − 1, , in the order of O(1), the density and momentum equations yield (4.12) where from (3.17) and (4.8), it follows ∇p (j) E = 0 for j = 1, · · · , i − 1. Using (4.10), we get (4.14) E Therefore, from (4.8) and (4.14), we get Now from (4.14) and (3.9), it follows for the energy equation

Numerical tests
In this section, we perform some numerical experiments to validate the high order accuracy, AP, AA, divergence-free, and good performances of our proposed scheme for the MHD system with all sonic Mach numbers. The third order SI IMEX-RK with an SA property from [8,30] is used for time discretizations. The fifth order finite difference WENO reconstruction [8,30,31,53] for the first order spatial derivatives, and a fourth order compact central difference discretization [8,10] for the second order derivatives will be applied for spatial discretizations. Our scheme is referred as the "IMEX" scheme. For some of the following examples, we will compare our results to reference solutions, which are produced by a fifth-order finite difference WENO scheme with a third-order strong stability preserving RK method. It is referred as an "ERK" scheme if without a CT method [32], while as an "ERKC" scheme if with a CT method [18]. To avoid excessive numerical dissipations for low sonic Mach numbers, for the explicit schemes "ERK" and "ERKC", here we also take a Lax-Friedrichs flux the same as in (3.11) and (3.12), and numerically we find that it leads to much better results as we can see in the following. Unless otherwise specified, the CFL number is taken as 0.25, and the gas constant is γ = 5/3. For 1D problems, the time step is with N + 1 computational grid points, and the definition ofĉ has been given in (3.12). For 2D problems, the time step is with (N x + 1) × (N y + 1) computational grid points. Here, the fast speeds in the x-direciton (ĉ (k,l) f,x ) and in the y-direciton (ĉ 5.1. One dimensional case. Since the divergence-free condition is satisfied automatically in 1D, therefore, the CT method will not need to be applied in the following 1D numerical simulations. with a periodic boundary condition. The exact solution to (5.1) propagates with a unit Alfvén speed (that is U (t, x) = U (0, x + t).). We set the time step to be ∆t = CFL∆x The sonic Mach number is ε = 1 and we compute up to a final time T = 1. The numerical L 1 , L 2 , and L ∞ errors and orders are shown in Table 5.1. We can observe that a fifth-order spatial accuracy is obtained.

5.2.
Two dimensional case. A CT method described in [18] will be used to maintain the divergencefree condition. In 2D, B = ∇ × A can be written as B x = ∂A z /∂y and B y = −∂A z /∂x. Therefore, we only need to update the third component A z of A, and ignore A x and A y . All solutions of the following tests maintain ∇ · B = 0 up to machine round off errors [16]. from Table 5.3 we can find around 5th-order accuracy for the momentum ρu, and 4th-order accuracy for the magnetic field component B x where a fourth-order central difference discretization is used for B = ∇ × A.
In this test, a loop of the magnetic field is advected. In Fig. 5.6, we show both the deviation of the density ρ away from 1 and the magnetic pressure B 2 /2 for the case of ε = 0.1 at T = √ 5. For the magnetic pressure, both schemes preserve the initial symmetric loop structure very well, and the IMEX scheme has almost the same results as the ERKC scheme. The one-dimensional cutting plots for the magnetic pressure along x = 0.2 show that both schemes slightly distort the profile at y = 0.15. For the deviation of the density, we can find the IMEX scheme is better than the ERKC scheme for preserving a symmetric structure, as compared with the results of the ERKC scheme on a uniform mesh of 384 × 192 cells. For this problem in the intermediate regime with ε = 0.1, we also compare the CPU cost of two schemes in Table 5.5. We can observe that the ERKC scheme costs less when ε is of O(1). However, the IMEX scheme is much more efficient in the intermediate and low sonic Mach regimes. The computational cost for the ERKC scheme increases very rapidly, as the sonic Mach number ε becomes smaller. This example demonstrates the main advantage of an AP scheme for such multiscale problems.  Example 5.9. (Magnetized Kelvin-Helmholtz instability) Finally, we consider a simulation of the magnetized Kelvin-Helmholtz instability problem [8,42]. The initial conditions are  Periodic boundary conditions are used for the conserved quantities, while a first-order extrapolation is utilized when the magnetic potential is outside of the computational domain. We run the solution up to T = 0.8 with γ = 1.4 and ε = 10 −6 on a mesh grid N x × N y = 256 × 128. We define the local sonic Mach number as M a = √ u 2 + v 2 / γp/ρ, and the ratio of the local sonic Mach number M ratio to the maximum value of M a as M ratio = M a / max(M a ), where the maximum is taken over all computational grid points. The vorticity ω = v x − u y , where v x and u y are discretized by the 4th order central difference, and M ratio are shown in Fig. 5.7. We observe that they are comparable to the results in [8,42], namely, our IMEX scheme can capture the incompressible MHD system in the low sonic Mach number regime very well, which verifies the AP property of our scheme.

Conclusion
In this paper, a high-order SI AP scheme with a discrete divergence-free property for the MHD equations has been developed. High-order accuracy in time has been obtained by SI temporal integrator based on an IMEX-RK framework. High-order accuracy in space has been achieved by finite difference WENO reconstructions. We have formally proved that the scheme is AP. We have also proved the AA property in the stiff limit as the sonic Mach number ε → 0, if the SI IMEX-RK scheme is SA. Numerical experiments in 1D and 2D have demonstrated the divergence-free for all ranges of sonic Mach numbers, AP and AA properties in the low sonic Mach limit. Compared to the explicit schemes as reference solutions, such as the ERK [32] and ERKC [18] schemes in 1D and 2D respectively, our IMEX scheme has almost the same performances as the ERK or ERKC scheme, but is much more efficient in the intermediate and low sonic Mach regimes.