A Six-moment Multi-fluid Plasma Model

TL;DR

This paper introduces a comprehensive six-moment multi-fluid plasma model that captures pressure anisotropy, Hall effects, and multiple charged fluids, providing a balanced approach between complexity and computational efficiency for plasma simulations.

Contribution

The paper develops a new six-moment multi-fluid plasma model that includes pressure anisotropy, Hall effects, and multiple charged fluids, bridging the gap between simpler and more complex models.

Findings

Model accurately captures pressure anisotropy and Hall effects.

Numerical scheme is robust and second-order accurate.

Suitable for coupling with particle-in-cell simulations.

Abstract

We present a six-moment multi-fluid model, which solves the governing equations for both ions and electrons, with pressure anisotropy along and perpendicular to the magnetic field direction, as well as the complete set of Maxwell equations. This set of equations includes the Hall effect, different temperatures for different species and pressure anisotropy. It is more comprehensive than the five-moment equations with isotropic pressures and significantly less expensive than the ten-moment equations with a full pressure tensors. Similarly to the five- and ten-moment equations, the wave speeds are naturally limited by the speed of light, which eliminates the issue of unlimited whistler wave speeds present in Hall magnetohydrodynamics (MHD). It is also possible to simulate multiple negatively charged fluids, which cannot be done in MHD models. The six-moment model is a reasonable description of the plasma outside magnetic reconnection regions and therefore well-suited to be coupled with an embedded particle-in-cell model that covers the reconnection region. Our numerical implementation uses a point-implicit scheme for the stiff source terms, and we use a second-order accurate Rusanov-type scheme with carefully selected wave speeds. For the plasma variables and the magnetic field the maximum wave speed is based on the fast magnetosonic speed of MHD with anisotropic pressures that we derive. For the electric field related variables the speed of light is used. The divergence of the magnetic field and Gauss's law are controlled with a hyperbolic-parabolic scheme. We present a number of numerical tests to demonstrate that this numerical model is robust without being excessively diffusive.

Tables

Table 1. Table 1: Initial conditions for the Brio-Wu Shock test.

	Left	Right
${\rho }_i$	1	0.125
${𝐮}_{𝐢}$	0	0
$p_i$	$5\times {10}^{-5}$	$5\times {10}^{-6}$
${\rho }_e$	1 $\cdot \frac{me}{mi}$	0.125$\cdot \frac{me}{mi}$
${𝐮}_{𝐞}$	0	0
$p_e$	$5\times {10}^{-5}$	$5\times {10}^{-6}$
$B_x$	$0.75\times {10}^{-2}$	$0.75\times {10}^{-2}$
$B_y$	$1\times {10}^{-2}$	$-1\times {10}^{-2}$
$B_z$	0	0
$𝐄$	0	0

Equations

$$\frac{\partial {\rho }_s}{\partial t}+\nabla \cdot ({\rho }_s{\mathbf{u}}_{\mathbf{s}})=0$$

$$\frac{\partial {\rho }_s{\mathbf{u}}_{\mathbf{s}}}{\partial t}+\nabla \cdot \left[{\rho }_s{\mathbf{u}}_{\mathbf{s}}{\mathbf{u}}_{\mathbf{s}}+p_{s\perp }\mathbf{I}+(p_{s\parallel }-p_{s\perp })\mathbf{bb}\right]=\frac{q_s}{m_s}{\rho }_s(\mathbf{E}+{\mathbf{u}}_{\mathbf{s}}\times \mathbf{B})$$

$$\frac{\partial p_{s\parallel }}{\partial t}+\nabla \cdot (p_{s\parallel }{\mathbf{u}}_{\mathbf{s}})=-2p_{s\parallel }\mathbf{b}\cdot (\mathbf{b}\cdot \nabla ){\mathbf{u}}_{\mathbf{s}}$$

$$\frac{\partial p_{s\perp }}{\partial t}+\nabla \cdot (p_{s\perp }{\mathbf{u}}_{\mathbf{s}})=-p_{s\perp }(\nabla \cdot {\mathbf{u}}_{\mathbf{s}})+p_{s\perp }\mathbf{b}\cdot (\mathbf{b}\cdot \nabla ){\mathbf{u}}_{\mathbf{s}}$$

$$\frac{\partial p_s}{\partial t}+\nabla \cdot (p_s{\mathbf{u}}_s)=(p_s-p_{s\parallel })\mathbf{b}\cdot (\mathbf{b}\cdot \nabla )\mathbf{u}-\left(p_s-\frac{p_{s\parallel }}{3}\right)\nabla \cdot {\mathbf{u}}_s$$

$$\frac{\partial e_s}{\partial t}+\nabla \cdot [{\mathbf{u}}_{\mathbf{s}}(e_s+p_s)+{\mathbf{u}}_{\mathbf{s}}\cdot (p_{s\parallel }-p_{s\perp })\mathbf{bb}]=\frac{q_s}{m_s}{\rho }_s{\mathbf{u}}_{\mathbf{s}}\cdot \mathbf{E}$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \times \mathbf{E}$$

$$\frac{\partial \mathbf{E}}{\partial t}-c^2\nabla \times \mathbf{B}$$

$$\nabla \cdot \mathbf{E}$$

$$\nabla \cdot \mathbf{B}$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \times \mathbf{E}+c_B\nabla \psi$$

$$\frac{\partial \mathbf{E}}{\partial t}-c^2\nabla \times \mathbf{B}+c_E\nabla \varphi$$

$$\frac{\partial \psi }{\partial t}+c_B\nabla \cdot \mathbf{B}$$

$$\frac{\partial \varphi }{\partial t}+c_E\nabla \cdot \mathbf{E}$$

$$\mathbf{E}=-{\mathbf{u}}_{\mathbf{e}}\times \mathbf{B}-\frac{1}{e{n}_e}\nabla \cdot [p_{e\perp }\mathbf{I}+(p_{e\parallel }-p_{e\perp })\mathbf{bb}]$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot [\rho \mathbf{uu}+p_{\perp }\mathbf{I}+(p_{\parallel }-p_{\perp })\mathbf{bb}]+\frac{\mathbf{B}}{{\mu }_0}\times (\nabla \times \mathbf{B})=0$$

$$\frac{\partial p_{\parallel }}{\partial t}+\nabla \cdot (p_{\parallel }\mathbf{u})+2p_{\parallel }\mathbf{b}\cdot (\mathbf{b}\cdot \nabla )\mathbf{u}=0$$

$$\frac{\partial p_{\perp }}{\partial t}+\nabla \cdot (p_{\perp }\mathbf{u})+p_{\perp }(\nabla \cdot \mathbf{u})-p_{\perp }\mathbf{b}\cdot (\mathbf{b}\cdot \nabla )\mathbf{u}=0$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{u}\times \mathbf{B})$$

