A new HLLD Riemann solver with Boris correction for reducing Alfv\'en speed

TL;DR

This paper introduces Boris-HLLD, a new Riemann solver for ideal MHD that incorporates Boris correction to reduce wave speeds, enabling larger timesteps and stable simulations of high Alfvén speed plasmas.

Contribution

The paper presents Boris-HLLD, a novel Riemann solver based on HLLD that applies Boris correction to improve stability and efficiency in MHD simulations with high Alfvén speeds.

Findings

Captures contact discontinuities sharply in shock tube tests.

Exhibits shock waves without overshoot using minmod limiter.

Stable for |u| up to ~0.5c at low Alfvén speeds and larger at high Alfvén speeds.

Abstract

A new Riemann solver is presented for the ideal magnetohydrodynamics (MHD) equations with the so-called Boris correction. The Boris correction is applied to reduce wave speeds, avoiding an extremely small timestep in MHD simulations. The proposed Riemann solver, Boris-HLLD, is based on the HLLD solver. As done by the original HLLD solver, (1) the Boris-HLLD solver has four intermediate states in the Riemann fan when left and right states are given, (2) it resolves the contact discontinuity, Alfv\'en waves, and fast waves, and (3) it satisfies all the jump conditions across shock waves and discontinuities except for slow shock waves. The results of a shock tube problem indicate that the scheme with the Boris-HLLD solver captures contact discontinuities sharply and it exhibits shock waves without any overshoot when using the minmod limiter. The stability tests show that the scheme is stable when $|u| \lesssim 0.5c$ for a low Alfv\'en speed ($V_A \lesssim c$), where $u$, $c$, and $V_A$ denote the gas velocity, speed of light, and Alfv\'en speed, respectively. For a high Alfv\'en speed ($V_A \gtrsim c$), where the plasma beta is relatively low in many cases, the stable region is large, $|u| \lesssim (0.6-1) c$. We discuss the effect of the Boris correction on physical quantities using several test problems. The Boris-HLLD scheme can be useful for problems with supersonic flows in which regions with a very low plasma beta appear in the computational domain.

Equations

$$\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \mathbb{F}=0,$$

$$\mathbf{U}=\left(\begin{array}[]{c}\rho\\ (1+V_{A}^{2}/c^{2})\rho\mathbf{u}\\ \mathbf{B}\\ e\end{array}\right),$$

$$\mathbb{F}=\left(\begin{array}[]{c}\rho\mathbf{u}\\ \rho\mathbf{u}\mathbf{u}+p_{T}-\mathbf{B}\mathbf{B}\\ \mathbf{u}\mathbf{B}-\mathbf{B}\mathbf{u}\\ (e+p_{T})\mathbf{u}-(\mathbf{B}\cdot\mathbf{u})\mathbf{B}\end{array}\right),$$

$$\mathbf{u}=(u,v,w{)}^T,$$

$$\mathbf{B}=(B_x,B_y,B_z{)}^T,$$

$$V_A^2=\frac{{\left|\mathbf{B}\right|}^2}{\rho },$$

$$e=\rho \frac{{\left|\mathbf{u}\right|}^2}{2}+\frac{p}{\gamma -1}+\frac{{\left|\mathbf{B}\right|}^2}{2},$$

$$p_T=p+\frac{{\left|\mathbf{B}\right|}^2}{2},$$

$${\rho }_A=\frac{|\mathbf{B}{|}^2}{c^2},$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x}=0$$

$$\mathbf{U}=\left(\begin{array}[]{c}\rho\\ (\rho+\rho_{A})u\\ (\rho+\rho_{A})v\\ (\rho+\rho_{A})w\\ B_{y}\\ B_{z}\\ e\end{array}\right),$$

$$\mathbf{F}=\left(\begin{array}[]{c}\rho u\\ \rho u^{2}+p_{T}-B_{x}^{2}\\ \rho uv-B_{x}B_{y}\\ \rho uw-B_{x}B_{z}\\ B_{y}u-B_{x}v\\ B_{z}u-B_{x}w\\ (e+p_{T})u-(\mathbf{B}\cdot\mathbf{u})B_{x}\end{array}\right).$$

$$S_{\alpha }{\mathbf{U}}_{\alpha }^{\ast }-{\mathbf{F}}_{\alpha }^{\ast }=S_{\alpha }{\mathbf{U}}_{\alpha }-{\mathbf{F}}_{\alpha },$$

$$S_{\alpha }^{\ast }{\mathbf{U}}_{\alpha }^{\ast \ast }-{\mathbf{F}}_{\alpha }^{\ast \ast }=S_{\alpha }^{\ast }{\mathbf{U}}_{\alpha }^{\ast }-{\mathbf{F}}_{\alpha }^{\ast },$$

$$S_M{\mathbf{U}}_R^{\ast \ast }-{\mathbf{F}}_R^{\ast \ast }=S_M{\mathbf{U}}_L^{\ast \ast }-{\mathbf{F}}_L^{\ast \ast },$$

$$u_L^{\ast }=u_L^{\ast \ast }=u_R^{\ast }=u_R^{\ast \ast }=S_M,$$

$$p_{TL}^{\ast }=p_{TL}^{\ast \ast }=p_{TR}^{\ast }=p_{TR}^{\ast \ast }.$$

$${\rho }_{\alpha }^{\ast \ast }={\rho }_{\alpha }^{\ast },$$

$${\left|{\mathbf{B}}_L^{\ast }\right|}^2={\left|{\mathbf{B}}_R^{\ast }\right|}^2={\left|{\mathbf{B}}_L^{\ast \ast }\right|}^2={\left|{\mathbf{B}}_R^{\ast \ast }\right|}^2={\rho }_A{c}^2,$$

