A high-order weighted finite difference scheme with a multi-state approximate Riemann solver for divergence-free magnetohydrodynamic simulations

TL;DR

This paper introduces a high-order finite difference scheme for ideal magnetohydrodynamics that ensures shock-capturing, divergence-free magnetic fields, and high accuracy through advanced interpolation, Riemann solvers, and a novel constrained transport method.

Contribution

It presents a new high-order finite difference scheme with a multi-state Riemann solver and a modified constrained transport method for divergence-free magnetic fields in MHD simulations.

Findings

The scheme achieves high-order accuracy and shock-capturing capabilities.

It maintains divergence-free magnetic fields with a new CT variant.

Benchmark tests demonstrate improved solution quality.

Abstract

We design a conservative finite difference scheme for ideal magnetohydrodynamic simulations that attains high-order accuracy, shock-capturing, and divergence-free condition of the magnetic field. The scheme interpolates pointwise physical variables from computational nodes to midpoints through a high-order nonlinear weighted average. The numerical flux is evaluated at the midpoint by a multi-state approximate Riemann solver for correct upwinding, and its spatial derivative is approximated by a high-order linear central difference to update the variables with designed order of accuracy and conservation. The magnetic and electric fields are defined at staggered grid points employed in the Constrained Transport (CT) method by Evans & Hawley (1988). We propose a new CT variant, in which the staggered electric field is evaluated so as to be consistent with the base one-dimensional Riemann solver and the staggered magnetic field is updated to be divergence-free as designed high-order finite difference representation. We demonstrate various benchmark tests to measure the performance of the present scheme. We discuss the effect of the choice of interpolation methods, Riemann solvers, and the treatment for the divergence-free condition on the quality of numerical solutions in detail.

Figures

Tables

Table 1. Table 1: L 1 subscript 𝐿 1 L_{1} errors (orders) of the Alfvén wave propagation problem.

$B_{⟂}$	$N=16$	32	64	128	256	512
First	4.80e-03 (-)	4.44e-03 ( 0.11)	3.27e-03 ( 0.44)	2.06e-03 ( 0.67)	1.17e-03 ( 0.82)	6.27e-04 ( 0.90)
MUSCL	4.76e-03 (-)	3.95e-04 ( 3.59)	1.07e-04 ( 1.88)	2.85e-05 ( 1.91)	7.27e-06 ( 1.97)	1.82e-06 ( 2.00)
W3I4D3	1.57e-03 (-)	5.55e-04 ( 1.50)	1.26e-04 ( 2.14)	2.67e-05 ( 2.23)	5.32e-06 ( 2.33)	1.05e-06 ( 2.34)
W4I4D5	1.84e-04 (-)	4.85e-06 ( 5.24)	1.39e-07 ( 5.12)	5.51e-09 ( 4.66)	2.90e-10 ( 4.25)	1.74e-11 ( 4.06)
W5I4D4	1.79e-04 (-)	5.01e-06 ( 5.16)	1.69e-07 ( 4.89)	8.30e-09 ( 4.35)	4.86e-10 ( 4.09)	3.00e-11 ( 4.02)
$B_z$	$N=16$	32	64	128	256	512
First	2.48e-03 (-)	3.32e-03 (-0.42)	2.48e-03 ( 0.42)	1.54e-03 ( 0.68)	8.68e-04 ( 0.83)	4.61e-04 ( 0.91)
MUSCL	1.78e-03 (-)	4.05e-04 ( 2.14)	1.05e-04 ( 1.94)	2.46e-05 ( 2.10)	5.30e-06 ( 2.21)	1.22e-06 ( 2.12)
W3I4D3	1.20e-03 (-)	4.32e-04 ( 1.48)	9.15e-05 ( 2.24)	1.76e-05 ( 2.38)	3.23e-06 ( 2.45)	5.83e-07 ( 2.47)
W4I4D5	1.59e-04 (-)	3.83e-06 ( 5.37)	9.48e-08 ( 5.34)	2.57e-09 ( 5.21)	7.41e-11 ( 5.12)	2.23e-12 ( 5.05)
W5I4D4	1.55e-04 (-)	3.76e-06 ( 5.36)	9.79e-08 ( 5.26)	3.29e-09 ( 4.90)	1.54e-10 ( 4.42)	8.82e-12 ( 4.12)

Table 2. Table 2: Relative differences (convergence rates) of the Orszag-Tang vortex problem from the reference solution. The solution is smooth at t = 0.2 ​ π 𝑡 0.2 𝜋 t=0.2\pi , but becomes discontinuous at t = π 𝑡 𝜋 t=\pi .

$t=0.2\pi$	$N=100$	200	400	800
MUSCL	1.17e-03 (-)	2.65e-04 ( 2.15)	6.05e-05 ( 2.13)	1.43e-05 ( 2.08)
HYBRID	7.74e-04 (-)	1.52e-04 ( 2.34)	2.87e-05 ( 2.41)	5.28e-06 ( 2.44)
W4I4D5	1.25e-04 (-)	8.14e-06 ( 3.94)	5.35e-07 ( 3.93)	3.71e-08 ( 3.85)
$t=\pi$	$N=100$	200	400	800
MUSCL	3.42e-02 (-)	1.65e-02 ( 1.05)	7.80e-03 ( 1.08)	3.23e-03 ( 1.27)
HYBRID	3.17e-02 (-)	1.49e-02 ( 1.09)	6.68e-03 ( 1.16)	2.52e-03 ( 1.41)
W4I4D5	3.01e-02 (-)	1.38e-02 ( 1.12)	6.19e-03 ( 1.16)	2.25e-03 ( 1.46)

Equations

$$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho\mbox{\boldmath$u$}\right)=0,$$

$$\displaystyle\frac{\partial\left(\rho\mbox{\boldmath$u$}\right)}{\partial t}+\nabla\cdot\left[\rho\mbox{\boldmath$u$}\mbox{\boldmath$u$}+\left(P+\frac{|\mbox{\boldmath$B$}|^{2}}{2}\right)\mbox{\boldmath$\rm I$}-\mbox{\boldmath$B$}\mbox{\boldmath$B$}\right]=0,$$

$$\displaystyle\frac{\partial\mbox{\boldmath$B$}}{\partial t}+\nabla\times\mbox{\boldmath$E$}=0,$$

$$\displaystyle\frac{\partial e}{\partial t}+\nabla\cdot\left[\left(\frac{\rho u^{2}}{2}+\frac{\gamma P}{\gamma-1}\right)\mbox{\boldmath$u$}+\mbox{\boldmath$E$}\times\mbox{\boldmath$B$}\right]=0,$$

$$\displaystyle\mbox{\boldmath$E$}=-\mbox{\boldmath$u$}\times\mbox{\boldmath$B$},$$

$$\displaystyle P=\left(\gamma-1\right)\left(e-\frac{\rho|\mbox{\boldmath$u$}|^{2}}{2}-\frac{|\mbox{\boldmath$B$}|^{2}}{2}\right),$$

$$\displaystyle\nabla\cdot{\mbox{\boldmath$B$}}=0.$$

$$\displaystyle\frac{\partial\mbox{\boldmath$U$}}{\partial t}+\frac{\partial\mbox{\boldmath$F$}}{\partial x}=0,$$

$$\displaystyle\mbox{\boldmath$U$}=\left[\begin{array}[]{c}\rho\\ {\rho u}\\ {\rho v}\\ {\rho w}\\ {B_{y}}\\ {B_{z}}\\ e\\ \end{array}\right],\mbox{\boldmath$F$}=\left[\begin{array}[]{c}\rho u\\ \rho u^{2}+P+\left(B_{y}^{2}+B_{z}^{2}-B_{x}^{2}\right)/2\\ \rho vu-B_{x}B_{y}\\ \rho wu-B_{x}B_{z}\\ B_{y}u-B_{x}v\\ B_{z}u-B_{x}w\\ \left(e+P+B^{2}/2\right)u-B_{x}\left(\mbox{\boldmath$u$}\cdot\mbox{\boldmath$B$}\right)\\ \end{array}\right],$$

$$\displaystyle\frac{d\mbox{\boldmath$U$}^{n}_{i}}{dt}=-\left.\frac{\partial\mbox{\boldmath$F$}}{\partial x}\right|_{i},$$

$$\displaystyle\mbox{\boldmath$U$}^{*}$$

$$\displaystyle\mbox{\boldmath$U$}^{**}$$

$$\displaystyle\mbox{\boldmath$U$}^{n+1}$$

$$\displaystyle\left.\frac{\partial\mbox{\boldmath$F$}}{\partial x}\right|_{i}=\frac{\mbox{\boldmath$\hat{F}$}_{i+1/2}-\mbox{\boldmath$\hat{F}$}_{i-1/2}}{\Delta x}+O(\Delta x^{q}),$$

$$\displaystyle\left.\frac{\partial\mbox{\boldmath$F$}}{\partial x}\right|_{i}=\sum_{k\geq 0}d_{k}\frac{\mbox{\boldmath$\tilde{F}$}_{i+k+1/2}-\mbox{\boldmath$\tilde{F}$}_{i-k-1/2}}{\Delta x}+O(\Delta x^{p})+O(\Delta x^{q}),$$

$$\displaystyle{\cal D}_{4}(\mbox{\boldmath$\tilde{F}$})_{i}$$

$$\displaystyle{\cal D}_{4}(\mbox{\boldmath$\tilde{F}$})_{1}$$

$${\tilde{U}}_{i+1/2}^L$$

$${\mathcal{D}}_4({\tilde{U}}^L{)}_i={\frac{\partial U}{\partial x}|}_i+{\frac{1}{16}\frac{{\partial }^{4}U}{\partial x^4}|}_i\Delta x^3-{\frac{9}{320}\frac{{\partial }^{5}U}{\partial x^5}|}_i\Delta x^4+O(\Delta x^5).$$

$${\tilde{U}}_{i+1/2}^L$$

$${\mathcal{D}}_4({\tilde{U}}^L{)}_i={\frac{\partial U}{\partial x}|}_i-{\frac{3}{640}\frac{{\partial }^{5}U}{\partial x^5}|}_i\Delta x^4-{\frac{3}{256}\frac{{\partial }^{6}U}{\partial x^6}|}_i\Delta x^5+O(\Delta x^6).$$

$${\tilde{U}}_{i+1/2}^L$$

$${\mathcal{D}}_4(\tilde{U}{)}_i={\frac{\partial U}{\partial x}|}_i-{\frac{9}{640}\frac{{\partial }^{6}U}{\partial x^6}|}_i\Delta x^5+O(\Delta x^6).$$

$${\omega }^{\ast }=-j\exp (-j\omega x){\mathcal{D}}_4(U).$$

$${\tilde{U}}_{i+1/2}^L=\sum _k{w}_k{\tilde{U}}_{i+1/2}^{L,k},w_s=\frac{{\alpha }_s}{{\sum }_k{\alpha }_k}$$

$${\tilde{U}}_{i+1/2}^{L,0}$$

$${\tilde{U}}_{i+1/2}^{L,1}$$

$${\alpha }_s=c_s\left(1+\frac{\tau }{I{S}_s+ϵ}\right),$$

$$I{S}_0={\left(U_i-U_{i-1}\right)}^2,I{S}_1={\left(U_{i+1}-U_i\right)}^2,\tau ={\left(U_{i+1}-2U_i+U_{i-1}\right)}^2,$$

$$c_0=\frac{1}{4},c_1=\frac{3}{4},$$

$${\tilde{U}}_{i+1/2}^{L,0}$$

$${\tilde{U}}_{i+1/2}^{L,1}$$

$${\tilde{U}}_{i+1/2}^{L,2}$$

$${\alpha }_s=\frac{c_s}{{\left(I{S}_s+ϵ\right)}^2},$$

$$I{S}_0=\frac{13}{12}{\left(U_{i-2}-2U_{i-1}+U_i\right)}^2+\frac{1}{4}{\left(U_{i-2}-4U_{i-1}+3U_i\right)}^2,$$

$$I{S}_1=\frac{13}{12}{\left(U_{i-1}-2U_i+U_{i+1}\right)}^2+\frac{1}{4}{\left(U_{i-1}-U_{i+1}\right)}^2,$$

$$I{S}_2=\frac{13}{12}{\left(U_i-2U_{i+1}+U_{i+2}\right)}^2+\frac{1}{4}{\left(3U_i-4U_{i+1}+U_{i+2}\right)}^2,$$