$$\begin{array}{rl}\rho \frac{\partial \mathbf{u}}{\partial t}& +{\gamma }_A^2(\mathbf{I}+\frac{V_A^2}{c^2}\mathbf{bb})\cdot \{\rho (\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p_{\perp }+\nabla \cdot [(p_{\parallel }-p_{\perp })\mathbf{bb}]\}\\ & +\frac{{\gamma }_A^2}{{\mu }_0}\mathbf{B}\times [\nabla \times \mathbf{B}-\frac{1}{c^2}\mathbf{u}\times (\nabla \times \mathbf{E})-\frac{1}{c_0^2}\mathbf{u}\nabla \cdot \mathbf{E}]=0\end{array}$$

$${\gamma }_A=\frac{1}{\sqrt{1+\frac{V_A^2}{c^2}}}$$

$$(\lambda -u_x{)}^2{\mathrm{P}}_2(\lambda ){\mathrm{P}}_4(\lambda )=0$$

$${\mathrm{P}}_2=$$

$$\begin{array}{rl}{\mathrm{P}}_4=& (\lambda -u_x{)}^4-(\frac{2p_{\perp }}{\rho }+\frac{2p_{\parallel }-p_{\perp }}{\rho }{b}_x^2)(\lambda -u_x{)}^2-(c^2-{\lambda }^2)\frac{V_A^2}{c^2}[(\lambda -u_x{)}^2-\frac{3p_{\parallel }}{\rho }{b}_x^2]\\ & -[\frac{p_{\perp }^2}{{\rho }^2}(1-b_x^2)-\frac{3p_{\parallel }{p}_{\perp }}{{\rho }^2}(2-b_x^2)+\frac{3p_{\parallel }^2}{{\rho }^2}{b}_x^2]b_x^2\end{array}$$

$$\begin{split}{\rm P_{4}}=&(\lambda-u_{x})^{4}-(a^{2}+\frac{2p_{\perp}-3p_{\parallel}}{\rho}+\frac{2p_{\parallel}-p_{\perp}}{\rho}b_{x}^{2})(\lambda-u_{x})^{2}-(c^{2}-\lambda^{2})\frac{V_{A}^{2}}{c^{2}}[(\lambda-u_{x})^{2}-a^{2}b_{x}^{2}]\\ &-[\frac{p_{\perp}^{2}-3p_{\perp}p_{\parallel}}{\rho^{2}}(1-b_{x}^{2})+\frac{3p_{\parallel}^{2}}{\rho^{2}}b_{x}^{2}{\color[rgb]{1,0,0}+}\frac{5p_{e}}{3\rho}(\frac{4p_{\parallel}-p_{\perp}}{\rho}b_{x}^{2}-\frac{3p_{\parallel}}{\rho})\color[rgb]{1,0,0}{-\frac{3p_{\parallel}p_{\perp}}{\rho^{2}}}]b_{x}^{2}\end{split}$$

$${\lambda }_{1,2}=u_x$$

$$\begin{array}{rl}{\lambda }_{4,5}=& \frac{1}{2}{\gamma }_A^2[u_x-\frac{V_A^2}{c^2}(\mathbf{u}\cdot \mathbf{b})b_x]+\frac{u_x}{2}\\ & \pm \sqrt{{\gamma }_A^2\left(V_{A,x}^2+\frac{p_{\perp }-p_{\parallel }}{\rho }{b}_x^2\right)+{\left[\frac{1}{2}{\gamma }_A^2(u_x-\frac{V_A^2}{c^2}(\mathbf{u}\cdot \mathbf{b})b_x)+\frac{u_x}{2}\right]}^2}\end{array}$$

$${\lambda }_{3,4}=u_x\pm \sqrt{\frac{B_x^2}{{\mu }_0\rho }+\frac{p_{\perp }-p_{\parallel }}{\rho }{b}_x^2}$$

$$\begin{array}{rl}{\mathrm{P}}_4=& (\lambda -u_x{)}^4-\left(V_A^2+\frac{2p_{\perp }}{\rho }+\frac{2p_{\parallel }-p_{\perp }}{\rho }{b}_x^2\right)(\lambda -u_x{)}^2\\ & -\left[\frac{p_{\perp }^2}{{\rho }^2}(1-b_x^2)-\frac{3p_{\parallel }{p}_{\perp }}{{\rho }^2}(2-b_x^2)+\frac{3p_{\parallel }^2}{{\rho }^2}{b}_x^2-\frac{3p_{\parallel }}{\rho }{V}_A^2\right]b_x^2\end{array}$$

$$\begin{array}{rl}{\lambda }_{5,6,7,8}=u_x& \pm \frac{1}{\sqrt{2\rho }}\{(\frac{B^2}{{\mu }_0}+2p_{\perp }+(2p_{\parallel }-p_{\perp })b_x^2)\pm [(\frac{B^2}{{\mu }_0}+2p_{\perp }+(2p_{\parallel }-p_{\perp })b_x^2{)}^2\\ & +4(p_{\perp }^2{b}_x^2(1-b_x^2)-3p_{\parallel }{p}_{\perp }{b}_x^2(2-b_x^2)+3p_{\parallel }^2{b}_x^4-3p_{\parallel }\frac{B_x^2}{{\mu }_0}){]}^{1/2}{\}}^{1/2}\end{array}$$

$$\begin{array}{rl}{\lambda }_{5,6,7,8}=u_x& \pm \frac{1}{\sqrt{2\rho }}\{(\frac{B^2}{{\mu }_0}+2p_{\perp }+\frac{5}{3}{p}_e+(2p_{\parallel }-p_{\perp })b_x^2)\pm [(\frac{B^2}{{\mu }_0}+2p_{\perp }+\frac{5}{3}{p}_e+(2p_{\parallel }-p_{\perp })b_x^2{)}^2\\ & +4(p_{\perp }^2{b}_x^2(1-b_x^2)-3p_{\parallel }{p}_{\perp }{b}_x^2(2-b_x^2)+3p_{\parallel }^2{b}_x^4+\frac{5}{3}{p}_e(4p_{\parallel }{b}_x^2-p_{\perp }{b}_x^2-3p_{\parallel })b_x^2-3(p_{\parallel }+\frac{5}{3}{p}_e)\frac{B_x^2}{{\mu }_0}){]}^{1/2}{\}}^{1/2}\end{array}$$

$${\tilde{\lambda }}_{5,6}=$$

$${\tilde{\lambda }}_{7,8}=$$

$$\mathbf{U}_{impl}=\left(\begin{array}[]{c}\rho_{s}\mathbf{u_{s}}\\ \mathbf{E}\\ \end{array}\right)$$

$$\mathbf{S}_{impl}(\mathbf{U}_{impl})=\left(\begin{array}[]{c}\frac{q_{s}}{m_{s}}(\rho_{s}\mathbf{E}+\rho_{s}\mathbf{u_{s}}\times\mathbf{B})\\ -c^{2}\mu_{0}\sum\frac{q_{s}}{m_{s}}\rho_{s}\mathbf{u_{s}}\\ \end{array}\right)$$

$$\frac{\partial {\mathbf{U}}_{\mathbf{impl}}}{\partial t}={\mathbf{R}}_{expl}+{\mathbf{S}}_{impl}$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Magnetohydrodynamics (MHD) simulations have been widely carried out to understand the mechanisms behind different phenomena in plasma physics. MHD models assume that the Larmor radius (gyro radius) is much smaller than the characteristic length scale and the particle distribution function can be described by the fluid equations (continuity, momentum and pressure/energy equations). Magnetic ($\mathbf{B}$) and electric ($\mathbf{E}$) fields are needed to solve the governing equations for the ions and electrons. In the MHD approximation, the mass of the electrons is neglected so that we can obtain the approximate electric field from the electron momentum equation. In such an approximation, the magnetic field is frozen into the electron fluid. Further simplifications include ignoring the velocity difference between the electrons and ions in the induction equation and assuming equal temperatures and Maxwellian distributions for both electrons and ions, which leads to the ideal MHD approximation. In ideal MHD, the magnetic field lines are frozen into the plasma consisting of co-moving ions and electrons. The frozen-in condition can be relaxed by taking into account resistivity, however, resistivity is negligible for collisionless plasmas found in space and astrophysics, for example. Another improvement is to include the velocity difference between ions and electrons in the magnetic induction equation. The resulting Hall MHD model includes some of the ion physics. A further step towards a kinetic description is to allow for different electron and ion pressures and allow for pressure anisotropy (for example [1, 2]).

The MHD (including ideal, resistive, Hall and MHD with anisotropic pressure) description neglects electron inertia due to the finite mass of electrons and assumes perfect charge neutrality. The five-moment equations of Shumlak and Loverich [3] remove the assumption of the massless electrons and solve the full set of Maxwell equations to obtain the electric and magnetic fields. In addition, the full set of hydrodynamic equations with the Lorentz force on the right-hand-side are solved separately for the electron and ion fluids, so that the electron mass is taken into account and charge separation is allowed. Further development [4, 5, 6] lead to the ten-moment two-fluid plasma model that solves for full pressure tensors for both the electron and ion fluids. Wang et al. [5] compare their five- and ten-moment two-fluid plasma models with a Particle-In-Cell (PIC) kinetic model and show that their five- and ten-moment models can reproduce many important kinetic features observed in the PIC simulation. Alvarez-Laguna et al. [7] recently proposed a new numerical method, which contains implicit time integration to handle the stiffness of the system and properly scaling of the numerical dissipation from the electromagnetic field solver to the plasma flow solver, to simulate the multi-fluid plasma system.

In this manuscript, we propose another approach, a six-moment multi-fluid plasma model, which is in-between the five- and ten-moment models, by introducing pressure anisotropy for both ions and electrons along and perpendicular to the magnetic field direction. The six-moment approximation requires only two pressure components (parallel and perpendicular) per fluid instead of the six independent components of the full pressure tensors, and the six-moment equations are significantly simpler than the ten-moment equations. This means that the six-moment equations are less expensive to solve. In addition, the six-moment equations are likely to be valid in the vast majority of the plasma system where the electrons and ions are both magnetized, so the off-diagonal terms of the pressure tensor are negligible. Near reconnection regions the six-moment approximation is not valid, but even the ten-moment approximation has a difficult time to reproduce all aspects of kinetic reconnection, although there has been some promising progress [5]. An alternative approach, that we plan to employ in the future, is to use an embedded particle-in-cell (PIC) model [8, 9] to cover the reconnection region. We expect the six-moment equations to provide a good fluid model that can be coupled effectively to the embedded PIC model covering the reconnection site. We expect that using the six-moment model allows reducing the size of the PIC domain compared to the case when the fluid model is simpler (MHD or Hall MHD). An additional feature of the six-moment (also true for the five- and ten-moment) equations is that one can allow for multiple electron fluids or a mixture of electrons and negatively charged ions. These situations cannot be handled with the usual MHD models, as the densities and velocities of the multiple electron fluids cannot be determined from the charge neutrality and electric current.

In the following section we present the six-moment equations, then we derive the characteristic speeds in section 3 that are used in the discretization. The numerical scheme employing a reduced numerical dissipation in combination with a point-implicit discretization of the stiff source terms is discussed in section 4. We present several numerical tests in section 5 to demonstrate the capabilities of our six-moment model, and conclude with section 6.

2 Model Equations

The six-moment equations are an extension of the five-moment equations [3, 10] by introducing pressure anisotropy for both ions and electrons [1, 2]. Under this assumption, the pressure tensor can be approximated with $\mathbf{P}=p_{\perp}\mathbf{I}+(p_{\parallel}-p_{\perp})\mathbf{bb}$, where $\mathbf{I}$ is the identity matrix, $\mathbf{b}$ is the unit vector along the magnetic field direction, $p_{\parallel}$ is the pressure along the parallel direction of the magnetic field and $p_{\perp}$ is the pressure in the perpendicular direction. For monatomic gases, the six-moment equations for all charged fluids (indexed by $s$) can be written as:

[TABLE]

where $\rho$ and $\mathbf{u}$ denote the mass density and the velocity vector, respectively, and $q$ and $m$ are the charges and masses of the particles. For convenience of implementation, we solve the average pressure $p=\frac{2p_{\perp}+p_{\parallel}}{3}$ instead of the perpendicular pressure $p_{\perp}$ because $p$ is already a primitive variable solved by our MHD code BATS-R-US. The equation for $p$ can be obtained from combining equations 1c and 1d:

[TABLE]

Alternatively, we can solve for the hydrodynamic energy density $e=\frac{\rho\mathbf{u}^{2}}{2}+\frac{3}{2}p$ for each species:

[TABLE]

which can be beneficial to get better jump conditions across shock waves. Note, however, that the parallel pressure equation is still solved with the adiabatic assumption, so non-adiabatic heating is not properly captured. In addition, the magnetic energy is not included into the energy density, so the jump conditions are only approximate. In general, there can be many more source terms on the right hand sides of the above equations corresponding to gravity, charge exchange, chemical reactions, collisions, etc.

The electric field ($\mathbf{E}$) and magnetic field ($\mathbf{B}$) are obtained from the Maxwell equations:

[TABLE]

where $\varepsilon_{0}$ is the vacuum permittivity, $\mu_{0}$ is the vacuum permeability, $c=1/\sqrt{\varepsilon_{0}\mu_{0}}$ is the speed of light, $\rho_{c}=\sum_{s}(q_{s}/m_{s})\rho_{s}$ is the total charge density and $\mathbf{j}=\sum_{s}(q_{s}/m_{s})\rho_{s}\mathbf{u}_{s}$ is the current density.

Equations 4c and 4d are constraints on the initial conditions and analytically these conditions are preserved. Numerically, however, this is not guaranteed to hold. We use the hyperbolic/parabolic cleaning method [11, 12, 13] to control the numerical errors in these equations. We introduce the scalars $\psi$ and $\phi$ as additional independent variables and solve the following modified form of the Maxwell equations:

[TABLE]

where $c_{B}$ and $c_{E}$ are the hyperbolic propagation speeds, while $d_{B}$ and $d_{E}$ are the parabolic decay rates. To make the paper more self-contained, we provide a brief derivation in the Appendix to show how the hyperbolic/parabolic cleaning works for the six-moment equations. On the other hand without using a cleaning method, Balsara et al. [14] solved the magnetic and electric fields in plasma on a facially-collocated Yee-type mesh and proved that magnetic field is reconstructed in a divergence-free fashion and the electric field is reconstructed in a form that is consistent with the Gauss’ law. Balsara et al. [15, 16] further extended their Yee-type mesh algorithm to simulate the electrodynamics in material media.

3 Characteristic Wave Speeds

The fastest wave speed in the six-moment (also five- and ten-moment) equations is the speed of light $c$. Using $c$ in the numerical fluxes, however, makes the scheme rather diffusive. To reduce diffusion while maintaining stability, we use a point-implicit evaluation of the stiff source terms following Shumlak et al.[3], who proposed to ignore the Lorentz force terms and consequently the interactions between the charged fluids and the electromagnetic fields while calculating the characteristic wave speeds. Using this approach, the characteristic speed for each fluid will simply be its sound wave speed. We tried this approach, but found that it gives unsatisfactory results in several applications.

We take an alternative approach by considering the wave speeds of MHD with anisotropic electron and ion pressures instead. This takes into account fast magnetosonic waves, which is the proper wave speed in the MHD limit of the six-moment equations. On the other hand the electron sound speed and the whistler wave speed are not included, which reduces the numerical diffusivity, and similarly to Shumlak et al.[3] we rely on the point-implicit scheme to provide numerical stability. In the following discussion, we limit our derivation to a single ion fluid and a single electron fluid, and we use the subscript $i$ to denote the ion fluid while the subscript $e$ is for the electron fluid. The proper generalization to arbitrary number of fluids is left for future work (we currently employ some heuristic formulas that work reasonably in most cases but may not be valid in general).

It is important to note that the following equations in this section are only used to derive the characteristic speeds, but not used in the six-moment model at all. In the MHD approximation, the electric field is obtained from the electron momentum equation by ignoring the electron inertial terms, which gives

[TABLE]

The electron number density $n_{e}$ can be obtained from charge neutrality as $n_{e}=n_{i}q_{i}/e$ (or simply $n_{i}$ for singly charged ions). For the momentum equation, the electron velocity is expressed from the current density as $\mathbf{u_{e}}=\mathbf{u_{i}}-\mathbf{j}/(en_{e})$ resulting in the usual MHD Lorentz force $\mathbf{j}\times\mathbf{B}$ in the ion momentum equation. The current density is obtained from Ampere’s law (after dropping the displacement current) as $\mathbf{j}=\nabla\times\mathbf{B}/\mu_{0}$ as usual in the MHD approximation. In all the other equations we take $\mathbf{u_{e}}=\mathbf{u_{i}}$, so the governing equations for the ions become

[TABLE]

where $p_{\perp}=p_{i\perp}+p_{e\perp}$, $p_{\parallel}=p_{i\parallel}+p_{e\parallel}$, $\rho=\rho_{i}$ and $\mathbf{u}=\mathbf{u_{i}}$.

The magnetic field can be obtained from the classical ideal MHD induction equation ignoring the Hall terms, which can be written as

[TABLE]

In a six-moment simulation, the speed of light is usually reduced to speed up the simulation, in which case the reduced speed of light need to be properly set to make sure that it must be larger than any of the characteristic speeds. In such a system, the characteristic speeds may not be much smaller than the reduced speed of light, so the semi-relativistic situation needs to be considered. In the semi-relativistic case, we only need to modify the momentum equation from the classical limit. The non-conservative form of the momentum equation (Equation 7b) can be written as

[TABLE]

where

[TABLE]

is the Alfvén factor, $c_{0}$ is the true value of the speed of light and $c$ is the artificially reduced speed of light. The term $c_{0}^{-2}\mathbf{u}\nabla\cdot\mathbf{E}$ can be dropped because $u$ is much smaller than $c_{0}$ for the semi-relativistic limit and this term is much smaller than $\nabla\times\mathbf{B}$.

We want to obtain the characteristic wave speed of Equations (7) and (8) with the ion momentum equation replaced by Equation 9). First we write the one dimensional (along the $x$ direction) equations in the form $\frac{\partial\mathbf{U}}{\partial t}+\mathbf{M_{x}}\frac{\partial\mathbf{U}}{\partial x}=0$ where $\mathbf{M_{x}}$ is the characteristic matrix. In the MHD approximation the variable array reduces to $\mathbf{U}=(\rho,\mathbf{u},\mathbf{B},p_{\parallel},p_{\perp})=(\rho,u_{x},u_{y},u_{z},B_{x},B_{y},B_{z},p_{\parallel},p_{\perp})$. In 1D $B_{x}$ is a constant , so the variable array can be further reduced to $\mathbf{U}=(\rho,\mathbf{u},B_{y},B_{z},p_{\parallel},p_{\perp})$. These variables only depend on $x$ and $t$. We further simplify the problem by rotating the coordinate system such that the magnetic field is in the $x-y$ plane so that $B_{z}=0$. The characteristic matrix $\mathbf{M_{x}}$ is obtained by the Mathematica software:

\mathbf{M_{x}}=\left(\begin{array}[]{cccccccc}u_{x}&\rho&0&0&0&0&0&0\\ 0&\gamma_{A}^{2}u_{x}+\chi_{11}&\chi_{12}&\chi_{13}&\kappa_{1}&0&\eta_{12}&\eta_{12}\\ 0&\chi_{21}&\gamma_{A}^{2}{u_{x}}+\chi_{22}&\chi_{23}&\kappa_{2}&0&\eta_{21}&\eta_{22}\\ 0&\chi_{31}&\chi_{32}&\gamma_{A}^{2}{u_{x}}+\chi_{33}&\kappa_{3}&\nu&0&0\\ 0&B_{y}&-{B_{x}}&0&{u_{x}}&0&0&0\\ 0&0&0&-{B_{x}}&0&{u_{x}}&0&0\\ 0&p_{\parallel}(2{b_{x}}^{2}+1)&2{p_{\parallel}}{b_{x}}{b_{y}}&0&0&0&{u_{x}}&0\\ 0&p_{\perp}(2-{b_{x}}^{2})&-{p_{\perp}}{b_{x}}{b_{y}}&0&0&0&0&{u_{x}}\\ \end{array}\right)

where

\chi=\frac{\gamma_{A}^{2}}{\mu_{0}\rho c^{2}}\left(\begin{array}[]{ccc}(B_{x}^{2}-B_{y}^{2})u_{x}&2B_{x}B_{y}u_{x}&0\\ 2B_{x}B_{y}u_{x}&(B_{y}^{2}-B_{x}^{2})u_{x}&0\\ -B_{y}^{2}u_{z}&B_{x}B_{y}u_{z}&-B_{x}^{2}u_{x}-B_{x}B_{y}u_{y}\\ \end{array}\right)