$${\rho }_{\alpha }^{\ast }={\rho }_{\alpha }\frac{S_{\alpha }-u_{\alpha }}{S_{\alpha }-S_M}.$$

$$S_M=\frac{(S_R^{\prime }-u_R){\rho }_R{u}_R-(S_L^{\prime }-u_L){\rho }_L{u}_L-p_{TR}+p_{TL}}{(S_R^{\prime }-u_R){\rho }_R-(S_L^{\prime }-u_L){\rho }_L},$$

$$p_{T\alpha }^{\ast }=\frac{(S_R^{\prime }-u_R){\rho }_R{p}_{TL}-(S_L^{\prime }-u_L){\rho }_L{p}_{TR}+{\rho }_L{\rho }_R(S_R^{\prime }-u_R)(S_L^{\prime }-u_L)(u_R-u_L)}{(S_R^{\prime }-u_R){\rho }_R-(S_L^{\prime }-u_L){\rho }_L},$$

$$S_{\alpha }^{\prime }=\left(1+\frac{V_A^2}{c^2}\right)S_{\alpha }=\left(1+\frac{{\rho }_A}{{\rho }_{\alpha }}\right)S_{\alpha }=\left(1+\frac{|{\mathbf{B}}^{\ast }{|}^2}{{\rho }_{\alpha }{c}^2}\right)S_{\alpha }.$$

$$v_{\alpha }^{\ast }=v_{\alpha }-B_{y\alpha }{B}_x\frac{S_M-u_{\alpha }}{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-S_M)-B_x^2},$$

$$B_{y\alpha }^{\ast }=B_{y\alpha }\frac{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-u_{\alpha })-B_x^2}{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-S_M)-B_x^2}.$$

$$w_{\alpha }^{\ast }=w_{\alpha }-B_{z\alpha }{B}_x\frac{S_M-u_{\alpha }}{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-S_M)-B_x^2},$$

$$B_{z\alpha }^{\ast }=B_{z\alpha }\frac{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-u_{\alpha })-B_x^2}{{\rho }_{\alpha }(S_{\alpha }^{\prime }-u_{\alpha })(S_{\alpha }-S_M)-B_x^2}.$$

$$e_{\alpha }^{\ast }=\frac{(S_{\alpha }-u_{\alpha })e_{\alpha }-p_{T\alpha }{u}_{\alpha }+p_{T\alpha }^{\ast }{S}_M+B_x({\mathbf{u}}_{\alpha }\cdot {\mathbf{B}}_{\alpha }-{\mathbf{u}}_{\alpha }^{\ast }\cdot {\mathbf{B}}_{\alpha }^{\ast })}{S_{\alpha }-S_M}$$

$$S_{\alpha }^{\ast }=\frac{1}{2}\left[\left(1+{\gamma }_A^2\right)S_M\mp \sqrt{{\left(1-{\gamma }_A^2\right)}^2{S}_M^2+4{\gamma }_A^2{V}_{A,x}^2}\right],$$

$$V_{A,x}^2=\frac{B_x^2}{{\rho }_{\alpha }^{\ast }},$$

$${\gamma }_A=\frac{1}{\sqrt{1+V_A^2/c^2}}=\frac{1}{\sqrt{1+|{\mathbf{B}}^{\ast }{|}^2/({\rho }_{\alpha }^{\ast }{c}^2)}}=\frac{1}{\sqrt{1+{\rho }_A/{\rho }_{\alpha }^{\ast }}},$$

$$e_{\alpha }^{\ast \ast }=e_{\alpha }^{\ast }+\frac{B_x({\mathbf{u}}_{\alpha }^{\ast }\cdot {\mathbf{B}}_{\alpha }^{\ast }-{\mathbf{u}}_{\alpha }^{\ast \ast }\cdot {\mathbf{B}}_{\alpha }^{\ast \ast })}{S_{\alpha }^{\ast }-S_M}$$

$$v_L^{\ast \ast }=v_R^{\ast \ast }=v^{\ast \ast },$$

Full text

magnetic fields — MHD — methods: numerical — plasmas — shock waves

††journal: ApJ††software: SFUMATO (Matsumoto, 2007), Athena++ (Stone et al., 2019)

1 Introduction

Magnetohydrodynamics (MHD) simulations are widely used in astrophysics and space sciences. When MHD equations are solved using explicit methods, a high Alfven speed often arises in a region where the plasma beta is low, and it requires a small timestep, which makes long-term calculations very difficult. The Alfvén speed becomes high when the magnetic field is strong or the gas density is low. An extremely low density requires an extremely small timestep, which greatly increases the number of timesteps and thus the computational burden. Such difficulties often arise when gravity is taken into account in MHD simulations (e.g., Matsumoto & Tomisaka, 2004). In MHD simulations of the magnetospheres of strongly magnetized planets, the fast Alfvén wave due to a planetary dipole field causes the same difficulties (e.g., Tóth et al., 2012).

In order to avoid a very small timestep due to an extremely large Alfvén speed, two approaches are commonly adopted. The first approach is to artificially reduce the Lorentz force (Rempel et al., 2009). This method directly weakens the magnetic effects. The second approach is to impose a variable inertia so that the time rate of change of the flow velocity is reduced as the magnetic field strengthens, leading to slow Alfvén and magnetosonic waves. The latter approach is called the Boris correction (Boris, 1970). With this correction, the speeds of the Alfvén and magnetosonic waves are bounded by the speed of light, $c$, which can be set to an artificially low value. Moreover, using the simplified version of the Boris correction,the steady-state solutions are independent of $c$. Both approaches are included in the semi-relativistic MHD equations (Gombosi et al., 2002).