Full text

magnetic fields — magnetohydrodynamics (MHD) — methods: numerical

1 Introduction

The magnetohydrodynamic (MHD) modeling has been extensively applied to various macroscopic dynamics in space and astrophysical plasmas such as flares and coronal mass ejections in the solar corona, evolution of the solar wind through the interplanetary space, auroral substorms and magnetic storms in the planetary magnetosphere, formation of the heliosphere through the interaction between the solar wind plasma and the interstellar medium, turbulence and dynamo in magnetized accretion disks, and so on. The MHD simulation is an indispensable tool to study such plasmas since the system of equations is highly nonlinear and a variety of phenomena are coupled with each other. Rapid development in computational technology and science enables us to conduct high-resolution and long-term numerical simulations for global dynamics as well as local physics.

A common strategy to build an MHD simulation code for space and astrophysical plasmas would be based on an upwind-type shock-capturing method, the so-called Godunov’s method. The Godunov-type scheme is based on a finite volume scheme that updates volume-averaged variables defined in computational cells through an integral form of hyperbolic conservation laws in space and time (Godunov, 1959). Godunov’s method evaluates a numerical flux at a cell interface through an exact or approximate solution of the Riemann problem for left- and right-side variables as an initial state. The solution of the Riemann problem retains upwind property for eigenmodes in a hyperbolic system, thus implicitly introduces numerical diffusion enough to stabilize a simulation. The method is suited for the plasma that frequently contains supersonic flows, shocks and discontinuities. Godunov’s method has been extensively succeeded in hydrodynamic simulations, and utilized for modern MHD simulation codes as well (Gardiner & Stone, 2005; Mignone et al., 2007; Gardiner & Stone, 2008; Stone et al., 2008; Lee & Deane, 2009; Mignone & Tzeferacos, 2010; Lee, 2013; Matsumoto et al., 2016; Lee et al., 2017). The choice of Riemann solvers impacts the quality of numerical solutions.

Since the original Godunov’s method assumes variables as piecewise constant in a cell, a numerical error of the solution is reduced only as first order. Improving the accuracy of the method is necessary for practical simulations within a reasonable computational cost. A straightforward way is increasing the order of accuracy of a scheme. For Godunov-type schemes, it is achieved by reconstructing Riemann states with a desired order of accuracy. Piecewise linear interpolation is one of conventional methods to obtain left and right Riemann states in second-order accuracy, and the intermediate state through the Riemann solver. However, high-order linear interpolation methods generally suffer from numerical oscillation and then tend to crash the simulation of nonlinear hyperbolic equations. The robustness as well as the accuracy of the simulation rely on interpolation methods. A nonlinear interpolation method is designed to preserve high order of accuracy in smooth regions but degrade to first order at discontinuous regions so as to avoid the oscillation. The second-order MUSCL scheme (van Leer, 1979) is one of standard nonlinear interpolation methods, and has been widely adopted for practical MHD simulations. To obtain a second-order accurate solution, one may combine a spatially second-order scheme with a multi-step time integration method. Many nonlinear second-order schemes are built to be total variation diminishing (TVD) that inevitably degrade the accuracy at a profile extremum. To preserve the profile extremum, some third- and higher-order schemes have been implemented to MHD simulation codes (Jiang & Wu, 1999; Balsara & Shu, 2000; Matsumoto & Seki, 2008; Li, 2010; Mignone et al., 2010; Matsumoto et al., 2016; Lee et al., 2017). Another approach is to statically or adaptively refine mesh at regions to be resolved finely. The Adaptive Mesh Refinement (AMR) technique has been implemented to the MHD simulation code to tackle problems including extreme dynamic range in space and time, for example, star formation and heliosphere (Powell et al., 1999; Balsara, 2001; Fromang et al., 2006; Matsumoto, 2007; Ziegler, 2008; Stone et al., 2008; Mignone et al., 2012). The finite volume scheme is usually combined with the AMR since the scheme is flexible about the mesh structure. In AMR, single step time integration methods are better suited for the reduction of a computational cost of communications between different mesh levels (Gardiner & Stone, 2005, 2008; Lee & Deane, 2009; Mignone & Tzeferacos, 2010; Lee, 2013).

The Godunov-type scheme can be extended to multi-dimension by a dimension-by-dimension reconstruction. However, a multidimensional finite volume Godunov-type scheme is computationally expensive to ensure an exact order of accuracy better than two. This is because one needs to reconstruct the numerical flux along orthogonal directions for face average: this is not identical to a solution of the Riemann problem for face-averaged variables (Zhang et al., 2011; Büchmuller & Helzel, 2014; Lee et al., 2017). Alternatively, a finite difference scheme is designed for hyperbolic conservation laws to retain high-order of accuracy, conservation, and upwind property (Liu et al., 1994; Jiang & Shu, 1996; Deng & Zhang, 2000). Extension of the high-order finite difference scheme to multi-dimension is straightforward because it updates pointwise variables at computational nodes rather than volume-averaged ones. The calculation of source terms is also straightforward for the finite difference scheme, while the finite volume scheme needs their high-order reconstruction in cells to retain the spatial accuracy better than two. Note that the two schemes are essentially identical at the second-order level.

Multidimensional MHD simulations should take special care of the divergence-free condition for the magnetic field. The Godunov-type scheme discretizes hyperbolic conservation laws into a divergence form. Generally in the baseline Godunov scheme without any special divergence control, the magnetic field is not necessarily divergence-free and a finite divergence error may lead to an unphysical solution, even though the field is updated in the divergence form (Brackbill & Barnes, 1980). Recipes against this problem are classified into two types, divergence-cleaning and divergence-free methods. The divergence-cleaning method keeps a divergence error within a tolerable level by correcting the magnetic field through the Poisson equation (Brackbill & Barnes, 1980; Ryu et al., 1995; Crockett et al., 2005), additional source terms (Powell et al., 1999), or additional equations (Dedner et al., 2002). The divergence-free method utilizes staggered grid spacing so as to discretize the induction equation to be consistent with the Faraday’s law, termed as the Constrained Transport (CT) method (Evans & Hawley, 1988). The CT method defines the in-plane magnetic field at cell interfaces and the out-of-plane electric field at cell corners to preserve the divergence-free condition in a discretized manner within a machine precision (e.g., Tóth, 2000). Since the staggered electric field is not directly obtained from a one-dimensional Riemann solver, a variety of methods have been proposed to calculate it with retaining the upwind property (Hawley & Stone, 1995; Dai & Woodward, 1998; Ryu et al., 1998; Balsara & Spicer, 1999; Londrillo & Del Zanna, 2000). Balsara & Spicer (1999) calculated the electric field through an arithmetic average of the neighboring numerical flux for the in-plane magnetic field obtained from a one-dimensional Riemann solver. Their method has been employed as a baseline for subsequent schemes (Gardiner & Stone, 2005; Lee & Deane, 2009; Lee, 2013). Meanwhile, the solution of a multidimensional Riemann solver was derived for the staggered electric field (Londrillo & Del Zanna, 2004; Balsara, 2010, 2012; Amano, 2015).

Comparative studies between divergence-cleaning and divergence-free methods have been reported by Tóth (2000), Balsara & Kim (2004), and Miyoshi & Kusano (2011).

Numerous studies have been devoted to develop novel numerical schemes so as to attain accurate and robust MHD simulations. A high-order scheme combined with a shock-capturing method improves the resolution of shocks and discontinuities, smooth structures such as vortices and waves, and their interaction driven by a high speed flow at a high Reynolds number (Shu, 2009). A divergence-cleaning method is combined with a shock-capturing method in a straightforward manner. However, it might lead to a spurious solution for a stringent problem such as multiple shock interaction, in which a continuously-produced local divergence error affects globally through a cleaning method (Balsara & Kim, 2004). A recipe for combining a high-order (better than two) scheme, an MHD shock-capturing method, and a divergence-free method must be useful for space and astrophysical plasma simulations (Del Zanna et al., 2007; Balsara et al., 2009; Li, 2010; Felker & Stone, 2018).

The purpose of this paper is to design a numerical scheme for ideal MHD simulations, specifically focusing on the high-order interpolation and the treatment for the divergence-free condition. We build the scheme based on the finite difference scheme with several criteria: (1) the order of accuracy can be improved to an arbitrary level (in principle) through a dimension-by-dimension interpolation, (2) a variety of Riemann solvers are available, and (3) the multidimensional scheme is divergence-free and reduces to the base one-dimensional scheme for one-dimensional problems. Section 2 presents the scheme design. Details of the finite difference scheme for hyperbolic conservation laws are introduced in Section 2.1. The high-order interpolation methods are described in Sections 2.2-2.3, the Riemann solver is briefly introduced in Section 2.4, and the treatment for the divergence-free condition is presented in Section 2.5. Section 3 details numerical simulation results of various benchmark tests so as to measure the performance of the present scheme. We compare a variety of schemes to assess the effect of numerical techniques on the quality of solutions. Finally, Section 4 summarizes the paper.

2 High-order finite difference schemes for ideal MHD

The normalized, fully compressible ideal MHD equations are

[TABLE]

where \rho,\mbox{\boldmathu}=(u,v,w),\mbox{\boldmathB}=(B_{x},B_{y},B_{z}),\mbox{\boldmathE}=(E_{x},E_{y},E_{z}),e and $P$ are the fluid mass density, the velocity, the magnetic field, the electric field, the total energy density, and the gas pressure, respectively. The electric field and the gas pressure are determined from the Ohm’s law and the equation of state,

[TABLE]

where $\gamma$ is the specific heat ratio. The magnetic field should satisfy the divergence-free condition,

[TABLE]

To construct a finite difference scheme in a dimension-by-dimension fashion, we firstly consider the one-dimensional MHD equations in the conservative form as

[TABLE]

where

[TABLE]

and $B_{x}=\rm{constant}$ is given as a constraint in one dimension.

We discretize the physical variables $U$ at a time $t_{n}$ on a computational node $i$ as pointwise representation, \mbox{\boldmathU}^{n}_{i}=\mbox{\boldmathU}(x_{i},t_{n}). Eq. (8) is discretized as

[TABLE]

and is integrated in time with the third-order Strong Stability Preserving Runge-Kutta (SSPRK) method that is popular for the time integration of partial differential equations (Shu & Osher, 1988; Gottlieb & Shu, 1998),

[TABLE]

where $\cal D$ stands for a difference operator and $\Delta t$ is a time step. The quality of numerical solutions relies on the evaluation of the flux derivative as long as a temporal error is relatively small.

Conventional high-order finite difference schemes for hyperbolic conservation laws approximate the flux derivative as a two-point difference (Liu et al., 1994; Jiang & Shu, 1996; Jiang & Wu, 1999; Balsara & Shu, 2000; Mignone et al., 2010; Christlieb et al., 2014),

[TABLE]

where \mbox{\boldmath\hat{F}}_{i+1/2} is the pointwise numerical flux at a midpoint $i+1/2$ and $\Delta x$ is a grid size. The numerical flux $\hat{F}$ is evaluated so that its two-point difference is an approximation of the flux derivative with the desired order of accuracy $q$; it is not identical to the $q$th-order approximation of the physical flux \mbox{\boldmathF}(\mbox{\boldmathU}) itself for $q>2$. As is shown by e.g. Mignone et al. (2010), \mbox{\boldmath\hat{F}}_{i+1/2} is reconstructed from the pointwise physical flux \mbox{\boldmathF}_{i}=\mbox{\boldmathF}(\mbox{\boldmathU}_{i}) by the same reconstruction method used for finite volume schemes.