\kappa=\frac{\gamma_{A}^{2}}{\mu_{0}\rho c^{2}}\left(\begin{array}[]{c}(c^{2}-u_{x}^{2})B_{y}+(2\mu_{0}\rho c^{2}B^{-2}+1)b_{x}^{2}B_{y}(p_{\perp}-p_{\parallel})\rho^{-1}\\ (u_{x}^{2}-c^{2})B_{x}+(b_{y}^{2}+(b_{y}^{2}-b_{x}^{2})\mu_{0}\rho c^{2}B^{-2})B_{x}(p_{\perp}-p_{\parallel})\rho^{-1}\\ -B_{y}u_{x}u_{z}\\ \end{array}\right)

$\nu=\frac{\gamma_{A}^{2}}{\mu_{0}\rho c^{2}}[(u_{x}^{2}-c^{2})B_{x}+B_{y}u_{x}u_{y}-\mu_{0}c^{2}B_{x}B^{-2}(p_{\perp}-p_{\parallel})]$

\eta=\frac{1}{\rho}\left(\begin{array}[]{cc}b_{x}^{2}&\gamma_{A}^{2}b_{y}^{2}\\ b_{x}b_{y}&-\gamma_{A}^{2}b_{x}b_{y}\\ \end{array}\right)

The matrix $\mathbf{M_{x}}$ is identical to the submatrix (the upper left $8\times 8$ elements) of the characteristic matrix that Meng et al. [2] (hereafter Paper I) obtained with pressure anisotropy in ions and isotropy in electrons. We note that here the parallel and perpendicular pressures are the sums of the ion and electron pressures. We also correct a typo in the second element of the matrix $\kappa$ in Paper I, where $B_{x}^{2}$ should be $B_{x}$. We use Mathematica to solve the characteristic equation $\det(\mathbf{M_{x}}-\lambda\mathbf{I})=0$ and after some tedious algebra, the characteristic equation can be written as

[TABLE]

where the wave speed $\lambda$ is one of the eigenvalues of $\mathbf{M_{x}}$ and $\rm P_{2}$ and $\rm P_{4}$ are second- and fourth-order polynomials, respectively:

[TABLE]

where $V_{A}^{2}=B^{2}/(\mu_{0}\rho)$ is the square of the classical Alfvén speed.

The $\rm P_{2}$ and $\rm P_{4}$ polynomials are identical to the $\rm P_{2}$ and $\rm P_{4}$ expressions in Paper I after substituing $p_{e}=0$ for the isotropic electron pressure in Paper I. There are, however, two typos in $\rm P_{4}$ in Paper I. The correct expression should be (with the corrections highlighted in red):

[TABLE]

We note that these typos in [2] are only in the published paper, and the equations used in Maple to derive the wave speeds and the wave speeds implemented into the code are all correct. There are eight eigenvalues for the characteristic equation and each one is associated with one characteristic wave. Two of the eigenvalues are straightforward:

[TABLE]

which are the two entropy waves related to $p_{\perp}$ and $p_{\parallel}$.

3.0.1 Alfvén wave

As the $\rm P_{2}$ polynomial is the same as in Paper I, the roots corresponding to the Alfvén wave speeds are the same too:

[TABLE]

where $V_{A,x}=V_{A}b_{x}=\sqrt{B_{x}^{2}/(\mu_{0}\rho)}$ is the classical Alfvén wave speed in the $x$ direction. It is important to point out that even though this formula looks the same as the solution in Paper I, the physical meaning is not the same because in our case, $p_{\perp}$ and $p_{\parallel}$ are the sum of parallel and perpendicular pressures of ions and electrons, which means that in our case, the electron pressure does contribute to the Alfvén wave speed.

In the classical limit ($V_{A}\ll c$ and $\gamma_{A}\rightarrow 1$), the solutions reduce to

[TABLE]

3.0.2 Fast and slow magnetosonic waves