Methods that utilize a reduced Alfvén speed have been used not only in steady-state simulations but also in time-dependent simulations of star formation (e.g., Allen et al., 2003), accretion disks (e.g., Miller & Stone, 2000; Parkin, 2014), and solar physics (e.g., Rempel, 2017), and in space sciences (e.g., Lyon et al., 2004; Tóth et al., 2012). Shock-capturing schemes usually employ the total variation diminishing (TVD) approach, but conventional Riemann solvers, such as the Lax-Friedrichs scheme and the Harten-Lax-van Leer (HLL) scheme (Harten et a., 1983), have also been adopted. A high-resolution, semi-relativistic Riemann solver that resolves many shocks and discontinuities is desirable. The HLLD Riemann solver (Miyoshi & Kusano, 2005) is one of the most widely used high-resolution schemes, being used in simulation codes such as SFUMATO (Matsumoto, 2007; Matsumoto et al., 2017), Athena++ (Stone et al., 2019; Takasao et al., 2018), and PLUTO (Mignone et al., 2012). It resolves the contact discontinuity, Alfvén waves, and fast waves. Although Parkin (2014) used the HLLD solver with the Boris correction, a simplified approach, in which the numerical flux is not compatible with the Boris correction, was adopted.

In this paper, we propose a scheme that incorporates the Boris correction into the HLLD Riemann solver. The rest of this paper is organized as follows. In Section 2, the HLLD solver with the Boris correction is derived. The results of numerical tests are presented in Section 3. Finally, a summary of the main results and a discussion are presented in Section 4.

2 Incorporation of Boris correction into HLLD solver

2.1 Governing equations

The governing equations of the semi-relativistic equation with the Boris simplification are given as follows (Gombosi et al., 2002),

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{U}$, $\mathbb{F}$, $\rho$, $\mathbf{u}$, $p$, $p_{T}$, $e$, $\mathbf{B}$, $V_{A}$, and $c$ are the state vector, flux, density, velocity, pressure, total pressure, total energy, magnetic field, Alfvén speed, and speed of light, respectively. The superscript $T$ denotes the transpose of a vector. Hereafter, we refer to this formulation as the Boris correction. The difference between this formulation and the original MHD equations is the factor $\left(1+V_{A}^{2}/c^{2}\right)$ in the three components of momentum (see Equation (2)). This factor indicates that inertia becomes large where the Alfvén speed is high compared to the reduced speed of light $c$. For later convenience, we define $\rho_{A}$ as,

[TABLE]

and write the three components of momentum as $\left(\rho+\rho_{A}\right)\mathbf{u}$, indicating that $\rho_{A}$ is the extra inertia caused by the magnetic fields.

For one-dimensional problems, the governing equations reduce to,

[TABLE]

[TABLE]

[TABLE]

2.2 Assumptions for Riemann solver

