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

###### Abstract.

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.

###### Key words and phrases:

all sonic Mach number; MHD equations; divergence-free; asymptotic preserving;
SI IMEX-RK; finite difference WENO

Wei Chen

School of Mathematical Sciences, Xiamen University

Xiamen, Fujian, 361005, P.R. China

Email<EMAIL_ADDRESS>
Kailiang Wu***This work is supported in part by NSFC grant 12171227

Department of Mathematics & SUSTech International Center for Mathematics

Southern University of Science and Technology

National Center for Applied Mathematics Shenzhen (NCAMS)

Shenzhen, Guangdong 518055, China

Email<EMAIL_ADDRESS>
Tao Xiong†††Corresponding author. The work of this author was partially
supported by NSFC grant No. 11971025, and the Strategic Priority Research
Program of the Chinese Academy of Sciences Grant No. XDA25010401.

School of Mathematical Sciences, Xiamen University

Fujian Provincial Key Laboratory of Mathematical Modeling and High-Performance
Scientific Computing

Xiamen, Fujian, 361005, P.R. China

Email<EMAIL_ADDRESS>

##  1\. Introduction

Ideal magnetohydrodynamic (MHD) equations are widely used in the modeling of
weather prediction, astrophysics, as well as laboratory plasma applications
such as flows in tokamaks and stellarators. Many shock capturing schemes with
explicit time discretizations have been developed for solving compressible
ideal MHD equations, including high order discontinuous Galerkin [43, 44, 60,
59], finite difference [49, 15, 18, 16, 17], and finite volume schemes [42,
63, 2, 54, 62], etc.

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 the computed solution, cause the loss
of pressure positivity [58, 61], and/or lead to the failure of the 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 [4, 1, 20, 21, 15, 16, 17, 18, 29] are followed.

Recently, the MHD system in the low sonic Mach limit has attracted a lot of
interests [45, 42, 25, 33, 34, 19, 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 [35, 37, 36, 38, 27, 24, 23] 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, 7, 6, 9, 10,
8, 40, 39, 50, 13, 47, 11]. 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 [10, 8], 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 $\varepsilon\rightarrow 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 $\varepsilon\rightarrow 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.

##  2\. Equations of compressible ideal MHD

The ideal MHD equations in a conservative form can be written as

(2.1) $\frac{\partial}{\partial t}\begin{bmatrix}\rho\\\ \rho\mathbf{u}\\\
\mathbf{B}\\\ E\end{bmatrix}+\nabla\cdot\begin{bmatrix}\rho\mathbf{u}\\\
\rho\mathbf{u}\otimes\mathbf{u}+(p+\frac{1}{2}{\|\mathbf{B}\|}^{2})\mathbf{I}-\mathbf{B}\otimes\mathbf{B}\\\
\mathbf{u}\otimes\mathbf{B}-\mathbf{B}\otimes\mathbf{u}\\\
(E+p+\frac{1}{2}{\|\mathbf{B}\|}^{2})\mathbf{u}-\mathbf{B}(\mathbf{u}\cdot\mathbf{B})\end{bmatrix}=\mathbf{0},$

where $\rho$ is the fluid mass density, $\mathbf{u}=(u,v,w)$ is the velocity,
$\mathbf{B}=(B_{x},B_{y},B_{z})$ is the magnetic field, $E$ is the total
energy density, $p$ is the gas pressure, and $\|\cdot\|$ is the Euclidean
vector norm. The total energy is given by

(2.2)
$E=\frac{p}{\gamma-1}+\frac{1}{2}(\rho{\|\mathbf{u}\|}^{2}+{\|\mathbf{B}\|}^{2}),$

where $\gamma$ is the specific heat ratio. The system is subject to a
divergence-free condition of the magnetic field

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

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 $\nabla\cdot{\bf u}=0$ constraint for the incompressible
Navier-Stokes equations), but rather an involution [29, 18, 16].

###  2.1. Eigenstructure of the MHD system

For (2.1), the eigenvalues of the Jacobian for the flux function along a
normal direction ${\bf n}$ are given as:

(2.4) $\lambda_{1,8}={\bf u}\cdot{\bf n}\mp c_{f},\quad\lambda_{2,7}={\bf
u}\cdot{\bf n}\mp c_{a},\quad\lambda_{3,6}={\bf u}\cdot{\bf n}\mp
c_{s},\quad\lambda_{4,5}={\bf u}\cdot{\bf n},\quad$