The exact solutions of $\rm P4$ are too complicated to obtain. We follow the approach suggested in Paper I to obtain the approximate fast and slow magnetosonic wave speeds. We first obtain the solutions in the classical limit, in which case $\rm P4$ simplifies to

[TABLE]

The solutions can be easily obtained as:

[TABLE]

The solutions look the same as the formula obtained in Paper I (when neglecting the electron pressure) and Baranov et al. (1970) [17]. We would like to correct another typo in $\lambda_{5,6,7,8}$ in Paper I. The first term in the second line of the expression should be $p_{\perp}^{2}b_{x}^{2}(1-b_{x}^{2})$. The complete correct expression is

[TABLE]

The next step is to extend the above solutions to the semi-relativistic case by considering some special cases (for example, $\mathbf{u}=0$ and $b_{x}=1$). The steps are the same as in Paper I (after setting the isotropic electron pressure to zero and adding the anisotropic electron pressure to the total pressure) so we do not repeat the procedure here. The final approximate formulas for the fast and slow wave speeds can be written as

[TABLE]

where $\overline{a}^{2}=a^{2}(1+\frac{V_{A,x}^{2}}{c^{2}})+\frac{2p_{\perp}-3p_{\parallel}}{\rho}+\frac{2p_{\parallel}-p_{\perp}}{\rho}b_{x}^{2}$, $b^{2}=\frac{b_{x}^{2}}{\rho^{2}}[3p_{\parallel}p_{\perp}(2-b_{x}^{2})-p_{\perp}^{2}(1-b_{x}^{2})-3p_{\parallel}^{2}b_{x}^{2}]$, $\overline{V}_{A}^{2}=V_{A}^{2}(1-\gamma_{A}^{2}\frac{u_{x}^{2}}{c^{2}})$ and $\overline{V}_{A,x}^{2}=V_{A,x}^{2}(1-\gamma_{A}^{2}\frac{u_{x}^{2}}{c^{2}})$. Paper I showed numerically that these approximate speeds are accurate in most of the practically important parameter regime.

4 Discretization

In the following subsections, we describe how we discretize the fluxes (the pure divergence terms on the left hand side) and the source terms on the right hand side of the six-moment equations.

The time step is limited by the Courant-Friedrichs-Lewy (CFL) condition based on the speed of light. In practice, we can reduce the speed of light to a value that is a factor of 2-3 faster than the fastest flow and fast wave speed obtained in the previous section to speed up the simulation.

4.1 Source terms

The stiff source terms are evaluated by a new point-implicit scheme. Only the momenta and the electric field are involved, so the implicit variables are

[TABLE]

In the momentum equations and the Maxwell equation for the electric field we split the various terms into two groups: the fluxes and non-stiff source terms $\mathbf{R}_{expl}$ and the stiff source terms $\mathbf{S}_{impl}$ containing the Lorentz force terms in the momentum equations and the $c^{2}\mu_{0}\mathbf{j}$ term in the Maxwell equations. The stiff source terms can be written as:

[TABLE]

which shows that $\mathbf{S}_{impl}$ is linear in $\mathbf{U}_{impl}$. The implicit variables $\rho_{s}\mathbf{u_{s}}$ (an independent variable and could be denoted as $\mathbf{m_{s}}$) and $\mathbf{E}$ are multiplied with explicit variables $\rho_{s}$ and $\mathbf{B}$. The point-implicit update is respective to the time level $n$, which is shown in Equation (24). With this notation the six-moment equations for the $\mathbf{U}_{impl}$ variables can be written as

[TABLE]

We use the following steps to update the point implicit variables (the rest of the variables are updated with a simple explicit step):

[TABLE]

where $0.5\leq\beta\leq 1$ is the time centering parameter and $\mathbf{R}_{expl}=\mathbf{R}_{expl}(\mathbf{U}_{impl})$ is restricted to the point implicit variables $\mathbf{U}_{impl}$. The three steps can be combined into a single update:

[TABLE]

where $I$ is the identity matrix and the matrix $({I}/{\Delta t}-\beta{\partial\mathbf{S}_{impl}}/{\partial\mathbf{U}_{impl}})$ is obtained analytically and inverted numerically. We note that $\mathbf{S}_{impl}$ is linear so it is very easy to calculate the partial derivatives. For example, for one ion and one electron fluids

[TABLE]

where $r_{i}=-c^{2}\mu_{0}{q_{i}}/{m_{i}}$, $r_{e}=-c^{2}\mu_{0}{q_{e}}/{m_{e}}$, $w_{i}=\rho_{i}{q_{i}}/{m_{i}}$ and $w_{e}=\rho_{e}{q_{e}}/{m_{e}}$.

This particular discretization of the point-implicit scheme has a very important property: it preserves steady state independent of the time step. If $\mathbf{R}_{expl}^{n}+\mathbf{S}_{impl}^{n}=0$, then $\mathbf{U}^{n+1}_{impl}=\mathbf{U}^{n}_{impl}$ independent of the value of $\Delta t$. This property is crucial when the source terms are very stiff, as is the case here. Alternative forms of the point-implicit scheme that do not have this property can produce incorrect solutions.

For a second-order in time scheme, we use the point-implicit update in both the predictor and corrector steps. In the predictor step the time step is $\Delta t/2$ and we set $\beta=1$, while in the corrector step we use the full time step $\Delta t$, the time centered value for $\mathbf{R}_{expl}^{n+1/2}$ and $\beta=1/2$ to get second order accuracy for the point-implicit term:

[TABLE]

where $\alpha_{1}=0$ or $1/2$ and $\alpha_{2}=1/2\ \text{or}\ 1$ depending on the time levels of the explicit variables being used. The $\alpha_{2}=1/2$ option will achieve second-order accuracy, while $\alpha_{2}=1$ is more robust but not perfectly second order accurate. Our current implementation uses $\alpha_{1}=1/2$ and $\alpha_{2}=1$, i.e. the already updated explicit variables.

This two-stage scheme also has the steady state conserving property. In steady state, the explicit update does not change the explicit variables, so $\mathbf{U}_{expl}^{n+\alpha_{1}}=\mathbf{U}^{n}_{expl}$. For the implicit variables, if $\mathbf{R}_{expl}^{n}+\mathbf{S}_{impl}(\mathbf{U}^{n}_{impl},\mathbf{U}^{n+\alpha_{1}}_{expl})=\mathbf{R}_{expl}^{n}+\mathbf{S}_{impl}^{n}=0$, then $\mathbf{U}_{impl}^{n+1/2}=\mathbf{U}^{n}_{impl}$ in the first stage. In the second stage the explicit update does not change the explicit variables, so $\mathbf{U}_{expl}^{n+\alpha_{2}}=\mathbf{U}^{n}_{expl}$ and $\mathbf{R}_{expl}^{n+1/2}+\mathbf{S}_{impl}(\mathbf{U}^{n}_{impl},\mathbf{U}^{n+\alpha_{2}}_{expl})=\mathbf{R}_{expl}^{n}+\mathbf{S}_{impl}^{n}=0$, and consequently $\mathbf{U}^{n+1}_{impl}=\mathbf{U}^{n}_{impl}$.

It is worth to mention that Balsara et al. [14] applied multiple stages with their Runge-Kutta implicit-explicit (IMEX) methods and could achieve a more accurate implicit-source discretization than our method. Abgrall and Kumar [18] used a very similar implicit source treatment with $\frac{\partial\mathbf{S}_{impl}}{\partial\mathbf{U}_{impl}}$ taken from the time level $n+1$ and showed that their point-implicit treatment could unconditionally preserve positivity.

4.2 Physical fluxes

The left hand sides of Equations (1 - 3) and (5) contain pure divergence terms and are obtained by the local Lax-Friedrichs or Rusanov scheme [19]:

[TABLE]