As is the case for the original HLLD solver, the HLLD solver with the Boris correction (hereafter Boris-HLLD solver) is assumed to have four intermediate states, $\mathbf{U}_{L}^{*}$, $\mathbf{U}_{L}^{**}$, $\mathbf{U}_{R}^{**}$, and $\mathbf{U}_{R}^{*}$, for the given left and right states, $\mathbf{U}_{L}$ and $\mathbf{U}_{R}$ (Figure 1). The states are separated by the five wave speeds $S_{L}$, $S_{L}^{*}$, $S_{M}$, $S_{R}^{*}$, and $S_{R}$. The wave speeds $S_{L}$ and $S_{R}$ are the highest wave speeds, usually the speeds of the fast waves, $S_{L}^{*}$ and $S_{R}^{*}$ are the speeds of the Alfvén waves, and $S_{M}$ is the speed of an entropy wave. Hereafter, we denote the physical variables in the intermediate states with the superscripts $*/**$ and subscripts $L/R$; for example, the density in the state $\mathbf{U}^{*}_{L}$ is $\rho^{*}_{L}$. Similarly, the fluxes, $\mathbf{F}_{L}$, $\mathbf{F}_{L}^{*}$, $\mathbf{F}_{L}^{**}$, $\mathbf{F}_{R}^{**}$, $\mathbf{F}_{R}^{*}$, and $\mathbf{F}_{R}$, are defined. For example, the flux in the state $\mathbf{U}^{*}_{L}$ is $\mathbf{F}_{L}^{*}=\mathbf{F}(\rho^{*}_{L},\mathbf{u}^{*}_{L},\mathbf{B}^{*}_{L},p_{TL}^{*},e^{*}_{L})$, where the function $\mathbf{F}(\rho,\mathbf{u},\mathbf{B},p_{T},e)$ is given by Equation (12). Throughout all the states, $B_{x}$ is constant because of the divergence-free feature of the magnetic field.

We construct the following jump condition for each wave,

[TABLE]

where $\alpha=L$ or $R$. Equations (13), (14), and (15) are the jump conditions across $S_{\alpha}$, $S^{*}_{\alpha}$, and $S_{M}$, respectively. The intermediate states are derived based on these jump conditions.

First, we assume that $u$ and $p_{T}$ are constant in the Riemann fan, yielding,

[TABLE]

and

[TABLE]

This assumption is the same as that in the original HLLD solver. Equation (16) and the mass conservation component of the jump condition across $S^{*}_{\alpha}$ (Equation (14)) give,

[TABLE]

for $S^{*}_{\alpha}\neq S_{M}$. This relation is the same as that for the original HLLD solver.

The $x$ momentum component of the three jump conditions across $S_{M}$ (Equation (15)) and $S^{*}_{\alpha}$ (Equation (14)) gives,

[TABLE]

indicating that $\rho_{A}$ is constant throughout the intermediate states.

2.3 Evaluation of $\mathbf{U}^{*}_{\alpha}$

Considering the jump conditions across $S_{\alpha}$ (Equation (13)), the mass conservation component yields,

[TABLE]

The $x$ momentum and mass conservation components give,

[TABLE]

[TABLE]

where

[TABLE]

The $y$ components of the momentum and induction equations of the jump conditions give,

[TABLE]

[TABLE]

Similarly, the $z$ components give,

[TABLE]

[TABLE]

From the energy conservation component,

[TABLE]

is derived.

2.4 Evaluation of $\mathbf{U}^{**}_{\alpha}$

Considering the jump conditions across $S^{*}_{\alpha}$ (Equation (14)), the $y$ and $z$ components of the momentum and induction equations are satisfied for arbitrary $v^{*}_{\alpha}$, $w^{*}_{\alpha}$, $v^{**}_{\alpha}$, $w^{**}_{\alpha}$, $B_{y\alpha}^{*}$, $B_{z\alpha}^{*}$, $B_{y\alpha}^{**}$, and $B_{z\alpha}^{**}$ when the following equations hold,

[TABLE]

[TABLE]

[TABLE]

where the sign $\mp$ corresponds to $\alpha=L$ and $R$. The subscript $\alpha$ is omitted in $V_{A,x}$ and $\gamma_{A}$ for simplicity, although they depend on the right and left states. For $V_{A}\ll c$, $S^{*}_{\alpha}\rightarrow S_{M}\mp V_{A,x}$ holds, which coincides with $S^{*}_{\alpha}$ in the original HLLD solver. For the special case of $S_{M}=0$ and $V_{A}\gg c$, $S^{*}_{\alpha}\rightarrow\mp c$, indicating that the wave speed is bounded by $c$. For the more general case of $V_{A}\gg c$, $S^{*}_{\alpha}\rightarrow S_{M}$ or 0. Equation (29) was also derived by Gombosi et al. (2002). From the energy component of the jump conditions,

[TABLE]

is obtained. The remaining components (mass conservation and $x$ component of the momentum) hold with those variables.

Considering the jump condition across $S_{M}$ (Equation (15)), the $y$ and $z$ components of the momentum and induction equations yield,

[TABLE]

where the same relationships hold in the original HLLD solver. These variables are also satisfied with the energy component of the jump condition.

The jump condition through the Riemann fan,

[TABLE]

and the jump condition across $S^{*}_{\alpha}$ (Equation (14)) are rewritten as,

[TABLE]

The $y$ and $z$ components of the momentum and induction equations of Equation (38) yield,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For $B_{x}=0$, the Alfvén wave propagating in the $x$ direction does not exist, and $\mathbf{U}^{*}_{L}$ and $\mathbf{U}^{*}_{R}$ are adopted as the intermediate states instead of $\mathbf{U}^{**}_{L}$ and $\mathbf{U}^{**}_{R}$, respectively, as in the original HLLD solver.

2.5 Approximation of $\rho_{A}$

Once $\rho_{A}$ is obtained, all the components of the intermediate states can be obtained. Formally, $\rho_{A}$ should be derived from Equations (19), (23), (25), and (27). However, the simultaneous equations are too complicated to obtain a simple formula for numerical computation. Here, we assume that $\rho_{A}$ is approximated by the so-called HLL average,

[TABLE]

The $B_{y}$ and $B_{z}$ components are,

[TABLE]

By using $\overline{B}_{y}$ and $\overline{B}_{z}$, we adopt

[TABLE]

2.6 $S_{L}$ and $S_{R}$

The wave speeds $S_{L}$ and $S_{R}$ are the highest wave speeds in the two directions, and specify the expansion of the Riemann fan. The following speed settings are usable in practice,

[TABLE]

where $c_{\mathrm{fast},L}^{B}$ and $c_{\mathrm{fast},R}^{B}$ are the speeds of the fast wave in the left and right states, respectively, for the Boris correction. For the governing equations (Equations (10)–(12)), Gombosi et al. (2002) solved the speed of the fast wave for the case of $\mathbf{u}=0$,

[TABLE]

For $V_{A}\ll c$, the wave speed coincides with that of the classical fast wave. For $V_{A}\rightarrow\infty$, $c_{\mathrm{fast}}^{B}\rightarrow c$, and it is bounded by $c$. We adopt Equation (48) for $c_{\mathrm{fast},L}^{B}$ and $c_{\mathrm{fast},R}^{B}$ in the numerical computations, substituting the primitive variables of the corresponding right and left states.

The appropriate order of the wave speeds, $S_{L}\leq S^{*}_{L}$ and $S^{*}_{R}\leq S_{R}$, is not guaranteed when Equations (46), (47), and (48) are adopted as $S_{L}$ and $S_{R}$. We therefore arrange the order of the wave speeds using the following correction,

[TABLE]

In addition, the denominators of Equations (24)–(27) can become zero. For very small denominators, e.g., ${\rho_{\alpha}(S^{\prime}_{\alpha}-u_{\alpha})(S_{\alpha}-S_{M})-B_{x}^{2}}<\epsilon p^{*}_{T\alpha}$, the wave speed $S_{L}$ is reduced and $S_{R}$ is increased by a small value in order to avoid a zero denominator, where $\epsilon\,(=10^{-6})$ denotes a small factor. This corresponds to an extension of the Riemann fan. The small value for reducing/increasing $S_{L/R}$ is set to $(S_{R}-S_{L})\epsilon$, where $S_{R}-S_{L}$ is the total width of the Riemann fan.

When $S_{L}$ and $S_{R}$ are corrected according to the rearrangement of the order of wave speeds and the correction for zero denominators, the dependent variables $\rho^{*}_{\alpha}$, $\rho_{A}$, $S_{M}$, $S^{*}_{\alpha}$, and $S^{\prime}_{\alpha}$ should be recalculated in principle. However, we only recalculate $S^{\prime}_{\alpha}$ because the recalculation of all the dependent variables does not guarantee an appropriate ordering of wave speeds and non-zero denominators. Other treatments may be possible depending on the implementation.

The correction factor $\gamma_{A}$ in Equation (48) produces an over-correction of the fast wave speed; it is reduced to $c_{\mathrm{fast}}^{B}$ even for the case of $c_{\mathrm{fast}}\lesssim c$, where $c_{\mathrm{fast}}$ $(=c_{\mathrm{fast}}^{B}\gamma_{A}^{-1})$ denotes the speed of the classical fast wave. The speed of the classical fast wave is therefore adopted if $c_{\mathrm{fast}}\leq c$; otherwise, Equation (48) is adopted. This switching of the fast wave is effective for the case where both the gas velocity and fast wave speed are low. The effects of switching are discussed in Section 3.4.

2.7 Numerical flux

Given the left and right states $\mathbf{U}_{L}$ and $\mathbf{U}_{R}$, the intermediate states $\mathbf{U}_{L}^{*}$, $\mathbf{U}_{L}^{**}$, $\mathbf{U}_{R}^{**}$, and $\mathbf{U}_{R}^{*}$ are obtained using the primitive variables derived above. When the fluxes in the right and left states $\mathbf{F}_{L}$ and $\mathbf{F}_{R}$ are given, the numerical fluxes are sequentially obtained with the jump conditions for all the intermediate states,

[TABLE]

Alternatively, the following method may be easier than Equation (51),

[TABLE]

where the function $\mathbf{F}(\rho,\mathbf{u},\mathbf{B},p_{T},e)$ is given by Equation (12). As in the original HLLD solver, the numerical flux is switched according to the wave speeds,

[TABLE]

For example, for the case shown in Figure 1, $\mathbf{F}_{L}^{**}$ is adopted as the numerical flux.

The intermediate states and waves derived in Sections 2.2 to 2.6 coincide with those in the original HLLD solver in the limit of $V_{A}\ll c$. The solutions are therefore expected to coincide with those obtained by the original HLLD solver in this limit.

3 Numerical tests

3.1 Implementation

For the test calculations, we use SFUMATO (Matsumoto, 2007), in which the Boris-HLLD solver is implemented. Adaptive mesh refinement is switched off (i.e., uniform grids are utilized). The scheme has second-order accuracy in time and space with the predictor-corrector method and the Monotonic Upwind Scheme for Conservation Laws (MUSCL), respectively. The minmod limiter is adopted as a slope limiter in the MUSCL unless explicitly mentioned. Hyperbolic divergence cleaning (Dedner et al., 2002) is adopted for the $\nabla\cdot\mathbf{B}$ treatment. The specific heat is set to $\gamma=5/3$ for all problems.

Incorporating the Boris-HLLD scheme is easy; the numerical flux of the original HLLD solver is replaced by that of the Boris-HLLD solver, given in Equation (53). The state vector $\mathbf{U}$ is modified according to Equation (2). The timestep is determined based on the Courant-Friedrichs-Lewy (CFL) condition, and the reduced fast wave speed is adopted for the three-dimensional case,

[TABLE]

where $C_{\mathrm{CFL}}$ denotes the CFL number (typically 0.7), the subscripts $(i,j,k)$ specify a cell, and $(\Delta x,\Delta y,\Delta z)$ denote the cell widths. The wave speed $c_{\mathrm{fast}}^{B,\mathrm{max}}$ is the maximum speed of the reduced fast wave. The reduced fast wave speed is given by Equation (48); its maximum value is given by,

[TABLE]

This CFL condition is the same as that for the classical MHD solver except for the factor $\gamma_{A}=\left(1+V_{A}^{2}/c^{2}\right)^{-1/2}$ in Equation (55). For the one-dimensional case, Equation (54) reduces to,

[TABLE]

As shown later, the proposed scheme becomes unstable in a certain situation. In order to confirm that the instability comes from the formulation of the Boris simplification, not from the discretization of the Boris-HLLD solver introduced in Section 2, we incorporate the Boris correction also into the HLL Riemann solver. The numerical flux of the HLL solver with the Boris correction (hereafter Boris-HLL solver) is given by

[TABLE]

where

[TABLE]

for $\alpha=L$ or $R$. The signal speed of the Alfvén wave is denoted by $\lambda_{A,\alpha}^{\pm}$. Equations (58) and (59) are the same as Equations (46) and (47), but they also arrange the order of the wave speeds in a way similar to Equation (50). This arrangement is necessary for the Alfvén and sound waves to be stable in stability tests in Section 3.4. The switching of the fast wave described in Section 2.6 is implemented because it is effective also in the Boris-HLL solver. The Boris-HLL solver is implemented in SFUMATO. The scheme is the same as the Boris-HLLD scheme except for the numerical flux, providing a fair comparison between the Boris-HLLD and Boris-HLL solvers.

In order to confirm that the instability does not arise from the implementation of the code, we use Athena++ (Stone et al., 2019) for comparison. The Athena++ code also adopts a scheme with second-order accuracy in time and space with the predictor-corrector method and MUSCL, respectively. The van Leer limiter is adopted as a slope limiter in the MUSCL. For the $\nabla\cdot\mathbf{B}$ treatment, the constraint transport method is adopted, conserving the initial $\nabla\cdot\mathbf{B}$ within machine accuracy (Stone & Gardiner, 2009).

3.2 Shock tube problem

The standard MHD shock tube problem proposed by Brio & Wu (1988) is solved. In the computational domain $x\in[-5,5]$ with 256 mesh points, the initial state is set to $\rho=1$, $B_{y}=1$, and $p=1$ for $x<0$, and $\rho=0.125$, $B_{y}=-1$, and $p=0.1$ for $x\geq 0$. Throughout the computational domain, $B_{x}=0.75$, $B_{z}=0$, and $\mathbf{u}=0$ are constant. The Alfvén speed is $V_{A}=1.25$ for $x<0$ and 3.5 for $x>0$. The wave speeds are reduced considerably for $x>0$ when the speed of light is set to $c=3$ or 2. For comparison, the original HLLD solver and the Boris-HLL solver are used in addition to the proposed Boris-HLLD scheme.

Figure 2 shows the distributions of physical variables in the shock tube problem. As shown, the shock waves and discontinuities are sharply resolved without any overshoot by the Boris-HLLD solver. Because of the Boris correction, the wave speeds are reduced as $c$ decreases. The fast rarefaction wave slows down, and the location of the wave front is in good agreement with that expected from Equation (48), as indicated by the arrows in the figure. Both velocity components ($u$ and $v$) decrease because of the extra inertia $\rho_{A}$, as can be seen for the fast rarefaction wave and the slow shock at $x\simeq 1.5$. The density, pressure, and magnetic field also change. The decrease in the velocity jump leads to a decrease in the density jump at the slow shock. In contrast, the pressure jump increases, resulting in a considerable increase in the temperature of the post-shock gas, as indicated by the $p/\rho$ distribution. This test problem demonstrates the influence of the Boris correction on the solution.

The solution is not considerably affected by the Boris correction for $x<0$. This is because the Alfvén speed is considerably lower than the speed of light in this region. The solution converges to that of the original HLLD solver as $c$ increases.

Setting $c$ to a low value decreases the number of timesteps. For the shock tube problem examined here, the original HLLD scheme requires 142 timesteps, while the Boris-HLLD scheme requires 105 and 94 timesteps for $c=3$ and 2, respectively.

Figure 2 also compares the solutions between the Boris-HLLD and Boris-HLL solvers in the case of $c=2$. The solutions obtained with the two solvers are in good agreement. However, the Boris-HLLD solver exhibits sharper profiles of the contact discontinuity in the distributions of $\rho$ and $p/\rho$ than the Boris-HLL solver does because the Boris-HLLD solver resolves a contact discontinuity. The difference in the sharpness between the two solvers is approximately the same as the difference between the original HLLD and HLL solvers. At the slow shock, the Boris-HLL solver produces small overshoots in the distributions of $u$ and $p$. The size of the overshoots depends on a slope limiter adopted in the MUSCL. When we use the van Leer limiter, which is a steeper limiter than the minmod limiter, both the Boris-HLLD and Boris-HLL scheme show overshoots at the slow shock, but the Boris-HLLD scheme exhibits a considerably smaller overshoot than the Boris-HLL scheme. The solution with the Boris-HLLD solver is more accurate than that with the Boris-HLL solver because it exhibits a sharp contact discontinuity and a slow shock with a small or no overshoot.

3.3 Linear Alfvén waves

We consider linear Alfvén waves propagating parallel to a uniform magnetic field $B_{0}$ with a uniform density $\rho_{0}$. The linear analysis of the one-dimensional governing equations (Equations (10)–(12)) leads to the following eigen mode of the Alfvén wave,

[TABLE]

where

[TABLE]

The initial condition was constructed according to Equation (63) with $t=0$, $\rho_{0}=1$, $B_{0}=1$, and $\delta_{\mathrm{pert}}=10^{-5}$. The classical Alfvén speed is therefore $V_{A}=1$. The initial condition also has $\rho=\rho_{0}$, $p=1$, $u=w=0$, $B_{x}=B_{0}$, and $B_{z}=0$. The wavelength is set to $L=1$ ($k=2\pi$). The computational domain is $x\in[-L/2,L/2]$ with a uniform grid with 128 mesh points. The periodic boundary condition is imposed at $x=-L/2$ and $L/2$. The calculation is terminated at $t_{\mathrm{last}}=1$, which is the wave crossing time for $c=\infty$.

Figure 3 (left panel) shows the profiles of $v$ for different $c$ values at $t=1$. The wave with $c=1000$ propagates a distance of one wavelength. The travel distance becomes shorter for a lower $c$. Figure 3 (right panel) shows the wave velocity measured in the calculations. In order to measure the velocity for each wave, the travel distance of the wave was evaluated with a phase offset of the first mode of the Fourier transform on the wave profile at $t=t_{\mathrm{last}}$. The measured wave velocities are in agreement with the theoretical values, which are shown by the solid line. The wave velocities are limited by $c$. Those with high $c$ asymptotically approach the classical Alfvén speed ($V_{A}=1$ in this case).

3.4 Stability of linear waves

The stability of the Alfvén wave, sound wave, and fast magnetosonic wave is investigated numerically. The propagation of a linear wave is calculated, and the increase in the amplitude of the wave is measured for each parameter. We use the amplification factor of the wave as an indicator of instability.

Alfvén waves propagating along a magnetic field are considered, as in Section 3.3, but the bulk motion of the gas is also taken into account. An Alfvén wave with a small amplitude in the moving gas has the following eigen mode for the governing equations (Equations (10)–(12)),

[TABLE]

where $\lambda_{A}$ is the signal speed of the Alfvén wave, given by,

[TABLE]

The initial condition was constructed according to Equation (67). We set $\rho=\rho_{0}=1$, $p=p_{0}$, $w=0$, $B_{x}=B_{0}$, $B_{z}=0$, and $\delta_{\mathrm{pert}}=10^{-4}$ at $t=0$. The wavelength is set to $L=1$ ($k=2\pi$). The computational domain is $x\in[-L/2,L/2]$ with a uniform grid with 128 mesh points. The periodic boundary condition is imposed at $x=-L/2$ and $L/2$. The calculation is terminated at $t_{\mathrm{last}}=1$.

We change $u_{0}$, $B_{0}$, and $p_{0}$ in the range of $V_{A}/c\in[0.1,3.0]$, $u_{0}/c\in[-2,2]$, and $a/c\in[0.1,3.0]$ in two-dimensional spaces of $V_{A}/c-u_{0}/c$ and $a/c-u_{0}/c$ (see Figure 4 for the parameter space). The speed of light is set to $c=10$. The numerical tests show that all the Alfvén waves are stable irrespective of $u_{0}$, $V_{A}$, and $a$ for both the Boris-HLLD and Boris-HLL schemes. This is because the wave speed of the eigen mode (Equation (68)) is always real.

In the semi-relativistic MHD formulations proposed by Gombosi et al. (2002), the Alfvén wave parallel to the magnetic field is unstable when $u_{0}>V_{A}$. The difference between our result and their result arises from the adopted form of the equations of motion. Gombosi et al. (2002) adopted the semi-relativistic equations of motion, which consider the off-diagonal terms in momentum, whereas we used the equations based on the Boris correction, neglecting the off-diagonal terms. Although the equations with off-diagonal terms are physically preferred, those with the Boris correction result in a more stable scheme.

Next, a sound wave propagating along a magnetic field is considered. For $a>V_{A}$, the sound wave corresponds to the fast wave, and the Alfvén wave degenerates to the slow wave. For $a<V_{A}$, the Alfvén wave degenerates to the fast wave, and the sound wave corresponds to the slow wave. The eigen mode of the sound wave for the MHD equations with the Boris correction (Equations (10)–(12)) is expressed as,

[TABLE]

where $\lambda_{s}$ is the signal speed of the sound wave, given by,

[TABLE]

The signal speed $\lambda_{s}$ is real because $\gamma>1$ and $\gamma_{A}<1$, indicating that the sound wave is always stable. Although this wave is a pure sound wave parallel to the magnetic field, the signal speed depends on the magnetic field through the factor $\gamma_{A}$. We examine the propagation of sound waves. The settings are the same as those for the Alfvén wave test except for the initial perturbation. The numerical tests show that the sound wave is stable irrespective of $a$, $V_{A}$, and $u_{0}$ for both the Boris-HLLD and Boris-HLL schemes.

Finally, we consider a magnetosonic wave propagating perpendicular to the magnetic field. This wave corresponds to the fast wave. Note that no Alfvén or sound waves exist perpendicular to the magnetic field. Since the eigen mode of the fast wave is complex for governing equations (Equations (10)–(12)), we simply extend the fast wave in the classical MHD, and the initial condition is given by,

[TABLE]

where $c_{\mathrm{fast}\perp}$ is the classical speed of the fast wave propagating perpendicular to the magnetic field, defined as,

[TABLE]

According to Equation (48), $\gamma_{A}c_{\mathrm{fast}\perp}$ is the asymptotic speed of the fast wave. We set $B_{x}=B_{z}=0$ and $v=w=0$ in the initial condition. The other parameters are the same as those in the previous tests of the Alfvén and sound waves.

Figure 4 (top panels) shows the distribution of the amplification factor of the magnetosonic wave at $t=1$ for the Boris-HLLD scheme. The scheme is stable (blue regions) for $|u_{0}|\lesssim(0.6-1)c$ when the magnetic field is relatively strong ($V_{A}\gtrsim c$), as shown in left and right panels. The boundary between the stable and unstable regions weakly depends on $V_{A}/c$ and $a/c$. In the relatively weak magnetic field case ($V_{A}\lesssim c$), the unstable region extends down to $|u_{0}|\sim 0.5c$ (left panel). In the stable (blue) regions, the amplitude is distributed around unity because the initial conditions are not pure eigen modes. The amplitude of the wave oscillates as time proceeds even for a stable wave. Note that the left edge of the left panel corresponds to the non-magnetized case ($V_{A}/c=0$), and the scheme is stable irrespective of $u_{0}$ there.

Figure 4 (bottom panels) shows the distribution of the amplification factor for the Boris-HLL scheme for comparison. Both the Boris-HLLD and Boris-HLL schemes are unstable for large $|u_{0}|/c$, indicating that the instability comes from the governing equations of the Boris correction rather than the discretization of the schemes. The Boris-HLLD scheme shows larger stable regions than the Boris-HLL scheme in the diagrams, indicating that the Boris-HLLD scheme is more stable than the Boris-HLL scheme for the magnetosonic wave.

Figure 4 shows the results in the case of the minmod limiter in the MUSCL. We confirm that the stable regions do not change even when we use the other limiter, e.g., the van Leer limiter, and when we adopt a scheme with a spatially first order accuracy without the MUSCL. These results indicate that the difference in the stable regions between the Boris-HLLD and Boris-HLL solvers is attributed to the difference between the solvers.

Figure 5 compares the stability of the magnetosonic wave between the schemes with and without the switching of the fast wave speed introduced in Section 2.6. Without the switching, instability arises for $c_{\mathrm{fast}}\lesssim u_{0}\lesssim\gamma_{A}c_{\mathrm{fast}}$ (the red narrow region in the right panel). The regions above and below the lines represent super- and sub-magnetosonic flows, respectively. In the unstable region, the numerical flux $\mathbf{F}_{L}$ is adopted, which corresponds to a fully upwind difference. This numerical test indicates that this unstable region is stabilized by adopting the numerical flux for sub-magnetosonic flow $\mathbf{F}_{L}^{*}$. For the case of the Boris-HLL solver, we observe the same results as the case of the Boris-HLLD solver.

3.5 Orszag-Tang vortex problem

The Orszag-Tang vortex problem (Orszag & Tang, 1979) is widely used as a two-dimensional test problem for MHD schemes. The computational domain is $x,y\in[0,2\pi]$ with $256^{2}$ mesh points. The periodic boundary conditions are imposed on the edges of the computational box, $x=0,2\pi$, and $y=0,2\pi$. The initial condition has distributions of $\rho=\gamma^{2}$, $p=\gamma$, $u=-\sin y$, $v=\sin x$, $v=0$, $B_{x}=-B_{0}\sin y$, $B_{y}=B_{0}\sin 2x$, $B_{z}=0$, and $B_{0}=1$. In this problem, we set a pressure floor ($p_{\mathrm{floor}}=10^{-2}$) in order to detect the onset of numerical instability.

Figure 6 shows the gas velocity and Alfvén speed distributions obtained using the original HLLD scheme and the Boris-HLLD scheme with different speeds of light. The solution obtained with $c=10$ is quite similar to that of the HLLD scheme. The time evolutions for the two schemes coincide, as shown in Figure 7 (compare orange and red lines). For $c=2$ (blue lines), the Alfvén speed exceeds $c$ at the later stages. Even when the gas velocity is lower than $c$, the maximum velocity becomes larger than that in the HLLD solution. The distributions of gas velocity and Alfvén speed are slightly affected by $c$ (bottom left panel in Figure 6). For $c=1$, the gas velocity is higher than $c$ in a considerable area, in which numerical instability is prominent (bottom right panel in Figure 6). Checkerboard-like instability appears in regions where the gas velocity exceeds $c$. The minimum pressure reaches the floor value in the early stage, as depicted by the arrow in Figure 7. The maximum gas velocity grows exponentially at this time.

Figure 8 shows the same models as those in Figure 6 but calculated using Athena++ for comparison. For $c=10$ and 2, the solutions obtained with the two codes are in agreement, indicating that the difference between the solutions with different $c$ values results from the Boris correction; it does not depend on the implementation of the code. As mentioned, SFUMATO adopts the hyperbolic divergence cleaning method and Athena++ adopts the constraint transport method for the treatment of $\nabla\cdot\mathbf{B}$. For $c=1$, checkerboard-like instability appears when the gas velocity exceeds $c$. Although numerical instability sometimes depends on the implementation of numerical schemes, checkerboard-like instability appeared in the solutions obtained with SFUAMTO as well as those obtained with Athena++. This indicates that the numerical instability here arises from the basic equations adopted.

4 Summary and discussion

We proposed a high-resolution scheme for the ideal MHD equations that incorporates the simplified version of the Boris correction (Gombosi et al., 2002) into the HLLD Riemann solver. The Boris correction introduces an extra inertia term in the momentum equation. The extra inertia is reflected by the magnetic field strength, and reduces the wave speeds. The wave speeds are bounded by the speed of light, which can be set to an artificially low value in order to avoid an extremely small timestep. As done by the original HLLD solver, the proposed scheme resolves four intermediate states separated by five waves: two fast waves, two Alfvén waves, and an entropy wave. In the limit of $V_{A}\ll c$, all the intermediate states and the numerical fluxes converge to those of the original HLLD.

Incorporating the Boris-HLLD scheme into existing code is simple. The numerical flux is replaced by that of the Boris-HLLD solver, the state vector is modified for the Boris correction, and the CFL condition is modified so that the wave speed is multiplied by the factor $\gamma_{A}$.

We performed a stability analysis and showed the parameter space in which the scheme is stable. The scheme is stable when $|u|\lesssim 0.5c$ for a low Alfvén speed ($V_{A}\lesssim c$). For a high Alfvén speed ($V_{A}\gtrsim c$), the stable region becomes large for $|u|\lesssim(0.6-1)c$. The Boris-HLLD scheme shows larger stable regions than the Boris-HLL scheme. The scheme can be unstable even when $V_{A}<c$, and the semi-relativistic treatment is not necessary there. In this case, one can switch the scheme to the original HLLD scheme or adopt a sufficiently high value of $c$ to avoid instability. Practically, setting the speed of light to several times higher than the maximum gas speed is an acceptable compromise (Rempel, 2017).

We showed the effects of the Boris correction on the solutions of non-steady-state problems (shock tube and the Orszag-Tang vortex problems). The Boris-HLLD scheme captures a contact discontinuity more sharply than the Boris-HLL scheme does. Although the semi-relativistic scheme including the Boris correction is powerful for stringent timestep problems, one has to check the impact of the modification, especially the dynamics in the region where $V_{A}\gtrsim c$. In other regions, the solution will be only weakly affected by the Boris correction. This scheme is therefore useful for avoiding an extremely high Alfvén speed in a relatively small volume in the computational domain. A conventional treatment for such a high Alfvén speed is to introduce a density floor (e.g., Bai & Stone, 2013). However, with a density floor, mass and energy are unphysically injected into the computational domain. When self-gravity is taken into account, the influence of a density floor is more serious because it can increase gravity. The Boris-HLLD solver is an alternative method that overcomes these difficulties.

Numerical computations were carried out in part on XC50 (ATERUI II) at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan. S.T. acknowledges support by the Research Fellowship of the Japan Society for the Promotion of Science (JSPS). This research was supported by JSPS KAKENHI Grant Numbers 15K04756, 16J02063, 17K05671, 17H02863, 17K05394, 18H04449, 18H05437, and 18K13579.