The example procedure of the conventional finite difference scheme is as follows: (i) the physical flux \mbox{\boldmathF}_{i} is calculated at the node, and is split into positive and negative fluxes for correct upwinding, \mbox{\boldmathF}_{i}=\mbox{\boldmathF}^{+}_{i}+\mbox{\boldmathF}^{-}_{i}, where {\partial\mbox{\boldmathF}^{+}}/{\partial\mbox{\boldmathU}}>0 and {\partial\mbox{\boldmathF}^{-}}/{\partial\mbox{\boldmathU}}<0, (ii) the numerical fluxes are reconstructed to the left side at $x_{i+1/2}$ and the right side at $x_{i-1/2}$, \mbox{\boldmath\hat{F}}^{+}_{i+1/2} and \mbox{\boldmath\hat{F}}^{-}_{i-1/2}, with the positive and negative fluxes \mbox{\boldmathF}^{\pm} from Step (i) in an upwind-biased stencil $(i-r,\dots,i+r)$. For example, the optimal fifth-order $(r=2)$ upstream central scheme gives \mbox{\boldmath\hat{F}}^{\pm}_{i\pm 1/2}=(2\mbox{\boldmathF}^{\pm}_{i\mp 2}-13\mbox{\boldmathF}^{\pm}_{i\mp 1}+47\mbox{\boldmathF}^{\pm}_{i}+27\mbox{\boldmathF}^{\pm}_{i\pm 1}-3\mbox{\boldmathF}^{\pm}_{i\pm 2})/60. To suppress numerical oscillation caused by linear reconstruction methods, one employs a nonlinear method such as the WENO scheme (Jiang & Shu, 1996), and (iii) the numerical flux is obtained at the midpoint, \mbox{\boldmath\hat{F}}_{i+1/2}=\mbox{\boldmath\hat{F}}^{+}_{i+1/2}+\mbox{\boldmath\hat{F}}^{-}_{i+1/2}, and then Eq. (24) is integrated through Eq. (26).

Steger & Warming (1981) developed a flux-vector splitting method for hydrodynamic simulations, but their method cannot be utilized for MHD simulations because the MHD equations are not homogeneous of degree one, \mbox{\boldmathF}\neq A\mbox{\boldmathU}, where A={\partial\mbox{\boldmathF}}/{\partial\mbox{\boldmathU}} is the Jacobian matrix. One may use the Lax-Friedrichs method, \mbox{\boldmathF}^{\pm}_{i}=\left(\mbox{\boldmathF}_{i}\pm|\lambda|_{\rm max}\mbox{\boldmathU}_{i}\right)/2, where $|\lambda|_{\rm max}$ is the maximum absolute speed among eigenmodes. Although the Lax-Friedrichs method is applicable to any hyperbolic conservation laws, a solution suffers from numerical diffusion with a time scale of $\Delta x/|\lambda|_{\rm max}$ that may be much faster than a dynamical time scale (Minoshima et al., 2015).

2.1 Differencing

We adopt another type of high-order finite difference schemes, termed as the Weighted Compact Nonlinear Scheme (WCNS; Deng & Zhang, 2000). This scheme is built based on the compact finite difference scheme in Lele (1992), which expresses a linear combination of flux derivatives around a node as a function of fluxes at midpoints (Eq. (4) in Deng & Zhang (2000)). For low computational cost and simplicity in parallel computation, we use an explicit form to approximate the flux derivative at the node $i$ as

[TABLE]

where we consider uniform grid spacing. Here, the pointwise numerical flux \mbox{\boldmath\tilde{F}}_{i+1/2} is a high-order approximation of the physical flux at the midpoint \mbox{\boldmathF}(\mbox{\boldmathU}(x_{i+1/2})) itself. The truncation error terms $O(\Delta x^{p})$ and $O(\Delta x^{q})$ come from the approximations of the numerical flux and the derivative, and the spatial accuracy depends on both. Conservative property is guaranteed as long as the coefficients $d_{k}$ are constant. We use the fourth-order linear central difference,

[TABLE]

Sixth-, eighth- and tenth-order differences are found in Zhang et al. (2008). At the leftmost interior point $i=1$, one can use the third-order boundary scheme,

[TABLE]

An expression for the rightmost interior point is symmetrically derived. The flux derivative is approximated by the numerical fluxes only similar to Deng & Zhang (2000) and Del Zanna et al. (2007). Meanwhile, one may employ a combination of the physical fluxes at the nodes and the numerical fluxes at the midpoints that makes a stencil smaller (Nonomura & Fujii, 2013; Chen et al., 2016).

This scheme allows us to interpolate the physical variables (Deng & Zhang, 2000) or the physical fluxes (Zhang et al., 2008) for the calculation of the numerical fluxes at the midpoints, whereas the conventional finite difference scheme directly reconstructs the numerical fluxes at the midpoints from the physical fluxes at the nodes. The present scheme adopts the following procedure: (i) the physical variables are interpolated to the left side at $x_{i+1/2}$ and the right side at $x_{i-1/2}$, \mbox{\boldmath\tilde{U}}^{L}_{i+1/2} and \mbox{\boldmath\tilde{U}}^{R}_{i-1/2}, in an upwind-biased stencil $(i-r,\dots,i+r)$. They should be a high-order approximation of \mbox{\boldmathU}(x_{i\pm 1/2}), (ii) using \mbox{\boldmath\tilde{U}}^{L,R}_{i+1/2} from Step (i) as the initial left and right Riemann states, the numerical flux is evaluated at the midpoint by a Riemann solver, \mbox{\boldmath\tilde{F}}_{i+1/2}, which is a high-order approximation of \mbox{\boldmathF}(\mbox{\boldmathU}(x_{i+1/2})) itself, and (iii) Eq. (24) is integrated through d\mbox{\boldmathU}^{n}_{i}/dt=-{\cal D}_{4}(\mbox{\boldmath\tilde{F}})_{i}.

Since the present scheme interpolates the physical variables to the midpoints, one can utilize a variety of Riemann solvers, some of which resolve multiple eigenmodes in the Riemann fan that are smeared in the Lax-Friedrichs method. This is a main advantage against the conventional finite difference scheme for MHD simulations. A disadvantage of the scheme is the increase of the number of grid points necessary to update. The explicit WCNS scheme achieves the $(2r+1)$th order of accuracy by combining a $(2r+1)$th-order interpolation and a $(2r+2)$th-order difference. The $(2r+2)$th-order flux difference at the node $i$ is approximated by the numerical fluxes at the midpoints $(i-1/2-r,\dots,i+1/2+r)$, and the numerical flux at the rightmost midpoint $(i+1/2+r)$ is reconstructed in the stencil $(i,\dots,i+2r+1)$. Therefore, the $(2r+1)$th-order explicit WCNS scheme needs $(4r+3)$ grid points to update the physical variables at the node $i$. The example grid spacing for $r=1$ is displayed in Figure 1. On the other hand, the $(2r+1)$th-order conventional finite difference scheme needs $(2r+3)$ grid points.

2.2 Interpolation

In this subsection, we consider a scalar linear advection case in which $F$ is proportional to $U$, e.g., $F=aU$ for simplicity. The scheme calculates the physical variable at the midpoint $i+1/2$ in an upwind-biased stencil $(i-r,\dots,i+r)$. The optimal third-order interpolation of ${\tilde{U}}^{L}_{i+1/2}$ from the point values ${U}_{i}$ in the stencil $(i-1,\dots,i+1)$ is

[TABLE]

The equation for ${\tilde{U}}^{R}_{i-1/2}$ is symmetric with respect to $i$. Assuming the advection velocity is $+1$, i.e., $F=U$, ${\tilde{F}}_{i+1/2}$ in Eq. (2.1) is replaced by ${\tilde{U}}^{L}_{i+1/2}$ and then Eq. (2.1) reduces to

[TABLE]

The combination of the third-order interpolation in Eq. (2.2) and the fourth-order difference in Eq. (2.1) has the third-order dissipation as a leading error. Similarly, the fifth-order interpolation in the stencil $(i-2,\dots,i+2)$ is

[TABLE]

and Eq. (2.1) reduces to

[TABLE]

The combination of the fifth-order interpolation in Eq. (2.2) and the fourth-order difference in Eq. (2.1) has the fourth-order dispersion as a leading error. One can improve the leading error to be fifth-order dissipation by using a sixth-order difference.

The fifth-order interpolation in Eq. (2.2) is not optimal for this scheme because the order of accuracy is four in spite of the five-point interpolation. Instead, we propose the following fourth-order interpolation that is appropriate to Eq. (2.1),

[TABLE]

Replacing ${\tilde{F}}_{i+1/2}$ in Eq. (2.1) by ${\tilde{U}}^{L}_{i+1/2}$, one find that the $O(\Delta x^{4})$ error cancels and the leading error is fifth-order dissipation,

[TABLE]

In Appendix A, we combine seventh- and sixth-order interpolation methods with a sixth-order difference, giving sixth- and seventh-order finite difference schemes.

In order to confirm the dissipative and dispersive properties of the above interpolation methods in Eqs. (2.2)-(35), we perform the Fourier analysis similar to Lele (1992). Using the periodic function for a scalar variable $U_{i}=\exp(j\omega x_{i})$ (here $j$ is the imaginary unit), we interpolate it to the midpoints and calculate ${\cal D}_{4}(U)$. Subsequently, we evaluate the modified wave number as

[TABLE]

Figure 2 shows the real and imaginary parts of the modified wave number as a function of the actual wave number. Evidently, the exact solution is ${\rm Re}(\omega^{*})=\omega,{\rm Im}(\omega^{*})=0$. The combinations of the third-, fourth-, and fifth-order interpolation methods (Eqs. (2.2), (2.2), and (2.2)) with the fourth-order difference (Eq. (2.1)) are referred to as 3I4D3, 4I4D5, and 5I4D4, where the first digit represents the order of interpolation, the second digit is the order of difference, and the last digit is the designed order of accuracy, respectively. The third- and fifth-order upstream central schemes (UP3 and UP5) are shown for comparison. The imaginary part is not zero (dissipative) because the interpolation methods are upwind biased. The numerical dispersion and dissipation are improved with increasing the order of accuracy of the interpolation. The 4I4D5 scheme improves/degrades the numerical dispersion/dissipation against the 5I4D4 scheme, which can be explained by the fact that the absolute value of the imaginary part in the 4I4D5 scheme is 1.2 times larger than the 5I4D4 scheme, as is found in Eqs. (33) and (35). In terms of the modified wave number, the schemes are slightly better than the same-order upstream central schemes as is mentioned in Deng & Zhang (2000).

In practical compressible hydrodynamic simulations, high-order linear interpolation methods do not work due to severe numerical oscillation around discontinuities. Therefore, the WCNS scheme adopts the same methodology as the WENO scheme to suppress the numerical oscillation. The interpolation is expressed as a convex combination of substencils,

[TABLE]

where ${\tilde{U}}^{L,k}_{i+1/2}$ is the interpolated value in the $k$th substencil and the nonlinear weight $w_{s}$ is a function of the smoothness of the profile measured in the substencil. The nonlinear weight is designed so as to be close to the optimal weight when the profile is sufficiently smooth in all substencils, but nearly zero in the substencil containing discontinuities. For the third-order case (Eq. (2.2)), the interpolations in two substencils $(i-1,i)$ and $(i,i+1)$ are

[TABLE]

We employ the weighted function from the third-order energy stable WENO scheme (Yamaleev & Carpenter, 2009),

[TABLE]

where $\epsilon=10^{-40}$ is a positive real number to avoid the denominator becoming zero. For the fifth- and fourth-order cases (Eqs. (2.2) and (2.2)), the interpolations in three substencils $(i-2,i-1,i)$, $(i-1,i,i+1)$, and $(i,i+1,i+2)$ are found in Zhang et al. (2008) (Eqs. (2-16) - (2-18) therein),

[TABLE]

The weighted function is the same as used in Zhang et al. (2008),

[TABLE]

where we use $\epsilon=10^{-40}$. The smoothness indicator $IS_{s}$ in Eq. (41) measures the average of the square of first- and second-order derivatives in a cell $[x_{i-1/2},x_{i+1/2}]$, whereas the indicator in Deng & Zhang (2000) measures it at a point $x_{i}$. Their difference is very small (factor $13/12$ is replaced by 1 in Deng & Zhang (2000)) and does not affect numerical simulation results in this paper. The accuracy may be further improved by adopting sophisticated weighted functions such as in Henrick et al. (2005), Borges et al. (2008), Ha et al. (2013), and Fan et al. (2014).

Figure 3 shows the numerical solutions of the one-dimensional advection $\partial U/\partial t+\partial U/\partial x=0$ over ten periods. We set the initial condition of a combination of rectangular and triangle waves, a Gaussian profile, and a half ellipse, and use 200 grid points and a CFL number of 0.4; this is the same condition used in Jiang & Shu (1996) and Suresh & Huynh (1997). From top to bottom, we employ a second-order scheme based on MUSCL reconstruction with Monotonized Central flux limiter (MUSCL-MC), the Weighted 3I4D3 scheme (W3I4D3, Eqs. (38) and (39)), and the Weighted 4I4D5 scheme (W4I4D5, Eqs. (40) and (41)). All schemes effectively suppress numerical oscillation that will arise around discontinuities in linear schemes. The MUSCL-MC scheme shows an asymmetric profile toward the direction of advection. The weighted high-order schemes presented here improve numerical dispersion, and then remedy the asymmetry. The W3I4D3 scheme considerably dissipates a profile extremum like TVD schemes (e.g., MUSCL-MC). Using the five-point interpolation, the W4I4D5 scheme succeeds in distinguishing the extremum from the discontinuity and preserving it similar to fifth-order schemes (e.g., Suresh & Huynh, 1997). Note that the solution of the W4I4D5 scheme is very close to that of the fifth-order WENO scheme (Fig. 4.2. in Suresh & Huynh (1997)).