where $U$ is one of the conservative variables, $F_{i+1/2}=\frac{F_{i+1/2}(U_{R})+F_{i+1/2}(U_{L})}{2}$ is the corresponding flux and $a$ is the fastest wave speed. The superscripts $n$ and $n+1$ refer to the time levels, while the subscript $i$ correspond to the spatial grid indexes. The subscript $i+1/2$ is the cell face between the cell centers $i$ and $i+1$, while the $R$ and $L$ subscripts correspond to the right and left extrapolated face values using some standard TVD type slope limiter. The second term on the right hand side is the numerical diffusion providing stability. In the standard Rusanov scheme the $a$ coefficients are set to the fastest wave speed (for the local state variables), which in principle should be the speed of light $c$. However, the use of speed of light to numerically diffuse all the variables greatly increases the numerical diffusion, and it requires a very fine grid and high computational cost to obtain an accurate solution.

As an alternate approach, we find that it is sufficient to use the fast magnetosonic speeds from Equation 20 and set $a_{i+1/2}=\max(|\lambda_{7}|,|\lambda_{8}|)$ for the plasma quantities $\rho_{s}$, $\mathbf{u}_{s}$, $p_{s}$ and $p_{\parallel,s}$, the magnetic field $\mathbf{B}$ and its hyperbolic cleaning variable $\psi$, while the electric field $\mathbf{E}$ and its hyperbolic cleaning variable $\phi$ need to be numerically diffused by the speed of light using $a_{i+1/2}=c$. These variable dependent wave speeds can substantially reduce the numerical diffusion for the plasma quantities and the magnetic field and improve the accuracy of the solution substantially. However, the numerical diffusion algorithm proposed here may not work for super thermal electrons, while their sound speeds can be very large compared to cold or warm electrons. In such a case, we provide an option to numerically diffuse a sub set of the variables with the speed of light, for example the electron density, momentum and pressure may be diffused with the speed of light, if necessary.

4.3 Relaxation towards isotropy