where $\lambda_{1,8}$, $\lambda_{2,7}$, $\lambda_{3,6}$, $\lambda_{4}$, and
$\lambda_{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

(2.5) $\left\\{\begin{aligned} a&:=\sqrt{\frac{\gamma p}{\rho}},\\\
c_{a}&:=\sqrt{\frac{{({\bf B}\cdot{\bf n})}^{2}}{\rho}},\\\
c_{f,s}&:={\left\\{\frac{1}{2}\left[a^{2}+\frac{{\|{\bf
B}\|}^{2}}{\rho}\pm\sqrt{{\left(a^{2}+\frac{{\|{\bf
B}\|}^{2}}{\rho}\right)}^{2}-4a^{2}{c^{2}_{a}}}\right]\right\\}}^{\frac{1}{2}}.\end{aligned}\right.$

A fixed ordering for the MHD eigenvalues is

$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\leq\lambda_{4}\leq\lambda_{5}\leq\lambda_{6}\leq\lambda_{7}\leq\lambda_{8}.$

Spectral decomposition for this system has already been described in [32, 12,
14, 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 $\rho_{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}\sqrt{\rho_{0}}/B_{0}=1$. Then we define the following dimensionless
variables:

$\hat{{{\bf x}}}=\frac{{{\bf
x}}}{x_{0}},\quad\hat{t}=\frac{t}{t_{0}},\quad\hat{\rho}=\frac{\rho}{\rho_{0}},\quad\hat{{\bf
u}}=\frac{{\bf u}}{u_{0}},\quad\hat{{\bf B}}=\frac{{\bf
B}}{B_{0}},\quad\hat{p}=\frac{p}{p_{0}}.$

Inserting these new variables into (2.1), we obtain the following non-
dimensionalised MHD equations[42]:

(2.6) $\frac{\partial}{\partial t}\begin{bmatrix}\rho\\\ \rho\mathbf{u}\\\
\mathbf{B}\\\ E\\\ \end{bmatrix}+\nabla\cdot\begin{bmatrix}\rho\mathbf{u}\\\
\rho\mathbf{u}\otimes\mathbf{u}+(\frac{p}{\varepsilon^{2}}+\frac{1}{2}{\|\mathbf{B}\|}^{2})\mathbf{I}-\mathbf{B}\otimes\mathbf{B}\\\
\mathbf{u}\otimes\mathbf{B}-\mathbf{B}\otimes\mathbf{u}\\\
(E+p+\frac{\varepsilon^{2}}{2}{\|\mathbf{B}\|}^{2})\mathbf{u}-\varepsilon^{2}\mathbf{B}(\mathbf{u}\cdot\mathbf{B})\end{bmatrix}=\mathbf{0},$

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 $\varepsilon$ are now given by

(2.7) $E:=\frac{p}{\gamma-1}+\frac{\varepsilon^{2}}{2}(\rho{\|{\bf
u}\|}^{2}+{\|{\bf B}\|}^{2}),\quad\nabla\cdot{\bf
B}=0,\quad\varepsilon:=u_{0}\sqrt{\frac{\rho_{0}}{p_{0}}}.$

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

$\left\\{\begin{aligned} a&:=\frac{1}{\varepsilon}\sqrt{\frac{\gamma
p}{\rho}},\\\ c_{a}&:=\sqrt{\frac{{({\bf B}\cdot{\bf n})}^{2}}{\rho}},\\\
c_{f,s}&:={\left\\{\frac{1}{2}\left[a^{2}+\frac{{\|{\bf
B}\|}^{2}}{\rho}\pm\sqrt{{\left(a^{2}+\frac{{\|{\bf
B}\|}^{2}}{\rho}\right)}^{2}-4a^{2}{c^{2}_{a}}}\right]\right\\}}^{\frac{1}{2}}.\end{aligned}\right.$

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 $\varepsilon=u_{0}\sqrt{\rho_{0}}/\sqrt{p_{0}}\ll
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
$\lambda_{1,8}$ to be much faster as compared to other waves. In this
situation, for explicit shock-capturing schemes, the time step subjects to

$\Delta t={\rm CFL}\frac{\Delta x}{\mathop{\max}\limits_{1\leq i\leq
8}|\lambda_{i}|}=\mathcal{O}(\varepsilon\Delta x),$

where $\Delta t$ is the time step size, $\Delta x$ is the mesh size, and ${\rm
CFL}$ is the time stability CFL number. On one hand, highly inaccurate
solutions will appear due to the excessive numerical viscosity (scales as
$\varepsilon^{-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]:

(2.8) $\left\\{\begin{aligned}
&p(\mathbf{x},t)=p_{0}(\mathbf{x},t)+\varepsilon
p_{1}(\mathbf{x},t)+\varepsilon^{2}p_{2}(\mathbf{x},t)+\cdots,\\\
&\rho(\mathbf{x},t)=\rho_{0}(\mathbf{x},t)+\varepsilon\rho_{1}(\mathbf{x},t)+\varepsilon^{2}\rho_{2}(\mathbf{x},t)+\cdots,\\\
&\mathbf{u}(\mathbf{x},t)=\mathbf{u}_{0}(\mathbf{x},t)+\varepsilon\mathbf{u}_{1}(\mathbf{x},t)+\varepsilon^{2}\mathbf{u}_{2}(\mathbf{x},t)+\cdots,\\\
&\mathbf{B}(\mathbf{x},t)=\mathbf{B}_{0}(\mathbf{x},t)+\varepsilon\mathbf{B}_{1}(\mathbf{x},t)+\varepsilon^{2}\mathbf{B}_{2}(\mathbf{x},t)+\cdots.\end{aligned}\right.$

Inserting (2.8) into (2.6) and (2.7), equating to zero for different orders of
$\varepsilon$, corresponding to the first three leading orders, we have

  * •

$\mathcal{O}(\varepsilon^{-2})$

$\nabla p_{0}=0,$

  * •

$\mathcal{O}(\varepsilon^{-1})$

$\nabla p_{1}=0,$

  * •

$\mathcal{O}(\varepsilon^{0})$

(2.9) $\frac{\partial}{\partial t}\begin{bmatrix}\rho_{0}\\\
\rho_{0}\mathbf{u}_{0}\\\ \mathbf{B}_{0}\\\ E_{0}\\\
\end{bmatrix}+\nabla\cdot\begin{bmatrix}\rho_{0}\mathbf{u}_{0}\\\
\rho_{0}\mathbf{u}_{0}\otimes\mathbf{u}_{0}+(p_{2}+\frac{1}{2}{\|\mathbf{B}_{0}\|}^{2})\mathbf{I}-\mathbf{B}_{0}\otimes\mathbf{B}_{0}\\\
\mathbf{u}_{0}\otimes\mathbf{B}_{0}-\mathbf{B}_{0}\otimes\mathbf{u}_{0}\\\
(E_{0}+p_{0})\mathbf{u}_{0}\end{bmatrix}=\mathbf{0},$

where we have assumed a similar expansion of $E$ as $p$ in (2.8), and take
$E_{0}=p_{0}/(\gamma-1)$. Denoting the material derivative ${\rm
d/d}t=\partial/\partial t+{\bf u}_{0}\cdot\nabla$, from the energy equation in
(2.9), $E_{0}=p_{0}/(\gamma-1)$, with $\nabla p_{0}=0$, we have

(2.10) $\nabla\cdot{\bf u}_{0}=-\frac{1}{p_{0}\gamma}\frac{{\rm d}p_{0}}{{\rm
d}t}.$

Taking into account periodic or no-slip ${\bf u}\cdot{\bf n}=0$ boundary
conditions on the spatial domain $\Omega$, and integrating (2.10) over
$\Omega$, we get $\int_{\Omega}\nabla\cdot{\bf u}_{0}d{\bf
x}=\int_{\partial\Omega}{\bf u}_{0}\cdot{\bf n}ds=0$, where ${\bf n}$ is the
unit outward normal vector along $\partial\Omega$. This implies $p_{0}$ is
constant not only in space but also in time, that is $p_{0}={\rm Const}$. A
direct result $\nabla\cdot{\bf u}_{0}=0$ can be obtained through (2.10) with
$p_{0}={\rm Const}$. Inserting $\nabla\cdot{\bf u}_{0}=0$ into the mass-
continuity equation in (2.9), we get ${\rm d}\rho_{0}/{\rm d}t=0$ which means
that the density is constant along the characteristic trajectories.

For the order $\mathcal{O}(\varepsilon)$ in (2.6), we have:

$\partial E_{1}/\partial t+\nabla\cdot[(E_{0}+p_{0}){\bf
u}_{1}+(E_{1}+p_{1}){\bf u}_{0}]=0.$

Due to $E_{1}=p_{1}/(\gamma-1)$ from the EOS in (2.7), similarly as above, we
get $p_{1}={\rm Const}$ and $\nabla\cdot{\bf u}_{1}=0$. In this case, we
include $\varepsilon\mathbf{u}_{1}$ into $\mathbf{u}_{0}$, and $\varepsilon
p_{1}$ into $p_{0}$, the asymptotic expansions of ${\bf u}$ and $p$ become

(2.11) $\mathbf{u}=\mathbf{u}_{0}+\varepsilon^{2}\mathbf{u}_{2}+\cdots,$

and

(2.12) $p=p_{0}+\varepsilon^{2}p_{2}+\cdots.$

Using such expansions for $\mathbf{u}$ and $p$, as $\varepsilon\rightarrow 0$,
the leading order of the system (2.6) yields

(2.13) $\left\\{\begin{aligned} &\frac{\partial}{\partial
t}\begin{bmatrix}\rho_{0}\mathbf{u}_{0}\\\ \mathbf{B}_{0}\\\
\end{bmatrix}+\nabla\cdot\begin{bmatrix}\rho_{0}\mathbf{u}_{0}\otimes\mathbf{u}_{0}+(p_{2}+\frac{1}{2}{\|\mathbf{B}_{0}\|}^{2})\mathbf{I}-\mathbf{B}_{0}\otimes\mathbf{B}_{0}\\\
\mathbf{u}_{0}\otimes\mathbf{B}_{0}-\mathbf{B}_{0}\otimes\mathbf{u}_{0}\end{bmatrix}=\mathbf{0},\\\
&\frac{{\rm d}\rho_{0}}{{\rm d}t}=0,\quad p_{0}={\rm
Const},\quad\nabla\cdot{\bf u}_{0}=0,\quad\nabla\cdot{\bf B}=0.\\\
\end{aligned}\right.$

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 ${\bf u}$ and $p$ from
(2.11) and (2.12). According to [22, 41], the well-prepared initial conditions
for (2.6) are given as

(2.14) $\left\\{\begin{aligned} &\rho(t=0,{\bf x})=\rho_{0}({\bf
x})+\mathcal{O}(\varepsilon)>0,\quad&&p(t=0,{\bf
x})=p_{0}+\varepsilon^{2}p_{2}({\bf x})+\cdots,\\\ &{\bf u}(t=0,{\bf x})={\bf
u}_{0}({\bf x})+\mathcal{O}(\varepsilon^{2}),\quad&&{\bf B}(t=0,{\bf x})={\bf
B}_{0}({\bf x})+\mathcal{O}(\varepsilon),\end{aligned}\right.$

where $\rho_{0}({\bf x})$ is a strictly positive function, $p_{0}={\rm
Const}$, $\nabla\cdot{\bf u}_{0}=0$, and $\nabla\cdot{\bf B}_{0}=0$.

###  2.4. CT for $\nabla\cdot{\bf B}=0$

Many CT methods have been proposed to satisfy the discrete divergence-free
condition $\nabla\cdot{\bf B}=0$ for the ideal MHD equations [18, 17, 16, 15].
Here we consider the one described in [18], in which the magnetic field is
written as the curl of a magnetic vector potential ${\bf A}$:

(2.15) ${\bf B}=\nabla\times{\bf A}.$

Due to $\nabla\cdot({\bf u}\otimes{\bf B}-{\bf B}\otimes{\bf
u})=\nabla\times({\bf B}\times{\bf u})$, the equation of the magnetic field
${\bf B}$ in the MHD system can be rewritten as

(2.16) $\frac{\partial{\bf B}}{\partial t}+\nabla\times({\bf B}\times{\bf
u})=0.$

Inserting (2.15) into (2.16), we further get

$\nabla\times\left(\frac{\partial{\bf A}}{\partial t}+(\nabla\times{\bf
A})\times{\bf u}\right)=0,$

which implies the existence of a scalar function $\psi$ such that

$\frac{\partial{\bf A}}{\partial t}+(\nabla\times{\bf A})\times{\bf
u}=-\nabla\psi.$

Stable solutions can be obtained by introducing a Weyl gauge, i.e. $\psi\equiv
0$ [29], and the evolution equation for the vector potential becomes

$\frac{\partial{\bf A}}{\partial t}+(\nabla\times{\bf A})\times{\bf u}=0.$

It is easy to check that $\nabla\cdot{\bf B}=\nabla\cdot(\nabla\times{\bf
A})=0$, i.e. the divergence-free condition will be satisfied for all times. In
summary, we get the following MHD system:

(2.17) $\left\\{\begin{aligned} &\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho{\bf u}\right)=0,\\\ &\frac{\partial\rho{\bf
u}}{\partial
t}+\nabla\cdot\left(\rho\mathbf{u}\otimes\mathbf{u}+\left(\frac{p}{\varepsilon^{2}}+\frac{1}{2}{\|\mathbf{B}\|}^{2}\right)\mathbf{I}-\mathbf{B}\otimes\mathbf{B}\right)=\mathbf{0},\\\
&\frac{\partial{\bf A}}{\partial t}+(\nabla\times{\bf A})\times{\bf
u}=\mathbf{0},\\\ &\frac{\partial E}{\partial
t}+\nabla\cdot\left(\left(E+p+\frac{\varepsilon^{2}}{2}{\|\mathbf{B}\|}^{2}\right)\mathbf{u}-\varepsilon^{2}\mathbf{B}\left(\mathbf{u}\cdot\mathbf{B}\right)\right)=0,\\\
&{\bf B}-\nabla\times{\bf A}=\mathbf{0}.\\\ \end{aligned}\right.$

Corresponding to (2.17), with well-prepared initial conditions (2.14), the low
sonic Mach limit incompressible MHD system becomes

(2.18) $\left\\{\begin{aligned} &\frac{d\rho_{0}}{dt}=0,\quad p_{0}={\rm
Const},\quad\nabla\cdot{\bf u}_{0}=0,\quad\nabla\cdot\mathbf{B}_{0}=0,\\\
&\frac{\partial\rho_{0}\mathbf{u}_{0}}{\partial
t}+\nabla\cdot\left(\rho_{0}\mathbf{u}_{0}\otimes\mathbf{u}_{0}+\left(p_{2}+\frac{1}{2}{\|\mathbf{B}_{0}\|}^{2}\right)\mathbf{I}-\mathbf{B}_{0}\otimes\mathbf{B}_{0}\right)=\mathbf{0},\\\
&\frac{\partial{\bf A}_{0}}{\partial t}+(\nabla\times{\bf A}_{0})\times{\bf
u}_{0}=\mathbf{0},\\\ &{\bf B}_{0}-\nabla\times{\bf
A}_{0}=\mathbf{0}.\end{aligned}\right.$

##  3\. 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 $\varepsilon$; 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
$\varepsilon$. 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.1. First order SI IMEX scheme

In the following, we start with a first-order SI IMEX time discretization for
(2.17), while keeping space continuous:

(3.1a) $\frac{\rho^{n+1}-\rho^{n}}{\Delta t}+\nabla\cdot\mathbf{q}^{n}=0,$
(3.1b) $\frac{\mathbf{q}^{n+1}-\mathbf{q}^{n}}{\Delta
t}+\nabla\cdot\left(\frac{\mathbf{q}^{n}\otimes\mathbf{q}^{n}}{\rho^{n}}-\mathbf{B}^{n}\otimes\mathbf{B}^{n}+\left(\frac{1}{2}{\|\mathbf{B}^{n}\|}^{2}+p^{n}\right)\mathbf{I}\right)+\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\nabla
p^{n+1}=\mathbf{0},$ (3.1c) $\frac{\mathbf{A}^{n+1}-\mathbf{A}^{n}}{\Delta
t}+\left(\nabla\times{\bf
A}^{n}\right)\times\frac{\mathbf{q}^{n}}{\rho^{n}}=\mathbf{0},$ (3.1d)
$\frac{E^{n+1}-E^{n}}{\Delta
t}+\nabla\cdot\left(\frac{E^{n}+p^{n}}{\rho^{n+1}}\mathbf{q}^{n+1}+\varepsilon^{2}\left(\frac{1}{2}{\|\mathbf{B}^{n}\|}^{2}\frac{\mathbf{q}^{n}}{\rho^{n}}-\left(\frac{\mathbf{q}^{n}}{\rho^{n}}\cdot\mathbf{B}^{n}\right)\mathbf{B}^{n}\right)\right)=0,$
(3.1e) ${\bf B}^{n+1}-\nabla\times{\bf A}^{n+1}=\mathbf{0},$

where $\mathbf{q}=\rho{\bf u}$. The updating flow chart based on the SI IMEX
scheme (3.1) is as follows:

  * •

Step 1. Update $\rho^{n+1}$ and ${\bf A}^{n+1}$ from (3.1a) and (3.1c)
respectively, then update ${\bf B}^{n+1}$ from (3.1e).

  * •

Step 2. We rewrite (3.1b) as

(3.2) $\mathbf{q}^{n+1}=\mathbf{q}^{*}-\Delta
t\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\nabla p^{n+1},$

with

$\mathbf{q}^{*}=\mathbf{q}^{n}-\Delta
t\nabla\cdot\left(\frac{\mathbf{q}^{n}\otimes\mathbf{q}^{n}}{\rho^{n}}-\mathbf{B}^{n}\otimes\mathbf{B}^{n}+\left(\frac{1}{2}{\|\mathbf{B}^{n}\|}^{2}+p^{n}\right)\mathbf{I}\right).$

Substituting $\mathbf{q}^{n+1}$ into (3.1d), we can further get

(3.3) $E^{n+1}=E^{*}+\Delta
t^{2}\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\nabla\cdot(H^{n}\nabla
p^{n+1}),$

where $H^{n}=(E^{n}+p^{n})/\rho^{n+1}$ and

$E^{*}=E^{n}-\Delta t\nabla\cdot(H^{n}\mathbf{q}^{*})-\Delta
t\varepsilon^{2}\nabla\cdot\left(\frac{1}{2}{\|\mathbf{B}^{n}\|}^{2}\frac{\mathbf{q}^{n}}{\rho^{n}}-\left(\frac{\mathbf{q}^{n}}{\rho^{n}}\cdot\mathbf{B}^{n}\right)\mathbf{B}^{n}\right).$

Corresponding to the pressure in the incompressible limit, we introduce a
pressure perturbation $p^{n+1}_{2}$, defined as

(3.4) $p^{n+1}_{2}=(p^{n+1}-\bar{p}^{n})/\varepsilon^{2},$

where $\bar{p}^{n}=\int_{\Omega}p^{n}d{\bf x}/|\Omega|$, and we replace
$E^{n+1}$ by
$p^{n+1}/(\gamma-1)+\varepsilon^{2}({\|\mathbf{q}^{n}\|}^{2}/(2\rho^{n})+{\|\mathbf{B}^{n}\|}^{2}/2)$
from (3.3). We rewrite (3.3) as

(3.5) $\frac{\varepsilon^{2}}{\gamma-1}p^{n+1}_{2}-\Delta
t^{2}(1-\varepsilon^{2})\nabla\cdot(H^{n}\nabla p^{n+1}_{2})=E^{**},$

where
$E^{**}=E^{*}-\bar{p}^{n}/(\gamma-1)-\varepsilon^{2}({\|\mathbf{q}^{n}\|}^{2}/(2\rho^{n})+{\|\mathbf{B}^{n}\|}^{2}/2)$.
We can now update $p_{2}^{n+1}$ from the elliptic equation (3.5).

  * •

Step 3. We update $\mathbf{q}^{n+1}$ from (3.2) by using (3.4), and then
$E^{n+1}$ from (3.1d).

As described in [8], we may set $\nabla\cdot(\alpha p^{n}{\bf I})$ and
$(1-\alpha\varepsilon^{2})\nabla p^{n+1}/\varepsilon^{2}$ to replace
$\nabla\cdot(p^{n}{\bf I})$ and $(1-\varepsilon^{2})\nabla
p^{n+1}/\varepsilon^{2}$ in (3.1b) with $\alpha=1$ for $\varepsilon<1$ and
$\alpha=1/\varepsilon^{2}$ for $\varepsilon\geq 1$. Note that if
$\varepsilon\geq 1$ the implicit pressure contribution in (3.1b) vanishes, so
the momentum $\mathbf{q}^{n+1}$ is evaluated explicitly. With updated
$\rho^{n+1}$, ${\bf A}^{n+1}$, ${\bf B}^{n+1}$, and $\mathbf{q}^{n+1}$,
$E^{n+1}$ can also be updated in an explicit way from (3.1d). In one
dimension, the divergence-free condition is satisfied automatically, so there
is no need to apply a CT method. We can replace (3.1c) and (3.1e) with

(3.6) $\frac{\mathbf{B}^{n+1}-\mathbf{B}^{n}}{\Delta
t}+\nabla\cdot\left(\frac{\mathbf{q}^{n}\otimes\mathbf{B}^{n}}{\rho^{n}}-\frac{\mathbf{B}^{n}\otimes\mathbf{q}^{n}}{\rho^{n}}\right)=0.$

The first step in the flow chart becomes: update $\rho^{n+1}$ and ${\bf
B}^{n+1}$ from (3.1a) and (3.6) respectively, and other steps remain
unchanged.

###  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:

$\left\\{\begin{aligned} &U_{t}=\mathcal{H}(U,U,{\bf B}),\\\ &{\bf
B}=\nabla\times{\bf A},\end{aligned}\right.$

where $U={(\rho,\mathbf{q},{\bf A},E)}^{T}$ and
$\mathcal{H}:\mathbb{R}^{n}\times\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{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 ${\bf B}$, and
subscript “I” to express an implicit treatment for the second $U$. For the
system (2.17), the function $\mathcal{H}(U_{E},U_{I},{\bf B}_{E})$ is defined
as

$\displaystyle\mathcal{H}(U_{E},U_{I},{\bf B}_{E})=$
$\displaystyle-\nabla\cdot\begin{pmatrix}\mathbf{q}_{E}\\\
\frac{\mathbf{q}_{E}\otimes\mathbf{q}_{E}}{\rho_{E}}-\mathbf{B}_{E}\otimes\mathbf{B}_{E}+(\frac{1}{2}{\|\mathbf{B}_{E}\|}^{2}+p_{E})\mathbf{I}\\\
\mathbf{0}\\\ 0\\\ \end{pmatrix}-\begin{pmatrix}0\\\ \mathbf{0}\\\
(\nabla\times{\bf A}_{E})\times\frac{\mathbf{q}_{E}}{\rho_{E}}\\\
0\end{pmatrix}$ $\displaystyle-\nabla\cdot\begin{pmatrix}0\\\ \mathbf{0}\\\
\mathbf{0}\\\
\varepsilon^{2}(\frac{1}{2}{\|\mathbf{B}_{E}\|}^{2}\frac{\mathbf{q}_{E}}{\rho_{E}}-(\frac{\mathbf{q}_{E}}{\rho_{E}}\cdot\mathbf{B}_{E})\mathbf{B}_{E})\end{pmatrix}-\nabla\cdot\begin{pmatrix}0\\\
(1-\varepsilon^{2})p_{I,2}\mathbf{I}\\\ \mathbf{0}\\\
H_{I}\mathbf{q}_{I}\end{pmatrix}$ $\displaystyle=$
$\displaystyle-\nabla\cdot\mathcal{F}_{E}^{1}-\mathcal{G}_{E}-\nabla\cdot\mathcal{F}_{E}^{2}-\nabla\cdot\mathcal{F}_{SI},$

with $H_{I}=(E_{E}+p_{E})/\rho_{I}$ and
$p_{I,2}=(p_{I}-\bar{p}_{E})/\varepsilon^{2}$.

To obtain high-order accuracy in time, we can apply a multi-stage IMEX-RK time
discretization with a double Butcher $tableau$,

$\begin{array}[]{c|c}\tilde{c}&\tilde{A}\\\ \hline\cr\vspace{-0.25cm}\hfil\\\
&\tilde{b}^{T}\end{array}\ \ \ \ \ \qquad\begin{array}[]{c|c}{c}&{A}\\\
\hline\cr\vspace{-0.25cm}\hfil\\\ &{b^{T}}\end{array},$

where $\tilde{A}=(\tilde{a}_{ij})$ is an $s\times s$ matrix for an explicit
scheme, with $\tilde{a}_{ij}=0$ for $j\geq i$, and $A=({a}_{ij})$ is an
$s\times 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 vectors
$\tilde{c}=(\tilde{c}_{1},...,\tilde{c}_{s})^{T}$,
$\tilde{b}=(\tilde{b}_{1},...,\tilde{b}_{s})^{T}$, $c=(c_{1},...,c_{s})^{T}$,
and $b=(b_{1},...,b_{s})^{T}$ complete the characterization of the scheme. The
coefficients $\tilde{c}$ and $c$ are given by the usual relation

$\displaystyle\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a}_{ij},\ \ \
c_{i}=\sum_{j=1}^{i}a_{ij}.$

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_{s}^{T}A$, with $e_{s}=(0,\cdots 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_{E}^{(0)}=U_{I}^{(0)}=U^{n}$ and ${\bf
B}_{E}^{(0)}={\bf B}^{n}$, for inner stages $i=1\text{ to }s$:

  * •

First update the solution $U_{E}^{(i)}$ for the explicit part

(3.7) $U_{E}^{(i)}=U^{n}+\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\mathcal{H}(U_{E}^{(j)},U_{I}^{(j)},{\bf
B}_{E}^{(j)}),\quad{\bf B}_{E}^{(i)}=\nabla\times{\bf A}_{E}^{(i)}.$

  * •

Update the known values for the implicit part $U_{*}^{(i)}$, where

(3.8) $U_{*}^{(i)}=U^{n}+\Delta
t\sum^{i-1}_{j=1}a_{ij}\mathcal{H}(U_{E}^{(j)},U_{I}^{(j)},{\bf
B}_{E}^{(j)}),$

and then solve

(3.9) $U_{I}^{(i)}=U_{*}^{(i)}+\Delta
ta_{ii}\mathcal{H}(U_{E}^{(i)},U_{I}^{(i)},{\bf B}_{E}^{(i)}).$

  * •

Finally, compute the solutions $U^{n+1}$ and ${\bf B}^{n+1}$ at time level
$t^{n+1}$ from

(3.10) $U^{n+1}=U^{n}+\Delta
t\sum_{i=1}^{s}b_{i}\mathcal{H}(U_{E}^{(i)},U_{I}^{(i)},{\bf
B}_{E}^{(i)}),\quad{\bf B}^{n+1}=\nabla\times{\bf A}^{n+1}.$

The first order scheme (3.1) is the same as applying the following Butcher
table to (2.17):

$\begin{array}[]{c|c}0&0\\\
\hline\cr&1\end{array}\qquad\qquad\begin{array}[]{c|c}1&1\\\
\hline\cr&1\end{array},$

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 $\mathcal{F}_{E}^{1}$ 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
$\mathcal{F}_{E}^{2}$ and $\mathcal{F}_{SI}$. As for $\mathcal{G}_{E}$, which
is a Hamiltonian, a WENO procedure described in [31] will be applied.

  * •

Component-wise WENO $\nabla_{W}$ for $\nabla\cdot\mathcal{F}_{E}^{2}$ and
$\nabla\cdot\mathcal{F}_{SI}$:
Let ${\bf v}=(\rho,\mathbf{q},{\bf B},E)^{T}$ and $\mathbf{f}=\mathbf{f}({\bf
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}<\cdots<x_{N}=b$. A conservative finite
difference WENO scheme for system (2.1) can be written in the following flux-
difference form:

$(\mathbf{v}_{i})_{t}=\frac{1}{\Delta
x}(\hat{\mathbf{f}}_{i+\frac{1}{2}}-\hat{\mathbf{f}}_{i-\frac{1}{2}}),$

where $\hat{\mathbf{f}}_{i\pm\frac{1}{2}}$ are the numerical fluxes, which can
be reconstructed by the following procedures:

    1. (1)

Compute the physical flux at each grid point: $\mathbf{f}_{i}=\mathbf{f}({\bf
v}_{i})$.

    2. (2)

Perform a Lax–Friedrichs flux vector splitting for each component of the
physical variables to get

(3.11)
$\mathbf{f}^{\pm}_{i}=\frac{1}{2}(\mathbf{f}_{i}\pm\alpha\mathbf{v}_{i}),$

with $\alpha=\mathop{\max}\limits_{0\leq k\leq N}(|{\bf
u}_{k}|+\hat{c}_{f}^{(k)})$, where

(3.12) $\left\\{\begin{aligned} &\hat{a}^{(k)}:={\rm
min}\left(\frac{1}{\varepsilon},1\right)\sqrt{\frac{\gamma
p_{k}}{\rho_{k}}},\quad\hat{c}_{a}^{(k)}:=\sqrt{\frac{{({\bf B}_{k}\cdot{\bf
n})}^{2}}{\rho_{k}}},\\\
&\hat{c}_{f,s}^{(k)}:={\left\\{\frac{1}{2}\left[\left(\hat{a}^{(k)}\right)^{2}+\frac{{\|{\bf
B}_{k}\|}^{2}}{\rho_{k}}\pm\sqrt{{\left(\left(\hat{a}^{(k)}\right)^{2}+\frac{{\|{\bf
B}_{k}\|}^{2}}{\rho_{k}}\right)}^{2}-\left(2{\hat{a}^{(k)}}{\hat{c}^{(k)}_{a}}\right)^{2}}\right]\right\\}}^{\frac{1}{2}}.\end{aligned}\right.$

    3. (3)

Perform a finite difference WENO reconstruction to obtain upwind and downwind
numerical fluxes:

$\hat{\mathbf{f}}^{+}_{i+\frac{1}{2}}=\Phi_{\rm
WENO5}(\mathbf{f}^{+}_{i-2},\mathbf{f}^{+}_{i-1},\mathbf{f}^{+}_{i},\mathbf{f}^{+}_{i+1},\mathbf{f}^{+}_{i+2}),\quad\hat{\mathbf{f}}^{-}_{i+\frac{1}{2}}=\Phi_{\rm
WENO5}(\mathbf{f}^{-}_{i+3},\mathbf{f}^{-}_{i+2},\mathbf{\mathbf{f}}^{-}_{i+1},\mathbf{f}^{-}_{i},\mathbf{f}^{-}_{i-1}).$

Here $\Phi_{\rm WENO5}$ denotes a 5th-order WENO reconstruction [31], while
other orders can also be used. Then we set

$\hat{\mathbf{f}}_{i+\frac{1}{2}}=\hat{\mathbf{f}}^{+}_{i+\frac{1}{2}}+\hat{\mathbf{f}}^{-}_{i+\frac{1}{2}}.$

  * •

Characteristic-wise WENO $\nabla_{CW}$ for $\nabla\cdot\mathcal{F}_{E}^{1}$:
The Jacobian matrix for the flux function
$\partial\mathbf{f}/\partial\mathbf{v}$, has a spectral decomposition of the
form $\partial\mathbf{f}/\partial\mathbf{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 $\nabla_{CW}$, we have the following procedures:

    1. (1)

Compute the physical flux at each grid point: $\mathbf{f}_{i}=\mathbf{f}({\bf
v}_{i})$.

    2. (2)

Compute the average state ${\bf v}_{i+\frac{1}{2}}$ by using an arithmetic
mean:

${\bf v}_{i+\frac{1}{2}}=\frac{1}{2}({\bf v}_{i}+{\bf v}_{i+1}).$

    3. (3)

Compute the right and left eigenvectors of the flux Jacobian matrix
$\partial\mathbf{f}/\partial\mathbf{v}$ by taking ${\bf v}={\bf
v}_{i+\frac{1}{2}}$ at $x=x_{i+\frac{1}{2}}$:

$R_{i+\frac{1}{2}}=R({\bf v}_{i+\frac{1}{2}}),\quad L_{i+\frac{1}{2}}=L({\bf
v}_{i+\frac{1}{2}}).$

    4. (4)

Project the solution and physical flux into the characteristic space:

$\mathbf{w}_{j}=L_{i+\frac{1}{2}}{\bf
v}_{j},\quad\mathbf{g}_{j}=L_{i+\frac{1}{2}}\mathbf{f}_{j},$

for all $j$ in the numerical stencil associated with $x=x_{i+\frac{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$.

    5. (5)

Perform a Lax–Friedrichs flux vector splitting for each component of the
characteristic variables to compute

$\mathbf{g}^{\pm}_{j}=\frac{1}{2}(\mathbf{g}_{j}\pm\alpha\mathbf{w}_{j}),$

where $\alpha$ has the same expression in the component-wise WENO (3.11).

    6. (6)

Perform a WENO reconstruction $\Phi_{\rm WENO5}$ on each component of
$\mathbf{g}^{\pm}_{j}$ to obtain their corresponding numerical fluxes, and
then set

$\hat{\mathbf{g}}_{i+\frac{1}{2}}=\hat{\mathbf{g}}^{+}_{i+\frac{1}{2}}+\hat{\mathbf{g}}^{-}_{i+\frac{1}{2}}.$

    7. (7)

Project the numerical flux $\hat{\mathbf{g}}_{i+\frac{1}{2}}$ back to the
physical space:

$\hat{\mathbf{f}}_{i+\frac{1}{2}}=R_{i+\frac{1}{2}}\hat{\mathbf{g}}_{i+\frac{1}{2}}.$

However, when the sonic Mach number $\varepsilon$ 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 $\varepsilon=\mathcal{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 $\varepsilon=1$,
which works well from our numerical results.

  * •

A Hamiltonian WENO reconstruction WENO-HJ $(\cdot)_{HJ}$ for
$\mathcal{G}_{E}$:
We take a 1D Hamilton–Jacobi equation as an example

(3.13) $v_{t}+\mathcal{HT}(t,x,v,v_{x})=0,$

where $\mathcal{HT}$ is the Hamiltonian. As described in [31], a semi-discrete
form for (3.13) is

$v_{t,i}=-\widehat{\mathcal{HT}}(t,x_{i},v_{i},v^{-}_{x,i},v^{+}_{x,i}),$

where $\widehat{\mathcal{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)
$\widehat{\mathcal{HT}}\left(t,x,v,u^{-},u^{+}\right)=\mathcal{HT}\left(t,x,v,\frac{u^{-}+u^{+}}{2}\right)-\beta\left(\frac{u^{+}-u^{-}}{2}\right),$

where

$\beta=\mathop{\max}\limits_{u\in I(u^{-},u^{+})}|\mathcal{HT}_{u}(t,x,v,u)|,$

and $I(u^{-},u^{+})$ is the interval between $u^{-}$ and $u^{+}$. The values
of $v^{+}_{x,i}$ and $v^{-}_{x,i}$ are obtained from

$\displaystyle v^{+}_{x,i}$ $\displaystyle=\Phi_{\rm
WENO5}\left(\frac{\Delta^{+}v_{i+2}}{\Delta x},\frac{\Delta^{+}v_{i+1}}{\Delta
x},\frac{\Delta^{+}v_{i}}{\Delta x},\frac{\Delta^{+}v_{i-1}}{\Delta
x},\frac{\Delta^{+}v_{i-2}}{\Delta x}\right),$ $\displaystyle v^{-}_{x,i}$
$\displaystyle=\Phi_{\rm WENO5}\left(\frac{\Delta^{+}v_{i-3}}{\Delta
x},\frac{\Delta^{+}v_{i-2}}{\Delta x},\frac{\Delta^{+}v_{i-1}}{\Delta
x},\frac{\Delta^{+}v_{i}}{\Delta x},\frac{\Delta^{+}v_{i+1}}{\Delta
x}\right),$

with $\Delta^{+}v_{j}=v_{j+1}-v_{j}$. Such a Hamiltonian WENO reconstruction
helps us to control unphysical oscillations both in $v_{x}$ and $v$. For
$\mathcal{G}_{E}$ in 2D, the Hamiltonian is

$\mathcal{HT}\left(t,x,y,\frac{\partial A_{z}}{\partial x},\frac{\partial
A_{z}}{\partial y}\right)=u(t,x,y)\frac{\partial A_{z}}{\partial
x}+v(t,x,y)\frac{\partial A_{z}}{\partial y},$

and the numerical Hamiltonian (3.14) can be applied along each direction for
$u(t,x,y)\frac{\partial A_{z}}{\partial x}$ and $v(t,x,y)\frac{\partial
A_{3}}{\partial y}$, respectively. Here $A_{z}$ is the third component of
${\bf 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 $i$
($1\leq i\leq s$), we first compute $U^{(i)}_{E}$ from (3.7)

$\rho^{(i)}_{E}=\rho^{n}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\nabla_{CW}\cdot\mathbf{q}^{(j)}_{E},$
$\mathbf{q}^{(i)}_{E}=\mathbf{q}^{n}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\left(\nabla_{CW}\cdot\left(\frac{\mathbf{q}^{(j)}_{E}\otimes\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}-\mathbf{B}^{(j)}_{E}\otimes\mathbf{B}^{(j)}_{E}+\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{E}\|}^{2}+p^{(j)}_{E}\right)\mathbf{I}\right)+(1-\varepsilon^{2})\nabla_{W}p^{(j)}_{I,2}\right),$
$\mathbf{A}^{(i)}_{E}=\mathbf{A}^{n}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\left(\left(\nabla\times{\bf
A}^{(j)}_{E}\right)\times\frac{\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}\right)_{HJ},$
$E^{(i)}_{E}=E^{n}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\left(\nabla_{W}\cdot\left(\varepsilon^{2}\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{E}\|}^{2}\frac{\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}-\left(\frac{\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}\cdot\mathbf{B}^{(j)}_{E}\right)\mathbf{B}^{(j)}_{E}\right)\right)+\nabla_{W}\cdot\left(\bar{H}^{(j)}\mathbf{q}^{(j)}_{I}\right)\right),$

with $\bar{H}^{(j)}=(E_{E}^{(j)}+p_{E}^{(j)})/\rho_{I}^{(j)}$. Then we obtain

${\bf B}_{E}^{(i)}=\nabla\times{\bf A}_{E}^{(i)},$

with a 4th-order central difference discretization.

  * •

Step 2. We compute $U^{(i)}_{*}$ in (3.8)

$\rho^{(i)}_{*}=\rho^{n}-\Delta
t\sum^{i-1}_{j=1}a_{ij}\nabla_{CW}\cdot\mathbf{q}^{(j)}_{E},$
$\mathbf{q}^{(i)}_{*}=\mathbf{q}^{n}-\Delta
t\sum^{i-1}_{j=1}a_{ij}\left(\nabla_{CW}\cdot\left(\frac{\mathbf{q}^{(j)}_{E}\otimes\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}-\mathbf{B}^{(j)}_{E}\otimes\mathbf{B}^{(j)}_{E}+\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{E}\|}^{2}+p^{(j)}_{E}\right)\mathbf{I}\right)+(1-\varepsilon^{2})\nabla_{W}p^{(j)}_{I,2}\right),$
$E^{(i)}_{*}=E^{n}-\Delta
t\sum^{i-1}_{j=1}a_{ij}\left(\nabla_{W}\cdot\left(\varepsilon^{2}\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{E}\|}^{2}\frac{\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}-\left(\frac{\mathbf{q}^{(j)}_{E}}{\rho^{(j)}_{E}}\cdot\mathbf{B}^{(j)}_{E}\right)\mathbf{B}^{(j)}_{E}\right)\right)+\nabla_{W}\cdot\left(\bar{H}^{(j)}\mathbf{q}^{(j)}_{I}\right)\right).$

  * •

Step 3. Solve $U_{I}^{(i)}$ from (3.9).

    1. (1)

In components, $U_{I}^{(i)}$ satisfies

(3.15a) $\rho^{(i)}_{I}=\rho^{(i)}_{*}-\Delta
ta_{ii}\nabla_{CW}\cdot\mathbf{q}^{(i)}_{E},$ (3.15b)
$\mathbf{q}^{(i)}_{I}=\mathbf{q}^{(i)}_{**}-\Delta
ta_{ii}(\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\nabla p^{(i)}_{I}),$ (3.15c)
$E^{(i)}_{I}=E^{(i)}_{**}-\Delta
ta_{ii}\nabla\cdot(\bar{H}^{(i)}\mathbf{q}^{(i)}_{I}),$

where

$\mathbf{q}^{(i)}_{**}=\mathbf{q}^{(i)}_{*}-\Delta
ta_{ii}\nabla_{CW}\cdot\left(\frac{\mathbf{q}^{(i)}_{E}\otimes\mathbf{q}^{(i)}_{E}}{\rho^{(i)}_{E}}-\mathbf{B}^{(i)}_{E}\otimes\mathbf{B}^{(i)}_{E}+\left(\frac{1}{2}{\|\mathbf{B}^{(i)}_{E}\|}^{2}+p^{(i)}_{E}\right)\mathbf{I}\right),$
$E^{(i)}_{**}=E^{(i)}_{*}-\Delta
ta_{ii}\nabla_{W}\cdot\left(\varepsilon^{2}\left(\frac{1}{2}{\|\mathbf{B}^{(i)}_{E}\|}^{2}\frac{\mathbf{q}^{(i)}_{E}}{\rho^{(i)}_{E}}-\left(\frac{\mathbf{q}^{(i)}_{E}}{\rho^{(i)}_{E}}\cdot\mathbf{B}^{(i)}_{E}\right)\mathbf{B}^{(i)}_{E}\right)\right).$

    2. (2)

We can obtain $\rho_{I}^{(i)}$ from (3.15a) directly. To solve (3.15b) and
(3.15c), we substitute $\mathbf{q}^{(i)}_{I}$ of (3.15b) into (3.15c), and
obtain

(3.16) $E^{(i)}_{I}=E^{\circ}_{I}+\left(1-\varepsilon^{2}\right)\Delta
t^{2}a^{2}_{ii}\nabla\cdot\left(\bar{H}^{(i)}\left(\frac{\nabla
p^{(i)}_{I}}{\varepsilon^{2}}\right)\right),$

with $E_{I}^{(i)}$ following the EOS
$E_{I}^{(i)}=p_{I}^{(i)}/(\gamma-1)+\varepsilon^{2}({\|\mathbf{q}_{E}^{(i)}\|}^{2}/(2\rho_{E}^{(i)})+{\|\mathbf{B}_{E}^{(i)}\|}^{2}/2)$
and $E^{\circ}_{I}=E^{(i)}_{**}-\Delta
ta_{ii}\nabla_{W}\cdot(\bar{H}^{(i)}\mathbf{q}^{(i)}_{**})$. By substituting

(3.17) $p^{(i)}_{I}=\bar{p}^{(i)}_{E}+\varepsilon^{2}p^{(i)}_{I,2}$

into (3.16), where $\bar{p}^{(i)}_{E}=\int_{\Omega}p_{E}^{(i)}d{\bf
x}/|\Omega|$, we solve

(3.18)
$\frac{\varepsilon^{2}}{\gamma-1}p^{(i)}_{I,2}-\left(1-\varepsilon^{2}\right)\Delta
t^{2}a^{2}_{ii}\nabla\cdot\left(\bar{H}^{(i)}\left(\nabla
p^{(i)}_{I,2}\right)\right)=E^{\circ\circ}_{I},$

with
$E^{\circ\circ}_{I}=E^{\circ}_{I}-\bar{p}^{(i)}_{E}/(\gamma-1)-\varepsilon^{2}({\|\mathbf{q}^{(i)}_{E}\|}^{2}/(2\rho^{(i)}_{E})+{\|\mathbf{B}^{(i)}_{E}\|}^{2}/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 $\nabla\cdot(\bar{H}^{(i)}\nabla p^{(i)}_{I,2})$. The
second-order spatial derivative is then discretized by a compact
discretization as proposed in [10].

    3. (3)

With $p^{(i)}_{I,2}$ solved from (3.18), we can update $\mathbf{q}^{(i)}_{I}$
from

$\mathbf{q}^{(i)}_{I}=\mathbf{q}^{(i)}_{**}-\Delta
ta_{ii}\left(\left(1-\varepsilon^{2}\right)\nabla_{W}p^{(i)}_{I,2}\right),$

and successively update $E^{(i)}_{I}$ from

$E^{(i)}_{I}=E^{(i)}_{**}-\Delta
ta_{ii}\nabla_{W}\cdot\left(\bar{H}^{(i)}\mathbf{q}^{(i)}_{I}\right).$

  * •

Step 4. Finally, we update the numerical solution $U^{n+1}=U^{(s)}_{I}$, and
${\bf B}^{n+1}=\nabla\times{\bf 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 (LABEL:S2_E26).

###### Proof.

We assume the following expansions of the solutions at time $t^{n}$:

(4.1) $\left\\{\begin{aligned}
&p^{n}(\mathbf{x})=p^{n}_{0}+\varepsilon^{2}p^{n}_{2}(\mathbf{x}),\quad\rho^{n}(\mathbf{x})=\rho^{n}_{0}(\mathbf{x})+\mathcal{O}(\varepsilon),\\\
&\mathbf{u}^{n}(\mathbf{x})=\mathbf{u}^{n}_{0}(\mathbf{x})+\mathcal{O}(\varepsilon^{2}),\quad\mathbf{A}^{n}(\mathbf{x})=\mathbf{A}^{n}_{0}(\mathbf{x})+\mathcal{O}(\varepsilon),\\\
&\mathbf{B}^{n}(\mathbf{x})=\mathbf{B}^{n}_{0}(\mathbf{x})+\mathcal{O}(\varepsilon),&\end{aligned}\right.$

with periodic or no-slip boundary conditions, and well-prepared conditions
$p_{0}^{n}=(\gamma-1)E^{n}_{0}={\rm Const}$,
$\nabla\cdot\mathbf{u}^{n}_{0}(\mathbf{x})=0$ and $\nabla\cdot{\bf B}^{n}({\bf
x})=0$. No matter what value ${\bf A}^{n+1}$ takes, we can obtain
$\nabla\cdot{\bf B}^{n+1}=\nabla\cdot(\nabla\times{\bf A}^{n+1})=0$ from
(3.1e). We plug (4.1) into the semi-discrete scheme (3.1) and EOS (3.3). For
the leading order $\mathcal{O}(\varepsilon^{-2})$, we obtain

$\nabla p_{0}^{n+1}=0,$

i.e. $p_{0}^{n+1}$ is constant in space. Equating to zero the
$\mathcal{O}(\varepsilon^{0})$ terms, we have the following equations:

(4.2a) $\frac{\rho_{0}^{n+1}-\rho_{0}^{n}}{\Delta
t}+\nabla\cdot(\rho_{0}^{n}{\bf u}_{0}^{n})=0,$ (4.2b)
$\frac{\rho_{0}^{n+1}\mathbf{u}_{0}^{n+1}-\rho_{0}^{n}\mathbf{u}_{0}^{n}}{\Delta
t}+\nabla\cdot\left(\rho_{0}^{n}{\bf u}_{0}^{n}\otimes{\bf
u}_{0}^{n}-\mathbf{B}_{0}^{n}\otimes\mathbf{B}_{0}^{n}+\frac{1}{2}{\|\mathbf{B}_{0}^{n}\|}^{2}\mathbf{I}\right)+\nabla
p^{n+1}_{2}=\mathbf{0},$ (4.2c)
$\frac{\mathbf{A}_{0}^{n+1}-\mathbf{A}_{0}^{n}}{\Delta t}+(\nabla\times{\bf
A}_{0}^{n})\times{\bf u}_{0}^{n}=\mathbf{0},$ (4.2d)
$\frac{E_{0}^{n+1}-E_{0}^{n}}{\Delta
t}+\nabla\cdot\left(\bar{H}_{0}^{n}\left(\rho_{0}^{n+1}{\bf
u}_{0}^{n+1}\right)\right)=0,$ (4.2e) ${\bf B}_{0}^{n+1}-\nabla\times{\bf
A}_{0}^{n+1}=\mathbf{0},$ (4.2f) $E_{0}^{n+1}=\frac{p_{0}^{n+1}}{\gamma-1},$

with

(4.3)
$\bar{H}_{0}^{n}=\frac{p_{0}^{n}+E_{0}^{n}}{\rho_{0}^{n+1}}=\frac{\gamma}{\gamma-1}\frac{p_{0}^{n}}{\rho_{0}^{n+1}}.$

Due to (4.2f) and (4.3), it is easy to check that (4.2d) is equivalent to

(4.4) $\frac{p_{0}^{n+1}-p_{0}^{n}}{\Delta t}+\gamma p_{0}^{n}\nabla\cdot{\bf
u}_{0}^{n+1}=0.$

Taking into account a spatial domain $\Omega$ with periodic or no-slip
boundary conditions and integrating (4.4) over $\Omega$, we further get
$\int_{\Omega}\nabla\cdot{\bf u}_{0}^{n+1}d{\bf x}=\int_{\partial\Omega}{\bf
u}_{0}^{n+1}\cdot{\bf n}ds=0$, where ${\bf n}$ is the unit outward normal
vector along $\partial\Omega$. This implies $p_{0}^{n+1}=p^{n}_{0}={\rm
Const}$. A direct result $\nabla\cdot{\bf u}_{0}^{n+1}=0$ can be obtained
through combing (4.4) with $p_{0}^{n+1}=p^{n}_{0}$. The equations
(4.2a)-(4.2c) and (4.2e) correspond to consistent discretizations for other
equations in (LABEL:S2_E26), so that the AP property is obtained. ∎

###  4.2. 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:

######  Theorem 4.2.

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 $\Omega\subset\mathbb{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 $(\rho^{0}({\bf
x}),\rho^{0}({\bf x}){\bf u}^{0}({\bf x}),p^{0},$ ${\bf A}^{0}({\bf x}),{\bf
B}^{0}({\bf x}))^{T}$ are well-prepared, namely in the form of (4.1). If we
denote $(\rho^{1}({\bf x};\varepsilon),\rho^{1}({\bf x};\varepsilon){\bf
u}^{1}({\bf x};\varepsilon),$ $p^{1}({\bf x};\varepsilon),{\bf A}^{1}({\bf
x};\varepsilon),{\bf B}^{1}({\bf x};\varepsilon))^{T}$ as the numerical
solution after one-time step, then, with $p_{*}$ being a constant, we have:

(4.5) $\lim_{\varepsilon\to 0}p^{1}({\bf
x};\varepsilon)=p_{*},\quad\lim_{\varepsilon\to 0}\nabla\cdot{\bf u}^{1}({\bf
x};\varepsilon)=0,\quad\nabla\cdot{\bf B}^{1}({\bf x};\varepsilon)=0.$

Furthermore, let ${\bf V}_{inc}({\bf x};t)=(\rho_{inc}({\bf
x};t),\rho_{inc}({\bf x};t){\bf u}_{inc}({\bf x};t),p_{inc}({\bf x};t),{\bf
A}_{inc}({\bf x};t),{\bf B}_{inc}({\bf x};t))^{T}$ be the exact solution of
the incompressible MHD equations (LABEL:S2_E26) with the same initial data,
one has the following one-step error estimate

(4.6) $\lim_{\varepsilon\to 0}{\bf V}^{1}({\bf x};\varepsilon)={\bf
V}_{inc}({\bf x};\Delta t)+\mathcal{O}(\Delta t^{r+1}),$

i.e. the scheme is AA.

###### Proof.

We consider the first step from $t^{0}=0$ to $t^{1}=\Delta 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):

$\left\\{\begin{aligned}
&p^{0}(\mathbf{x})=p_{*}+\varepsilon^{2}p^{0}_{2}(\mathbf{x}),\quad\rho^{0}(\mathbf{x})=\rho^{0}_{inc}+\mathcal{O}(\varepsilon),\\\
&\mathbf{u}^{0}(\mathbf{x})=\mathbf{u}^{0}_{inc}+\mathcal{O}(\varepsilon^{2}),\quad\mathbf{A}^{0}(\mathbf{x})=\mathbf{A}^{0}_{inc}+\mathcal{O}(\varepsilon),\quad\mathbf{B}^{0}(\mathbf{x})=\mathbf{B}^{0}_{inc}+\mathcal{O}(\varepsilon),\end{aligned}\right.$

where

$\displaystyle\rho^{0}_{inc}:=\rho_{inc}({\bf
x},0),\quad\mathbf{u}^{0}_{inc}:=\mathbf{u}_{inc}({\bf
x},0),\quad\mathbf{A}^{0}_{inc}:=\mathbf{A}_{inc}({\bf
x},0),\quad\mathbf{B}^{0}_{inc}:=\mathbf{B}_{inc}({\bf x},0).$

By well-prepared assumptions we have: $p_{inc}({\bf x},0):=p_{*}$ is a
constant independent of time and space, $\nabla\cdot{\bf u}_{inc}^{0}=0$ and
$\nabla\cdot{\bf B}^{0}({\bf x})=0$.

Now we consider a formal $\varepsilon$-expansion of the quantities:

$U_{I}^{(i)}=(\rho_{I}^{(i)},\mathbf{q}_{I}^{(i)},{\bf
A}_{I}^{(i)},E_{I}^{(i)}),\quad
U_{E}^{(i)}=(\rho_{E}^{(i)},\mathbf{q}_{E}^{(i)},{\bf
A}_{E}^{(i)},E_{E}^{(i)}),$

and ${\bf B}_{E}^{(i)}$, with $\mathbf{q}_{I}^{(i)}=\rho_{I}^{(i)}{\bf
u}_{I}^{(i)}$ and $\mathbf{q}_{E}^{(i)}=\rho_{E}^{(i)}{\bf u}_{E}^{(i)}$. For
example, the expansions of the density and pressure are

(4.7) $\left\\{\begin{aligned}
&\rho_{I}^{(i)}=\rho_{0,I}^{(i)}+\varepsilon\rho_{1,I}^{(i)}+\cdots,\quad&\rho_{E}^{(i)}&=\rho_{0,E}^{(i)}+\varepsilon\rho_{1,E}^{(i)}+\cdots,\\\
&p_{I}^{(i)}=p_{0,I}^{(i)}+\varepsilon^{2}p_{2,I}^{(i)}+\cdots,\quad&p_{E}^{(i)}&=p_{0,E}^{(i)}+\varepsilon^{2}p_{2,E}^{(i)}+\cdots.\end{aligned}\right.$

To prove the theorem, we use mathematical induction.

  * •

AA property for internal stages $i=1,\cdots,s.$ When $i=1$, it leads to the
same AP analysis for the scheme (3.1) with $\Delta t$ replaced by
$a_{11}\Delta t$. To prove the result for $i>1$, we make use of the induction
hypothesis, assuming the property holds for $j\leq i-1$, and prove that it
still holds for $j=i$. For $j=1,\cdots,i-1$, we have

(4.8) $p_{0,E}^{(j)}=p_{0,I}^{(j)}=p_{*},\quad
E_{0,E}^{(j)}=E_{0,I}^{(j)}=\frac{p_{*}}{\gamma-1},\quad\nabla\cdot{\bf
u}_{0,I}^{(j)}=0,\quad\nabla\cdot{\bf B}_{E}^{(j)}=0.$

Now we insert the expansions (4.7) into the explicit step in (3.7), and from
the energy equation we get

(4.9) $E_{0,E}^{(i)}=E_{inc}^{0}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\nabla\cdot\left(\bar{H}_{0}^{(j)}\mathbf{q}_{0,I}^{(j)}\right),$

with $E_{inc}^{0}=p_{*}/(\gamma-1)$ and for $j=1,\cdots,i-1$,

(4.10)
$\bar{H}_{0}^{(j)}=\frac{E_{0,E}^{(j)}+p_{0,E}^{(j)}}{\rho_{0,I}^{(j)}}=\frac{\gamma}{\gamma-1}\frac{p_{*}}{\rho_{0,I}^{(j)}}.$

Now substituting (4.8), (4.10), and $E_{0,E}^{(i)}=p_{0,E}^{(i)}/(\gamma-1)$
into (4.9), we obtain

$p_{0,E}^{(i)}=p_{*}-\Delta t\gamma
p_{*}\sum^{i-1}_{j=1}\tilde{a}_{ij}\nabla\cdot{\bf u}_{0,I}^{(j)}=p_{*},$

and $E_{0,E}^{(i)}=p_{*}/(\gamma-1)$ is also a constant for the stage $i$.
From (3.7), we have

(4.11) $\left\\{\begin{aligned} &\mathbf{A}^{(i)}_{E}=\mathbf{A}^{n}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\left(\left(\nabla\times{\bf
A}^{(j)}_{E}\right)\times\mathbf{u}^{(j)}_{E}\right),\\\ &{\bf
B}_{E}^{(i)}=\nabla\times{\bf A}_{E}^{(i)}.\end{aligned}\right.$

From (LABEL:416), we get $\nabla\cdot{\bf
B}_{E}^{(i)}=\nabla\cdot\left(\nabla\times{\bf A}_{E}^{(i)}\right)=0$. From
(3.7), in the order of $\mathcal{O}(1)$, the density and momentum equations
yield

(4.12) $\left\\{\begin{aligned} &\rho_{0,E}^{(i)}=\rho_{inc}^{0}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\nabla\cdot\mathbf{q}_{0,E}^{(j)},\\\
&\mathbf{q}^{(i)}_{0,E}=\mathbf{q}^{0}_{inc}-\Delta
t\sum^{i-1}_{j=1}\tilde{a}_{ij}\left(\nabla\cdot\left(\rho_{0,E}^{(j)}\mathbf{u}_{0,E}^{(j)}\otimes\mathbf{u}_{0,E}^{(j)}-\mathbf{B}^{(j)}_{0,E}\otimes\mathbf{B}^{(j)}_{0,E}+\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{0,E}\|}^{2}\right)\mathbf{I}\right)+\nabla
p^{(j)}_{2,I}\right),\end{aligned}\right.$

with $\mathbf{q}^{0}_{inc}=(\rho{\bf u})^{0}_{inc}$ and $\nabla
p_{0,E}^{(j)}=0$ for $j=1,\cdots,i-1$.
Similarly, inserting expansions (4.7) into (3.8), in the order of
$\mathcal{O}(1)$, we get

(4.13) $\left\\{\begin{aligned} &\rho_{0,*}^{(i)}=\rho_{inc}^{0}-\Delta
t\sum^{i-1}_{j=1}{a}_{ij}\nabla\cdot\mathbf{q}_{0,E}^{(j)},\\\
&\mathbf{q}^{(i)}_{0,*}=\mathbf{q}^{0}_{inc}-\Delta
t\sum^{i-1}_{j=1}{a}_{ij}\left(\nabla\cdot\left(\rho_{0,E}^{(j)}\mathbf{u}_{0,E}^{(j)}\otimes\mathbf{u}_{0,E}^{(j)}-\mathbf{B}^{(j)}_{0,E}\otimes\mathbf{B}^{(j)}_{0,E}+\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{0,E}\|}^{2}\right)\mathbf{I}\right)+\nabla
p^{(j)}_{2,I}\right),\end{aligned}\right.$

where from (3.17) and (4.8), it follows $\nabla\bar{p}^{(j)}_{E}=0$ for
$j=1,\cdots,i-1$. Using (4.10), we get

(4.14) $E_{0,*}^{(i)}=E_{inc}^{0}-\Delta
t\sum^{i-1}_{j=1}{a}_{ij}\nabla\cdot\left(\bar{H}_{0}^{(j)}\mathbf{q}_{0,I}^{(j)}\right)=\frac{p_{*}}{\gamma-1}-\frac{\gamma
p_{*}}{\gamma-1}\Delta t\sum^{i-1}_{j=1}{a}_{ij}\nabla\cdot{\bf
u}_{0,I}^{(j)}.$

Therefore, from (4.8) and (4.14), we get

$E_{0,*}^{(i)}=\frac{p_{*}}{\gamma-1}.$

Now from (4.14) and (3.9), it follows for the energy equation

$E^{(i)}_{0,I}=E^{(i)}_{0,*}-\Delta
ta_{ii}\nabla\cdot(\bar{H}^{(i)}_{0}\mathbf{q}^{(i)}_{0,I})=\frac{p_{*}}{\gamma-1}-\frac{\gamma
p_{*}}{\gamma-1}\Delta t{a}_{ii}\nabla\cdot{\bf u}_{0,I}^{(i)}$

Considering the EOS
$E_{I}^{(i)}=p_{I}^{(i)}/(\gamma-1)+\varepsilon^{2}({\|\mathbf{q}_{E}^{(i)}\|}^{2}/(2\rho_{E}^{(i)})+{\|\mathbf{B}_{E}^{(i)}\|}^{2}/2)$,
to zeroth order in $\varepsilon$, we get
$E_{0,I}^{(i)}=p_{0,I}^{(i)}/(\gamma-1)$, and we obtain for the pressure

$p^{(i)}_{0,I}=p_{*}-\gamma p_{*}\Delta t{a}_{ii}\nabla\cdot{\bf
u}_{0,I}^{(i)}.$

Integrating it over the spatial bounded domain $\Omega$, and assuming periodic
or compact support boundary conditions, we first obtain $p^{(i)}_{0,I}=p_{*}$,
and then we get $\nabla\cdot{\bf u}_{0,I}^{(i)}=0$ at the stage $i$.

Finally, from (3.9), considering (4.13), we get for the density and momentum:

(4.15) $\left\\{\begin{aligned} &\rho_{0,I}^{(i)}=\rho_{inc}^{0}-\Delta
t\sum^{i}_{j=1}{a}_{ij}\nabla\cdot\mathbf{q}_{0,E}^{(j)},\\\
&\mathbf{q}^{(i)}_{0,I}=\mathbf{q}^{0}_{inc}-\Delta
t\sum^{i}_{j=1}{a}_{ij}\left(\nabla\cdot\left(\rho_{0,E}^{(j)}\mathbf{u}_{0,E}^{(j)}\otimes\mathbf{u}_{0,E}^{(j)}-\mathbf{B}^{(j)}_{0,E}\otimes\mathbf{B}^{(j)}_{0,E}+\left(\frac{1}{2}{\|\mathbf{B}^{(j)}_{0,E}\|}^{2}\right)\mathbf{I}\right)+\nabla
p^{(j)}_{2,I}\right).\end{aligned}\right.$

Then equations (LABEL:416), (4.12), and (4.15), with $\nabla\cdot{\bf
B}_{E}^{(i)}=0$, $p_{0,E}^{(i)}=p_{0,I}^{(i)}=p_{*}$, and the divergence-free
leading order velocity, i.e. $\nabla\cdot{\bf u}_{0,I}^{(i)}=0$, provide a
multi-stage SI IMEX-RK discretization of the limiting system (LABEL:S2_E26)
for the internal stage $i$.

  * •

AA property for the updated numerical solution. Assuming that SI IMEX-RK
scheme in Section 3.2 is SA, then the numerical solution coincides with the
last internal stage $s$, and by setting $i=s$, we get

(4.16) $p_{0}^{1}=p_{0,I}^{(s)}=p_{*},\quad\nabla\cdot{\bf
u}_{0}^{1}=\nabla\cdot{\bf u}_{0,I}^{(s)}=0,$

and for the magnetic field

(4.17) $\nabla\cdot{\bf B}_{E}^{1}=\nabla\cdot(\nabla\times{\bf
A}_{E}^{(1)})=0,$

namely we have (4.5). Now if we denote

${\bf V}_{inc}({\bf x};t)=(\rho_{inc}({\bf x};t),\rho_{inc}({\bf x};t){\bf
u}_{inc}({\bf x};t),p_{inc}({\bf x};t),{\bf A}_{inc}({\bf x};t),{\bf
B}_{inc}({\bf x};t))^{T}$

as the exact solution of the incompressible MHD equations (LABEL:S2_E26), with
initial data ${\bf V}_{inc}({\bf x};0)=(\rho^{0}({\bf x}),\rho^{0}({\bf
x}){\bf u}^{0}({\bf x}),p^{0},{\bf A}^{0}({\bf x}),{\bf B}^{0}({\bf x}))^{T}$,
from equations (LABEL:416), (4.12), (4.15), (4.16), and (4.17), one gets in
the low sonic Mach limit where $\varepsilon=0$, a high order SI IMEX-RK scheme
of order $r$ for the numerical solutions of (LABEL:S2_E26), that is, the
scheme (3.7)-(3.10) is AA, and (4.6) holds.

∎

##  5\. 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, 53, 31] 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 ${\rm CFL}$ number is taken as $0.25$, and the gas
constant is $\gamma=5/3$. For 1D problems, the time step is

$\Delta t=\frac{{\rm CFL}\Delta x}{\mathop{\max}\limits_{0\leq k\leq
N}\left(|u_{k}|+\hat{c}_{f}^{(k)}\right)},$

with $N+1$ computational grid points, and the definition of
$\hat{c}_{f}^{(k)}$ has been given in (3.12). For 2D problems, the time step
is

$\Delta t={\rm CFL}/\left(\frac{\mathop{\max}\limits_{0\leq k\leq N_{x},0\leq
l\leq N_{y}}\left(|u_{k,l}|+\hat{c}_{f,x}^{(k,l)}\right)}{\Delta
x}+\frac{\mathop{\max}\limits_{0\leq k\leq N_{x},0\leq l\leq
N_{y}}\left(|v_{k,l}|+\hat{c}_{f,y}^{(k,l)}\right)}{\Delta y}\right),$

with $(N_{x}+1)\times(N_{y}+1)$ computational grid points. Here, the fast
speeds in the $x$-direciton ($\hat{c}_{f,x}^{(k,l)}$) and in the $y$-direciton
($\hat{c}_{f,y}^{(k,l)}$) are defined as

$\left\\{\begin{aligned}
\hat{c}_{f,x}^{(k,l)}={\left\\{\frac{1}{2}\left[\left(\hat{a}^{(k,l)}\right)^{2}+\frac{{\|{\bf
B}_{k,l}\|}^{2}}{\rho_{k,l}}\pm\sqrt{{\left(\left(\hat{a}^{({k,l})}\right)^{2}+\frac{{\|{\bf
B}_{k,l}\|}^{2}}{\rho_{k,l}}\right)}^{2}-\left(2{\hat{a}^{({k,l})}}\sqrt{\frac{B_{x;{k,l}}^{2}}{\rho_{k,l}}}\right)^{2}}\right]\right\\}}^{\frac{1}{2}},\\\
\hat{c}_{f,y}^{(k,l)}={\left\\{\frac{1}{2}\left[\left(\hat{a}^{(k,l)}\right)^{2}+\frac{{\|{\bf
B}_{k,l}\|}^{2}}{\rho_{k,l}}\pm\sqrt{{\left(\left(\hat{a}^{({k,l})}\right)^{2}+\frac{{\|{\bf
B}_{k,l}\|}^{2}}{\rho_{k,l}}\right)}^{2}-\left(2{\hat{a}^{({k,l})}}\sqrt{\frac{B_{y;{k,l}}^{2}}{\rho_{k,l}}}\right)^{2}}\right]\right\\}}^{\frac{1}{2}},\end{aligned}\right.$

where

$\hat{a}^{({k,l})}=\min\left(\frac{1}{\varepsilon},1\right)\sqrt{\frac{\gamma
p_{k,l}}{\rho_{k,l}}}.$

###  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.

######  Example 5.1.

(Accuracy test) To assess the convergence of our proposed scheme, we simulate
the nonlinear circularly polarized Alfvén wave problem in 1D. The smooth
initial conditions are given as [16]:

(5.1) $\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(0,x)$
$\displaystyle\quad=(1,0,0.1\sin(2\pi x),0.1\cos(2\pi x),1,0.1\sin(2\pi
x),0.1\cos(2\pi x),0.1)\quad x\in[0,1],$

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

$\Delta t=\frac{{\rm CFL}\Delta x^{\frac{5}{3}}}{\mathop{\max}\limits_{0\leq
k\leq N}\left(|u_{k}|+\hat{c}_{f}^{(k)}\right)}.$

The sonic Mach number is $\varepsilon=1$ and we compute up to a final time
$T=1$. The numerical $L_{1}$, $L_{2}$, and $L_{\infty}$ errors and orders are
shown in Table 5.1. We can observe that a fifth-order spatial accuracy is
obtained.

Table 5.1.  Example 5.1. The $L_{1}$, $L_{2}$, and $L_{\infty}$ errors and orders for $\rho v$ with $\varepsilon=1$ at $T=1$.  $N$ | 10 | 20 | 40 | 80 | 160
---|---|---|---|---|---
$L_{1}$ error | 6.24E-04 | 2.03E-05 | 6.41E-07 | 2.01E-08 | 6.27E-10
order | – | 4.94 | 4.99 | 5.00 | 5.00
$L_{2}$ error | 6.83E-04 | 2.25E-05 | 7.11E-07 | 2.23E-08 | 6.97E-10
order | – | 4.93 | 4.98 | 5.00 | 5.00
$L_{\infty}$ error | 9.65E-04 | 3.15E-05 | 1.00E-06 | 3.15E-08 | 9.85E-10
order | – | 4.94 | 4.97 | 4.99 | 5.00

######  Example 5.2.

(Accuracy test for a range of $\varepsilon$) We take a well-prepared initial
data which is $\varepsilon-$dependent:

$\left\\{\begin{aligned} \rho(x,0)&=1+\varepsilon^{2}\sin^{2}(2\pi
x),\quad&p(x,0)&=(1+\varepsilon^{2}\sin^{2}(2\pi x))^{\gamma},\\\
u(x,0)&=\varepsilon^{2}\sin(2\pi x),\quad&B_{x}(x,0)&=0.5,\\\
v(x,0)&=\sin(2\pi x)+\varepsilon^{2}\cos(2\pi
x),\quad&B_{y}(x,0)&=(1+\varepsilon^{2})\sin(2\pi x),\\\
w(x,0)&=0,\quad&B_{z}(x,0)&=(1+\varepsilon^{2})\cos(2\pi x),\\\
\end{aligned}\right.$

on the domain $\Omega=[0,1]$ with $\gamma=1.4$ and periodic boundary
conditions. Reference solutions are computed with $N=320$. Four different
sonic Mach numbers $\varepsilon=1,10^{-2},10^{-6}$, and $0$ are taken, with a
final time $T=0.05$. Numerical errors and orders are shown in Table 5.2. For
the intermediate value of $\varepsilon=10^{-2}$, order reduction is observed,
which is a typical behavior of high-order IMEX schemes for these multi-scale
problems [8]. Around 5th-order accuracy for $\varepsilon=1,10^{-6}$, and $0$
is observed.

Table 5.2.  Example 5.2. The $L_{1}$ errors and orders for $\rho v$ with $\varepsilon=1,10^{-2},10^{-6}$, and $0$.  $N$ $\varepsilon$ |  $\varepsilon=1$ | $\varepsilon=10^{-2}$ | $\varepsilon=10^{-6}$ | $\varepsilon=0$
---|---|---|---|---
$L_{1}$ error | order |  $L_{1}$ error | order |  $L_{1}$ error | order |  $L_{1}$ error | order
10 | 3.08E-02 | – | 4.93E-04 | – | 4.82E-04 | – | 4.82E-04 | –
20 | 6.47E-03 | 2.25 | 7.18E-05 | 2.78 | 1.61E-05 | 4.90 | 1.61E-05 | 4.90
40 | 8.32E-04 | 2.96 | 3.33E-05 | 1.11 | 5.22E-07 | 4.95 | 5.22E-07 | 4.95
80 | 3.74E-05 | 4.48 | 1.01E-05 | 1.72 | 1.76E-08 | 4.89 | 1.76E-08 | 4.89
160 | 1.22E-06 | 4.94 | 3.34E-06 | 1.60 | 6.65E-10 | 4.72 | 6.65E-10 | 4.72

######  Example 5.3.

(MHD shock tube) In this example, we consider a 1D MHD shock tube in the
compressible regime where the sonic Mach number is of $\mathcal{O}(1)$. The
initial conditions are as follows

$\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(0,x)$
$\displaystyle\quad=\left\\{\begin{aligned}
&(1,0,0,0,0.75,1,0,1)\quad&x&\in[0,0.5),\\\
&(0.125,0,0,0,0.75,-1,0,0.1)\quad&x&\in[0.5,1],\end{aligned}\right.$

with a final time $T=0.1$. The specific heat ratio is set to be $\gamma=2$
with $\varepsilon=1$ on the domain $\Omega=[0,1]$. Reflective boundary
conditions are considered, and we take $N=200$. The results are shown in Fig.
5.1, and the solutions obtained by the ERK scheme with $N=2400$ are displayed
as references. From these results, it can be seen that our IMEX scheme can
capture the MHD waves very well.

Figure 5.1.  MHD shock tube solutions for Example 5.3. Left: density $\rho$;
Right: magnetic field $B_{y}$. Every two points are plotted for the IMEX
scheme.

###  5.2. Two dimensional case

A CT method described in [18] will be used to maintain the divergence-free
condition. In 2D, ${\bf B}=\nabla\times{\bf A}$ can be written as
$B_{x}=\partial A_{z}/\partial y$ and $B_{y}=-\partial A_{z}/\partial x$.
Therefore, we only need to update the third component $A_{z}$ of ${\bf A}$,
and ignore $A_{x}$ and $A_{y}$. All solutions of the following tests maintain
$\nabla\cdot{\bf B}=0$ up to machine round off errors[16].

######  Example 5.4.

(Accuracy test) The 2D version of the smooth Alfvén wave problem is obtained
by rotating the direction of propagation by an angle of $\theta$, so that the
wave now propagates along the direction ${\bf
n}=(-\cos{\theta},-\sin{\theta})$ on the domain
$\Omega=[0,1/\cos{\theta}]\times[0,1/\sin{\theta}]$ [16]. In this numerical
test, periodic boundary conditions are applied on all both directions with
$\varepsilon=1,\theta=\pi/4$, and $T=1$. Similarly, if we take

$\Delta t={\rm CFL}/\left(\frac{\mathop{\max}\limits_{0\leq k\leq N_{x},0\leq
l\leq N_{y}}\left(|u_{k,l}|+\hat{c}_{f,x}^{(k,l)}\right)}{\Delta
x^{\frac{5}{3}}}+\frac{\mathop{\max}\limits_{0\leq k\leq N_{x},0\leq l\leq
N_{y}}\left(|v_{k,l}|+\hat{c}_{f,y}^{(k,l)}\right)}{\Delta
y^{\frac{5}{3}}}\right),$

from Table 5.3 we can find around 5th-order accuracy for the momentum $\rho
u$, and 4th-order accuracy for the magnetic field component $B_{x}$ where a
fourth-order central difference discretization is used for ${\bf
B}=\nabla\times{\bf A}$.

Table 5.3.  Example 5.4. The $L_{1}$, $L_{2}$, and $L_{\infty}$ errors and orders for $\rho u$ and $B_{x}$ with $\varepsilon=1$ and $T=1$.  | $N_{x}\times N_{y}$ | $8\times 8$ | $16\times 16$ | $32\times 32$ | $64\times 64$ | $128\times 128$
---|---|---|---|---|---|---
$\rho u$ |  $L_{1}$ error | 1.54E-03 | 5.73E-05 | 1.86E-06 | 5.88E-08 | 1.84E-09
order | – | 4.75 | 4.94 | 4.99 | 5.00
$L_{2}$ error | 1.80E-03 | 6.47E-05 | 2.08E-06 | 6.55E-08 | 2.07E-09
order | – | 4.80 | 4.96 | 4.99 | 4.99
$L_{\infty}$ error | 2.55E-03 | 9.26E-05 | 2.98E-06 | 9.44E-08 | 3.00E-09
order | – | 4.79 | 4.96 | 4.98 | 4.98
$B_{x}$ |  $L_{1}$ error | 9.29E-04 | 6.36E-05 | 3.17E-06 | 1.69E-07 | 9.64E-09
order | – | 3.87 | 4.33 | 4.23 | 4.13
$L_{2}$ error | 1.08E-03 | 7.02E-05 | 3.51E-06 | 1.88E-07 | 1.07E-08
order | – | 3.94 | 4.32 | 4.23 | 4.13
$L_{\infty}$ error | 1.29E-03 | 9.69E-05 | 5.06E-06 | 2.70E-07 | 1.52E-08
order | – | 3.73 | 4.26 | 4.23 | 4.15

######  Example 5.5.

(Accuracy test for a range of $\varepsilon$) In this 2D case, we set a well-
prepared initial condition as follows

$\left\\{\begin{aligned}
\rho(x,y,0)&=1+\varepsilon^{2}\sin^{2}(2\pi(x+y)),\quad&p(x,y,0)&=(1+\varepsilon^{2}\sin^{2}(2\pi(x+y)))^{\gamma},\\\
u(x,y,0)&=\sin(2\pi(x-y))+\varepsilon^{2}\sin(2\pi(x+y)),\quad&B_{x}(x,y,0)&=-\frac{1}{\sqrt{2}}\sin(2\pi(x+y)),\\\
v(x,y,0)&=\sin(2\pi(x-y))+\varepsilon^{2}\cos(2\pi(x+y)),\quad&B_{y}(x,y,0)&=\frac{1}{\sqrt{2}}\sin(2\pi(x+y)),\\\
w(x,y,0)&=0,\quad&B_{z}(x,y,0)&=\cos(2\pi(x+y)),\\\ \end{aligned}\right.$

with the initial magnetic potential
$A_{z}(x,y,0)=\cos(2\pi(x+y))/(2\sqrt{2}\pi)$. As in Example 5.2, four
different sonic Mach numbers are considered with periodic boundary conditions
on the domain $\Omega=[0,1]\times[0,1]$. We compute the solution up to a final
time $T=0.01$ on mesh grid points of $N^{2}$. The $L_{1}$ errors and orders of
the accuracy are shown in Table 5.4. From this Table, we can see high-order
accuracy can be obtained for $\varepsilon=1,10^{-6}$, and $0$. However, under
the current mesh sizes, order degeneracy can also be found for
$\varepsilon=10^{-2}$, which is similar to the 1D case.

Table 5.4.  Example 5.5. The $L_{1}$ errors and orders for $\rho u$ with $\varepsilon=1,10^{-2},10^{-6}$, and $0$.  $N$ $\varepsilon$ |  $\varepsilon=1$ | $\varepsilon=10^{-2}$ | $\varepsilon=10^{-6}$ | $\varepsilon=0$
---|---|---|---|---
$L_{1}$ error | order |  $L_{1}$ error | order |  $L_{1}$ error | order |  $L_{1}$ error | order
8 | 8.98E-02 | – | 5.18E-03 | – | 7.65E-03 | – | 7.65E-03 | –
16 | 1.39E-02 | 2.69 | 2.09E-03 | 1.31 | 8.74E-04 | 3.13 | 8.74E-04 | 3.13
32 | 6.76E-04 | 4.36 | 2.09E-03 | – | 1.68E-05 | 5.70 | 1.68E-05 | 5.70
64 | 2.28E-05 | 4.89 | 2.41E-03 | – | 7.42E-07 | 4.50 | 7.43E-07 | 4.50
128 | 7.15E-07 | 5.00 | 6.41E-04 | 1.91 | 1.42E-08 | 5.70 | 1.40E-08 | 5.73

######  Example 5.6.

(Orszag-Tang vortex) Next we consider the Orszag–Tang vortex problem, which is
widely considered as a standard test for MHD [18, 16, 56, 3]. The problem has
smooth initial conditions

$\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(x,y,0)$ $\displaystyle=$
$\displaystyle(\gamma^{2},-\sin{(y)},\sin{(x)},0,-\sin{(y)},\sin{(2x)},0,\gamma),$

with the initial magnetic potential: $A_{z}=0.5\cos{(2x)}+\cos{(y)}$ on the
domain $\Omega=[0,2\pi]\times[0,2\pi]$. Periodic boundary conditions are
imposed on all boundaries. As time evolves, the solution forms several shock
waves and a vortex structure in the middle of the computational domain. We set
$\varepsilon=1$, $N_{x}=N_{y}=192$ and present the density at $T=0.5,T=2,T=3$,
and $T=4$ in Fig. 5.2. A slice of the pressure at $y=0.625\pi$ and $T=3$ is
shown on the right panel of Fig. 5.3. We find our IMEX scheme can successfully
capture the shocks. We do not observe significant oscillations in any of the
conserved quantities, and our results are in good agreement with those given
in [18, 16].

Figure 5.2.  Density for Orszag-Tang problem of Example 5.6 on the mesh grid
$192^{2}$.

Figure 5.3.  Example 5.6. Thermal Pressure for Orszag-Tang problem on the mesh
grid $192^{2}$. Every two points are plotted for the IMEX scheme for the
cutting plot.

######  Example 5.7.

(MHD blast wave) MHD blast wave has been commonly used to test numerical
methods for the MHD system [15, 28]. Here, we set periodic boundary conditions
on the domain $\Omega=[-0.5,0.5]\times[-0.5,0.5]$. The initial conditions are

$\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(x,y,0)$
$\displaystyle\quad=\left\\{\begin{aligned}
&(1,0,0,0,10\sin{\theta},10\cos{\theta},0,100)\quad&r&\leq 0.125,\\\
&(1,0,0,0,10\sin{\theta},10\cos{\theta},0,10)\quad&r&>0.125,\end{aligned}\right.$

with $\theta=\pi/4$ and $r=\sqrt{x^{2}+y^{2}}$. The initial magnetic potential
is given by $A_{z}(x,y,0)=5\sqrt{2}(-x+y)$ with a first order extrapolation if
outside the boundary. In Fig. 5.4, we show the density and magnetic pressure
obtained by the IMEX and ERKC schemes on a $200\times 200$ mesh with
$\varepsilon=0.9$ at $T=0.02$. The slices of the density and the magnetic
pressure along $y=-x+0.5$ are also shown in Fig. 5.5. We see that the solution
of the IMEX scheme is close to the reference ERKC solutions and no obvious
oscillations are observed. In our test, the maximum error of discrete
divergence of the magnetic field is $2.87\times 10^{-11}$ for the IMEX scheme
and $2.77\times 10^{-11}$ for the ERKC scheme, both are close to machine
precision.

Figure 5.4.  Example 5.7. MHD blast wave problem on the mesh grid $200^{2}$.
Top: IMEX; Bottom: ERKC. 35 equally spaced contours are used.

Figure 5.5.  Example 5.7 MHD blast wave problem. 1D cutting plots along
$y=-x+0.5$. Every two points are plotted for the IMEX scheme for the cutting
plots.

######  Example 5.8.

(Field loop advection) Now, we show the 2D field loop advection test [45, 49]
to assess the capability of capturing a tangential discontinuity for a
multidimensional flow. The computational domain is
$\Omega=[-1,1]\times[-0.5,0.5]$ with periodic boundary conditions. The domain
is divided with $256\times 128$ uniform cells. The initial conditions are set
as

$\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(r,0)$
$\displaystyle\quad=(1,v_{0}\cos{(\theta_{0})},v_{0}\sin{(\theta_{0})},0,\partial
A_{z}/\partial y,-\partial A_{z}/\partial x,0,1)$

with
$v_{0}=1,\cos{(\theta_{0})}=-2/\sqrt{5},\sin{(\theta_{0})}=1/\sqrt{5},r=\sqrt{x^{2}+y^{2}}$,
and

$\displaystyle A_{z}(r)=\left\\{\begin{aligned} &10^{-3}(0.3-r)\quad&r\leq
0.3,&\\\ &0\quad&otherwise.&\end{aligned}\right.$

In this test, a loop of the magnetic field is advected. In Fig. 5.6, we show
both the deviation of the density $\rho$ away from 1 and the magnetic pressure
$\|{\bf B}\|^{2}/2$ for the case of $\varepsilon=0.1$ at $T=\sqrt{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\times 192$ cells. For this problem in the
intermediate regime with $\varepsilon=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
$\varepsilon$ is of $\mathcal{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
$\varepsilon$ becomes smaller. This example demonstrates the main advantage of
an AP scheme for such multiscale problems.

Figure 5.6.  Example 5.8. Field loop advection on the mesh grid $256\times 128$. Top left: magnetic pressure (IMEX); Middle left: magnetic pressure (ERKC); Bottom left: 1D cutting plot for the magnetic pressure; Top right: deviation of density (IMEX); Middle right: deviation of density (ERKC); Bottom right: 1D cutting plot for the deviation of density. Table 5.5.  Example 5.8. The CPU cost (seconds) for the IMEX and ERKC schemes on the mesh grid $256\times 128$ at $T=\sqrt{5}$.  $\varepsilon$ | 0.5 | 0.1 | 0.05
---|---|---|---
IMEX | 7634.880 | 5687.467 | 8015.959
ERKC | 4736.466 | 19719.964 | 38666.192

######  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

$\displaystyle(\rho,u,v,w,B_{x},B_{y},B_{z},p)(x,y,0)$
$\displaystyle\quad=(\gamma,1-2\eta(x),0.1\sin{(2\pi x)},0,0.1,0,0,1),$

with a magnetic potential $A_{z}(x,y,0)=0.1y$ on the domain
$\Omega=[0,2]\times[-0.5,0.5]$, and

$\displaystyle\eta(x)=\left\\{\begin{aligned}
&\frac{1}{2}\left(1+\sin\left(16\pi\left(y+\frac{1}{4}\right)\right)\right)\quad&y\in[-\frac{9}{32},-\frac{7}{32}],&\\\
&1\quad&y\in[-\frac{7}{32},\enspace\;\frac{7}{32}],&\\\
&\frac{1}{2}\left(1-\sin\left(16\pi\left(y-\frac{1}{4}\right)\right)\right)\quad&y\in[\enspace\;\frac{7}{32},\enspace\;\frac{9}{32}],&\\\
&0\quad&otherwise.&\end{aligned}\right.$

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
$\gamma=1.4$ and $\varepsilon=10^{-6}$ on a mesh grid $N_{x}\times
N_{y}=256\times 128$. We define the local sonic Mach number as
$M_{a}=\sqrt{u^{2}+v^{2}}/\sqrt{\gamma p/\rho}$, 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 $\omega=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 [42, 8], 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.

Figure 5.7.  Magnetized Kelvin–Helmholtz instability problem for Example 5.9
on the mesh grid $256\times 128$. Left: Vorticity $\omega=v_{x}-u_{y}$; Right:
The ratio of the local sonic Mach number $M_{ratio}$.

##  6\. 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
$\varepsilon\rightarrow 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.

## References

  * [1] D. S. Balsara.  Second-order-accurate schemes for magnetohydrodynamics with divergence-free reconstruction.  The Astrophysical Journal Supplement Series, 151(1):149, 2004.
  * [2] D. S. Balsara and M. Dumbser.  Divergence-free MHD on unstructured meshes using high order finite volume schemes based on multidimensional Riemann solvers.  Journal of Computational Physics, 299:687–715, 2015.
  * [3] D. S. Balsara, M. Dumbser, and R. Abgrall.  Multidimensional HLLC Riemann solver for unstructured meshes–with application to Euler and MHD flows.  Journal of Computational Physics, 261:172–208, 2014.
  * [4] D. S. Balsara and D. S. Spicer.  A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations.  Journal of Computational Physics, 149(2):270–292, 1999.
  * [5] S. Boscarino.  Error analysis of IMEX Runge–Kutta methods derived from differential-algebraic systems.  SIAM Journal on Numerical Analysis, 45(4):1600–1621, 2007.
  * [6] S. Boscarino, F. Filbet, and G. Russo.  High order semi-implicit schemes for time dependent partial differential equations.  Journal of Scientific Computing, 68(3):975–1001, 2016.
  * [7] S. Boscarino, L. Pareschi, and G. Russo.  Implicit-explicit Runge–Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit.  SIAM Journal on Scientific Computing, 35(1):A22–A51, 2013.
  * [8] S. Boscarino, J. Qiu, G. Russo, and T. Xiong.  High Order Semi-implicit WENO Schemes for All-Mach Full Euler System of Gas Dynamics.  SIAM Journal on Scientific Computing, 44(2):B368–B394, 2022.
  * [9] S. Boscarino, J.-M. Qiu, and G. Russo.  Implicit-explicit integral deferred correction methods for stiff problems.  SIAM Journal on Scientific Computing, 40(2):A787–A816, 2018.
  * [10] 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.
  * [11] 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.
  * [12] M. Brio and C. C. Wu.  An upwind differencing scheme for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics, 75(2):400–422, 1988.
  * [13] W. Cao, F. Zeng, Z. Zhang, and G. E. Karniadakis.  Implicit-explicit difference schemes for nonlinear fractional differential equations with nonsmooth solutions.  SIAM Journal on Scientific Computing, 38(5):A3070–A3093, 2016.
  * [14] P. Cargo and G. Gallice.  Roe matrices for ideal MHD and systematic construction of Roe matrices for systems of conservation laws.  Journal of Computational Physics, 136(2):446–466, 1997.
  * [15] A. J. Christlieb, X. Feng, Y. Jiang, and Q. Tang.  A high-order finite difference WENO scheme for ideal magnetohydrodynamics on curvilinear meshes.  SIAM Journal on Scientific Computing, 40(4):A2631–A2666, 2018.
  * [16] A. J. Christlieb, X. Feng, D. C. Seal, and Q. Tang.  A high-order positivity-preserving single-stage single-step method for the ideal magnetohydrodynamic equations.  Journal of Computational Physics, 316:218–242, 2016.
  * [17] A. J. Christlieb, Y. Liu, Q. Tang, and Z. Xu.  Positivity-preserving finite difference weighted ENO schemes with constrained transport for ideal magnetohydrodynamic equations.  SIAM Journal on Scientific Computing, 37(4):A1825–A1845, 2015.
  * [18] A. J. Christlieb, J. A. Rossmanith, and Q. Tang.  Finite difference weighted essentially non-oscillatory schemes with constrained transport for ideal magnetohydrodynamics.  Journal of Computational Physics, 268:302–325, 2014.
  * [19] W. Cui, Y. Ou, and D. Ren.  Incompressible limit of full compressible magnetohydrodynamic equations with well-prepared data in 3-D bounded domains.  Journal of Mathematical Analysis and Applications, 427(1):263–288, 2015.
  * [20] W. Dai and P. R. Woodward.  A simple finite difference scheme for multidimensional magnetohydrodynamical equations.  Journal of Computational Physics, 142(2):331–369, 1998.
  * [21] W. Dai and P. R. Woodward.  On the divergence-free condition and conservation laws in numerical simulations for supersonic magnetohydrodynamical flows.  The Astrophysical Journal, 494(1):317, 1998.
  * [22] S. Dellacherie.  Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number.  Journal of Computational Physics, 229(4):978–1016, 2010.
  * [23] G. Dimarco, R. Loubère, and M.-H. Vignal.  Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit.  SIAM journal on Scientific Computing, 39(5):A2099–A2128, 2017.
  * [24] G. Dimarco and L. Pareschi.  High order asymptotic-preserving schemes for the Boltzmann equation.  Comptes Rendus Mathematique, 350(9-10):481–486, 2012.
  * [25] M. Dumbser, D. S. Balsara, M. Tavelli, and F. Fambri.  A divergence-free semi-implicit finite volume scheme for ideal, viscous, and resistive magnetohydrodynamics.  International Journal for Numerical Methods in Fluids, 89(1-2):16–42, 2019.
  * [26] C. R. Evans and J. F. Hawley.  Simulation of magnetohydrodynamic flows: A constrained transport method.  The Astrophysical Journal, 332:659–677, 1988.
  * [27] F. Filbet and S. Jin.  A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources.  Journal of Computational Physics, 229(20):7625–7648, 2010.
  * [28] J. Han and H. Tang.  An adaptive moving mesh method for two-dimensional ideal magnetohydrodynamics.  Journal of Computational Physics, 220(2):791–812, 2007.
  * [29] C. Helzel, J. A. Rossmanith, and B. Taetz.  An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations.  Journal of Computational Physics, 230(10):3803–3829, 2011.
  * [30] G. Huang, Y. Xing, and T. Xiong.  High order well-balanced asymptotic preserving finite difference WENO schemes for the shallow water equations in all Froude numbers.  Journal of Computational Physics, 463:111255, 2022.
  * [31] G.-S. Jiang and D. Peng.  Weighted ENO schemes for Hamilton–Jacobi equations.  SIAM Journal on Scientific Computing, 21(6):2126–2143, 2000.
  * [32] G.-S. Jiang and C.-C. Wu.  A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics, 150(2):561–594, 1999.
  * [33] S. Jiang, Q. Ju, and F. Li.  Incompressible limit of the compressible magnetohydrodynamic equations with periodic boundary conditions.  Communications in Mathematical Physics, 297(2):371–400, 2010.
  * [34] S. Jiang, Q. Ju, and F. Li.  Low Mach number limit for the multi-dimensional full magnetohydrodynamic equations.  Nonlinearity, 25(5):1351, 2012.
  * [35] S. Jin.  Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations.  SIAM Journal on Scientific Computing, 21:441–454, 1999.
  * [36] S. Jin.  Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review.  Lecture notes for summer school on methods and models of kinetic theory (M$\&$MKT), Porto Ercole (Grosseto, Italy), pages 177–216, 2010.
  * [37] S. Jin.  Asymptotic-preserving schemes for multiscale physical problems.  Acta Numerica, 31:415–489, 2022.
  * [38] S. Jin, M. Tang, and X. Zhang.  A spatial-temporal asymptotic preserving scheme for radiation magnetohydrodynamics in the equilibrium and non-equilibrium diffusion limit.  Journal of Computational Physics, 452:110895, 2022.
  * [39] S. Kadioglu, D. Knoll, M. Sussman, and R. Martineau.  A second order JFNK-based IMEX method for single and multi-phase flows.  In Computational Fluid Dynamics 2010, pages 549–554. Springer, 2011\.
  * [40] A. Kanevsky, M. H. Carpenter, D. Gottlieb, and J. S. Hesthaven.  Application of implicit–explicit high order Runge–Kutta methods to discontinuous-Galerkin schemes.  Journal of Computational Physics, 225(2):1753–1781, 2007.
  * [41] R. Klein.  Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics I: One-dimensional flow.  Journal of Computational Physics, 121(2):213–237, 1995.
  * [42] G. Leidi, C. Birke, R. Andrassy, J. Higl, P. Edelmann, G. Wiest, C. Klingenberg, and F. Röpke.  A finite-volume scheme for modelling compressible MHD flows at low Mach numbers in stellar interiors.  https://arxiv.org/abs/2210.01641, 2022.
  * [43] F. Li and C.-W. Shu.  Locally divergence-free discontinuous Galerkin methods for MHD equations.  Journal of Scientific Computing, 22(1):413–442, 2005.
  * [44] F. Li and L. Xu.  Arbitrary order exactly divergence-free central discontinuous Galerkin methods for ideal MHD equations.  Journal of Computational Physics, 231(6):2655–2675, 2012.
  * [45] T. Mamashita, K. Kitamura, and T. Minoshima.  SLAU2-HLLD numerical flux with wiggle-sensor for stable low mach Magnetohydrodynamics simulations.  Computers $\&$ Fluids, 231:105165, 2021.
  * [46] W. H. Matthaeus and M. R. Brown.  Nearly incompressible magnetohydrodynamics at low Mach number.  The Physics of Fluids, 31(12):3634–3644, 1988.
  * [47] F. Meng, J. Banks, W. Henshaw, and D. Schwendeman.  Fourth-order accurate fractional-step IMEX schemes for the incompressible Navier–Stokes equations on moving overlapping grids.  Computer Methods in Applied Mechanics and Engineering, 366:113040, 2020.
  * [48] T. Minoshima and T. Miyoshi.  A low-dissipation HLLD approximate Riemann solver for a very wide range of Mach numbers.  Journal of Computational Physics, 446:110639, 2021.
  * [49] T. Minoshima, T. Miyoshi, and Y. Matsumoto.  A high-order weighted finite difference scheme with a multistate approximate Riemann solver for divergence-free magnetohydrodynamic simulations.  The Astrophysical Journal Supplement Series, 242(2):14, 2019.
  * [50] O. O’Reilly, E. M. Dunham, and J. Nordström.  Simulation of wave propagation along fluid-filled cracks using high-order summation-by-parts operators and implicit-explicit time stepping.  SIAM Journal on Scientific Computing, 39(4):B675–B702, 2017.
  * [51] K. G. Powell.  An approximate Riemann solver for magnetohydrodynamics.  In Upwind and High-Resolution Schemes, pages 570–583. Springer, 1997.
  * [52] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw.  A solution-adaptive upwind scheme for ideal magnetohydrodynamics.  Journal of Computational Physics, 154(2):284–309, 1999.
  * [53] C.-W. Shu.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, pages 325–432, 1998.
  * [54] A. Susanto, L. Ivan, H. De Sterck, and C. P. Groth.  High-order central ENO finite-volume scheme for ideal MHD.  Journal of Computational Physics, 250:141–164, 2013.
  * [55] K. Takahashi and S. Yamada.  Regular and non-regular solutions of the Riemann problem in ideal magnetohydrodynamics.  Journal of Plasma Physics, 79(3):335–356, 2013.
  * [56] H. Tang and K. Xu.  A high-order gas-kinetic method for multidimensional ideal magnetohydrodynamics.  Journal of Computational Physics, 165(1):69–88, 2000.
  * [57] M. Tavelli and M. Dumbser.  A pressure-based semi-implicit space–time discontinuous Galerkin method on staggered unstructured meshes for the solution of the compressible Navier–Stokes equations at all Mach numbers.  Journal of Computational Physics, 341:341–376, 2017.
  * [58] K. Wu.  Positivity-preserving analysis of numerical schemes for ideal magnetohydrodynamics.  SIAM Journal on Numerical Analysis, 56(4):2124–2147, 2018.
  * [59] K. Wu, H. Jiang, and C.-W. Shu.  Provably positive central DG schemes via geometric quasilinearization for ideal MHD equations.  SIAM Journal on Numerical Analysis, in press, 2022.
  * [60] K. Wu and C.-W. Shu.  A provably positive discontinuous Galerkin method for multidimensional ideal magnetohydrodynamics.  SIAM Journal on Scientific Computing, 40(5):B1302–B1329, 2018.
  * [61] K. Wu and C.-W. Shu.  Provably positive high-order schemes for ideal magnetohydrodynamics: analysis on general meshes.  Numerische Mathematik, 142(4):995–1047, 2019.
  * [62] K. Wu, D. Xiu, and X. Zhong.  A WENO-Based stochastic Galerkin scheme for ideal MHD equations with random inputs.  Communications in Computational Physics, 30(2):423–447, 2021.
  * [63] Z. Xu, D. S. Balsara, and H. Du.  Divergence-free WENO reconstruction-based finite volume scheme for solving ideal MHD equations on triangular meshes.  Communications in Computational Physics, 19(4):841–880, 2016.