2.3 Characteristic decomposition

Upwind schemes are developed based on the property of the linear advection equation. When one applies the scheme to nonlinear hyperbolic conservation laws, it is more appropriate to interpolate characteristic variables rather than conservative or primitive variables (Qiu & Shu, 2002; Shu, 2009). For the primitive variables in MHD, \mbox{\boldmathV}=(\rho,u,v,w,B_{y},B_{z},P)^{T}, the non-conservative form of Eq. (8)

[TABLE]

is diagonalized by left and right matrices $L$ and $R$, and then the equations reduce to nonlinear advection equations for characteristic variables d\mbox{\boldmathQ}=\mbox{\boldmathL}\cdot d\mbox{\boldmathV}. The matrices consist of left and right eigenvectors of $k$th characteristic variables \mbox{\boldmathl}^{k} and \mbox{\boldmathr}^{k},

[TABLE]

Using the eigenvectors given by Stone et al. (2008), the interpolation is carried out by the following steps: (i) $V$ is calculated from $U$, and then \mbox{\boldmathL},\mbox{\boldmathR} are calculated from $V$, (ii) $V$ is converted to $Q$ in each stencil through \mbox{\boldmathQ}_{i+s}=\mbox{\boldmathL}_{i}\cdot\mbox{\boldmathV}_{i+s} where $s=-r,\dots,r$ covers the stencil, (iii) the characteristic variables at the midpoints, \mbox{\boldmath\tilde{Q}}^{L}_{i+1/2},\mbox{\boldmath\tilde{Q}}^{R}_{i-1/2}, are interpolated from \mbox{\boldmathQ}_{i+s}, and (iv) $Q$ is converted to $V$ through \mbox{\boldmath\tilde{V}}^{L}_{i+1/2}=\mbox{\boldmathR}_{i}\cdot\mbox{\boldmath\tilde{Q}}^{L}_{i+1/2}, \mbox{\boldmath\tilde{V}}^{R}_{i-1/2}=\mbox{\boldmathR}_{i}\cdot\mbox{\boldmath\tilde{Q}}^{R}_{i-1/2}, and then $V$ to $U$.

2.4 Riemann solver

The physical variables are prepared at the midpoints to evaluate the numerical fluxes. A common strategy for the calculation of the numerical flux in hyperbolic conservation laws is to solve the Riemann problem for \mbox{\boldmath\tilde{U}}^{L,R}_{i+1/2} as the initial left and right states. A number of linearized and nonlinear approximate Riemann solvers have been proposed and succeeded in practical MHD simulation studies (Brio & Wu, 1988; Balsara, 1998; Li, 2005; Miyoshi & Kusano, 2005) . Among them, we adopt the HLL-type approximate Riemann solvers because of their simplicity, good performance, and widespread acceptance (Kritsuk et al., 2011). The HLL Riemann solver is one of the simplest solver, in which the solution of the Riemann problem is approximated as a single state bounded by the fastest and slowest waves (Harten et al., 1983). For MHD simulations, we employ the HLLD Riemann solver developed by Miyoshi & Kusano (2005) to improve the accuracy of Alfvén and entropy waves against the HLL Riemann solver.

2.5 Multi-dimension

The finite difference scheme designed to solve hyperbolic conservation laws in a divergence form is extended to multi-dimension by a dimension-by-dimension fashion. Applying the scheme to the MHD equations, the induction equation (3) may be written into the divergence form, \partial\mbox{\boldmathB}/\partial t+\nabla\cdot(\mbox{\boldmathu}\mbox{\boldmathB}-\mbox{\boldmathB}\mbox{\boldmathu})=0. However, this form does not numerically guarantee \nabla\cdot\mbox{\boldmathB}=0 in multi-dimension. The error is accumulated, and could finally crash a simulation run. Evans & Hawley (1988) proposed the Constrained Transport (CT) method to avoid this problem, in which the induction equation is discretized on staggered grid points and is solved in the curl form to preserve \nabla\cdot\mbox{\boldmathB} within a machine precision. The present two-dimensional scheme employs the same grid spacing as the CT method in Figure 4(a) (extension to three dimension is straightforward). The in-plane magnetic field and the out-of-plane electric field components are defined at the staggered grid points as pointwise representation (not face-averaged or line-averaged ones),

[TABLE]

A conventional two-point difference formula for the induction equation preserves \nabla\cdot\mbox{\boldmathB} in second-order accuracy. To ensure \nabla\cdot\mbox{\boldmathB}=0 in the spatial accuracy higher than two, we adopt a similar methodology in Londrillo & Del Zanna (2004) and Del Zanna et al. (2007) to the present scheme. As they mentioned, we should measure the derivative in \nabla\cdot\mbox{\boldmathB} by the same algorithm used for the flux derivative. Therefore, the present scheme measures \nabla\cdot\mbox{\boldmathB} from the staggered magnetic field as

[TABLE]

where ${\cal D}_{4,x}$ and ${\cal D}_{4,y}$ stand for the fourth-order difference (Eq. (2.1)) in the $x$ and $y$ directions. It is obvious that the above \nabla\cdot\mbox{\boldmathB} is preserved in a discretized manner when the magnetic field is updated through

[TABLE]

This is an explicit formula and is essentially same as the method used in Del Zanna et al. (2007). On the other hand, the finite difference scheme in Londrillo & Del Zanna (2004) solves the induction equation for “primary” magnetic field $\hat{B}$, not $B$ itself, which approximates the first derivative as a two-point difference irrespective of the order of accuracy. Since they are related implicitly, their scheme needs a matrix inversion to obtain $B$ from $\hat{B}$. See their Appendix A for detail.

A critical issue is how to calculate the staggered electric field $E_{z,i+1/2,j+1/2}$, being consistent with the solution of the Riemann problem. A straightforward strategy is to approximate it from the solutions of the one-dimensional Riemann problem at the neighboring midpoints. The one-dimensional procedure (Sections 2.2-2.4) returns the numerical fluxes in the $x$ and $y$ directions, \tilde{\mbox{\boldmathF}}_{i+1/2,j} and \bar{\mbox{\boldmathG}}_{i,j+1/2}, where the tilde (bar) notation means that a variable is interpolated in the $x$ $(y)$ direction. A simple calculation method for the staggered electric field is an arithmetic average proposed by Balsara & Spicer (1999), termed as the FluxCT method. A second-order accurate solution is

[TABLE]

where $\tilde{F}^{B_{y}}$ and $\bar{G}^{B_{x}}$ stand for the numerical fluxes of $B_{y}$ and $B_{x}$ in the $x$ and $y$ directions, respectively. A high-order extension of the FluxCT method is provided by Lee & Deane (2009).

The numerical fluxes obtained from the Riemann solver $F,G$ would be expressed as a sum of the non-dissipative terms ${\cal F},{\cal G}$, and the upwind-correction terms ${\Phi},{\Gamma}$ in the $x$ and $y$ directions that contain numerical diffusion to stabilize a simulation,

[TABLE]

Substituting Eq. (61) into Eq. (58), one obtain

[TABLE]

A solution of one-dimensional problems in two dimension with the above two-dimensional method is not equivalent to one with the base one-dimensional scheme. Following Gardiner & Stone (2005), let us consider a one-dimensional problem aligned in the $x$-direction that ignores $j$-dependence, $\bar{\cal G}^{Bx}_{i,j+1/2}={\cal G}^{Bx}_{i,j}=-{F}^{B_{y}}_{i,j}$, and the numerical diffusion in the $y$ direction, $\bar{\Gamma}^{B_{x}}=0$. Using a piecewise-constant approximation, e.g., $\tilde{\cal F}^{B_{y}}_{i+1/2}=\left({F}^{B_{y}}_{i}+{F}^{B_{y}}_{i+1}\right)/2$, Eq. (62) becomes

[TABLE]

Owing to the arithmetic average in Eq. (58), the upwind-correction terms in Eq. (63) are halved from the base numerical flux in Eq. (61). The reduction of upwind property might lead to a solution more oscillatory.

Variants of the FluxCT method have been proposed to take into account for the full upwind property. Balsara & Spicer (1999) themselves considered a directional biasing with respect to the gas pressure gradient rather than the arithmetic average in Eq. (58). Gardiner & Stone (2005) modified Eq. (58) by introducing the first-order derivative of the electric field. With proper approximation of the derivative, their scheme reduces to the base one-dimensional scheme for one-dimensional problems. Lee (2013) incorporated an upwind biasing with respect to the entropy mode into the high-order FluxCT method.

In this paper, we present a simple modification of the high-order FluxCT method so as to be consistent with the base one-dimensional scheme. We explicitly calculate the non-dissipative terms in the numerical flux of the magnetic field (Eq. (61)) through a linear interpolation,

[TABLE]

where the superscripts $B,T$ denote the bottom and top sides at the midpoint, and the left arrow stands for the linear interpolation method (e.g., Eq. (2.2)) with the same order of accuracy as the nonlinear interpolation method used before. The upwind-correction terms are readily obtained from Eq. (61). Subsequently, we interpolate the non-dissipative and upwind-correction terms toward the edge along the orthogonal direction (shown in Figure 4(b)),

[TABLE]

where the interpolation method is same as used in Eq. (2.5). The staggered electric field is calculated through their arithmetic average,

[TABLE]

The first term on the right-hand side of Eq. (66) corresponds to the non-dissipative approximation of the electric field. The second and third terms are the upwind-correction terms in the $x$ and $y$ directions, which are twice as large as the FluxCT method in Eq. (62) and thus give the same amount of numerical diffusion as the base numerical flux in Eq. (61).

Lee (2013) pointed out that a simple arithmetic average for the staggered electric field might lead to spurious oscillations around discontinuous regions, even though the proper amount of numerical diffusion is recovered like Eq. (66). This is observed, for example, when a two-dimensional object flows along the direction almost parallel to one of the coordinate axes. In order to deal with this problem, we modify Eq. (66) by taking a direction-biased average with respect to the Alfvén mode,

[TABLE]

where

[TABLE]

and $\delta=10^{-6}$ is a small positive number to avoid the denominator becoming zero. We take a second-order arithmetic average for the primitive variables to calculate Eq. (68) at the edge for simplicity, e.g., $\rho_{i+1/2,j+1/2}=(\rho_{i,j}+\rho_{i+1,j}+\rho_{i,j+1}+\rho_{i+1,j+1})/4$.

The present method has two main differences from the previous CT variants proposed by Gardiner & Stone (2005) and Lee (2013): the method is able to extend to an arbitrary high order of accuracy whereas the order of accuracy is limited to be two in Gardiner & Stone (2005) (see also Felker & Stone (2018)), and the method reduces to the base one-dimensional scheme for one-dimensional problems whereas it is not always satisfied in Lee (2013). Hereafter, we refer to Eq. (67) as the Central Upwind CT (CUCT) method.

2.6 Summary of numerical procedure

Here, we summarize the numerical procedure for the present scheme. Initially, we define the fluid conservative variables $(\rho,\rho u,\rho v,\rho w,e)$ at each node $(i,j,k)$, and the staggered magnetic field at the midpoint $(B_{x,i+1/2,j,k},B_{y,i,j+1/2,k},B_{z,i,j,k+1/2})$ as pointwise representation. The magnetic field satisfies the divergence-free condition in a discretized manner, (\nabla\cdot\mbox{\boldmathB})={\cal D}_{4,x}(B_{x})+{\cal D}_{4,y}(B_{y})+{\cal D}_{4,z}(B_{z})=0. The magnetic field is calculated at the node from the midpoints by a central interpolation with the same order of accuracy of the scheme. For example, the fourth-order interpolation gives (e.g., Lele, 1992)

[TABLE]

which can improve the accuracy of $B_{x,i}$ against the two-point average. Then, we perform the following steps in each substep of the third-order SSPRK time integration (Eq. (25)).