Paper I discussed three kinds of instabilities (fire hose, mirror and proton cyclotron instabilities), which will push the pressure tensor towards isotropy in the context of a single anisotropic ion fluid. As our model contain multiple fluids, the stability criteria become much more complicated. For sake of simplicity, we implemented a simple exponential decay term for each fluid (a right-hand-side source term in Equation (1c), similar to that suggested in Paper I:

[TABLE]

where $\tau_{s}$ is the relaxation time which relaxes the $p_{\parallel}$ towards $p$. In the extreme case when $\tau$ is extremely small the anisotropy will relax to isotropy immediately, in which case the six-moment simulation becomes a five-moment simulation.

The source term is applied in the same way numerically as discussed in Paper I, in a split manner at the end of the time step and discretized point-implicitly for the sake of numerical stability:

[TABLE]

where $\Delta t$ is the stable time step, $*$ and $n+1$ are the incomplete and final time levels.

5 Numerical Tests

We perform a number of numerical tests to verify the robustness of the six-moment solver, including the light wave test to check the magnetic and electric field solver, the fast wave test to check the propagation of the fast magnetosonic wave. We also perform the GEM reconnection challenge [20] to test the applicability of the six-moment model to simulate the magnetic reconnection. All the tests are performed in normalized units.

5.1 Light Wave

The light wave test is set up on a 1-D grid between $x=\pm 5\times 10^{-4}$ with periodic boundary conditions. The ion and electron masses are $m_{i}=1$ and $m_{e}=0.01$, respectively. The initially uniform fluid states are set to $\rho_{i}=1$, $\mathbf{u_{i}}=0$, $p_{i,\perp}=p_{i,\parallel}=5\times 10^{-6}$, $\rho_{e}=0.01$, $\mathbf{u_{e}}=0$, $p_{e,\perp}=p_{e,\parallel}=5\times 10^{-6}$. We set the speed of light to $c=10$ and the magnetic and electric fields are perturbed with sinusoidal waves as $\delta B_{y}=-0.01\cos(2000\pi x)$ and $\delta E_{z}=0.01\cos(2000\pi x)$. Figure 1 shows the simulation results with 800 grid cells after the wave propagating one full period showing the expected solution.

We did a grid convergence study with $n_{x}=$ 100, 200, 400 and 800 grid cells using the 2nd order Rusanov scheme. The errors are calculated as the $L_{1}$ norm of the difference of the solution after 1 period relative to the initial condition. Figure 2 shows the grid convergence rate for $B_{y}$, which is very close to the 2nd order convergence rate, as expected.

5.2 Brio-Wu Shock

We carry out the Brio-Wu shock [21] test with isotropic pressure for both ions and electrons (five-moment) in this session. The five-moment equations are chosen so the solution can be compared with published results [10], and also because the exact solution of the Brio-Wu shock for the six-moment system is unknown and depends on the pressure anisotropy behind the shock, which is not determined by conservation laws. At this time we have no physics based relaxation of the pressure anisotropy implemented for the six-moment equations.

The test is set up on a 1-D grid between $x=\pm 0.5$ with $10^{4}$ cells and open boundary conditions. As suggested by Hakim et al. [10], the ion inertial length plays an important role when the electron fluid is taken into account. So we set the ion mass $m_{i}$ to 1, 0.1, and 0.001, respectively, while the ion mass to electron mass ratio is fixed (${m_{i}}/{m_{e}}=1836$). The initial conditions are listed in Table 1. Figure 3 shows the simulation results at $t=10$. As expected, these results are very close to the results published in Hakim et al. [10], except that there is a spike in the ion density at about $x=0.1$ for $m_{i}=1$ and $m_{i}=0.1$, and the oscillations do not start next to the shock at about $x=0.05$ for $m_{i}=0.001$. The small differences may come from the different schemes and/or the point-implicit evaluations.

5.3 Firehose instability

We perform a test of the firehose instability on a 1-D grid between $x=\pm 6$ with $10^{4}$ cells and periodic boundary conditions. We apply similar parameters as suggested in Paper I, which is $\rho_{i}=1$, $\mathbf{u_{i}}=\mathbf{u_{e}}=0$, $B_{x}=10$, $B_{y}=B_{z}=0$, $\mathbf{E}=0$, $p_{i,\parallel}=p_{e,\parallel}=52$, $p_{i,\perp}=p_{e,\perp}=55/3$. Due to the relatively small characteristic length, we set $m_{i}=0.001$, in which case the ion inertial length is much smaller than the characteristic length so that the result is close to the classical MHD limit. The ion mass to electron mass ratio ${m_{i}}/{m_{e}}$ is set to 1000 and $n_{e}=n_{i}$ due to quasi-neutrality. We impose small perturbations on the background as $\delta u_{i,y}=\delta u_{e,y}=0.01\cos(k_{A}x)=0.01\cos(2\pi x/6)$, $\delta B_{y}=0.1\cos(2\pi x/6+\pi/2)$ and the electric field is perturbed with the relation $\mathbf{E}=-\mathbf{u_{e}}\times\mathbf{B}$. Because the Alfvén speed ($v_{A}$) is $\sqrt{-1}$, which is obtained from Equation (16), so the perturbations will not propagate but start to grow exponentially with $\exp(|{v_{A}}|k_{A}t)$.

Figure 4 plots $u_{i,y}$ and $B_{y}$ at three different times, which shows the growth patten of the firehose instability. There are fast growing oscillations near the local extrema and sign change, which are short wave length perturbations caused by the numerical errors. In the six-moment model, ion and electron kinetics as well as point-implicit source terms are involved, which make the six-moment model more complicated than the ideal anisotropic MHD. So it is not unexpected that the short wave length perturbations appear much sooner than in the ideal anisotropic MHD test presented in Paper I. Figure 5 shows the agreement of the simulated growth rate of the average kinetic energy $E_{k,y}=\rho_{i}u_{i,y}^{2}/2$ and the analytical growth rate ($2|{v_{A}}|k_{A}$, as the perturbations grow exponentially with $\exp(|{v_{A}}|k_{A}t)$) is good until $t=0.085$. After $t=0.085$ the short wavelength perturbation becomes significant and the growth rate of the kinetic energy deviates from the theoretical expectation.

5.4 Fast Magnetosonic Wave

Daldorff et al. [8] used the initial conditions

[TABLE]

to set up a fast magnetosonic wave for the MHD equations with anisotropic ion pressure . Here $c_{f}=\frac{\omega}{k}=\sqrt{\frac{B_{0}^{2}+2p_{0}}{\rho_{i,0}}}$ is the fast wave propagation speed moving perpendicular relative to the magnetic field direction and $\gamma=\frac{5}{3}$ is the adiabatic index. This is an exact solution for infinitesimal perturbation $\delta$.

We follow Daldorff et al. [8] to create a test for the six-moment model. When deriving the characteristic wave speeds in the classical MHD limit, we showed that the ion and electron pressures can be combined when pressure anisotropy exists in both ions and electrons. This means that we can adopt the above initial conditions by simply splitting the pressure evenly between ions and electrons.

[TABLE]

The simulation domain is on a 2-D grid bounded between $-80/3<x<80/3$ and $-20<y<20$. A single full wave with a rotation $\phi=\tan^{-1}(4/3)$ relative to the $x$ axis is used, which means that the wavelength is $\lambda=32$, $k=2\pi/\lambda\approx 0.1964$ and $T\ \text{(period)}=\lambda/c_{f}=640$ . We set $n_{0}=1$, $m_{i}=1$, $m_{e}=0.01$, $p_{0}=4.5\times 10^{-4}$, $B_{0}=0.04$, which gives the fast magnetosonic speed as $c_{f}=0.05$. The perturbation amplitude $\delta$ is set to 0.1, which is moderately non-linear, and the speed of light is $c=1$ (much larger than $c_{f}$ so that the system is in the classical limit). Figure 6 and 7 show the propagation of the simulated fast magnetosonic wave on the 2-D grid with 512 cells along the $x$ direction and 384 cells along the $y$ direction. The simulation results are consistent with Daldorff et al. [8].

We carry out a grid convergence study for the fast wave test on the 2-D grid with $n_{x}=$ 64, 128, 256, 512 and 1024 cells along the $x$ direction and $n_{y}=$48, 96, 192, 384 and 768 cells along the $y$ direction, respectively, using the 2nd order Rusanov scheme. As we don’t have an analytical solution for this non-linear fast wave, we use the simulation on the grid with $n_{x}=1024$ and $n_{y}=768$ as a reference solution and calculate the relative error of $E_{z}$. Figure 8 shows the grid convergence rate, which is very close to the 2nd order, as expected. We note that this test exercises the complete set of equations and their discretization, including the point-implicit scheme.

5.5 GEM Reconnection

The Geospace Environmental Modeling (GEM) reconnection challenge [20] has been widely applied to test a physics model’s capability to simulate the reconnection process. Many models, including resistive MHD [20], Hall MHD [22, 20, 23], hybrid models [20] and Particl-In-Cell (PIC) models [24, 25], have been compared. The general conclusion of [20] was that including the Hall physics is the minimum requirement to correctly capture the fast magnetic reconnection rate.

The five-, six- and ten-moment multi-fluid models all include the Hall effect automatically by solving the exact Maxwell equations and allowing for different electron and ion velocities. Hakim et al. [10] carried out the GEM reconnection challenge and obtained some complex flow features in the electron fluid. Wang et al. [5] compared their five-moment and ten-moment two-fluid plasma model with a Particle-In-Cell (PIC) model and showed that their five- and ten-moment models can reasonably reproduce some of the important electron kinetic features observed in the PIC simulation during magnetic reconnections. In this section, we perform six-moment simulations for the GEM reconnection challenge.

The classical GEM reconnection challenge is based on the Harris current sheet equilibrium model, and such equilibrium only occurs when the inertial terms from the electrons are neglected. The classical Harris current sheet is not in equilibrium state for multi-fluid with electron fluid and particle-in-cell models. Eventhough we can obtain similar results as Hakim et al. with our five-moment model, their initial conditions are obtained from the classical MHD limit and it is not applicable to the multi-fluid model because their initial conditions without perturbations are not in equilibrium state.

We obtain the initial conditions starting from the oppositely directed magnetic fields

[TABLE]

where $B_{0}$ is the background magnetic field and $\lambda$ is the width of current sheet. The current density is then given by

[TABLE]

Multiple equilibrium states can exist with different plasma conditions. We choose a uniform ion fluid background ($n_{i}$ and $p_{i}$ are uniform) with $u_{i}=0$. There is no charge separation initially, so $n_{e}=n_{i}$ and $\mathbf{E}=0$. The ion fluid is in equilibrium in the unperturbed system. The current is carried by the electrons with the velocity

[TABLE]

Substituting $n_{e}$, $\mathbf{u_{e}}$, $\mathbf{E}$ and $\mathbf{B}$ into the electron momentum equation, the electron fluid is in equilibrium when $p_{e,\parallel}=p_{e,\perp}=p_{e,0}+0.5\cdot(B_{0}^{2}-B_{x}^{2})$, where $p_{e,0}$ is the background electron pressure.

We set the background plasma parameters as $B_{0}=0.07$, $\lambda=0.5$, $m_{i}=1$, $m_{e}=1/25$, $n_{i}=n_{e}=1.225$, $\mathbf{u_{i}}=0$ $p_{i,\perp}=p_{i,\parallel}=2.45\times 10^{-3}$ and $p_{e,0}=4.9\times 10^{-3}$. The speed of light is set to $c=1$. The simulation domain is a 2-D grid between $\pm L_{x}=\pm 25.6$ in the $x$ direction with 512 cells and $\pm L_{y}=\pm 12.8$ in the $y$ direction with 256 cells giving the grid resolution $\Delta x=\Delta y=0.05$. Both the ion inertial length $d_{i}=\frac{1}{e}\sqrt{\frac{m_{i}}{\mu_{0}n_{i}}}\approx 0.90$ and the electron skin depth $d_{e}=\frac{1}{e}\sqrt{\frac{m_{e}}{\mu_{0}n_{e}}}\approx 0.18$ are reasonably well resolved. A periodic boundary condition is applied in the $x$ direction, while a reflecting boundary condition is applied in the $y$ direction. We apply the same form of perturbation to the magnetic field $\delta\mathbf{B}=\mathbf{e_{z}}\times\nabla\chi$ as suggested by Birn et al. [20], which is $\chi=\chi_{0}\cos(2\pi x/L_{x})\cos(\pi x/L_{y})$ and $\chi_{0}=0.1B_{0}$.

In order to understand how pressure anisotropy affects the magnetic reconnection, we investigate four different scenarios: 1. the relaxation time $\tau=1\times 10^{-5}$ for both ions and electrons, which is an approximate five-moment simulation; 2. $\tau=1\times 10^{-5}$ for ions and no relaxation constrains for electrons, which is isotropic ions and anisotropic electrons; 3. $\tau=1\times 10^{-5}$ for electrons and no relaxation constrains for ions, which is anisotropic ions and isotropic electrons; and 4. no relaxation for either ions or electrons, which is a six-moment simulation.

Figures 9, 11 and 10 show the magnetic field, ion velocity and electron velocity magnitudes at $t=495.01$ for the four cases, respectively. The approximate five-moment simulation shows the same features as previous five-moment simulations performed in the literature [10, 4, 5]. The most surprising finding is that a pure six-moment simulation provides a completely different solution than the five-moment simulation, despite the fact that the pressure anisotropy remains within 10% of being isotropic for both ions and electrons during the whole simulation. Wang et al. [5] showed that with the full pressure tensor, their ten-moment simulation can represent the magnetic reconnection process reasonably well compared to a PIC simulation. This suggests that both the isotropic (five-moment) and full pressure tensor (ten-moment) equations give fast reconnection, but the anisotropic pressure (six-moment) does not. From these four cases, we also find that the isotropic electrons case is closer to the approximate five-moment case, which implies that the ion pressure anisotropy has more impact than the electron anisotropy.

The reconnected magnetic flux $F=\int_{0}^{L_{x}/2}|B_{y}|dx$ is typically used to measure how fast the magnetic reconnection occurs. In the approximate five-moment simulation, a magnetic island exists in the middle of the domain, so simply using this formula would provide a larger magnetic flux than the simulations that do not have a magnetic island. The reconnected magnetic flux formula is modified to $F=\int_{0}^{L_{x}/2}\max(0,B_{y})dx$, where only $B_{y}>0$ is taken into account, so that the magnetic island does not contribute to the reconnected flux.

Figure 12 shows the reconnected flux for all cases. We find that the original reconnected flux formula provides a faster reconnection rate for the approximate five-moment simulation than the PIC simulation while the revised formula suggests a different result, which can be explained by the fact that the magnetic island in the five-moment simulation contributes to the reconnected flux in the original formula while it would not in the revised formula and no magnetic islands exist in the PIC simulation. We arrive at the conclusion that the revised formula provides a better description of the reconnected flux if magnetic islands exist. Figure 12 also shows that the pure six-moment simulation provides a very low reconnection rate, as expected from Figure 9 - 10. Overall, the approximate five-moment simulation has the closest reconnection rate to the PIC simulation. The isotropic electrons with anisotropic ions case has a faster reconnection rate than the anisotropic electrons with isotropic electrons case, but they are both slower than the approximate five-moment simulation. We also checked that if we start from the approximate five-moment simulation, then turn off the relaxation towards isotropy for both ions and electrons, then the already ongoing reconnection gets suppressed due to the developing pressure anisotropies. This suggests that the six-moment equations produce slow reconnection even if the simulation is started from a fast reconnection scenario.

6 Conclusions

In this manuscript, we have developed a new model, the six-moment multi-fluid plasma model, to simulate both ions and electrons with pressure anisotropy in a system when they can be described by the fluid equations. The new six-moment model solves for the full set of the electron continuity, momentum and pressure equations, as well as the exact Maxwell equations. The six-moment model can simulate the light wave, Langmuir wave and MHD waves accurately if the grid is fine enough to resolve the corresponding wave length.

We use a steady-state conserving point-implicit time discretization for the stiff source terms associated with the Lorentz force terms in the momentum equations and the $c^{2}\mu_{0}\mathbf{j}$ term in the Maxwell equations. The point-implicit time integration is combined with a spatial discretization based on a Lax-Friedrichs type numerical flux employing the fast magnetosonic speed for the plasma and magnetic field variables and the light speed to the electric field related variables. Using the magnetosonic speed instead of the light speed for the majority of the variables greatly reduces the numerical dissipation. Our numerical tests show that this approach is sufficient to maintain stability in most circumstances. Our implementation provides an option to use the speed of light for the Lax-Friedrichs flux for an arbitrary subset of variables if needed.

A surprising discovery is that the six-moment model cannot provide a good description of the magnetic reconnection process. However, our goal is not to use the six-moment model to study the reconnection process. We plan to couple the six-moment model with an embedded PIC code [26, 9], similar as Daldorff et al. [8]. In such an approach, the magnetic reconnection region will be simulated with the PIC code and other regions will be simulated with the six-moment model. We expect that the six-moment equations are able to describe the vast majority of the plasma system outside magnetic reconnection regions where off-diagonal elements of the pressure tensor are negligible and with such an approach, the size of the PIC domain can be reduced.

An additional feature of the six-moment (also five- and ten-moment) equations is that it allows the use of multiple electron populations. The densities and velocities of multiple electron fluids cannot be approximated from charge neutrality and the current density, so an MHD approximation is not possible. This means that the six-moment equations can be applied to plasmas with thermal and super-thermal populations, or counter streaming populations, or populations of different origins (for example solar wind and ionospheric).

Finally, the six-moment equations provide a reasonable description of typical collisionless plasma conditions where the random motions along the field lines and the gyration perpendicular to the field naturally result in independent parallel and perpendicular pressures, but the off-diagonal terms are usually negligible in the vast majority of the simulation domain. The six-moment equations are only moderately more complicated than the 5-moment equations, but much simpler than the ten-moment equations, which results in lower computational cost and simpler implementation in comparison with the latter.

Acknowledgements

This work was supported by the INSPIRE NSF grant PHY-1513379, the NSF PREEVENTS grant 1663800, the NSF strategic capability grant AGS-1322543, the NASA grant NNX14AE75G, the US Rosetta Project with the JPL subcontract 1266313 under the Rosetta NASA grant NMO710889.

The authors would like to acknowledge the following high-performance computing resources: the Blue Waters super computer by the NSF PRAC grant ACI-1640510, the Pleiades computer by NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, and Yellowstone (ark:/85065/d7wd3xhc) and Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation.

The six-moment solver is built within BATS-R-US and publicly available through the csem.engin.umich.edu/tools/swmf website after registration.

Appendix A Hyperbolic/parabolic cleaning for the Maxwell equations

In this appendix, we briefly derive how the hyperbolic/parabolic cleaning works for the Maxwell equations. We start from the modified Maxwell equations:

[TABLE]

Taking ($\nabla\cdot$) of Equation (36a) gives

[TABLE]

With the expression $\nabla\cdot\mathbf{B}=-\frac{1}{c_{B}}\frac{\partial\psi}{\partial t}-\frac{d_{B}}{c_{B}}\psi$ obtained from Equation (36c), Equation (37) can be written as

[TABLE]

which is the damped wave equation, so the hyperbolic/parabolic cleaning variable $\psi$ propagates isotropically with speed $c_{B}$ and decays at a rate $d_{B}$.

Taking $(\frac{1}{c_{B}}\frac{\partial}{\partial t}+\frac{d_{B}}{c_{B}})$ of Equation (37) gives

[TABLE]

Taking the Laplace operator of the relationship $\nabla\cdot\mathbf{B}=-\frac{1}{c_{B}}\frac{\partial\psi}{\partial t}-\frac{d_{B}}{c_{B}}\psi$ from Equation (36c) gives

[TABLE]

so that Equation (39) becomes:

[TABLE]

which is the same damped wave equation as was obtained for $\psi$ showing that $(\nabla\cdot\mathbf{B})$ will also propagate with speed $c_{B}$ and decay at a rate $d_{B}$.

In a similar fashion, we take $\nabla\cdot$ of Equation (36b) and obtain

[TABLE]

With the expression $\nabla\cdot\mathbf{E}=\frac{1}{\varepsilon_{0}}\rho_{c}-\frac{d_{E}}{c_{E}}\phi-\frac{1}{c_{E}}\frac{\partial\phi}{\partial t}$ obtained from Equation (36d) and $\frac{\partial\rho_{c}}{\partial t}+\nabla\cdot\mathbf{j}=0$ from the ion and electron continuity equations, Equation (42) can be written as

[TABLE]

which shows that the hyperbolic/parabolic cleaning variable $\phi$ for the electric field $\mathbf{E}$ propagates with $c_{E}$ and has a decay rate of $d_{E}$.

Taking $(\frac{1}{c_{E}}\frac{\partial}{\partial t}+\frac{d_{E}}{c_{E}})$ of Equation (42) gives

[TABLE]

Taking the Laplace operator of $\nabla\cdot\mathbf{E}=-\frac{1}{c_{E}}\frac{\partial\phi}{\partial t}-\frac{d_{E}}{c_{E}}\phi+\frac{1}{\varepsilon_{0}}\rho_{c}$ from Equation (36d) gives $\nabla^{2}(\nabla\cdot\mathbf{E})=-(\frac{1}{c_{E}}\frac{\partial}{\partial t}+\frac{d_{E}}{c_{E}})\nabla^{2}\phi+\frac{1}{\varepsilon_{0}}\nabla^{2}\rho_{c}$, which can be substituted into Equation (44) and arrive at

[TABLE]

which shows that $(\nabla\cdot\mathbf{E}-\frac{\rho_{c}}{\varepsilon_{0}})$ satisfies the damped wave equation and has the same behavior as $\phi$.