Primitive variables \mbox{\boldmathV}=(\rho,u,v,w,B_{x},B_{y},B_{z},P)^{T} are calculated from the conservative variables $U$ at each node. 2. 2.

Characteristic variables \mbox{\boldmathQ}_{i,j,k} are converted from the primitive variables in the stencil $(i-r,\dots,i+r)$ along the $x$ direction (Section 2.3). They are interpolated to the left and right sides at the right and left midpoints by the selected method (Section 2.2), \mbox{\boldmath\tilde{Q}}^{L}_{i+1/2,j,k} and \mbox{\boldmath\tilde{Q}}^{R}_{i-1/2,j,k}, and then are converted to the primitive variables. 3. 3.

The approximate HLLD Riemann solver calculates the numerical flux in the $x$ direction, \tilde{\mbox{\boldmathF}}_{i+1/2,j,k} (Section 2.4). 4. 4.

The above one-dimensional procedure (Step 2-3) is performed along the $y$ and $z$ directions. 5. 5.

The staggered electric field $E_{z,i+1/2,j+1/2,k}$ is calculated by the selected CT method (Section 2.5). The $x$ and $y$ components, $E_{x,i,j+1/2,k+1/2}$ and $E_{y,i+1/2,j,k+1/2},$ are obtained in the same manner. 6. 6.

The flux derivative is approximated by the fourth-order central difference (Eq. (2.1)), and then the fluid conservative variables are updated in the divergence form (Eqs. (1), (2), (4)). The staggered magnetic field is updated in the curl form (Eq. (3)), and then it is interpolated to the node by Eq. (69).

A higher-order interpolation method tends to make a solution more oscillatory that might lead to negative density or pressure and crash a simulation run. In order to suppress the numerical oscillation around shocks, we employ a lower-order method (e.g., the third-order method in Eqs. (38) and (39) as opposed to the fifth-order method in Eqs. (40) and (41)) to the left and right state interpolations at the right and left midpoints (e.g., \mbox{\boldmath\tilde{Q}}^{L}_{i+1/2,j,k} and \mbox{\boldmath\bar{Q}}^{R}_{i-1/2,j,k}), when a node $(i,j,k)$ is judged as a shocked-node based on (\nabla\cdot\mbox{\boldmathu})_{i,j,k}<0. We apply this hybrid method to the Orszag-Tang vortex problem in Section 3.3 and a blast wave problem in Section 3.5.

3 Numerical tests

This section presents numerical simulation results of MHD test problems including propagation, shock, and fundamental plasma phenomena to discuss the capability of the present scheme. In the following benchmark tests, we employ the W4I4D5 scheme in Eqs. (40) and (41), the HLLD Riemann solver, and the CUCT method in Eq. (67) for the divergence-free treatment as a fiducial scheme, referred to as the W4I4D5-HLLD-CUCT scheme. In Sec. 3.6 and 3.8, we also use the CUCT-AA method that adopts the arithmetic average in Eq. (66) as opposed to the direction-biased average in Eq. (67). We solve two- or three-dimensional ideal MHD equations with a specific heat ratio $\gamma=5/3$ and a CFL number of 0.4 unless otherwise stated.

3.1 Circularly polarized Alfvén waves

The circularly polarized Alfvén wave is utilized to check the accuracy of the designed code since the wave is an analytic solution for the MHD equations. Improving the accuracy of the Alfvén wave solution is of great importance to tackle with MHD turbulent flows. We test a problem similar to Tóth (2000). The two-dimensional computational domain is periodic with $0\leq x<1/\cos\alpha$ and $0\leq y<1/\sin\alpha$, where $\alpha=30^{\circ}$ is the wave propagation angle relative to the $x$ axis. The domain is resolved by $N\times N$ grid points and the grid number increases as $N=16,32,64,128,256,512$. The initial condition is

[TABLE]

where the subscripts $\parallel,\perp$ stand for the parallel and perpendicular directions relative to the wave propagation. They are related to the variables in an original coordinate system through $x_{\parallel}=x\cos\alpha+y\sin\alpha,(u,B_{x})=(u_{\parallel},B_{\parallel})\cos\alpha-(u_{\perp},B_{\perp})\sin\alpha,$ and $(v,B_{y})=(u_{\parallel},B_{\parallel})\sin\alpha+(u_{\perp},B_{\perp})\cos\alpha$. The wave amplitude is $\delta V=\delta B=0.01$ so that the wave propagates in negative direction with the Alfvén speed $V_{A}=B_{\parallel}/\sqrt{\mathstrut{\rho}}=1$. The simulations are conducted by five schemes with different interpolation methods, the first order, MUSCL-MC, W3I4D3, W4I4D5, and W5I4D4. The first-order Euler time integration is employed for the first-order scheme, the second-order SSPRK time integration for the MUSCL-MC scheme, and the third-order SSPRK time integration for the remaining schemes. In order to assess the numerical error attributed to the spatial discretization, the CFL number is changed as ${\rm CFL}=16\sqrt{\mathstrut{3}}/5N$ $(\Delta t\propto\Delta x^{2})$ so that the temporal discretization errors in the five schemes are negligible compared to the spatial discretization errors.

We measure numerical errors of the in-plane and out-of-plane magnetic fields, $B_{\perp}$ and $B_{z}$, at $t=2$. Table 1 summarizes the $L_{1}$ errors and the orders in parenthesis. Except for the W3I4D3 scheme, the four schemes converge at the designed order of accuracy. The W3I4D3 scheme is slightly better than the MUSCL-MC scheme, but it does not achieve the optimal third order of accuracy even with the finest grid points. The nonlinear weight for this scheme (Eq. (39)) may not converge to the optimal weight, even though the unit wavelength is resolved by 512 grid points. The W3I4D3 scheme requires 1.5 times longer computational time (due to the order of the Runge-Kutta method) and a larger stencil than the MUSCL-MC scheme.

The errors in the W4I4D5 and W5I4D4 schemes are remarkably improved from the three lower-order schemes. The W4I4D5 scheme achieves fifth order of accuracy in the out-of-plane $(B_{z})$ case whereas W5I4D4 scheme converges at fourth order in both the in-plane $(B_{\perp})$ and out-of-plane cases, because the W4I4D5 scheme is designed so as to satisfy the fifth order for the linear advection problem (Eq. (35)). In the in-plane case, the scheme calculates the staggered electric field by averaging the fourth-order accurate solutions of the one-dimensional Riemann solver in the two-dimensional plane (Eq. (67)). This limits the order of accuracy of the W4I4D5 scheme to be four in the in-plane case (and in a three-dimensional case as well). Nevertheless, the error in the W4I4D5 scheme is slightly smaller than the W5I4D4 scheme at least in this problem. The computational times of the two schemes are almost equal, and about two times longer than the MUSCL-MC scheme.

3.2 Kelvin-Helmholtz instability

The boundary layer is a common feature in a variety of plasma environments. For example, the magnetopause in the planetary magnetosphere is a boundary between the magnetospheric and the solar wind plasmas. The Kelvin-Helmholtz instability is thought to contribute to the transport of the tailward-flowing solar wind plasma into the Earth magnetosphere (e.g., Hasegawa et al., 2004). We simulate the linear growth of the Kelvin-Helmholtz instability so as to assess the capability to resolve a shear flow. The two-dimensional computational domain ranging $0\leq x<L$ and $-L/2\leq y<L/2\;(L=14.0)$ is resolved by $N\times N$ grid points. The boundary condition is periodic and symmetric in the $x$ and $y$ directions. The initial condition has a velocity shear, $u=-(V_{0}/2)\tanh(y/\lambda)$, uniform density and pressure $\rho=\rho_{0}=1.0,P=P_{0}=500$, and a uniform out-of-plane magnetic field $B_{z}=B_{0}=1.0$, where $\lambda=1.0$ and $V_{0}=1.0$ is normalized by the Alfvén speed. We use $\gamma=2.0$. The system is essentially hydrodynamic because the in-plane magnetic field is kept zero. To initiate the instability, we impose a small $(1\%)$ perturbation to the $y$-component of the velocity around the boundary with a wavelength equal to $L$, which corresponds to the fastest growing mode under the adopted initial condition (Miura & Pritchett, 1982). The simulations are conducted by the MUSCL-MC and the W4I4D5 schemes.

Figure 5 shows the linear growth of the fastest growing mode from the simulations and the theory. The theoretical growth rate, $\omega_{i}=0.095V_{0}/\lambda$, is obtained by numerically solving the linearized equations for the velocity and the total pressure $P_{t}=P+B_{z}^{2}/2$ as an eigenvalue problem,

[TABLE]

where $P_{t0}=P_{0}+B_{0}^{2}/2=500.5$, $k_{x}=2\pi/L$, and $i$ is the imaginary unit.

The W4I4D5 scheme converges to the theory at resolutions $N\geq 64$, while the MUSCL-MC scheme requires $N=128$ for convergence. The computational time of the W4I4D5 scheme with $N=64$ is about 4 times faster than the MUSCL-MC scheme with $N=128$. This problem demonstrates that the present high-order scheme is highly cost effective for the shear flow against the low-order scheme.

3.3 Orszag-Tang vortex

We conduct the so-called Orszag-Tang vortex problem, which is a well-known two-dimensional MHD test so as to verify the capability to capture multiple interactions of shock waves and vortices (Orszag & Tang, 1979). The periodic computational domain with $0<x<2\pi$ and $0<y<2\pi$ is resolved by $N\times N$ grid points. The initial condition is $(\rho,u,v,w,B_{x},B_{y},B_{z},P)=(\gamma^{2},-\sin(y),\sin(x),0,-\sin(y),\sin(2x),0,\gamma)$. We solve the problem with the MUSCL-MC scheme, the hybrid W3I4D3-W4I4D5 scheme (see Section 2.6), and the W4I4D5 scheme at the resolution of $N=100,200,400,800$. The solution with the W4I4D5 scheme at $N=1600$ is employed as a reference.

Figure 6 (a)-(d) compares the two-dimensional profile of the temperature $T=2P/\rho$ at $t=\pi$ obtained from the three schemes with $N=200$ and the reference solution. The overall structure, such as cold and hot regions and the position of discontinuities, agrees well with previous works (e.g., Fig. 12 in Miyoshi & Kusano, 2005). The three schemes give almost identical solutions at this period, implying that the spatial order of accuracy is degraded and the quality of solution less depends on the designed optimal order of accuracy when the discontinuities become dominant. Panels (e) and (f) compare the one-dimensional profiles at $y=0.64\pi$ and $y=\pi$ among different resolutions with the W4I4D5 scheme. Compared to the reference, the profile at $N=100$ (green cross) is degraded especially around spike regions (e.g., $(x,y)=(0.4,0.64\pi)$ in panel (e)), and shows a spurious oscillation around $(x,y)=(0,\pi)$ in panel (f). These errors are found also in the MUSCL-MC and hybrid W3I4D3-W4I4D5 schemes, and are largely improved at resolutions $N\geq 200$.

We measure relative differences of the temperature profile at different resolutions from the reference one $T_{\rm ref}$, defined as $\delta T=\sum|T-T_{\rm ref}|/\sum T_{\rm ref}$. Table 2 summarizes the relative difference and the rate of convergence in parenthesis at two different periods $t=0.2\pi$ and $t=\pi$. The MUSCL-MC and W4I4D5 schemes converge at the designed order of accuracy when the profile is smooth at $t=0.2\pi$. In contrast, the rate of convergence of the hybrid W3I4D3-W4I4D5 scheme is worse than the optimal fourth order at $t=0.2\pi$. As is presented in Section 2.6, this hybrid method adopts the low-order interpolation when the flow is divergent irrespective of the degree of compressibility: this is a crude criterion. Actually, the flow is compressible in $38\%$ of the domain at this period that would limit the overall convergence. To sophisticate the method, for example, we may modify the criterion based on the strength of shocks. When the discontinuities are dominant at $t=\pi$ in Figure 6, their relative differences are close and convergence rate drops to less than two.

3.4 Oblique Shock

The shock tube problem is a standard test for hydrodynamic simulations so as to examine the accuracy and robustness for discontinuities and rarefaction waves. In addition, divergence-free preservation is a critical issue for multidimensional MHD shocks since it tends to be violated especially at shocks in non-divergence-free schemes. We test the familiar Brio-Wu shock tube problem (Brio & Wu, 1988) in two dimensions. The initial left and right states are $[\rho,u_{\parallel},u_{\perp},w,B_{\perp},B_{z},P]=[1,0,0,0,1,0,1]$ and $[0.125,0,0,0,-1,0,0.1]$, where the subscripts have the same meaning in Section 3.1. We use $B_{\parallel}=0.75$ and $\gamma=2.0$. The shock tube is tilted by angles $\alpha=\tan^{-1}(1/2)=26.6^{\circ}$ and $45^{\circ}$ relative to the $x$ axis so as to assess the divergence-free condition. The domain of the shock tube $0\leq x<\cos\alpha$ and $0\leq y<\sin\alpha$ is resolved by $N\times N/4$ (for $\alpha=26.6^{\circ}$) and $N\times N$ (for $\alpha=45^{\circ}$) grid points where $N=200$. In order to apply the periodic boundary condition, the computational domain is extended to $0\leq x<2/\cos\alpha$ and $0\leq y<2/\sin\alpha$, and the whole domain is resolved by $mN\times mN$ grid points where $m=2.5$ for $\alpha=26.6^{\circ}$ and $4$ for $\alpha=45^{\circ}$ (Kawai, 2013). The initial condition is set to be the left state in $0\leq x_{\parallel}<0.5,1.5\leq x_{\parallel}<2.5,3.5\leq x_{\parallel}<4$ ($x_{\parallel}=x\cos\alpha+y\sin\alpha$), and the right state otherwise. The initial jump condition is given as a hyperbolic tangent function with a finite thickness of $0.5\Delta x/\cos\alpha$.

Figure 7 shows the primitive variables at angles of $\alpha=26.6^{\circ}$ (blue) and $45^{\circ}$ (red) solved by the HLLD-CUCT scheme at $t=0.1$. The both cases capture the expected discontinuities and rarefaction waves without numerical instabilities (from left to right; the fast rarefaction wave, the slow compound wave, the contact discontinuity, the slow shock, and the fast rarefaction wave), and agree well with the one-dimensional reference solution with 2000 grid points (black). Note that the two-dimensional scheme for $\alpha=0^{\circ}$ reduces to the base one-dimensional scheme. The bottom left panel shows the deviation of the parallel magnetic field component from the initial value, $\Delta B_{\parallel}$. The parallel magnetic field is preserved extremely well within an error of $\sim 10^{-8}$ at this period. We find that the \nabla\cdot\mbox{\boldmathB} error is preserved within an error of $\sim 10^{-11}$ (not shown here).

In order to compare the divergence-free and the non-divergence-free schemes, we combine the HLLD scheme with the hyperbolic divergence cleaning method (GLM; Dedner et al., 2002). This scheme introduces a scalar variable $\psi$ to diminish the \nabla\cdot\mbox{\boldmathB} error from a simulation domain through propagation and diffusion. The in-plane magnetic field and $\psi$ obey the following equations,

[TABLE]

where $c_{h}$ is the propagation speed of the additional eigenmode and $c_{p}^{2}$ is the diffusion coefficient. In Figure 8, this scheme also captures discontinuities and rarefaction waves as the HLLD-CUCT scheme. The main difference is the parallel magnetic field (bottom left panel) with an error of $\sim 1\times 10^{-2}$ mainly at the front of the compound shock and the slow shock. The error comes from the equation of the parallel magnetic field,

[TABLE]

The scalar variable $\psi$ may have a finite value at \nabla\cdot\mbox{\boldmathB}\neq 0. Since the grid spacing $(\Delta x,\Delta y)=(2/mN\cos\alpha,2/mN\sin\alpha)$ is set so that the shock surface is uniform and parallel to the diagonal axis in the cell $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$, the variables are the left state at $(i-1/2,j),(i,j-1/2)$ and the right state at $(i+1/2,j),(i,j+1/2)$. Therefore, Eq. (90) is discretized as

[TABLE]

It may be finite around discontinuities. In Appendix B, we examine the dependence of the solution on the parameter in the hyperbolic divergence cleaning method ($c_{h}$ and $c_{p}$).

3.5 Blast wave in strongly magnetized medium

This test problem treats the propagation of strong MHD shocks in a two-dimensional strongly magnetized medium so as to assess the robustness of a scheme. Many similar problems have been tested as a simple example of astrophysical phenomena (Londrillo & Del Zanna, 2000; Balsara & Kim, 2004; Mignone et al., 2007; Stone et al., 2008). The periodic computational domain with $-2<x<2$ and $-2<y<2$ is resolved by $1024\times 1024$ grid points. We adopt the initial condition used in Londrillo & Del Zanna (2000): the ambient medium $(\rho,u,v,w,B_{x},B_{y},B_{z},P)=(1,0,0,0,B_{0}\sin(\theta),B_{0}\cos(\theta),0,1)$ where $B_{0}=10$, $\theta=45^{\circ}$, and then a high pressure cylinder imposed at the center of the domain, $P=100$ for $\sqrt{\mathstrut{x^{2}+y^{2}}}<0.125$. The plasma beta in the ambient medium is very low, \beta=2P/|\mbox{\boldmathB}|^{2}=0.02, thus this is a stringent condition for conservative schemes. We compare the CUCT method with the GLM method (Eqs. (88) and (89)) and the high-order FluxCT method (Eq. (66) with replacing the denominator in the second and third terms on the right-hand side by 4), so as to examine the effects of the finite \nabla\cdot\mbox{\boldmathB} error and the amount of numerical diffusion for the CT method on the quality of numerical solution in the presence of strong low-$\beta$ shocks. We use the hybrid W3I4D3-W4I4D5 scheme (Sec. 2.6), because the W4I4D5 scheme is found not to be enough to suppress the numerical oscillation in this stringent problem irrespective of the divergence-free treatment.

Figure 9 (a)-(c) shows the profile of the magnetic energy |\mbox{\boldmathB}|^{2}/2 at $t=0.02$ for the HLLD-GLM, the HLLD-FluxCT, and the HLLD-CUCT schemes. The three schemes capture complex shocks and rarefaction waves without numerical instabilities, and they are in good agreement with the previous work (Fig. 13 in Londrillo & Del Zanna (2000)), indicating that the present schemes well preserve the axial symmetry even though the ambient magnetic field is oblique relative to the coordinate axis. At this period, the three schemes give the similar result.

The difference among the schemes becomes prominent as the system evolves. Figure 9 (d)-(f) shows the gas pressure profile at $t=0.1$ for the three schemes. The gas pressure is susceptible to a numerical error of the magnetic field especially at low-$\beta$ shock regions, since the present schemes calculate it by subtracting the magnetic energy from the total energy density (Eq. (6)). The HLLD-GLM scheme (d) shows spurious profiles at fast shock fronts. In panel (g), we check the one-dimensional profile along the outermost fast shock front (denoted by dashed lines in panels (d)-(f)). The HLLD-GLM scheme suffers from the oscillation in both the downstream and upstream regions of the shock, at which a \nabla\cdot\mbox{\boldmathB} error of $O(1)$ is observed. The oscillatory profile in the gas pressure is attributed to the finite \nabla\cdot\mbox{\boldmathB} production at the shock and its propagation through the GLM method. The oscillation is seen at the downstream region also in the two CT schemes, but its amplitude is reduced especially in the HLLD-CUCT scheme compared to the HLLD-GLM scheme. We then suggest that the CT methods are better suited for low-$\beta$ shocks than the GLM method.

The HLLD-FluxCT scheme in Figure 9 (e) shows spurious profiles at the slow shock front (along the dot-dashed line). In Figure 9 (h), we check the one-dimensional profile along this shock front. The HLLD-FluxCT scheme suffers from an intolerable level of numerical oscillation that leads to a negative pressure and finally crashes the simulation. The HLLD-CUCT scheme largely suppresses this oscillation due to the recovery of the proper amount of numerical diffusion for the induction equation against the HLLD-FluxCT scheme. The \nabla\cdot\mbox{\boldmathB} error in the HLLD-GLM scheme is relatively small, $O(10^{-1})$, compared to that at the fast shock, thus it would give a less oscillatory solution as the HLLD-CUCT scheme. This problem demonstrates that the CUCT method can improve the robustness of the divergence-free MHD scheme for low-$\beta$ flows.

3.6 Field loop advection

We conduct a stringent problem of the advection of a weakly magnetized loop proposed by Gardiner & Stone (2005), so as to assess the capability to capture a tangential discontinuity in a multidimensional flow. The two-dimensional computational domain is periodic with $-1\leq x<1$ and $-0.5\leq y<0.5$, and is resolved by $256\times 128$ grid points. The initial condition is

[TABLE]

where $v_{0}=\sqrt{\mathstrut{5}},A_{0}=10^{-3},R=0.3,r=\sqrt{\mathstrut{x^{2}+y^{2}}}$. The in-plane magnetic field is initialized by numerically calculating the curl of the vector potential through Eq. (2.1) so as to satisfy \nabla\cdot\mbox{\boldmathB}=0 in a discretized form (Eq. (55)). The ambient gas pressure $P_{0}$ is set to be much higher than the magnetic pressure in the loop so that the loop behaves almost like a passive scalar. The loop flows toward the positive direction with an angle $\theta=\tan^{-1}(1/2)$ relative to the $x$ axis. The simulations are conducted by three different choices of the Riemann solver and the divergence-free treatment, the HLLD-GLM, the HLL-CUCT, and the HLLD-CUCT schemes.

Figure 10 shows the magnetic pressure |\mbox{\boldmathB}|^{2}/2 for the case of $P_{0}=1$ at $t=8$. The in-plane flow is supersonic, $v_{0}>c_{s}=\sqrt{\mathstrut{\gamma P/\rho}}=1.3$. All schemes successfully preserve the initial symmetric loop structure. Numerical errors are mainly observed around the leading edge of the loop at $(x,y)=(0.2,0.2)$. One-dimensional profile at $x=0.2$ in panel (d) shows that the HLLD-CUCT scheme slightly distorts the profile at $y=0.15$. The thickness of the leading edge is nearly equal among the three schemes, but it is broadened at the side edge of the loop $(x,y)=(0.2,-0.2)$ in the HLLD-GLM scheme.

Gardiner & Stone (2008) and Lee & Deane (2009) argued that the in-plane geometry could be distorted by the \nabla\cdot\mbox{\boldmathB} error through \partial B_{z}/\partial t\propto w\nabla\cdot\mbox{\boldmathB}. The HLL-CUCT and HLLD-CUCT schemes keep \nabla\cdot\mbox{\boldmathB} within $\sim 10^{-14}$, and the amplitude of the out-of-plane magnetic field within $\sim 10^{-10}$. On the other hand, the \nabla\cdot\mbox{\boldmathB} error and the amplitude of the out-of-plane magnetic field are around $1\times 10^{-4}$ and $2\times 10^{-5}$ (2 % of the in-plane magnetic field) at $t=32$, and continue to increase in the HLLD-GLM scheme.

The problem is more stringent for a higher plasma beta (weaker compressibility) case. Figure 11 shows the case of $P_{0}=10^{4}$ so that the in-plane flow is subsonic $(v_{0}\ll c_{s}=130)$. The loop structure is seriously distorted in the HLLD-GLM and the HLL-CUCT schemes owing to the overestimation of the amount of numerical diffusion. The time scale of the numerical diffusion is proportional to the transit time of the fast mode wave in the HLL-CUCT scheme (Minoshima et al., 2015). Therefore, the magnetic loop quickly diffuses outward and inward in the course of the flow, and the inward diffusion leads to the cancellation of the anti-parallel field at the center of the loop. As a result, the dilute ring-shaped profile is observed in the HLL-CUCT scheme (Figure 11 (b)).

The unphysical solution in the HLLD-GLM scheme (Figure 11 (a)) is attributed to the numerical flux for the additional equations (Eqs. (88) and (89)), which is evaluated by the Lax-Friedrichs flux splitting method (Eq. (42) in Dedner et al. (2002)). The method introduces the numerical diffusion of $B_{x}$ $(B_{y})$ in the $x$ $(y)$ direction, and its coefficient is proportional to $c_{h}$. The propagation speed $c_{h}$ is determined from the CFL condition so as to quickly diminish the \nabla\cdot\mbox{\boldmathB} error, and thus it is much faster than the flow speed in this problem that leads to rapid numerical diffusion of $B_{x}$ in the $x$ direction observed as the horizontal structure. Similarly, the vertical structure corresponds to the numerical diffusion of $B_{y}$ in the $y$ direction. To avoid this unphysical solution, one has to carefully switch the Lax-Friedrichs flux splitting to a simple average so that the spatial derivatives in Eqs. (88) and (89) reduce to the central difference.

In contrast to the above two schemes, the HLLD-CUCT scheme (Figure 11 (c)) can preserve the symmetric loop structure. The thickness of the discontinuity is nearly unchanged from the former case with $P_{0}=1$ (see panel (d) in Figures 10 and 11). Although the solutions of the HLLD-GLM and the HLL-CUCT schemes are slightly better in the supersonic case (Figure 10), the HLLD-CUCT scheme has a capability to capture the tangential discontinuity in the multidimensional flow irrespective of the plasma beta (compressibility) by virtue of the combination of the multi-state Riemann solver and the divergence-free method.

Additionally, we conduct a small angle advection problem for a three-dimensional loop proposed by Lee (2013). The three-dimensional periodic computational domain with $-0.5\leq x<0.5,-0.5\leq y<0.5$, and $-1\leq z<1$ is resolved by $100\times 100\times 200$ grid points. We adopt the initial condition of the magnetic field used in Gardiner & Stone (2008): the loop is tilted around the $y$ axis by an angle of $\omega=\tan^{-1}(1/2)$, and its magnetic field components are described by the curl of the vector potential in the coordinate system $(x_{1},x_{2},x_{3})=(x\cos\omega+z\sin\omega,y,-x\sin\omega+z\cos\omega)$,

[TABLE]

where $r=\sqrt{\mathstrut{x_{1}^{2}+x_{2}^{2}}}$. For fluid variables, we set $\rho=P=1$, $(u,v,w)=(\cos\theta,\sin\theta,2)$ and $\theta=\tan^{-1}(0.01)$ so that the loop flows almost in the $x$-$z$ plane. Figure 12 compares the result at $t=1.0$ with two different CUCT methods; (a) the arithmetic average (Eq. (66), the HLLD-CUCT-AA scheme), and (b) the direction-biased average (Eq. (67), the HLLD-CUCT scheme). Obviously, the HLLD-CUCT-AA scheme suffers from the numerical oscillations in the $x$-$y$ and $y$-$z$ planes that arise at the side edge of the loop and then propagate inward. As is pointed out by Lee (2013), this is because the scheme takes the electric field average across the discontinuity in the $y$ direction (first term in Eq. (66)), even though the numerical diffusion is biased in the flow direction ($x$ or $z$). The oscillation is largely suppressed by the HLLD-CUCT scheme, because the scheme biases the electric field with respect to the direction of the Alfvén mode in each two-dimensional plane.

3.7 Richtmyer-Meshkov instability

The impact of a shock against a corrugated contact discontinuity can induce a shear flow at the interface and drive the so-called Richtmyer-Meshkov instability (RMI, Richtmyer, 1960; Meshkov, 1972). Such a situation is expected in astrophysical environments, for example, the propagation of a supernova shock in the inhomogeneous interstellar medium (ISM). The RMI in magnetized plasmas has a potential to amplify the magnetic field beyond the compression by a single shock, and thus it is of great interest to astrophysics. Sano et al. (2012, 2013) conducted two-dimensional MHD simulations and argued that the RMI can be a powerful candidate for the mechanism of the magnetic field amplification at supernova remnants, at which strong field regions are localized and their maximum strength reaches $100$ times larger than the typical value in the ISM. For practical application of the present scheme, we simulate the nonlinear evolution of the RMI. The two-dimensional computational domain is taken in the rest frame of a shock at $y=0$. The initial state in the upstream region $y>0$ is $(\rho,u,v,w,B_{x},B_{y},B_{z},P)=(1.0,0,-1.0,0,0.00034641,0,0,0.006)$ so that the Mach number M=|\mbox{\boldmathu}|/c_{s}=10 and the plasma beta \beta=2P/|\mbox{\boldmathB}|^{2}=10^{5}. The initial state in the downstream region $y<0$ is calculated from the Rankine-Hugoniot condition for perpendicular MHD shocks, $(\rho,u,v,w,B_{x},B_{y},B_{z},P)=(3.8835,0,-0.25750,0,0.0013453,0,0,0.74850)$. A corrugated contact discontinuity is imposed in the upstream region, $y_{\rm cd}=Y_{0}+\psi_{0}\cos(2\pi x/\lambda)$, where $Y_{0}=1.0$, $\psi_{0}=0.1$ is a corrugation amplitude, and $\lambda=1.0$ is the wavelength. We consider that the density is decreased to $\rho=0.01$ behind the contact discontinuity $y>y_{\rm cd}$ (see Figure 13(a)). The domain is extended in the $y$ direction so as to ignore the boundary effect, $0\leq x<\lambda$ and $-80\lambda\leq y<80\lambda$, and is resolved by $N\times 160N$ grid points. We use $N=64$ for a fiducial run. The boundary condition is periodic in the $x$ direction, and is fixed to be the initial state in the $y$ direction.

The contact discontinuity flows in the negative $y$ direction and strikes the shock at $t=1$. Then, a transmitted shock and a reflected rarefaction wave propagate in the positive and negative $y$ directions. Due to the corrugation of the contact discontinuity, the surfaces of these waves are also rippled that induce the shear flow at the interface. Subsequently, the RMI grows and the interface develops nonlinearly. Figure 13(b) shows the density profile at $t=20$. Arrows represent the direction of the rotational velocity, which is obtained by decomposing the velocity field into rotational, compressible, and mean components, \mbox{\boldmathu}_{\rm R},\mbox{\boldmathu}_{\rm C},\mbox{\boldmathu}_{M}. They satisfy

[TABLE]

where an overline means spatial averaging. The mushroom-shaped spike develops with a height of $20$ times the initial corrugation amplitude. The primary vortex in $9.5<y<11$ forms the spike, and its velocity is approximately $10-20\%$ of the upstream velocity. In addition, the secondary vortex and the associated ring-shaped profile are visible in $10<y<10.5$. The rotational flow is observed only around the spike. Compared to the rotation, the compressible flow does not show a significant feature around the spike, and its velocity is approximately 10 times slower than the rotational one (not shown here).

The magnetic field can be amplified by stretching and compression (Sano et al., 2012). Figure 13(c) shows the profile of the magnetic field strength |\mbox{\boldmathB}| normalized by the upstream field strength. The strong field regions are observed as filamentary structures along the interface, e.g., $y=0.2,0.8$, and the field strength reaches 90 times the upstream value. They correlate well with the rotational flow, indicating that the field is amplified mainly by the stretching. The compression would less contribute to the field amplification because the compressible velocity is much slower than the rotational velocity around the spike. Figure 13(d) shows the time profile of the rotational energy \rho|\mbox{\boldmathu}_{R}|^{2}/2 and the magnetic energy increased by the RMI |\mbox{\boldmathB}|^{2}/2-|\mbox{\boldmathB}_{\rm ref}|^{2}/2 integrated over space, where \mbox{\boldmathB}_{\rm ref} is obtained from the run without the corrugation. Just after the impact of the shock against the contact discontinuity, the rotational energy increases impulsively, followed by the gradual decrease of the rotational energy and increase of the magnetic energy. The magnetic field behaves passively and is continuously amplified as long as its energy density is much smaller than the kinetic energy density. However, the magnetic energy saturates and then begins to decrease after it exceeds the rotational energy at $t=36$, which can be interpreted as the suppression by the amplified magnetic field itself when the Lorentz force becomes strong enough to prevent the flow (Cao et al., 2008).

In order to examine the resolution dependence of the RMI, we conduct simulation runs with different interpolation methods (MUSCL-MC and W4I4D5), different Riemann solvers (HLL and HLLD), and different divergence-free treatments (GLM and CUCT) at different resolutions $(N=32,64,128)$. Figure 14 compares the time profile of the maximum of the magnetic field strength |\mbox{\boldmathB}|_{\rm max} among different runs. Both the growth speed and the saturation level increase with increasing the resolution, as is reported by Sano et al. (2012) for weak initial field cases. The maximum field strength exceeds 100 times the upstream value in the high resolution run. The W4I4D5-HLLD-CUCT scheme exhibits higher saturation levels than the other schemes, indicating better accuracy of the solution. In particular, the HLLD scheme obtains a solution comparable to the HLL scheme with a half grid size in terms of the field amplification. The low saturation level in the HLL scheme is due to the overestimation of the amount of numerical diffusion in the high beta plasma (mentioned in Sec 3.6). The plasma beta value at the maximum field region is $\sim 10^{2}$ in this run. However, the value is decreased to $\sim 10$ and the discrepancy of the saturation level between the HLL and HLLD schemes is reduced in the run with an initial plasma beta of $10^{3}$ (not shown here).

3.8 Magnetic reconnection

Magnetic reconnection is a fundamental energy release process in space and terrestrial environments such as the solar corona and the planetary magnetosphere. Since it intrinsically contains a hierarchical structure, a variety of models have been utilized so as to understand its nature from the fully kinetic scale to the MHD scale (Birn et al., 2001, and collabolation papers). Among them, the resistive MHD is frequently adopted to model the macroscopic dynamics of the reconnection and its associated phenomena far beyond the kinetic scale. Using the ideal (not resistive) MHD equations, we simulate the evolution of an anti-parallel magnetic field configuration that is subject to reconnect numerically, so as to examine the effect of numerical diffusion and the robustness of a scheme similar to Gardiner & Stone (2005) and Lee & Deane (2009). The two-dimensional computational domain with $-L_{x}/2\leq x<L_{x}/2\;(L_{x}=512)$ and $-L_{y}/2\leq y\leq L_{y}/2\;(L_{y}=16)$ is resolved by $N\times N/32$ grid points. The boundary condition is periodic and symmetric in the $x$ and $y$ directions. The initial condition is a Harris current sheet configuration, \rho=\rho_{0},\mbox{\boldmathu}=0,\mbox{\boldmathB}=B_{0}\tanh(y/\lambda)\mbox{\boldmathe}_{x},P=[(1+\beta)B_{0}^{2}-|\mbox{\boldmathB}|^{2}]/2, where $\rho_{0}=B_{0}=1.0$, $\lambda=1.0$ is the current sheet thickness, and $\beta=0.2$ is the plasma beta in the upstream region, respectively. We impose a localized perturbation to the flux function $\delta A_{z}=0.05\lambda\exp[-(x^{2}+y^{2})/(2\lambda)^{2}]$ and a small $(1\%)$ uniform random perturbation to $v$ inside the current sheet. The localized perturbation transfers magnetic flux toward the positive and negative $x$ directions, and then forms a thin current sheet. The resistivity plays a critical role for the reconnection that is imposed numerically in this problem. The simulations are conducted by the HLL-CUCT, the HLLD-CUCT-AA, and the HLLD-CUCT schemes.

Figures 15 through 17 show the time evolution of the out-of-plane current density solved by the three schemes at $N=8192$. After the passage of the initial perturbation, the current sheet gets thin down to the grid scale until $t\sim 100$. Subsequently, the numerical resistivity switches on and the system evolves differently depending on the scheme. The current sheet evolves in the Y-shape structure, and large-scale plasmoids are attached there in the HLL-CUCT scheme (Figure 15). In contrast, the HLLD-CUCT-AA scheme shows the uniform current sheet elongated in the $x$ direction. In the current sheet, multiple plasmoids grow and then merge each other (Figure 16). The HLLD-CUCT scheme also shows the elongated current sheet, but the subsequent plasmoid formation is modest compared to the HLLD-CUCT-AA scheme (Figure 17). The growth of the plasmoid is considerably suppressed with decreasing the resolution to $N=4096$ in the HLLD-CUCT(-AA) schemes (not shown here).

Figure 18 shows the time profile of the spatially-integrated magnetic energy for the three schemes at $N=4096,8192$. The rate of dissipation of the magnetic energy is the fastest in the HLL-CUCT scheme, and the slowest in the HLLD-CUCT scheme. The magnetic energy dissipates faster in the higher resolution run with the HLLD-CUCT-AA scheme, implying that the plasmoid dynamics dominates the dissipation of the current sheet. Since the HLLD Riemann solver exactly captures an isolated and stationary tangential discontinuity without numerical diffusion (Miyoshi & Kusano, 2005), the rate of dissipation remains slow in the low resolution run when the plasmoid is absent and hence the current sheet is approximately static. The aspect ratio of the current sheet (length/thickness) increases with increasing the resolution that would increase the number of plasmoids as is demonstrated in the resistive MHD (Loureiro et al., 2007; Bhattacharjee et al., 2009). The HLLD-CUCT scheme shows the similar tendency as the HLLD-CUCT-AA scheme, although the rate of dissipation is slower owing to the weaker plasmoid activity. In contrast to the HLLD-CUCT(-AA) schemes, the rate of dissipation is faster in the lower resolution run with the HLL-CUCT scheme.

The overall dynamics is quite different between the HLL-CUCT and the HLLD-CUCT(-AA) schemes. The difference is already found in an early stage of the simulation run. Figure 19(a)-(b) compares the $x$-component of the velocity at $t=100$ with the HLL-CUCT and the HLLD-CUCT schemes. The HLLD-CUCT-AA scheme obtains an almost identical result to the HLLD-CUCT scheme until this period, thus it is not shown here. The outflow jet driven by the initial perturbation propagates rightward. The width of the outflow jet is unchanged from the center to the edge of the perturbation $(x=8)$ in the HLLD-CUCT scheme. On the other hand, the outflow jet broadens toward the downstream and forms a crab-claw structure in the HLL-CUCT scheme, and the maximum jet velocity $u=0.7$ is much faster than that in the HLLD-CUCT scheme $u=0.13$. This structure, along with the Y-shaped structure in Figure 15, is a resemblance to the Petscheck-type reconnection triggered by an ad-hoc localized resistivity (Zenitani & Miyoshi, 2011; Zenitani, 2015). This implies that the HLL-CUCT scheme tends to artificially localize the diffusion region under the adopted initial condition. Figure 19(c) compares the plasma beta value averaged over the current sheet $(-1<y<1)$ between the two schemes. The HLL-CUCT scheme shows that the plasma beta gradually increases to $\sim 100$ in the course of the flow from the center to $x=9$. As is discussed in Section 3.6, rapid numerical diffusion of the magnetic field occurs in the high beta region with the HLL-CUCT scheme, and consequently, it increases the gas pressure since we use the conservative scheme. The increased gas pressure helps to open the jet, and in turn, localize the diffusion region. This numerical artifact is not observed in the HLLD-CUCT(-AA) schemes because the HLLD Riemann solver is designed to capture a tangential discontinuity. The plasma beta is nearly uniform in the current sheet.

The different plasmoid evolution between the HLLD-CUCT-AA and the HLLD-CUCT schemes (Figures 16 and 17) may be attributed to the difference in the numerical Lundquist number that is proportional to the length of the current sheet. The current sheet is elongated to $|x|=30$ in the HLLD-CUCT-AA scheme, but it is shortened by $|x|=10$ in the HLLD-CUCT scheme at $t=300$ (middle panel in Figures 16 and 17). This discrepancy must be owing to the difference in the calculation of the electric field. While the CUCT-AA method takes the arithmetic average in Eq. (66), the CUCT method takes the directional biasing with respect to the Alfvén mode in Eq. (67), and its nonlinear coefficient $\Theta$ depends on the position (Eq. (68)). This nonlinear effect may impose a numerical resistivity varying on the $x$ position, and eventually, it may discourage the elongation of the current sheet. The nonlinearity of the numerical resistivity in the CUCT-AA method would be relatively weak, and it results in the elongated current sheet and the subsequent plasmoid formation that is a resemblance to the tearing instability for a uniform resistivity (Biskamp, 1986). In order to test this hypothesis, we apply the schemes to the visco-resistive MHD equations (Minoshima et al., 2016), and solve the same problem. Kinematic viscosity and resistivity coefficients, $\nu$ and $\eta$, are assumed constant and uniform, and the corresponding kinetic and magnetic Reynolds numbers are $V_{A}\lambda/\nu=V_{A}\lambda/\eta=10^{3}$. The coefficients employed here are larger than numerical ones in this problem. Figure 20 confirms the formation of multiple plasmoids and their nonlinear evolution in the visco-resistive MHD simulation with the HLLD-CUCT scheme (and a similar result is obtained with the HLLD-CUCT-AA scheme).

4 Summary

We have designed a conservative finite difference scheme for ideal MHD simulations based on the Weighted Compact Nonlinear Scheme (WCNS; Deng & Zhang, 2000; Zhang et al., 2008), which is easier to increase the order of accuracy in multi-dimension than finite volume schemes. The scheme reconstructs a high-order accurate numerical flux and approximates its spatial derivative by a high-order linear central difference. We have proposed the fourth-order upwind-biased nonlinear interpolation method (Eqs. (40) and (41)) to combine with the fourth-order central difference (Eq. (2.1)), which achieves optimal fifth order of accuracy for linear problems. The scheme can use the physical variables for the high-order interpolation and the calculation of the numerical flux at the midpoint, allowing us to employ a multi-state Riemann solver for MHD. Especially, the scheme employs the HLLD Riemann solver in order to robustly capture rotational, tangential, and contact discontinuities as well as fast mode shocks (Miyoshi & Kusano, 2005). We have proposed a new and simple variant of the Constrained Transport (CT) method, termed as the Central Upwind CT method (Eq. (67)), which is consistent with the base one-dimensional scheme, attains a designed order of accuracy in a multidimensional flow, and preserves the divergence-free condition of the magnetic field as designed high-order finite difference representation. Furthermore, this method implements a two-dimensional directional biasing so as to avoid the numerical oscillation in a specific coordinate axis pointed out by Lee (2013).

The performance of the present scheme has been measured through linear and nonlinear problems. The scheme achieves the designed order of accuracy, and it is more cost effective than lower-order schemes in multidimensional problems (Sec. 3.1-3.2). MHD shocks and discontinuities, smooth flows, and their interactions are resolved without generating unphysical divergence errors (Sec. 3.3-3.5). These problems clearly demonstrate that the CUCT method can improve the robustness of the MHD scheme for strong shocks. The HLLD-CUCT scheme correctly captures the tangential discontinuity in the multidimensional flow irrespective of the plasma beta by virtue of the combination of the multi-state Riemann solver and the divergence-free method (Sec. 3.6). One has to be cautious to use a single-state Riemann solver for high beta plasmas, in which the transit time of the fastest mode may be much faster than a dynamical time scale that seriously overestimates the amount of numerical diffusion. The scheme successfully resolves the magnetic field amplification associated with the Richtmyer-Meshkov instability (Sec. 3.7). The saturation level of the maximum field strength depends on a scheme design as well as a grid resolution: a high saturation level is obtained by a high resolution run and/or a high accurate scheme. The evolution of the current sheet turns to be susceptible to the calculation method of the electric field (Sec. 3.8). The HLLD-CUCT scheme tends to form an elongated current sheet and subsequent multiple plasmoids. The evolution of the plasmoid considerably differs between the arithmetic average (Eq. (66)) and the direction-biased average (Eq. (67)). On the other hand, the HLL-CUCT scheme tends to localize the current density and trigger the magnetic reconnection there.

Jiang & Wu (1999) pointed out that ideal MHD simulations with different Godunov-type schemes may not lead to identical or convergent results in the presence of magnetic field rotation, because the amount of numerical diffusion for the velocity, magnetic field, and temperature vary according to a scheme design. Minoshima et al. (2015) reported that the saturation level of the magnetorotational instability is susceptible to the numerical magnetic Prandtl number, which relies on the choice of the Riemann solver and the treatment for the divergence-free condition. It would be preferable to involve explicit diffusion terms, viscosity, resistivity, and thermal conductivity for practical MHD simulations rather than the ideal approximation. High accurate schemes allow us to use small diffusion coefficients for application to space and astrophysical plasma phenomena.

The CUCT method (Eq. (67)) has upwind-correction terms twice as large as the baseline FluxCT method (Eq. (62)), thus it would give a more robust solution as is demonstrated in Section 3.5. We also test the magnetic reconnection problem in Section 3.8 with the W4I4D5-HLLD-FluxCT method, but numerical oscillation grows to an intolerable level at an early stage and prevents the system to evolve. We then adopt the CUCT method rather than the FluxCT method to combine with the WCNS scheme and the HLLD Riemann solver. In contrast, Tóth (2000) conducted various benchmark tests and concluded that the FluxCT method is superior to other selected CT methods, one of which recovers the proper amount of numerical diffusion against the FluxCT method (Ryu et al., 1998). The difference between this paper and Tóth (2000) might be partly owing to the order of accuracy in space: the reduction of the upwind property in the FluxCT method could be more problematic for higher-order schemes. The CUCT method would be applicable to radiation and relativistic MHD equations because they solve the same Ohm’s law. The combination of the high-order WCNS scheme, the HLLD Riemann solver, and the CUCT method possesses good capabilities in accuracy, robustness, and computation cost. We expect that the present scheme is useful for practical MHD simulations of global space and astrophysical objects such as planetary magnetosphere, heliosphere, and accretion disks. They contain a wide variety of waves, shocks, and magnetic field inhomogeneities that would be faithfully resolved by the high-order, shock-capturing, and divergence-free scheme.

We thank the anonymous referee for carefully reviewing our manuscript and giving insightful comments that greatly improved the manuscript. This research is supported by JSPS KAKENHI Grant Numbers JP15K05369 (T. Minoshima) and JP15K04756 (T. Miyoshi). Numerical simulations were in part carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

Appendix A Sixth and seventh-order finite difference schemes

The sixth-order linear central difference can be written as

[TABLE]

If we apply the seventh-order left-side biased interpolation to $\tilde{U}$,

[TABLE]

the resulting scheme has the sixth-order dispersion as a leading error,

[TABLE]

Similar to Eq. (2.2), we derive the following sixth-order left-side biased interpolation that is appropriate to Eq. (A),

[TABLE]

Substituting Eq. (A) into Eq. (A) cancels the $O(\Delta x^{6})$ error, and the resulting scheme has the seventh-order dissipation as a leading error,

[TABLE]

Comparing Eqs. (A3) and (A5), one expects that the seventh-order scheme improves the dispersion error but the amount of numerical dissipation is increased by a factor of $2048/1792=1.14$ against the sixth-order scheme. This is readily confirmed from the Fourier analysis performed in Section 2.2.

Appendix B Parameter dependence of the rotated shock tube problem with the hyperbolic divergence cleaning method

The hyperbolic divergence cleaning method (Eqs. (88) and (89)) has two free parameters, $c_{h}$ and $c_{p}$, to diminish the \nabla\cdot\mbox{\boldmathB} error through propagation and diffusion. In Figure 21, we examine the $c_{h}$ dependence of \nabla\cdot\mbox{\boldmathB} and density in the Brio-Wu rotated shock tube problem (Sec. 3.4) solved by the W4I4D5-HLLD-GLM scheme. The shock tube is tilted by an angle of $45^{\circ}$ relative to the $x$ axis. The value of $c_{h}$ is set to the corresponding CFL numbers ${\rm CFL}(c_{h})=c_{h}\Delta t/\Delta x=0.4,0.5,0.6,0.7,0.8$ (note that the CFL number for the MHD part is fixed to be 0.4). The value of $c_{p}$ is fixed to be an optimal parameter proposed by Dedner et al. (2002), $c_{p}^{2}=0.18c_{h}$. A \nabla\cdot\mbox{\boldmathB} error of $\sim 2$ is generated at the slow compound wave and the slow shock irrespective of the value of $c_{h}$. Subsequently, the error propagates and diffuses toward the upstream and downstream directions. Grid-scale oscillations seen in ${\rm CFL}(c_{h})\leq 0.5$ is largely reduced for ${\rm CFL}(c_{h})\geq 0.6$. Although the density profile in panel (b) is in good agreement with the W4I4D5-HLLD-CUCT scheme for ${\rm CFL}(c_{h})\leq 0.5$, it gets worse for ${\rm CFL}(c_{h})\geq 0.6$ and this makes it hard to conclude that a solution with higher $c_{h}$ is always better than with lower $c_{h}$. The solution becomes more diffusive around the slow compound wave, the contact discontinuity, and the slow shock for higher $c_{h}$ value. The thickness of the slow shock for ${\rm CFL}(c_{h})=0.8$ is 2.5 times broader toward upstream than for ${\rm CFL}(c_{h})=0.4$. The solution is almost equivalent to one using a half CFL number for the MHD part, and the error is also observed with a tilt angle of $26.6^{\circ}$, indicating the acceptable upper limit of $c_{h}$.