# A stochastic kinetic scheme for multi-scale plasma transport with
uncertainty quantification

Tianbai Xiao<EMAIL_ADDRESS>Martin Frank<EMAIL_ADDRESS>Karlsruhe Institute of Technology, Karlsruhe, Germany

###### Abstract

Plasmas present a diverse set of behaviors in different regimes. Given the
intrinsic multi-scale nature of plasma dynamics, classical theoretical and
numerical methods are often employed at separate scales with corresponding
assumptions and approximations. Clearly, the coarse-grained modeling may
introduce considerable uncertainties between the field solutions of flow and
electromagnetic variables, and the real plasma physics. To study the
emergence, propagation and evolution of randomness from gyrations of charged
particles to magnetohydrodynamics poses great opportunities and challenges to
develop both sound theories and reliable numerical algorithms. In this paper,
a physics-oriented stochastic kinetic scheme will be developed that includes
random inputs from both flow and electromagnetic fields via a hybridization of
stochastic Galerkin and collocation methods. Based on the BGK-type relaxation
model of the multi-component Boltzmann equation, a scale-dependent kinetic
central-upwind flux function is designed in both physical and particle
velocity space, and the governing equations in the discrete temporal-spatial-
random domain are constructed. By solving Maxwell’s equations with the wave-
propagation method, the evolutions of ions, electrons and electromagnetic
field are coupled throughout the simulation. We prove that the scheme is
formally asymptotic-preserving in the Vlasov, magnetohydrodynamical, and
neutral Euler regimes with the inclusion of random variables. Therefore, it
can be used for the study of multi-scale and multi-physics plasma system under
the effects of uncertainties, and provide scale-adaptive physical solutions
under different ratios among numerical cell size, particle mean free path and
gyroradius (or time step, local particle collision time and plasma period).
Numerical experiments including one-dimensional Landau Damping, the two-stream
instability and the Brio-Wu shock tube problem with one- to three-dimensional
velocity settings, and each under stochastic initial conditions with one-
dimensional uncertainty, will be presented to validate the scheme.

###### keywords:

Boltzmann equation, plasma physics, kinetic theory, uncertainty
quantification, multi-scale methods

##  1 Introduction

Plasma applications cover an extremely wide range of density $n$ from $10^{6}$
to $10^{34}\ \rm m^{-3}$ and temperature $k_{B}T$ from $10^{-1}$ to $10^{10}\
\rm eV$ [1]. As F. Chen wrote in his famous monograph [2], ”What makes plasmas
particularly difficult to analyze is the fact that the densities fall in an
intermediate range. They behave sometimes like fluids, and sometimes like a
collection of individual particles.” A qualitative classification of typical
plasma regimes with respect to density and temperature is presented in Fig. 1.
In this paper, we confine ourselves to the study of classical, non-
relativistic and weakly coupled plasmas, e.g. the magnetosphere, whose
behaviors can be well described by the kinetic theory of gases.

Figure 1: An illustrative demonstration of different plasma regimes and
related phenomena [1].

Given the intrinsic multi-scale nature in plasma physics, classical theories
are devoted to different governing equations at different hierarchies. For
example, at a particle mean free path which is larger than the Debye length
$\lambda_{D}=(\varepsilon_{0}k_{B}T_{ele}/ne^{2})^{1/2}$, and a mean collision
time which is larger than the reciprocal of the plasma frequency
$\omega_{p}=(ne^{2}/\varepsilon_{0}m_{ion})^{1/2}$, the motions of ions and
electrons can be depicted statistically through kinetic equations, e.g. the
Vlasov equation and the Fokker-Planck-Landau equation. On the other hand, at
the macroscopic level with intensive intermolecular collisions, the fluid
dynamic equations are routinely used to model the collective behaviors of
charged particles, i.e. the magnetohydrodynamics (MHD) equations.

Since the 1950s, rapid development has been made in deterministic numerical
methods for plasma simulations [3, 4, 5, 6, 7, 8, 9, 10]. However, given the
coarse-grained approximation in the field theories of plasmas and errors
inherited from numerical simulations, considerable uncertainties may be
introduced inevitably. One typical example are the uncertain inputs of the
initial/boundary value problem. Furthermore, for the evaluation of collision
kernel in the kinetic equations, the phenomenological model parameters often
need to be calibrated by experiments to reproduce correct transport
coefficients, which introduce errors into the simulations. Another example
goes to the vacuum permeability employed in the Maxwell’s equations for
electromagnetic fields. In the SI system which has gone into force in 2019
[11], this value is measured experimentally as $\mu_{0}=8.8541878128\times
10^{-12}\rm F\cdot m^{-1}$, with a relative standard deviation being
$1.5\times 10^{-10}$.

Uncertainty quantification (UQ) is a thriving subject that quantifies one’s
lack of knowledge concerning a physical reality. It applies itself to answer
the challenging questions, e.g. how predictive are the simulation results from
the idealized models, and how can one explicitly assess the effects of
uncertainties on the quality of model predictions. Depending on the
methodology to model the random variables, the methods for UQ study can be
classified into intrusive and non-intrusive ones. In the former case, a series
of realizations of random inputs are generated based on a prescribed
probability distribution. Each realization is solved by a deterministic
solver, and then a post-processing is employed to estimated the uncertainties.
In contrast, intrusive methods work in a way such that we reformulate the
original deterministic system.

One prevalent intrusive strategy is the stochastic Galerkin (SG) method [12],
in which the solutions are expressed into orthogonal polynomials of the input
random parameters. It promises spectral convergence in random space when the
solution depends smoothly on the stochastic parameters [13]. In a nonlinear
Galerkin system, all the expansion coefficients are essentially coupled, which
becomes cumbersome in massive computations. The stochastic collocation (SC)
method [14], although a non-intrusive method, can be seen as a middle way. It
combines the strengths of non-intrusive sampling and SG by evaluating the
generalized polynomial chaos (gPC) expansions on quadrature points in random
space. As a result, a set of decoupled equations can be derived and solved
with deterministic solvers on each quadrature point. Provided the solutions
posses sufficient smoothness over random space, the SC methods maintain
similar convergence as SG, but suffers from aliasing errors due to limited
number of quadrature points.

Although the UQ field has undergone rapid development over the past few years,
its applications on plasma physics mainly focus on the two limits of Vlasov
[15, 16, 17, 18] and MHD [19, 20] with standard stochastic settings. To the
best of the authors’ knowledge, limited work has been conducted on the
evolutionary process of uncertainty in multi-scale physics. Given the
nonlinear system including intermolecular collisions, initial inputs, fluid-
surface interactions and geometric complexities, uncertainties may emerge from
molecular-level nature, develop upwards, affect macroscopic collective
behaviors, and vice versa. To study the emergence, propagation and evolution
of uncertainty poses great opportunities and challenges to develop both sound
theories and reliable multi-scale numerical algorithms.

It is noticeable that tracking the evolution of stochastic variables with
either polynomial chaos or quadrature rules is similar in spirit to solving
kinetic equations in phase space with moment or discrete velocity methods. The
advantages of SG and SC methods can be combined when the integrals that are
necessary for SG inside the algorithm are computed numerically using SC. In
this paper, we follow the strategy proposed in [21] and develop a stochastic
kinetic scheme for multi-scale plasma transport and we couple it to Maxwell’s
equations. Based on the Bhatnagar-Gross-Krook (BGK) type relaxation model for
multi-component plasmas, a scale-dependent central-upwind flux function is
constructed in both physical and particle velocity space, which considers
simultaneously the individual particle transports and their collective
behaviors. The update of source terms of plasma and electromagnetic fields are
solved in a coupled way implicitly. We thus combine the advantages of SG and
SC methods with the construction principle of kinetic schemes, and obtain an
efficient and accurate scheme for cross-scale BGK-Maxwell system with
uncertainties. The randomly initial inputs of both flow and electromagnetic
fields are considered.

The rest of this paper is organized as follows. Sec. 2 is a brief introduction
of kinetic theory of plasma and its stochastic formulation. Sec. 3 presents
the numerical implementation of the current scheme and detailed solution
algorithm. Sec. 4 includes numerical experiments to demonstrate the
performance of the current scheme and analyze some new physical observations.
The last section is the conclusion.

##  2 Deterministic and stochastic theories

###  2.1 Kinetic theory of plasmas

The gas kinetic theory describes the time-space evolution of particle
distribution function. With a separate modeling of particle transport and
collision processes, the evolution equation of monatomic plasmas writes as

$\frac{\partial f_{\alpha}}{\partial
t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}f_{\alpha}+\frac{q_{\alpha}}{m_{\alpha}}\left(\mathbf{E}+\mathbf{u}\times\mathbf{B}\right)\cdot\nabla_{\mathbf{u}}f_{\alpha}=Q_{\alpha},$
(1)

where $\alpha=ion,ele$ denotes a specific ion or electron,
$(q_{\alpha},m_{\alpha})$ are particle charge and mass, and $(\bf E,B)$ are
electric and magnetic fields respectively. For the Coulomb collisions between
charged particles, the limiting case of Boltzmann collision integral leads to
the Fokker-Planck-Landau operator,

$Q_{\alpha}=\sum_{\beta}^{N}\left\\{\nabla_{\mathbf{u}}\cdot\int\mathbf{\Phi}(\mathbf{u}-\mathbf{u}^{\prime})\left[f_{\beta}(\mathbf{u}^{\prime})\nabla_{\mathbf{u}}f_{\alpha}(\mathbf{u})-f_{\alpha}(\mathbf{u})\nabla_{\mathbf{u}}f_{\beta}(\mathbf{u}^{\prime})\right]d\mathbf{u}^{\prime}\right\\},$
(2)

where $(\mathbf{u},\ \mathbf{u}^{\prime})$ are the velocities of two classes
of particles and
$\mathbf{\Phi}(\mathbf{u})=(|\mathbf{u}|^{2}\mathbf{I}_{3}-\mathbf{u}\otimes\mathbf{u})/|\mathbf{u}|^{3}$
is a $3\times 3$ matrix.

Instead of solving the above collision integral directly, here we employ a
multi-component BGK-type relaxation model proposed by Andries, Aoki and
Perthame [22] in the current work to mimic the collision process,

$Q_{\alpha}=\nu_{\alpha}(\mathcal{M}_{\alpha}-f_{\alpha}),$ (3)

where $\nu_{\alpha}$ is the collision frequency. The equilibrium distribution
is defined based on the local modified macroscopic variables, i.e.,

$\mathcal{M}_{\alpha}=n_{\alpha}\left(\frac{m_{\alpha}}{2\pi
k_{B}\bar{T}_{\alpha}}\right)^{\frac{3}{2}}\exp\left(-\frac{m_{\alpha}}{2k_{B}\bar{T}_{\alpha}}(\mathbf{u}-\mathbf{\bar{U}}_{\alpha})^{2}\right),$
(4)

where $n_{\alpha}$ is number density, $m_{\alpha}$ is molecular mass and
$k_{B}$ is the Boltzmann constant. The design of modified temperature
$\bar{T}_{\alpha}$ and velocity $\mathbf{\bar{U}}_{\alpha}$ is based on the
idea that the macroscopic transfer rates in the moments equations derived from
BGK model should be consistent with that from multi-component Boltzmann
equation. For elastic scattering, the evaluation of modified variables for
Maxwell and hard sphere molecules can be written as,

$\displaystyle\mathbf{\bar{U}}_{\alpha}=\mathbf{U}_{\alpha}+\tau_{\alpha}\sum_{r}\frac{2m_{r}}{m_{\alpha}+m_{r}}\nu_{\alpha
r}(\mathbf{U}_{r}-\mathbf{U}_{\alpha}),$ (5)
$\displaystyle\frac{3}{2}k_{B}\bar{T}_{\alpha}=\frac{3}{2}k_{B}T_{\alpha}-\frac{m_{\alpha}}{2}(\mathbf{\bar{U}}_{\alpha}-\mathbf{U}_{\alpha})^{2}$
$\displaystyle+\tau_{\alpha}\sum_{r}\frac{4m_{\alpha}m_{r}}{(m_{\alpha}+m_{r})^{2}}\nu_{\alpha
r}\left[\frac{3}{2}k_{B}T_{r}-\frac{3}{2}k_{B}T_{\alpha}+\frac{m_{r}}{2}(\mathbf{U}_{r}-\mathbf{U}_{\alpha})^{2}\right],$

where $\nu_{\alpha r}$ is the frequency of intermolecular interactions which
can be derived through specific molecule models [23], and it determines the
relaxation time by $\tau_{\alpha}=1/\sum_{r}^{N}\nu_{\alpha r}$. Here we adopt
the hard-sphere molecules, i.e.,

$\nu_{\alpha
r}=\frac{4\sqrt{\pi}n_{r}}{3}\left(\frac{2k_{B}T_{\alpha}}{m_{\alpha}}+\frac{2k_{B}T_{r}}{m_{r}}\right)^{1/2}\left(\frac{d_{\alpha}+d_{r}}{2}\right)^{2},$
(6)

where $d$ is the kinetic molecule diameter.

Macroscopic conservative flow variables are related to the moments of the
particle distribution function,

$\mathbf{W}_{\alpha}=\left(\begin{matrix}\rho_{\alpha}\\\
\rho_{\alpha}\mathbf{U}_{\alpha}\\\
\rho_{\alpha}\mathscr{E}_{\alpha}\end{matrix}\right)=\int
m_{\alpha}f_{\alpha}\varpi d\mathbf{u},$ (7)

where
$\mathscr{E}_{\alpha}=(\mathbf{U}_{\alpha})^{2}/2+3k_{B}T_{\alpha}/2m_{\alpha}$
is total energy density, and
$\varpi=\left(1,\mathbf{u},\frac{1}{2}\mathbf{u}^{2}\right)^{T}$ is the vector
of collision invariants. Hence, macroscopic transport equations can be derived
by taking moments of the kinetic equation with respect to the collision
invariants, i.e.,

$\displaystyle\frac{\partial n_{\alpha}}{\partial
t}+\nabla_{\mathbf{x}}\cdot(n_{\alpha}\mathbf{U}_{\alpha})=0,$ (8)
$\displaystyle\frac{\partial(\rho_{\alpha}\mathbf{U}_{\alpha})}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho_{\alpha}\mathbf{U}_{\alpha}\mathbf{U}_{\alpha})=\nabla_{\mathbf{x}}\cdot\mathbf{P}_{\alpha}+n_{\alpha}q_{\alpha}(\mathbf{E}+\mathbf{U}_{\alpha}\times\mathbf{B})+\mathbf{R}_{\alpha},$
$\displaystyle\frac{\partial(\rho_{\alpha}\mathscr{E}_{\alpha})}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho_{\alpha}\mathscr{E}_{\alpha}\mathbf{U}_{\alpha})=\nabla_{\mathbf{x}}\cdot(\mathbf{P}_{\alpha}\mathbf{U}_{\alpha})-\nabla_{\mathbf{x}}\cdot\mathbf{q}_{\alpha}+n_{\alpha}q_{\alpha}\mathbf{U}_{\alpha}\cdot\mathbf{E}+H_{\alpha},$

where the source terms $\mathbf{R}_{\alpha}$ and $H_{\alpha}$ in the balance
laws come from the moments of collision term respectively,

$\displaystyle\mathbf{R}_{\alpha}$
$\displaystyle=\int\mathbf{u}m_{\alpha}{\nu_{\alpha}}(\mathcal{M}_{\alpha}-f_{\alpha})d\mathbf{u}=\sum_{r}\frac{2m_{\alpha}m_{r}}{m_{\alpha}+m_{r}}n_{\alpha}v_{\alpha
r}\left(\mathbf{U}_{r}-\mathbf{U}_{\alpha}\right),$ (9) $\displaystyle
H_{\alpha}$
$\displaystyle=\int\frac{1}{2}(\mathbf{u}-\mathbf{U})^{2}m_{\alpha}{\nu_{\alpha}}(\mathcal{M}_{\alpha}-f_{\alpha})d\mathbf{u}$
$\displaystyle=\sum_{r}\frac{4m_{\alpha}m_{r}}{\left(m_{\alpha}+m_{r}\right)^{2}}n_{\alpha}v_{\alpha
r}\left[\frac{3}{2}k_{B}T_{r}-\frac{3}{2}k_{B}T_{\alpha}+\frac{m_{r}}{2}\left(\mathbf{U}_{r}-\mathbf{U}_{\alpha}\right)^{2}\right].$

Eq.(8) is consistent with the well-known Braginskii’s two-fluid model [24].

###  2.2 Generalized polynomial chaos approximation of kinetic equation with
uncertainties

Several sources of uncertainty can be considered in the BGK equation. Here we
consider uncertain initial and boundary conditions, which turn the
deterministic system into stochastic case. We employ the generalized
polynomial chaos (gPC) expansion of particle distribution with degree $N$,
i.e.,

$f_{\alpha}(t,\mathbf{x},\mathbf{u},\mathbf{z})\simeq f_{\alpha
N}=\sum_{|i|=0}^{N}\hat{f}_{\alpha
i}(t,\mathbf{x},\mathbf{u})\Phi_{i}(\mathbf{z})=\hat{\bm{f}}_{\alpha}^{T}\bm{\Phi},$
(10)

where i could be a scalar or a $K$-dimensional vector
$i=(i_{1},i_{2},\cdots,i_{K})$ with $|i|=i_{1}+i_{2}+\cdots+i_{K}$. The
$\hat{f}_{\alpha\mathbf{i}}$ is the coefficient of $i$-th polynomial chaos
expansion, and the basis functions used are orthogonal polynomials
{$\Phi_{i}(\mathbf{z})$} satisfying the following constraints,

$\mathbb{E}[\Phi_{j}(\mathbf{z})\Phi_{k}(\mathbf{z})]=\gamma_{k}\delta_{jk},\quad
0\leq|j|,|k|\leq N,$ (11)

where

$\gamma_{k}=\mathbb{E}[\Phi_{k}^{2}(\mathbf{z})],\quad 0\leq|\mathbf{k}|\leq
N,$ (12)

are the normalization factors. The expectation value defines a scalar product,

$\mathbb{E}[\Phi_{j}(\mathbf{z})\Phi_{k}(\mathbf{z})]=\int_{I_{\mathbf{z}}}\Phi_{j}(\mathbf{z})\Phi_{k}(\mathbf{z})\varrho(\mathbf{z})d\mathbf{z},$
(13)

where $\varrho(\mathbf{z})$ is the probability density function. In practice,
it can be evaluated theoretically or with numerical quadrature rule, i.e.,

$\mathbb{E}[\Phi_{j}(\mathbf{z})\Phi_{k}(\mathbf{z})]=\sum_{i}\Phi_{j}(\mathbf{z}_{i})\Phi_{k}(\mathbf{z}_{i})w(\mathbf{z}_{i}),$
(14)

where $w(\mathbf{z}_{i})$ is the corresponding quadrature weight function in
random space. In the following we adopt a uniform notation
$\langle\Phi_{j}\Phi_{k}\rangle$ to denote the integrals over random space
from Eq.(13) and (14).

Given the correspondence between macroscopic and mesoscopic variables, from
Eq.(7) we can derive,

$\displaystyle\mathbf{W}_{\alpha}$ $\displaystyle\simeq\int f_{\alpha N}\varpi
d\mathbf{u}=\int\sum_{i}^{N}\hat{f}_{\alpha
i}(t,\mathbf{x},\mathbf{u})\Phi_{i}(\mathbf{z})\varpi
d\mathbf{u}=\sum_{i}\left(\int\hat{f}_{\alpha i}\varpi
d\mathbf{u}\right)\Phi_{i}$ (15) $\displaystyle\simeq\mathbf{W}_{\alpha
N}=\sum_{i}^{N}\hat{w}_{\alpha i}\Phi_{i}=\hat{\bm{w}}_{\alpha}^{T}\bm{\Phi}.$

After substituting the Eq.(10) into the kinetic equation (1) and (3), and
performing a Galerkin projection, we then obtain

$\frac{\partial\hat{\bm{f}}_{\alpha}}{\partial
t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}\hat{\bm{f}}_{\alpha}+\hat{\bm{G}}_{\alpha}=\hat{\bm{Q}}_{\alpha}=\nu_{\alpha}(\hat{\bm{m}}_{\alpha}-\hat{\bm{f}}_{\alpha}),$
(16)

where $\hat{\bm{Q}}_{\alpha}$ is the gPC coefficient vector of the projection
from collision operator to the polynomial basis,

$\hat{\bm{Q}}_{\alpha}=\nu_{\alpha}(\hat{\bm{m}}_{\alpha}-\hat{\bm{f}}_{\alpha}),$
(17)

with $\hat{\mathbf{m}}_{\alpha}$ being the vector of gPC coefficients of
Maxwellian distribution (which depends implicitly on the
$\hat{\bm{f}}_{\alpha}$) and $\nu_{\alpha}$ being a deterministic collision
frequency. The electromagnetic forcing term $\hat{\bm{G}}_{\alpha}$ is

$\hat{\bm{G}}_{\alpha}=\hat{\bm{G}}_{\alpha}^{T}\bm{\Phi}=\sum_{i}^{N}\hat{G}_{i}\Phi_{i},\
\hat{G}_{{i}}=\frac{q_{\alpha}}{m_{\alpha}}\frac{\sum_{{j}}^{N}\sum_{{k}}^{N}\left(\hat{\mathbf{E}}_{j}+\mathbf{u}\times\hat{\mathbf{B}}_{j}\right)\nabla_{\mathbf{u}}\hat{f}_{\alpha
k}\left\langle\Phi_{{j}}\Phi_{{k}}\Phi_{{i}}\right\rangle}{\left\langle\Phi_{{i}}^{2}\right\rangle}.$
(18)

Notice that the gPC coefficients of both $\hat{\bm{G}}_{\alpha}$ and
$\hat{\bm{Q}}_{\alpha}$ are nonlinear functions of the state variables (cf.
Eq.(5)).

###  2.3 Maxwellian distribution in generalized polynomial chaos

For a deterministic system, the evaluation of the Maxwellian distribution
given in Eq.(4) is straight-forward. However, given a generalized polynomial
chaos (gPC) system, the multiplication and division can’t be operated directly
on the stochastic moments without modifying the orthogonal basis. Starting
from a known particle distribution function in Eq.(10), here we draw a brief
outline to approximately evaluate the Maxwellian distribution function in the
gPC expansion.

1\. Derive the macroscopic conservative variables from particle distribution
function with gPC expansion,

$\mathbf{W}_{\alpha N}=\left(\begin{matrix}\rho_{\alpha N}\\\
(\rho_{\alpha}\mathbf{U}_{\alpha})_{N}\\\
(\rho_{\alpha}\mathscr{E}_{\alpha})_{N}\end{matrix}\right)=\sum_{i}^{N}\left(\int\hat{f}_{\alpha
i}\psi d\mathbf{u}\right)\Phi_{i};$ (19)

2\. Locate conservative variables on quadrature points $\mathbf{z}_{j}$ of
random space and calculate primitive variables, e.g. flow velocity

$\mathbf{U}_{\alpha}(\mathbf{z}_{j})=\frac{(\rho_{\alpha}\mathbf{U}_{\alpha})_{N}(\mathbf{z}_{j})}{\rho_{\alpha
N}(\mathbf{z}_{j})},$ (20)

and

$T_{\alpha}(\mathbf{z}_{j})=\frac{(\rho_{\alpha}E_{\alpha})_{N}(\mathbf{z}_{j})-(\rho_{\alpha}U_{\alpha})_{N}^{2}(\mathbf{z}_{j})/2\rho_{\alpha
N}(\mathbf{z}_{j})}{3k_{B}\rho_{\alpha N}/2m_{\alpha}},$ (21)

and then calculate the modified velocity and temperature via Eq.(5).

3\. Calculate Maxwellian distribution on quadrature points

$\mathcal{M}_{\alpha}(\mathbf{z}_{j})=n_{\alpha}(\mathbf{z}_{j})\left(\frac{m_{\alpha}}{2\pi
k_{B}\bar{T}_{\alpha}(\mathbf{z}_{j})}\right)^{\frac{3}{2}}e^{-\bar{\lambda}_{\alpha}(\mathbf{z}_{j})(\mathbf{u}-\bar{\mathbf{U}}_{\alpha}(\mathbf{z}_{j}))^{2}},$
(22)

and decompose it into a gPC expansion

$\mathcal{M}_{\alpha N}=\sum_{i}^{N}\hat{m}_{\alpha i}\Phi_{i},$ (23)

with each coefficient in the expansion being given by a quadrature rule

$\hat{m}_{\alpha
i}=\frac{\langle\mathcal{M}_{\alpha},\Phi_{i}\rangle}{\langle\Phi_{i}^{2}\rangle}=\frac{\sum_{j}\mathcal{M}_{\alpha}(\mathbf{z}_{j})\Phi_{\mathbf{i}}(\mathbf{z}_{j})\varrho(\mathbf{z}_{j})}{\int_{I_{\mathbf{z}}}(\Phi_{i}(\mathbf{z}))^{2}\varrho(\mathbf{z})d\mathbf{z}}.$
(24)

###  2.4 Maxwell’s equations

For the self-consistent problem of plasma dynamics, the evolutions of electric
and magnetic fields $(\mathbf{E},\mathbf{B})$ are coupled with the motions of
charged particles, which can be described by the linear Maxwell’s equations in
vacuum,

$\displaystyle\frac{\partial\mathbf{E}}{\partial
t}-c^{2}\nabla_{\mathbf{x}}\times\mathbf{B}=-\frac{1}{\varepsilon_{0}}\mathbf{J},$
(25) $\displaystyle\frac{\partial\mathbf{B}}{\partial
t}+\nabla_{\mathbf{x}}\times\mathbf{E}=0,$
$\displaystyle\nabla_{\mathbf{x}}\cdot\mathbf{E}=\frac{\sigma}{\varepsilon_{0}},$
$\displaystyle\nabla_{\mathbf{x}}\cdot\mathbf{B}=0.$

Here $\sigma=e(n_{i}-n_{e})$ is the net charge density, $\mathbf{J}$ is the
current, and the speed of light is related to the permeability and
permittivity of vacuum with $c=(\mu_{0}\varepsilon_{0})^{-1/2}$. To ensure the
divergence constraints in numerical simulations, some techniques can be used
in solving Maxwell’s equations. Here we employ the perfectly hyperbolic
Maxwell’s equations (PHM) [25],

$\displaystyle\frac{\partial\mathbf{E}}{\partial
t}-c^{2}\nabla_{\mathbf{x}}\times\mathbf{B}+\chi
c^{2}\nabla_{\mathbf{x}}\phi=-\frac{1}{\varepsilon_{0}}\mathbf{J},$ (26)
$\displaystyle\frac{\partial\mathbf{B}}{\partial
t}+\nabla_{\mathbf{x}}\times\mathbf{E}+\gamma\nabla_{\mathbf{x}}\psi=0,$
$\displaystyle\frac{1}{\chi}\frac{\partial\phi}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathbf{E}=\frac{\sigma}{\varepsilon_{0}},$
$\displaystyle\frac{\varepsilon_{0}\mu_{0}}{\gamma}\frac{\partial\psi}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathbf{B}=0,$

where $\phi,\psi$ are two additional correction potentials, and the
propagation speed of errors for the divergence of magnetic and electric fields
are $\gamma c$ and $\chi c$ correspondingly. With the stochastic Galerkin
formulation, the PHM system can be rewritten as

$\displaystyle\frac{\partial\hat{\bm{E}}}{\partial
t}-c^{2}\nabla_{\mathbf{x}}\times\hat{\bm{B}}+\chi
c^{2}\nabla_{\mathbf{x}}\hat{\bm{\phi}}=-(\hat{\bm{J}/\varepsilon_{0}}),$ (27)
$\displaystyle\frac{\partial\hat{\bm{B}}}{\partial
t}+\nabla_{\mathbf{x}}\times\hat{\bm{E}}+\gamma\nabla_{\mathbf{x}}\hat{\bm{\psi}}=0,$
$\displaystyle\frac{1}{\chi}\frac{\partial\hat{\bm{\phi}}}{\partial
t}+\nabla_{\mathbf{x}}\cdot\hat{\bm{E}}=(\hat{\sigma/\varepsilon_{0}}),$
$\displaystyle\frac{\varepsilon_{0}\mu_{0}}{\gamma}\frac{\partial\hat{\bm{\psi}}}{\partial
t}+\nabla_{\mathbf{x}}\cdot\hat{\bm{B}}=0.$

The parameters in Maxwell’s equations are assumed to be deterministic.

##  3 Solution algorithm

###  3.1 Update algorithm

The current numerical algorithm is implemented within the finite volume
framework. We adopt the notation of cell averaged macroscopic conservative
variables and particle distribution function in a control volume,

$\displaystyle\mathbf{W}_{\alpha}({t^{n},\mathbf{x}_{i},\mathbf{z}_{k}})=(\mathbf{W}_{\alpha})_{i,k}^{n}=\frac{1}{\Omega_{i}(\mathbf{x})\Omega_{k}(\mathbf{z})}\int_{\Omega_{i}}\int_{\Omega_{k}}\mathbf{W}_{\alpha}(t^{n},\mathbf{x},\mathbf{z})d\mathbf{x}d\mathbf{z},$
$\displaystyle
f_{\alpha}({t^{n},\mathbf{x}_{i},\mathbf{u}_{j},\mathbf{z}_{k}})=(f_{\alpha})_{i,j,k}^{n}=\frac{1}{\Omega_{i}(\mathbf{x})\Omega_{j}(\mathbf{u})\Omega_{k}(\mathbf{z})}\int_{\Omega_{i}}\int_{\Omega_{j}}\int_{\Omega_{k}}f_{\alpha}(t^{n},\mathbf{x},\mathbf{u},\mathbf{z})d\mathbf{x}d\mathbf{u}d\mathbf{z},$

along with the coefficient vector in the gPC expansions,

$\displaystyle\hat{\bm{W}}_{\alpha}(t^{n},\mathbf{x}_{i})=(\hat{\bm{W}}_{\alpha})_{i}^{n}=\frac{1}{\Omega_{i}(\mathbf{x})}\int_{\Omega_{i}}\hat{\bm{W}}_{\alpha}(t^{n},\mathbf{x})d\mathbf{x},$
$\displaystyle\hat{\bm{f}}_{\alpha}(t^{n},\mathbf{x}_{i},\mathbf{u}_{j})=(\hat{\bm{f}}_{\alpha})_{i,j}^{n}=\frac{1}{\Omega_{i}(\mathbf{x})\Omega_{j}(\mathbf{u})}\int_{\Omega_{i}}\int_{\Omega_{j}}\hat{\bm{f}}_{\alpha}(t^{n},\mathbf{x},\mathbf{u})d\mathbf{x}d\mathbf{u},$

where $\Omega_{i}$, $\Omega_{j}$ and $\Omega_{k}$ are the cell area in the
discretized physical, velocity and random space.

The update of the stochastic Galerkin coefficients for particle distribution
function can be written as,

$\displaystyle(\hat{\bm{f}}_{\alpha})_{i,j}^{n+1}=$
$\displaystyle(\hat{\bm{f}}_{\alpha})_{i,j}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{S_{r}\in\partial\Omega_{i}}S_{r}(\hat{\bm{F}}_{\alpha})_{r,j}^{f}dt$
(28)
$\displaystyle+\frac{1}{\Omega_{j}}\int_{t^{n}}^{t^{n+1}}\sum_{S_{r}\in\partial\Omega_{j}}S_{r}(\hat{\bm{F}}_{\alpha})_{i,r}^{f}dt+\int_{t^{n}}^{t^{n+1}}(\hat{\bm{Q}}_{\alpha})_{i,j}^{f}dt.$

where $(\hat{\bm{F}}_{\alpha})^{f}_{r}$ is the time-dependent fluxes for
distribution function at interface $r$ in physical and velocity space, $S_{r}$
is the interface area, and $\hat{\bm{Q}}_{\alpha}^{f}$ is the collision term.
Taking velocity moments of Eq.(28), we obtain the corresponding macroscopic
system,

$\displaystyle(\hat{\bm{W}}_{\alpha})_{i}^{n+1}=$
$\displaystyle(\hat{\bm{W}}_{\alpha})_{i}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{\mathbf{S}_{r}\in\partial\Omega_{i}}\mathbf{S}_{r}\cdot(\hat{\bm{F}}_{\alpha})_{r}^{W}dt$
(29)
$\displaystyle+\int_{t^{n}}^{t^{n+1}}(\hat{\bm{G}}_{\alpha})_{i}^{W}dt+\int_{t^{n}}^{t^{n+1}}(\hat{\bm{Q}}_{\alpha})_{i}^{W}dt,$

where $(\hat{\bm{F}}_{\alpha})^{W}_{r}$ is the flux functions for macroscopic
conservative variables, $\mathbf{S}_{r}=\bm{n}S_{r}$ is the interface area
vector, and $\hat{\bm{G}}_{\alpha}^{W}$ is the external force terms related to
Eq.(8).

Notice that the gPC coefficients of different orders are coupled through the
external force term $\hat{\bm{G}}_{\alpha}$ (17) and the collision term
$\hat{\bm{Q}}_{\alpha}$ (18). Instead of solving this large nonlinear system,
we can locate the gPC system onto the collocation points in random space. It
results a decoupled system that can be solved efficiently. Therefore, we
combine the advantages of stochastic Galerkin and collcation methods, which is
one of the novelties of this paper. To make use of it, in the solution
algorithm, we first update the gPC coefficients to an intermediate step
$t^{*}$,

$\displaystyle(\hat{\bm{W}}_{\alpha})_{i}^{*}=(\hat{\bm{W}}_{\alpha})_{i}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{\mathbf{S}_{r}\in\partial\Omega_{i}}\mathbf{S}_{r}\cdot(\hat{\bm{F}}_{\alpha})_{r}^{W}dt,$
(30)
$\displaystyle(\hat{\bm{f}}_{\alpha})_{i,j}^{*}=(\hat{\bm{f}}_{\alpha})_{i,j}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{S_{r}\in\partial\Omega_{i}}S_{r}(\hat{\bm{F}}_{\alpha})_{r,j}^{f}dt,$

which are then evaluated on random quadrature cell $\Omega_{k}$,

$\displaystyle\mathbf{W}_{i,k}^{*}=\mathbf{W}_{Ni}^{*}\left(\mathbf{z}_{k}\right)=\sum_{m}^{N}\hat{w}_{i,m}^{*}\left(\mathbf{z}_{k}\right)\Phi_{m}\left(\mathbf{z}_{k}\right),$
(31) $\displaystyle
f_{i,j,k}^{*}=f_{Ni,j}^{*}\left(\mathbf{z}_{k}\right)=\sum_{m}^{N}\hat{f}_{i,j,m}^{*}\left(\mathbf{z}_{k}\right)\Phi_{m}\left(\mathbf{z}_{k}\right).$

Afterwards, the collision and forcing term are evaluated on collocation points
via

$\left(\mathbf{W}_{\alpha}\right)_{i,k}^{n+1}=\left(\mathbf{W}_{\alpha}\right)_{i,k}^{*}+\int_{t^{n}}^{t^{n+1}}\left(\mathbf{G}_{\alpha}\right)_{i,k}^{W}dt+\int_{t^{n}}^{t^{n+1}}\left(\mathbf{Q}_{\alpha}\right)_{i,k}^{W}dt,$
(32)

and

$\left(f_{\alpha}\right)_{i,j,k}^{n+1}=\left(f_{\alpha}\right)_{i,j,k}^{*}+\frac{1}{\Omega_{j}}\int_{t^{n}}^{t^{n+1}}\sum_{S_{r}\in\partial\Omega_{j}}S_{r}\left(F_{\alpha}\right)_{i,r,k}^{f}dt+\int_{t^{n}}^{t^{n+1}}\left(Q_{\alpha}\right)_{i,j,k}^{f}dt,$
(33)

where $(F_{\alpha})^{f}_{r}$ is the numerical flux at interface $r$ in
velocity space. In the solution algorithm loop, Eq.(32) can be solved first,
and then the updated variables at $t^{n+1}$ can be employed to evaluate the
Maxwellian distribution in Eq.(33) implicitly.

For plasma transport, the evolution of electromagnetic field should be solved
in a coupled way with flow field. The hybrid Galerkin-collocation method is
employed as well, where the gPC coefficients are stepped to the intermediate
state first,

$\hat{\bm{M}}_{i}^{*}=\hat{\bm{M}}_{i}^{n}+\frac{1}{\Omega_{i}}\int_{t^{n}}^{t^{n+1}}\sum_{r}\Delta\mathbf{S}_{r}\cdot\hat{\bm{F}}_{r}^{M}dt,$
(34)

where $\hat{\bm{F}}^{M}$ is the flux functions for the electromagnetic fields
$(\hat{\bm{E}},\hat{\bm{B}},\hat{\bm{\phi}},\hat{\bm{\psi}})$. Then the source
terms are solved via,

$\displaystyle{\mathbf{E}}_{i}^{n+1}={\mathbf{E}}_{i}^{*}-\frac{e}{\varepsilon_{0}}\int_{t^{n}}^{t^{n+1}}\left({(n\mathbf{U})}_{ion}-{(n\mathbf{U})}_{ele}\right)dt,$
(35)
$\displaystyle{{\phi}}_{i}^{n+1}={{\phi}}_{i}^{*}+\frac{e}{\varepsilon_{0}}\int_{t^{n}}^{t^{n+1}}\left(n_{ion}-n_{ele}\right)dt.$

###  3.2 Fluxes computed using Stochastic Galerkin

####  3.2.1 Plasma flux

Based on the finite volume framework, a scale-dependent interface flux
function is needed in multi-scale modeling and simulation. Different from a
purely upwind flux which loses efficiency in the collisional limit, we here
develop a kinetic central-upwind flux function based on an integral solution
of the kinetic model equation. The integral solution originates from Kogan’s
monograph on rarefied gas dynamics [26] and has been employed by a series of
gas-kinetic schemes [27, 28, 29, 30, 31]. Let us rewrite the stochastic BGK
equation (16) for the gPC coefficients vector along the characteristics,

$\frac{D\hat{\bm{f}}_{\alpha}}{Dt}+\nu_{\alpha}\hat{\bm{f}}_{\alpha}=\nu_{\alpha}\hat{\bm{m}}_{\alpha}.$
(36)

We assume that the collision frequency $\nu_{\alpha}$ is kept fixed at the
value that can be computed from the macroscopic variables at the previous
timestep. Then the following integral solution holds along the
characteristics,

$\hat{\bm{f}}_{\alpha}(t,\mathbf{x},\mathbf{u})=\nu_{\alpha}\int_{0}^{t}\hat{\bm{m}}_{\alpha}(t^{\prime},\mathbf{x}^{\prime},\mathbf{u}^{\prime})e^{-\nu_{\alpha}(t-t^{\prime})}dt^{\prime}+e^{-\nu_{\alpha}t}\hat{\bm{f}}_{\alpha}(0,\mathbf{x}^{0},\mathbf{u}^{0}),$
(37)

where
$\mathbf{x}^{\prime}=\mathbf{x}-\mathbf{u}^{\prime}(t-t^{\prime})-\frac{1}{2}\mathbf{a}(t-t^{\prime})^{2}$
is the particle trajectory, and
$\mathbf{x}^{0}=\mathbf{x}-\mathbf{u}^{\prime}t-\frac{1}{2}\mathbf{a}t^{2}$ is
the location at initial time $t=0$. The above solution indicates a self-
conditioned mechanism for multi-scale gas dynamics. For example, when the
evolving time $t$ is much less than the mean collision time
$\tau=1/\nu_{\alpha}$, the latter term in Eq.(37) dominates and describes the
free transport of particles. And if $t$ is much larger than $\tau$, the second
term approaches to zero, and then the distribution function will be an
accumulation of equilibrium state along the characteristic lines, which
provides the underlying wave-propagation mechanism for hydrodynamic solutions.
In the following, we present a detailed strategy for the construction of
numerical fluxes. Since the electromagnetic force term $\mathbf{a}$ is a
stochastic variable, we do an operator splitting to evaluate fluxes in
physical and velocity space.

With the simplified notations of physical cell interface
$\mathbf{x}_{i+1/2}=0$ and initial time $t^{n}=0$, Eq.(37) along physical
trajectories of particles can be rewritten into the following form,

$\hat{\bm{f}}_{\alpha}(t,0,\mathbf{u}_{j})=\nu_{\alpha}\int_{0}^{t}\hat{\bm{m}}_{\alpha}(t^{\prime},\mathbf{x}^{\prime},\mathbf{u}_{j})e^{-\nu_{\alpha}(t-t^{\prime})}dt^{\prime}+e^{-\nu_{\alpha}t}\hat{\bm{f}}_{\alpha}(0,-\mathbf{u}_{j}t,\mathbf{u}_{j}),$
(38)

where $\hat{\bm{f}}_{\alpha}(0,-\mathbf{u}_{j}t,\mathbf{u}_{j})$ is the
initial distribution at each time step.

In the numerical scheme, the initial distribution function at cell interface
can be obtained through reconstruction, i.e.,

$(\hat{\bm{f}}_{\alpha})(0,\pm 0,\mathbf{u}_{j})=\left\\{\begin{aligned}
&(\hat{\bm{f}}_{\alpha})_{i+1/2,j}^{L},\quad x=0^{-},\\\
&(\hat{\bm{f}}_{\alpha})_{i+1/2,j}^{R},\quad x=0^{+}.\end{aligned}\right.$
(39)

The initial distributions $(\hat{\bm{f}}_{\alpha})_{i+1/2,j}^{L,R}$ at the
left and right hand sides of a cell interface are obtained through the van-
Leer limiter.

The macroscopic conservative variables in the gPC expansions at the interface
can be evaluated by taking moments over velocity space,

$\hat{\bm{w}}=\sum_{u_{j}>0}\hat{\bm{f}}_{i+1/2,j}^{L}\psi\Delta
u_{j}+\sum_{u_{j}<0}\hat{\bm{f}}_{i+1/2,j}^{R}\psi\Delta u_{j}.$

The equilibrium distribution at interface can be determined as illustrated in
Sec. 2.3.

After all coefficients are obtained, the interface distribution function
becomes

$\displaystyle\hat{\bm{f}}_{\alpha}(t,0,u_{j})=$
$\displaystyle\left(1-e^{-\nu_{\alpha}t}\right)(\hat{\bm{m}}_{\alpha})_{j}$
(40)
$\displaystyle+e^{-\nu_{\alpha}t}\left[(\hat{\bm{f}}_{\alpha})_{i+1/2,j}^{L}H\left[\mathbf{u}_{j}\right]+(\hat{\bm{f}}_{\alpha})_{i+1/2,j}^{R}(1-H\left[\mathbf{u}_{j}\right])\right],$

where $H(\mathbf{u})$ is the heaviside step function. The above interface
distribution function can be regarded as a combination of central difference
and upwind methods. With the variation of the ratio between evolving time $t$
(i.e., the time step in the computation) and collision time
$\tau_{\alpha}=1/\nu_{\alpha}$, the factor $e^{-\nu_{\alpha}t}$ plays as a
modulator and provides a self-adaptive solution between equilibrium and non-
equilibrium physics.

After the coefficients of distribution function at all orders are determined,
the corresponding gPC expansion can be expressed as,

$f_{\alpha
N}(t,0,\mathbf{u}_{j})=\sum_{m=0}^{N}\hat{f}_{m}(t,0,\mathbf{u}_{j})\Phi_{m}(\mathbf{z}),$
(41)

and the corresponding fluxes of particle distribution function and
conservative flow variables can be evaluated via

$\displaystyle{{F}}_{\alpha
N}^{f}(t,0,\mathbf{u}_{j},\mathbf{z})=\mathbf{u}_{j}f_{\alpha
N}(t,0,\mathbf{u}_{j},\mathbf{z}),$ (42) $\displaystyle{\mathbf{F}}_{\alpha
N}^{W}(t,0,\mathbf{z})=\int\mathbf{u}f_{\alpha
N}(t,0,\mathbf{u},\mathbf{z})\varpi d\mathbf{u}\simeq\sum
w_{j}\mathbf{u}_{j}f_{\alpha N}(t,0,\mathbf{u}_{j},\mathbf{z})\varpi_{j},$

where $\mathbf{u}_{j}$ denotes quadrature points in particle velocity space,
and $w_{j}$ is its integral weight in velocity space, and the time-integrated
fluxes in Eq.(29) and (28) can be evaluated with respect to time in Eq.(40).

####  3.2.2 Electromagnetic flux

Besides the flow variables, the numerical fluxes of electromagnetic fields in
Eq.(27) are calculated by the wave-propagation method developed by Hakim et
al. [32].

###  3.3 Fluxes and sources computed using stochastic collocation

Besides the construction of the interface flux, the source terms need to be
evaluated inside each control volume within each time step. In this part, we
show the detailed update algorithm for the collision and external force term
with stochastic collocation method.

####  3.3.1 Macroscopic source terms

Given the intermediate macroscopic variables in the discrete cell
$(\Omega_{i},\Omega_{j},\Omega_{k})$, the macroscopic system writes as
follows,

$\displaystyle(\rho_{\alpha}\mathbf{U}_{\alpha})^{n+1}_{i,k}=(\rho_{\alpha}\mathbf{U}_{\alpha})^{*}_{i,k}+\Delta
tq_{\alpha}\left(n_{\alpha}(\mathbf{E}+\mathbf{U}_{\alpha}\times\mathbf{B})\right)_{i,k}^{n+1}+\int_{t^{n}}^{t^{n+1}}\left(\nu_{\alpha}\rho_{\alpha}(\overline{\mathbf{U}}^{*}_{\alpha}-\mathbf{U}^{*}_{\alpha})\right)_{i,k}dt,$
(43)
$\displaystyle(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{n+1}=(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{*}+\Delta
tq_{\alpha}(n_{\alpha}\mathbf{E}\cdot\mathbf{U}_{\alpha})_{i,k}^{n+1}+\int_{t^{n}}^{t^{n+1}}\left(\nu_{\alpha}(\rho_{\alpha}\overline{\mathscr{E}}^{*}_{\alpha}-\rho_{\alpha}\mathscr{E}^{*}_{\alpha})\right)_{i,k}dt,$
$\displaystyle{\mathbf{E}}_{i,k}^{n+1}={\mathbf{E}}_{i,k}^{*}-\frac{e\Delta
t}{\varepsilon_{0}}\left({(n\mathbf{U})}_{ion}-{(n\mathbf{U})}_{ele}\right)^{n+1}_{i,k},$
$\displaystyle{{\phi}}_{i,k}^{n+1}={{\phi}}_{i,k}^{*}+\frac{e\Delta
t}{\varepsilon_{0}}\left({n}_{ion}-{n}_{ele}\right)^{n+1}_{i,k}.$

We split the above system into two parts. First, the flow variables are
evolved with respect to mixture source terms,

$\displaystyle(\rho_{\alpha}\mathbf{U}_{\alpha})^{**}_{i,k}=(\rho_{\alpha}\mathbf{U}_{\alpha})^{*}_{i,k}+\int_{t^{n}}^{t^{n+1}}\left(\nu_{\alpha}\rho^{n+1}_{\alpha}(\overline{\mathbf{U}}^{*}_{\alpha}-\mathbf{U}^{*}_{\alpha})\right)_{i,k}dt,$
(44)
$\displaystyle(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{**}=(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{*}+\int_{t^{n}}^{t^{n+1}}\left(\nu_{\alpha}\rho_{\alpha}^{n+1}\left(\overline{\mathscr{E}}_{\alpha}^{*}-\mathscr{E}_{\alpha}^{*}\right)\right)_{i,k}dt.$

The relaxation integrals are evaluated through the Rosenbrock method [33] to
overcome possible stiffness.

Then we solve the electromagnetic sources implicitly, in which a linear system
can be solved, i.e.,

$\displaystyle(\rho_{\alpha}\mathbf{U}_{\alpha})^{n+1}_{i,k}=(\rho_{\alpha}\mathbf{U}_{\alpha})^{**}_{i,k}+\Delta
tq_{\alpha}\left(n_{\alpha}(\mathbf{E}+\mathbf{U}_{\alpha}\times\mathbf{B})\right)_{i,k}^{n+1},$
(45)
$\displaystyle(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{n+1}=(\rho_{\alpha}\mathscr{E}_{\alpha})_{i,k}^{**}+\Delta
tq_{\alpha}(n_{\alpha}\mathbf{E}\cdot\mathbf{U}_{\alpha})_{i,k}^{n+1},$
$\displaystyle{\mathbf{E}}_{i,k}^{n+1}={\mathbf{E}}_{i,k}^{*}-\frac{e\Delta
t}{\varepsilon_{0}}\left({(n\mathbf{U})}_{ion}-{(n\mathbf{U})}_{ele}\right)^{n+1}_{i,k},$
$\displaystyle{{\phi}}_{i,k}^{n+1}={{\phi}}_{i,k}^{*}+\frac{e\Delta
t}{\varepsilon_{0}}\left({n}_{ion}-{n}_{ele}\right)^{n+1}_{i,k}.$

####  3.3.2 Particle distribution function

The updated macroscopic variables can be used to step the particle
distribution functions of ion and electron implicitly. Let us consider the
kinetic equation,

$\frac{\partial f_{\alpha}}{\partial
t}+\mathbf{a}_{\alpha}\cdot\nabla_{\mathbf{u}}f_{\alpha}=\nu_{\alpha}(\mathcal{M}_{\alpha}-f_{\alpha}),$
(46)

where the electromagnetic force is

$\mathbf{a}_{\alpha}=\frac{q_{\alpha}}{m_{\alpha}}\left(\mathbf{E}^{n+1}+\mathbf{u}_{\alpha}\times\mathbf{B}^{n+1}\right).$

Making the simplified notations of velocity cell interface
$\mathbf{u}_{j+1/2}=0$ and initial time $t^{n}=0$ again, we write the integral
solution of Eq.(46) as,

$\displaystyle f_{\alpha}(t,\mathbf{x}_{i},0,\mathbf{z}_{k})$
$\displaystyle=\nu_{\alpha}\int_{0}^{t}\mathcal{M}^{*}_{\alpha}(t^{\prime},\mathbf{x}_{i},\mathbf{u}^{\prime},\mathbf{z}_{k})e^{-\nu_{\alpha}(t-t^{\prime})}dt^{\prime}+e^{-\nu_{\alpha}t}f^{*}_{\alpha}(0,\mathbf{x}_{i},-\mathbf{a}_{\alpha}t,\mathbf{z}_{k})$
(47)
$\displaystyle=\mathcal{M}^{*}_{\alpha}(0,\mathbf{x},\mathbf{u}-\mathbf{a}_{\alpha}t,\mathbf{z}_{k})\left(1-\mathbf{e}^{-\nu_{\alpha}t}\right)+f^{*}_{\alpha}(0,\mathbf{x}_{i},\mathbf{u}-\mathbf{a}_{\alpha}t,\mathbf{z}_{k})\mathbf{e}^{-\nu_{\alpha}t}.$

Similar to physical space, the above interface distribution function can also
be regarded as a combination of central difference and upwind methods in
particle velocity space. With the variation of the ratio between evolving time
$t$ (i.e., the time step in the computation) and collision time
$\tau_{\alpha}=1/\nu_{\alpha}$, it provides a self-adaptive solution from
equilibrium to non-equilibrium. Based the above solution, the interface flux
in particle velocity space can be constructed as

${{F}}_{\alpha}^{f}(t,\mathbf{x}_{i},0,\mathbf{z}_{k})=\mathbf{a}_{\alpha}f_{\alpha}(t,\mathbf{x}_{i},0,\mathbf{z}_{k}),$
(48)

and the time-integrated flux can be evaluated directly with respect to $t$.

For the collision term, with the updated flow variables at $t^{n+1}$, an
implicit update can be arranged. Therefore, the update algorithm of particle
distribution function can be written as,

$\left(f_{\alpha}\right)_{i,j,k}^{n+1}=\left(f_{\alpha}\right)_{i,j,k}^{*}+\frac{1}{\Omega_{j}}\int_{t^{n}}^{t^{n+1}}\sum_{S_{r}\in\partial\Omega_{j}}S_{r}\left(F_{\alpha}\right)_{i,r,k}^{f}dt+\Delta
t\nu_{\alpha}^{n+1}\left((\mathcal{M}_{\alpha})_{i,j,k}^{n+1}-(f_{\alpha})_{i,j,k}^{n+1}\right).$
(49)

###  3.4 Time step

In the current scheme, the time step is determined by the Courant-Friedrichs-
Lewy condition in phase space,

$\Delta
t=\mathcal{C}\min\left(\frac{\min(|\Delta\mathbf{x}|)}{\max(|\mathbf{u}|)+\max(|\mathbf{U}|)},\frac{\min(|\Delta\mathbf{x}|)}{c},\frac{\min(|\Delta\mathbf{u}|)}{\max(|\mathbf{a}|)}\right),$
(50)

where $\mathcal{C}$ is the CFL number,
$\mathbf{u}=(\mathbf{u}_{ion},\mathbf{u}_{ele})$ is particle velocity,
$\mathbf{U}=(\mathbf{U}_{ion},\mathbf{U}_{ele})$ is fluid velocity, and $c$ is
speed of light. In the computation, $c$ usually takes a pseudo value which is
less than $3.0\times 10^{8}m/s$ but a few orders larger than particle
velocity.

###  3.5 Asymptotic analysis

In this part we will present theoretical analysis on the current kinetic model
and numerical algorithm, with special focus on their asymptotic limits.

####  3.5.1 Asymptotic limits of the BGK-Maxwell system

Let us return to the BGK equation (1). From Eq.(6), we see that the collision
frequency is positively correlated with plasma density and temperature. When
the plasma gets rarefied, the intensity of collision term decreases
correspondingly. As the collision frequency $\nu_{\alpha}$ goes to zero, the
BGK equation automatically reduces to the Vlasov equation, i.e.,

$\frac{\partial f_{\alpha}}{\partial
t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}f_{\alpha}+\frac{q_{\alpha}}{m_{\alpha}}\left(\mathbf{E}+\mathbf{u}_{\alpha}\times\mathbf{B}\right)\cdot\nabla_{\mathbf{u}}f_{\alpha}=0.$
(51)

On the other hand, as the collision frequency increases, then the two-fluid
system (8) is equivalent with the BGK equation, which describes the motions of
plasma as fluid. Under the quasineutral assumption $n_{i}\simeq n_{e}$ with
$\lambda_{D}\ll L_{0}$ ($L_{0}$ is the characteristic length of system),
$m_{i}\gg m_{e}$ and $\nu_{i},\ \nu_{e}\gg 0$, it can be further degenerated
to the single-fluid Hall-MHD equations,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathbf{U})=0,$ (52)
$\displaystyle\frac{\partial(\rho\mathbf{U})}{\partial
t}+\nabla_{\mathbf{x}}\cdot\left(\rho\mathbf{U}\mathbf{U}\right)=-\nabla_{\mathbf{x}}p+\mathbf{J}\times\mathbf{B},$
$\displaystyle\frac{\partial(\rho\mathscr{E})}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathscr{E}\mathbf{U})=\nabla_{\mathbf{x}}\cdot(p\mathbf{U})+\mathbf{J}\cdot\mathbf{E},$
$\displaystyle\mathbf{E}+\mathbf{U}\times\mathbf{B}=\eta\mathbf{J}+\frac{1}{en}(\mathbf{J}\times\mathbf{B}+\nabla_{\mathbf{x}}p),$
$\displaystyle\frac{\partial\sigma}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathbf{J}=0.$

where $\sigma$ is electric charge density and $\mathbf{J}$ is net current.
Together with Maxwell’s equations, this set is able to describe the
equilibrium state of the plasma.

The specific resistivity $\eta$ is proportional to the interaction frequency
between ions and electrons, i.e.,

$\eta=\frac{m_{e}}{ne^{2}}\nu_{ie}.$ (53)

Therefore, given the fully conductive conditions with $\nu_{ie}\to 0$ and
ignore the Hall current term and pressure gradient in the generalized Ohm’s
law, we get the ideal MHD equations,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathbf{U})=0,$ (54)
$\displaystyle\frac{\partial(\rho\mathbf{U})}{\partial
t}+\nabla_{\mathbf{x}}\cdot\left(\rho\mathbf{U}\mathbf{U}\right)=-\nabla_{\mathbf{x}}p+\frac{\left(\mathbf{B}\cdot\nabla_{\mathbf{x}}\right)\mathbf{B}}{\mu_{0}}-\nabla_{\mathbf{x}}\left(\frac{\mathbf{B}^{2}}{2\mu_{0}}\right),$
$\displaystyle\frac{\partial(\rho\mathscr{E})}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathscr{E}\mathbf{U})=\nabla_{\mathbf{x}}\cdot(p\mathbf{U})+\frac{1}{\mu_{0}}\rho\mathbf{U}\cdot\left(\nabla_{\mathbf{x}}\times\mathbf{B}\times\mathbf{B}\right),$
$\displaystyle\frac{\partial\mathbf{B}}{\partial
t}+\nabla_{\mathbf{x}}\times(\mathbf{U}\times\mathbf{B})=0.$

In contrast to the ideal MHD regime, we can also consider the non-conductive
limit where the interspecies molecular interaction is intensive with
$\nu_{ie}\to\infty$. As a result, the plasma now behaves like dielectric
material, and the two-fluid system deduces to the Euler equations,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathbf{U})=0,$ (55)
$\displaystyle\frac{\partial(\rho\mathbf{U})}{\partial
t}+\nabla_{\mathbf{x}}\cdot\left(\rho\mathbf{U}\mathbf{U}\right)=-\nabla_{\mathbf{x}}p,$
$\displaystyle\frac{\partial(\rho\mathscr{E})}{\partial
t}+\nabla_{\mathbf{x}}\cdot(\rho\mathscr{E}\mathbf{U})=\nabla_{\mathbf{x}}\cdot(p\mathbf{U}).$

####  3.5.2 Interface fluxes

Besides the theoretical modeling, the asymptotic preserving property for the
limiting solutions is also a preferred nature for the kinetic scheme. Let us
consider the interface fluxes constructed in Sec. 3.2 and 3.3 first. In the
collisionless limit where $\nu_{\alpha}$ approaches zero, the relation
$\nu_{\alpha}\Delta t\ll 1$ holds naturally, and the collision term in Eq.(1)
disappears. In this case, the interface distribution function in both Eq.(40)
and (47) go to the non-equilibrium distribution parts. For brevity, we write
them down into collocation form, i.e.,

$\displaystyle(f_{\alpha})_{i+1/2,j,k}=(f_{\alpha})_{i+1/2,j}^{L}H(\mathbf{u}_{j})+(f_{\alpha})_{i+1/2,j}^{R}(1-H(\mathbf{u}_{j})),$
(56)
$\displaystyle(f_{\alpha})_{i,j+1/2,k}=f^{*}_{\alpha}(\mathbf{x}_{i},\mathbf{u}-\mathbf{a}_{\alpha}t,\mathbf{z}_{k}),$

which constitute a fully upwind scheme for the Vlasov equation.

On the other hand, in the hydrodynamic regime with intensive collisions, the
interface distribution becomes,

$\displaystyle(f_{\alpha})_{i+1/2,j,k}=(\mathcal{M}_{\alpha})_{i+1/2,j},$ (57)
$\displaystyle(f_{\alpha})_{i,j+1/2,k}=\mathcal{M}^{*}_{\alpha}(\mathbf{x}_{i},\mathbf{u}-\mathbf{a}_{\alpha}t,\mathbf{z}_{k}),$

leading to near-equilibrium MHD or Euler solutions.

####  3.5.3 Source terms

Besides the flux functions, the treatment of source terms plays an important
role for the asymptotic property of scheme. Let us rewrite the electromagnetic
sources in Eq.(45) as follows,

$\displaystyle(\rho_{ion}\mathbf{U}_{ion})^{n+1}_{i,k}=(\rho_{ion}\mathbf{U}_{ion})^{**}_{i,k}+\Delta
te\left(n_{ion}(\mathbf{E}+\mathbf{U}_{ion}\times\mathbf{B})\right)_{i,k}^{n+1},$
(58)
$\displaystyle(\rho_{ele}\mathbf{U}_{ele})^{n+1}_{i,k}=(\rho_{ele}\mathbf{U}_{ele})^{**}_{i,k}-\Delta
te\left(n_{ele}(\mathbf{E}+\mathbf{U}_{ele}\times\mathbf{B})\right)_{i,k}^{n+1},$
$\displaystyle(\rho_{ion}\mathscr{E}_{ion})_{i,k}^{n+1}=(\rho_{ion}\mathscr{E}_{ion})_{i,k}^{**}+\Delta
te(n_{ion}\mathbf{E}\cdot\mathbf{U}_{ion})_{i,k}^{n+1},$
$\displaystyle(\rho_{ele}\mathscr{E}_{ele})_{i,k}^{n+1}=(\rho_{ele}\mathscr{E}_{ele})_{i,k}^{**}-\Delta
te(n_{ele}\mathbf{E}\cdot\mathbf{U}_{ele})_{i,k}^{n+1},$
$\displaystyle{\mathbf{E}}_{i,k}^{n+1}={\mathbf{E}}_{i,k}^{*}-\frac{e\Delta
t}{\varepsilon_{0}}\left({n}_{ion}\mathbf{U}_{ion}-n_{ele}{\mathbf{U}}_{ele}\right)^{n+1}_{i,k},$
$\displaystyle{{\phi}}_{i,k}^{n+1}={{\phi}}_{i,k}^{*}+\frac{e\Delta
t}{\varepsilon_{0}}\left({n}_{ion}-{n}_{ele}\right)^{n+1}_{i,k}.$

Adding the momentum and energy equations of ion and electron, we get

$\displaystyle(\rho\mathbf{U})^{n+1}_{i,k}-(\rho\mathbf{U})^{**}_{i,k}=\Delta
te\left((n_{ion}-n_{ele})\mathbf{E}\right)^{n+1}_{i,k}+\Delta
t\left((\mathbf{J}_{ion}+\mathbf{J}_{ele})\times\mathbf{B}\right)^{n+1}_{i,k},$
(59)
$\displaystyle(\rho\mathscr{E})^{n+1}_{i,k}-(\rho\mathscr{E})^{**}_{i,k}=\Delta
t\left((\mathbf{J}_{ion}+\mathbf{J}_{ele})\cdot\mathbf{E}\right)^{n+1}_{i,k}.$

With the quasineutral assumption, it can be simplified as

$\displaystyle(\rho\mathbf{U})^{n+1}_{i,k}-(\rho\mathbf{U})^{**}_{i,k}=\Delta
t\left(\mathbf{J}\times\mathbf{B}\right)^{n+1}_{i,k},$ (60)
$\displaystyle(\rho\mathscr{E})^{n+1}_{i,k}-(\rho\mathscr{E})^{**}_{i,k}=\Delta
t\left(\mathbf{J}\cdot\mathbf{E}\right)^{n+1}_{i,k}.$

At the same time, let us multiply the first equation in Eq.(58) by $m_{ele}$
and second by $m_{ion}$ and subtract the latter from the former, which
results,

$\displaystyle
m_{ion}m_{ele}\left(n\left(\mathbf{U}_{ion}-\mathbf{U}_{ele}\right)\right)^{n+1}_{i,k}-m_{ion}m_{ele}\left(n\left(\mathbf{U}_{ion}-\mathbf{U}_{ele}\right)\right)^{**}_{i,k}$
(61) $\displaystyle=\Delta
te\left(m_{ion}+m_{ele}\right)(n\mathbf{E})^{n+1}_{i,k}+\Delta
te\left(n\left(m_{ele}\mathbf{U}_{ion}+m_{ion}\mathbf{U}_{ele}\right)\times\mathbf{B}\right)^{n+1}_{i,k}.$

Given $m_{ion}\gg m_{ele}$, the above the equation can be degenerated to

$\mathbf{E}_{i,k}^{n+1}+\mathbf{U}_{i,k}^{n+1}\times\mathbf{B}_{i,k}^{n+1}=\frac{1}{e\rho_{i,k}^{n+1}}\left(\frac{m_{ion}m_{ele}}{e}(\mathbf{J}^{n+1}_{i,k}-\mathbf{J}^{**}_{i,k})+m_{ion}\mathbf{J}^{n+1}_{i,k}\times\mathbf{B}^{n+1}_{i,k}\right),$
(62)

which is the generalized Ohm’s law. At a time scale where the inertial effects
(i.e., cyclotron frequency) are unimportant, the time difference of current is
ignorable. The second term is the Hall current term, which becomes unimportant
in the ideal MHD limit. Let $\nu_{ie}\to 0$ with fully conductive assumption,
then we get the ideal Ohm’s law,

$\mathbf{E}_{i,k}^{n+1}+\mathbf{U}_{i,k}^{n+1}\times\mathbf{B}_{i,k}^{n+1}=0.$
(63)

###  3.6 Summary

The flowchart of the current solution algorithm is summarized in Fig. 2.

Figure 2: Flowchart of solution algorithm.

##  4 Numerical experiments

In this section, we will present some numerical results. The goal of numerical
experiments is not simply to validate the performance of the current scheme,
but also to present and analyze new physical observations. In order to
demonstrate the multi-scale nature of the algorithm, simulations from Vlasov
to magnetohydrodynamics (MHD) regimes are presented. The following
dimensionless flow variables are introduced in the simulations,

$\displaystyle\tilde{\mathbf{x}}=\frac{\mathbf{x}}{L_{0}},\tilde{t}=\frac{t}{L_{0}/U_{0}},\tilde{m}=\frac{m}{m_{ion}},\tilde{n}=\frac{n}{n_{0}},\tilde{\mathbf{U}}=\frac{\mathbf{U}}{U_{0}},\tilde{T}=\frac{T}{T_{0}},$
$\displaystyle\tilde{\mathbf{P}}=\frac{\mathbf{P}}{m_{ion}n_{0}U_{0}^{2}},\tilde{\mathbf{q}}=\frac{\mathbf{q}}{m_{ion}n_{0}U_{0}^{3}},\tilde{f}=\frac{f}{n_{0}U_{0}^{3}},\tilde{\mathbf{u}}=\frac{\mathbf{u}}{U_{0}},$
$\displaystyle\tilde{\mathbf{B}}=\frac{\mathbf{B}}{B_{0}},\tilde{\mathbf{E}}=\frac{\mathbf{E}}{B_{0}U_{0}},\tilde{\sigma}=\frac{\sigma}{en_{0}},\tilde{\mathbf{J}}=\frac{\mathbf{J}}{en_{0}U_{0}},\tilde{\lambda_{D}}=\frac{\lambda_{D}}{r_{g}},$

where $U_{0}=\sqrt{k_{B}T_{0}/m_{ion}}$ is the thermal velocity of ions,
$\lambda_{D}=\sqrt{\varepsilon_{0}k_{B}T_{0}/n_{0}e^{2}}$ is the Debye length,
and $r_{g}=m_{ion}U_{0}/eB_{0}$ is the gyroradius of ion in the reference
state. For brevity, the tilde notation for dimensionless variables will be
removed henceforth. Therefore, the dimensionless BGK-Maxwell system becomes,

$\displaystyle\frac{\partial f_{\alpha}}{\partial
t}+\mathbf{u}\cdot\nabla_{\mathbf{x}}f_{\alpha}+\frac{1}{r_{g}m_{\alpha}}(\mathbf{E}+\mathbf{u}\times\mathbf{B})\cdot\nabla_{\mathbf{u}}f_{\alpha}={\nu_{\alpha}}({\mathcal{M}_{\alpha}-f_{\alpha}}),$
(64) $\displaystyle\frac{\partial\mathbf{E}}{\partial
t}-c^{2}\nabla_{\mathbf{x}}\times\mathbf{B}=-\frac{1}{\lambda_{D}^{2}r_{g}}\mathbf{J},$
$\displaystyle\frac{\partial\mathbf{B}}{\partial
t}+\nabla_{\mathbf{x}}\times\mathbf{E}=0,$
$\displaystyle\frac{1}{\chi}\frac{\partial\phi}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathbf{E}=\frac{\sigma}{\lambda_{D}^{2}r_{g}},$
$\displaystyle\frac{1}{c^{2}\gamma}\frac{\partial\psi}{\partial
t}+\nabla_{\mathbf{x}}\cdot\mathbf{B}=0.$

###  4.1 Landau damping

The Landau damping is a physical phenomenon first predicted by Landau [34]
based on theoretical derivation, namely the exponential decay effect of
electromagnetic waves in collisionless plasmas. Here, we use the example of
Landau damping to verify the performance of the current stochastic scheme in
the Vlasov limit. For brevity, we consider the one-dimensional case first, and
thus the Maxwell’s equations for the electrostatic field degenerates into the
Poisson equation for the electric potential.

####  4.1.1 Linear case

Consider the following macroscopic system

$\left[\begin{array}[]{c}n_{ion}\\\ n_{ele}\\\ U\\\ T\\\
\end{array}\right]_{\text{t=0}}=\left[\begin{array}[]{c}1\\\
1+\alpha\xi\cos(kx)\\\ 0\\\ 1\\\ \end{array}\right],$

where $\alpha$ is the amplitude of electromagnetic wave, $k$ is the wave
number, and $\xi$ is a random parameter. The corresponding particle
distribution functions at $t=0$ are the Maxwellian,

$f_{\alpha}=n_{\alpha}\left(\frac{m_{\alpha}}{2\pi
T}\right)^{1/2}\exp\left(-\frac{m_{\alpha}}{2T}(u-U)^{2}\right).$

The computational setup is listed in Table 1. Given the minor amplitude
$\alpha$, the Vlasov-Poisson system can be regarded as the Maxwellian plus
lineazised perturbation near equilibrium. With the large mass ratio, we fix
the slow motions of ions as background, and solve the evolution of electrons.
The standard non-intrusive stochastic collocation method is also employed to
provide reference solution. Besides, the damping rate of electric field energy
can be derived by linear theory [5]. We hereby combine $\alpha$ and $\xi$ into
a new perturbation amplitude $\beta$, and derive the following relation as
benchmark,

$\mathcal{E}(x,t)\simeq 4\beta\times
0.3677e^{-0.1533t}\sin(kx)\cos(1.4156t-0.536245).$ (65)

Table 1: Computational setup of linear Landau damping. $t$ | $x$ | $N_{x}$ | $m_{ion}/m_{ele}$ | $u_{ele}$ | $N_{u}$ | $\alpha$
---|---|---|---|---|---|---
$(0,40]$ | $[-\pi/k,\pi/k]$ | 128 | $1836$ | $[-5,5]$ | 128 | 0.01
$k$ | $\xi$ | Polynomial |  $N_{z}$(gPC) |  $N_{z}$(quad) | Kn | cfl
0.5 | $\mathcal{U}(0,1)$ | Legendre | $5$ | $9$ | $100$ | 0.2

Fig. 3 presents the time evolution of electric energy. As is shown, the
current stochastic kinetic scheme provides equivalent numerical solution as
the standard collocation method. The expected energy damping rate is
consistent with the theoretical damping rate $-0.1533$, and the oscillation
frequency of the electromagnetic wave corresponds well to the theoretical
value $\omega=1.4156$.

Besides the evolution of expectation value, the stochastic scheme provides the
opportunity to study the uncertainty propagation modes simultaneously. As
shown in Fig. 3(b), the uncertainties travel along with the structure of
electromagnetic wave and present similar evolution patterns with mean field.
The magnitudes of expectation and variance are on the same order due to the
small perturbation strength.

Fig. 4 presents the time evolution of particle distribution function at the
domain center $x=0$. From the zoom-in expected distribution function around
$u=0$, we see that the low-speed particles resonate with the electromagnetic
wave and gradually absorb energy, resulting in decreased value of particle
distribution. In spite of the minor variation of expected value around the
Maxwellian for the linear damping case, the standard deviation of solutions
shows an increasingly symmetric oscillation in the velocity space during the
resonance process. It indicates a more significant sensitivity compared to the
mean field, and provides a quantitative description of the non-equilibrium
effects triggered by the particle-wave resonance.

####  4.1.2 Nonlinear case

The same initial conditions and computational setups as linear case are
followed, except for the enhanced amplitude $\alpha=0.5$. The evolution of
electric energy is presented in Fig. 5. As the amplitude of the
electromagnetic wave increases, nonlinear effects would emerge
correspondingly, resulting in a rise in energy after the initial damping.

Fig. 6 shows the time evolution of expected particle distribution over the
phase space $(x,u)$, and Fig. 7 picks out the expectation and variance of
particle distribution at physical domain center $x=0$. Given the increasing
intensity of radio field, here the particle distribution function is deformed
in phase space. Consistent with the linear case, the change of variance here
is more significant than expectation value, indicating a stronger sensibility.

###  4.2 Two-stream instability

The two-stream instability is another typical phenomenon in for collisionless
plasmas. To a certain extent, it can also be regarded as an inverse phenomenon
of Landau damping.

####  4.2.1 Linear case

Consider the following initial system,

$\left[\begin{array}[]{c}n_{ion}\\\ n_{ele}\\\ U\\\ T\\\ f_{ion}\\\ f_{ele}\\\
\end{array}\right]_{\text{t=0}}=\left[\begin{array}[]{c}1\\\
\frac{2}{7}\left[1+\alpha\xi\left(\frac{\cos(2kx)+\cos(3kx)}{1.2}+\cos(kx)\right)\right]\\\
0\\\ 1\\\ n_{ion}\left(\frac{m_{ion}}{2\pi
T}\right)^{1/2}\exp\left(-\frac{m_{ion}}{2T}(u-U)^{2}\right)\\\
n_{ele}\left(\frac{m_{ele}}{2\pi
T}\right)^{1/2}\left(1+5u^{2}\right)\exp\left(-\frac{m_{ele}}{2T}(u-U)^{2}\right)\\\
\end{array}\right],$

where the computational setup is listed in Table 2. Same as the Landau damping
case, we fix the slow motions of ions as background, and solve the evolution
of electrons.

Table 2: Computational setup of linear two-stream instability. $t$ | $x$ | $N_{x}$ | $m_{ion}/m_{ele}$ | $u_{ele}$ | $N_{u}$ | $\alpha$
---|---|---|---|---|---|---
$(0,70]$ | $[0,2\pi/k]$ | 128 | $1836$ | $[-5,5]$ | 128 | 0.001
$k$ | $\xi$ | Polynomial |  $N_{z}$(gPC) |  $N_{z}$(quad) | Kn | cfl
0.5 | $\mathcal{U}(0,1)$ | Legendre | $5$ | $9$ | $100$ | 0.2

The evolution of electric field energy is shown in Fig. 8, from which it can
be seen that the numerical and theoretical solutions [2] fit well, and the
uncertainties of electric energy propagate in the similar pattern as mean
field. Fig. 9 provides the particle distribution function over the phase space
$(x,u)$ at $t=70$. The swirling pattern is clearly identified in both
expectation and variance values. More fine structures can be seen in the
standard deviation, which provides a clearer way to quantify the stochastic
evolution of particles.

####  4.2.2 Nonlinear case

In the following let us consider the nonlinear case. The initial system is
given as,

$\left[\begin{array}[]{c}n_{ion}\\\ n_{ele}\\\ U\\\ T\\\ f_{ion}\\\ f_{ele}\\\
\end{array}\right]_{\text{t=0}}=\left[\begin{array}[]{c}1\\\
\frac{1}{2}\left(1+\alpha\xi*\cos(kx)\right)\\\ 0\\\ 0.09\\\
n_{ion}\left(\frac{m_{ion}}{2\pi
T}\right)^{1/2}\exp\left(-\frac{m_{ion}}{2T}(u-U)^{2}\right)\\\
\frac{1}{2}n_{ele}\left(\frac{m_{ele}}{2\pi
T}\right)^{1/2}\left(\exp\left(-\frac{m_{ele}}{2T}(u-0.99)^{2}\right)+\exp\left(-\frac{m_{ele}}{2T}(u+0.99)^{2}\right)\right)\\\
\end{array}\right].$

The computational setups adopted here is shown in Table 3.

Table 3: Computational setup of nonlinear two-stream instability. $t$ | $x$ | $N_{x}$ | $m_{ion}/m_{ele}$ | $u_{ele}$ | $N_{u}$ | $\alpha$
---|---|---|---|---|---|---
$(0,70]$ | $[0,2\pi/k]$ | 256 | $1836$ | $[-5,5]$ | 256 | 0.001
$k$ | $\xi$ | Polynomial |  $N_{z}$(gPC) |  $N_{z}$(quad) | Kn | cfl
0.5 | $\mathcal{U}(0,1)$ | Legendre | $5$ | $9$ | $100$ | 0.2

Fig. 10 shows the time evolution of expected particle distribution over the
phase space $(x,u)$. Due to the increasing kinetic energy from two particle
streams, the particle distribution function is stretched and warped in the
phase space. The pattern of standard deviation is similar as mean field, yet
presenting more fine-scale structures, indicating a stronger sensitivity with
respect to stochastic parameters.

###  4.3 Brio-Wu shock tube

After the validation of current scheme in the Vlasov limit, let us turn to
another limiting regime, i.e. the magnetohydrodynamic (MHD) transport problem
under continuum assumption. In this case, we employ the Brio-Wu shock tube [7]
as benchmark case for ideal MHD solutions.

In the simulation, the initial background of deterministic solutions is
consistent with the Sod problem, with an additional magnetic discontinuity,
i.e.,

$\left[\begin{array}[]{c}n_{ion}\\\ n_{ele}\\\ U\\\ p\\\ E_{x}\\\ E_{y}\\\
E_{z}\\\ B_{x}\\\ B_{y}\\\ B_{z}\\\ \phi\\\ \psi\\\
\end{array}\right]_{{L}}=\left[\begin{array}[]{c}1\\\ 1\\\ 0\\\ 1\\\ 0\\\ 0\\\
0\\\ 0.75\\\ 1\\\ 0\\\ 0\\\ 0\\\ \end{array}\right],$

for the left half, and

$\left[\begin{array}[]{c}n_{ion}\\\ n_{ele}\\\ U\\\ p\\\ E_{x}\\\ E_{y}\\\
E_{z}\\\ B_{x}\\\ B_{y}\\\ B_{z}\\\ \phi\\\ \psi\\\
\end{array}\right]_{{R}}=\left[\begin{array}[]{c}0.125\\\ 0.125\\\ 0\\\ 0.8\\\
0\\\ 0\\\ 0\\\ 0.75\\\ -1\\\ 0\\\ 0\\\ 0\\\ \end{array}\right],$

for the right half. The computational setup is listed in Table 4. To recover
the ideal MHD solutions, here the interspecies coefficients $\nu_{ie}$ is set
as zero.

Table 4: Computational setup of Brio-Wu shock tube. $t$ | $x$ | $N_{x}$ | $m_{ion}/m_{ele}$ | $u_{ion}$ | $u_{ele}$ | cfl
---|---|---|---|---|---|---
$(0,0.1]$ | $[0,1]$ | 400 | $1836$ | $[-5,5]$ | $[-5,5]\times\sqrt{\frac{m_{ion}}{m_{ele}}}$ | 0.3
$\xi$ | Polynomial |  $N_{z}$(gPC) |  $N_{z}$(quad) | Kn | $r_{g}$ | $\lambda_{D}$
$\mathcal{U}(0.95,1.05)$ | Legendre | $5$ | $9$ | $1.0\times 10^{-6}$ | $[0.003,100]$ | 0.01

####  4.3.1 Case 1: stochastic flow field

First, we consider the uncertainties of plasma density in the left half tube,
with $n_{ion,L}=n_{ele,L}=\xi$. The numerical expected solutions of
macroscopic variables under different reference gyroradius at $t=0.1$ are
shown in Fig. 11. It is known that the evolution of flow and electromagnetic
fields can be well described by MHD equations when $r_{g}$ is small. At
$r_{g}=0.003$, besides the typical wave structures in the classic Sod problem,
a compound wave emerges around the tube center. As illustrated by Brio and Wu
[7], with the nonconvexity of ideal MHD equation, this wave structure is
induced from initial magnetic discontinuity. From Fig. 11(a)-(c), we see the
current scheme from the BGK-Maxwell system provides equivalent numerical
solutions as ideal MHD equations, which identifies the asymptotic preserving
(AP) property of the current scheme in the MHD limit.

As the gyroradius $r_{g}$ increases, the electric and Lorentz force decreases
accordingly. The effects of charge separation are enhanced, and the omitted
terms, e.g. the Hall current term, become important. The plasma gradually
manifests itself as dielectric material. As a result, the magnetic diffusion
is enhanced, and the wave structures from discontinuous magnetic field $B_{y}$
are flattened. When $r_{g}=100$, at this point the gyroradius is much larger
than the characteristic length of the shock tube, and thus plasma behaves as
neutral gas, and the Brio-Wu shock tube degenerates into a standard Sod case.
In Fig. 11(m) and (n), we see the current scheme provides equivalent solutions
as Euler equations, indicating the AP property of the current scheme in the
Euler limit.

Fig. 12 presents the standard deviations at the same output instant.
Generally, the uncertainties travel along with the wave structure of
expectation values and present similar propagating patterns. Typical wave
structures serve as sources of local maximums of variance. Compared with the
expected value, its variance is more sensitive to physical discontinuities and
holds finer-scale structures due to the spectrum formulation in the random
space. As a result, the overshoots near contact discontinuity and shock, as
well as the oscillations around central compound wave are observed.

Fig. 13 presents the expectations and standard deviations of particle
distribution function at different reference gyrorarius. As is shown, the
discontinuities in macroscopic expectations and overshoots in standard
deviations come from the uncertainties contained in the particle distribution
function near the center of velocity space. From MHD to Euler regimes, the
randomness on particles get reduced and smoothed, resulting in gentle profiles
of macroscopic quantities.

####  4.3.2 Case 2: stochastic magnetic field

In the second case, we turn to the uncertainties from magnetic field in the
left half tube, with $B_{y,L}=\xi$. The numerical solutions of macroscopic
variables under different reference gyroradius at $t=0.1$ are shown in Fig. 14
and 15. For the expected values in Fig. 14, we see that different initial
stochastic conditions don’t significantly change the wave structures of mean
field. Both MHD and Euler limits are precisely preserved under stochastic
magnetic field.

Fig. 15 presents the corresponding standard deviations of macroscopic system.
As can be seen, the uncertainties manifest the consistent propagating patterns
of expected values. Different from the previous case, here the high-frequency
variance of density in compound wave region is much reduced, which indicates
this phenomenon comes from the uncertain densities. As $r_{g}$ increases, the
wave structures around tube center becomes even more complicated in the
transition regime $r_{g}=0.01$, and simplify again when it comes to
$r_{g}=0.1$, which demonstrates the complex nonlinearity from Hall current,
two-fluid effects, etc. The increasing gyroradius doesn’t lead to a simple
process of monotonic variation.

This test case clearly shows the consistency and distinction of propagation
modes between expectation value and variance. It also illustrates the capacity
of current scheme to simulate multi-scale and multi-physics plasma transports,
and capture the propagation of uncertainties in different regimes.

##  5 Conclusion

Plasma dynamics is associated with an intrinsic multi-scale nature due to the
large variations of particle density and temperature, as well as the
characteristic scales of the local structures. Based on the multi-component
BGK model and Maxwell’s equations, a stochastic kinetic scheme with hybrid
Galerkin-collocation strategy has been constructed in this paper, which allows
for a unified numerical simulation for multi-scale plasma physics. Based on
the cross-scale modeling, the solution algorithm is able to capture both
equilibrium magnetohyrodynamics and non-equilibrium gyrations of charged
particles simultaneously, and recover the scale-dependent plasma physics along
with the emergence, propagation, and evolution of randomness. The asymptotic-
preserving property of the scheme is validated through theoretical analysis
and numerical tests.

## Acknowledgement

The current research is funded by the Alexander von Humboldt Foundation.

## References

  * [1] National Research Council, Plasma Science Committee, Plasma 2010 Committee, et al.  Plasma science: advancing knowledge in the national interest, volume 3.  National Academies Press, 2008.
  * [2] Francis F Chen.  Introduction to plasma physics.  Springer Science & Business Media, 2012.
  * [3] Vahid Vahedi and Maheswaran Surendra.  A monte carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges.  Computer Physics Communications, 87(1-2):179–198, 1995.
  * [4] Francis Filbet, Eric Sonnendrücker, and Pierre Bertrand.  Conservative numerical schemes for the vlasov equation.  Journal of Computational Physics, 172(1):166–187, 2001.
  * [5] Nicolas Crouseilles, Thomas Respaud, and Eric Sonnendrücker.  A forward semi-lagrangian method for the numerical solution of the vlasov equation.  Computer Physics Communications, 180(10):1730–1745, 2009.
  * [6] Jing-Mei Qiu and Andrew Christlieb.  A conservative high order semi-lagrangian weno method for the vlasov equation.  Journal of Computational Physics, 229(4):1130–1149, 2010.
  * [7] Moysey Brio and Cheng Chin Wu.  An upwind differencing scheme for the equations of ideal magnetohydrodynamics.  Journal of computational physics, 75(2):400–422, 1988.
  * [8] Kenneth G Powell, Philip L Roe, Timur J Linde, Tamas I Gombosi, and Darren L De Zeeuw.  A solution-adaptive upwind scheme for ideal magnetohydrodynamics.  Journal of Computational Physics, 154(2):284–309, 1999.
  * [9] Pierre Degond, Fabrice Deluzet, Laurent Navoret, An-Bang Sun, and Marie-Hélène Vignal.  Asymptotic-preserving particle-in-cell method for the vlasov–poisson system near quasineutrality.  Journal of Computational Physics, 229(16):5630–5652, 2010.
  * [10] Pierre Degond and Fabrice Deluzet.  Asymptotic-preserving methods and multiscale models for plasma physics.  Journal of Computational Physics, 336:429–457, 2017.
  * [11] National Institute of Standards and Technology.  The NIST Reference on Constants, Units, and Uncertainty, 5 2019\.
  * [12] Dongbin Xiu.  Numerical methods for stochastic computations: a spectral method approach.  Princeton university press, 2010.
  * [13] Claudio Canuto and Alfio Quarteroni.  Approximation results for orthogonal polynomials in sobolev spaces.  Mathematics of Computation, 38(157):67–86, 1982.
  * [14] Dongbin Xiu and Jan S Hesthaven.  High-order collocation methods for differential equations with random inputs.  SIAM Journal on Scientific Computing, 27(3):1118–1139, 2005.
  * [15] Jingwei Hu, Shi Jin, and Ruiwen Shu.  A stochastic galerkin method for the fokker–planck–landau equation with random uncertainties.  In XVI International Conference on Hyperbolic Problems: Theory, Numerics, Applications, pages 1–19. Springer, 2016.
  * [16] Jingwei Hu and Shi Jin.  Uncertainty quantification for kinetic equations.  In Uncertainty Quantification for Hyperbolic and Kinetic Equations, pages 193–229. Springer, 2017.
  * [17] Shi Jin and Yuhua Zhu.  Hypocoercivity and uniform regularity for the vlasov–poisson–fokker–planck system with uncertainty and multiple scales.  SIAM Journal on Mathematical Analysis, 50(2):1790–1816, 2018.
  * [18] Zhiyan Ding and Shi Jin.  Random regularity of a nonlinear landau damping solution for the vlasov-poisson equations with random inputs.  International Journal for Uncertainty Quantification, 9(2), 2019\.
  * [19] Edward G Phillips and Howard C Elman.  A stochastic approach to uncertainty in the equations of mhd kinematics.  Journal of Computational Physics, 284:334–350, 2015.
  * [20] Kazuo Yamazaki.  Stochastic hall-magneto-hydrodynamics system in three and two and a half dimensions.  Journal of Statistical Physics, 166(2):368–397, 2017.
  * [21] Tianbai Xiao and Martin Frank.  A stochastic kinetic scheme for multi-scale flow transport with uncertainty quantification.  arXiv preprint arXiv:2002.00277, 2020.
  * [22] Pierre Andries, Kazuo Aoki, and Benoit Perthame.  A consistent bgk-type model for gas mixtures.  Journal of Statistical Physics, 106(5):993–1018, 2002.
  * [23] TF Morse.  Energy and momentum exchange between nonequipartition gases.  The Physics of Fluids, 6(10):1420–1427, 1963.
  * [24] SI Braginskii.  Transport processes in a plasma.  Reviews of plasma physics, 1, 1965.
  * [25] C-D Munz, Pascal Omnes, Rudolf Schneider, Eric Sonnendrücker, and Ursula Voss.  Divergence correction techniques for maxwell solvers based on a hyperbolic model.  Journal of Computational Physics, 161(2):484–511, 2000.
  * [26] Mikhail N Kogan.  Rarefied Gas Dynamics.  Plenum Press, 1969.
  * [27] Kun Xu and Juan-Chen Huang.  A unified gas-kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, 229(20):7747–7764, 2010.
  * [28] Tianbai Xiao, Qingdong Cai, and Kun Xu.  A well-balanced unified gas-kinetic scheme for multiscale flow transport under gravitational field.  Journal of Computational Physics, 332:475 – 491, 2017.
  * [29] Chang Liu and Kun Xu.  A unified gas kinetic scheme for continuum and rarefied flows v: Multiscale and multi-component plasma transport.  Communications in Computational Physics, 22(5):1175–1223, 2017\.
  * [30] Tianbai Xiao, Kun Xu, and Qingdong Cai.  A unified gas-kinetic scheme for multiscale and multicomponent flow transport.  Applied Mathematics and Mechanics, pages 1–18, 2019.
  * [31] Tianbai Xiao, Chang Liu, Kun Xu, and Qingdong Cai.  A velocity-space adaptive unified gas kinetic scheme for continuum and rarefied flows.  Journal of Computational Physics, page 109535, 2020.
  * [32] Ammar Hakim, John Loverich, and Uri Shumlak.  A high resolution wave propagation scheme for ideal two-fluid plasma equations.  Journal of Computational Physics, 219(1):418–442, 2006.
  * [33] Lawrence F Shampine.  Implementation of rosenbrock methods.  ACM Transactions on Mathematical Software (TOMS), 8(2):93–113, 1982\.
  * [34] L. Landau.  On the vibration of the electronic plasma.  Journal of Physics, 10(1):25 – 34, 1946.

Table 5: Nomenclature of stochastic kinetic scheme. $k_{B}$ | Boltzmann constant
---|---
$\varepsilon_{0}$ | vacuum permittivity
$\mu_{0}$ | vacuum permeability
$e$ | elementary charge
$c$ | speed of light
$q_{\alpha}$ | charge of species $\alpha$, with $q_{\alpha}=\pm e$ ($\alpha=ion,\ ele$)
$m_{\alpha}$ | particle mass of species $\alpha$
$n_{\alpha}$ | number density of species $\alpha$
$\rho_{\alpha}$ | density of species $\alpha$
$\mathbf{U}_{\alpha}$ | macroscopic velocity of species $\alpha$
$T_{\alpha}$ | temperature of species $\alpha$
$\mathbf{P}_{\alpha}$ | stress tensor of species $\alpha$
$p_{\alpha}$ | pressure of species $\alpha$
$\mathbf{q}_{\alpha}$ | heat flux of species $\alpha$
$\mathbf{W}_{\alpha}$ | macroscopic conservative variables of species $\alpha$ (density, momentum, energy)
$f_{\alpha}$ | particle distribution function of species $\alpha$
$\mathbf{u}$ | particle velocity
$\mathbf{a}_{\alpha}$ | electromagnetic force acting on species $\alpha$
$Q_{\alpha}$ | kinetic collision operator of species $\alpha$
$\mathcal{M}_{\alpha}$ | equilibrium distribution function of species $\alpha$
$\mathbf{E}$ | electric field
$\mathbf{B}$ | magnetic field
$\mathbf{\phi}$ | correction potential for $\mathbf{E}$ in perfectly hyperbolic Maxwell’s equations
$\mathbf{\psi}$ | correction potential for $\mathbf{B}$ in perfectly hyperbolic Maxwell’s equations
$\sigma$ | charge density
$\mathbf{J}$ | current density
$\mathbf{M}$ | abbreviation of electromagnetic variables ($\mathbf{E},\mathbf{B},\phi,\psi$)
$\varpi$ | vector of collision invariants
$\nu_{\alpha}$ | collision frequency of species $\alpha$
$\nu_{ie}$ | interaction frequency between ion and electron
$\tau_{\alpha}$ | Collision time of species $\alpha$ with $\tau_{\alpha}=1/\nu_{\alpha}$
$\mathbf{F}^{W}_{\alpha}$ | flux for macroscopic conservative variables of species $\alpha$
$\mathbf{F}^{f}_{\alpha}$ | flux for particle distribution function of species $\alpha$
$\mathbf{F}^{M}$ | flux for electromagnetic variables
${L}_{N}$ | generalized polynomial chaos expansion of any stochastic variable $L$ with degree $N$
$\hat{\bm{l}}$ | vector of generalized polynomial chaos coefficients of any stochastic variable $L$
$\Phi$ | orthogonal polynomials in random space $\mathbf{z}$
$\varrho$ | probability density function of random variable $\mathbf{z}$
$\rm Kn$ | Knudsen number
$r_{L}$ | gyroradius
$\lambda_{D}$ | Debye length
$H[x]$ | Heaviside step function

(a) Expectation value (b) Standard deviation

Figure 3: Expectation value and standard deviation of electric field energy
(logarithmic) in linear landau damping.

(a) Expectation value (b) Standard deviation

Figure 4: Time evolved expectation value and standard deviation of particle
distribution function at $x=0$ in linear landau damping.

(a) Expectation value (b) Standard deviation

Figure 5: Expectation value and standard deviation of electric field energy
(logarithmic) in nonlinear landau damping.

(a) $t=0$ (b) $t=5$ (c) $t=10$ (d) $t=20$

(e) $t=30$ (f) $t=50$

Figure 6: Time evolved expectation value of particle distribution function
over phase space $(x,u)$ in nonlinear landau damping.

(a) Expectation value (b) Standard deviation

Figure 7: Time evolved expectation value and standard deviation of particle
distribution function at $x=0$ in nonlinear landau damping.

(a) Expectation value (b) Standard deviation

Figure 8: Expectation value and standard deviation of electric field energy
(logarithmic) in linear two-stream instability.

(a) Expectation value (b) Standard deviation

Figure 9: Expectation value and standard deviation of particle distribution
function over phase space $(x,u)$ in linear two-stream instability.

(a) Expectation value (b) Standard deviation

Figure 10: Expectation value and standard deviation of particle distribution
function over phase space $(x,u)$ in nonlinear two-stream instability.

(a) $\mathbb{E}(N)$ (b) $\mathbb{E}(U)$ (c) $\mathbb{E}(B_{y})$ (d)
$\mathbb{E}(N)$

(e) $\mathbb{E}(U)$ (f) $\mathbb{E}(B_{y})$ (g) $\mathbb{E}(N)$ (h)
$\mathbb{E}(U)$

(i) $\mathbb{E}(B_{y})$ (j) $\mathbb{E}(N)$ (k) $\mathbb{E}(U)$ (l)
$\mathbb{E}(B_{y})$

(m) $\mathbb{E}(N)$ (n) $\mathbb{E}(U)$ (o) $\mathbb{E}(B_{y})$

Figure 11: Expectation values of $N$, $U$ and $B_{y}$ in Brio-Wu shock tube
with density uncertainty at $t=0.1$ (row 1: $r_{g}=0.003$, row 2:
$r_{g}=0.01$, row 3: $r_{g}=0.1$, row 4: $r_{g}=1$, row 5: $r_{g}=100$).

(a) $\mathbb{S}(N)$ (b) $\mathbb{S}(U)$ (c) $\mathbb{S}(B_{y})$ (d)
$\mathbb{S}(N)$

(e) $\mathbb{S}(U)$ (f) $\mathbb{S}(B_{y})$ (g) $\mathbb{S}(N)$ (h)
$\mathbb{S}(U)$

(i) $\mathbb{S}(B_{y})$ (j) $\mathbb{S}(N)$ (k) $\mathbb{S}(U)$ (l)
$\mathbb{S}(B_{y})$

(m) $\mathbb{S}(N)$ (n) $\mathbb{S}(U)$ (o) $\mathbb{S}(B_{y})$

Figure 12: Standard deviations of $N$, $U$ and $B_{y}$ in Brio-Wu shock tube
with density uncertainty at $t=0.1$ (row 1: $r_{g}=0.003$, row 2:
$r_{g}=0.01$, row 3: $r_{g}=0.1$, row 4: $r_{g}=1$, row 5: $r_{g}=100$).

(a) $\mathbb{E}(h_{ion})$ (b) $\mathbb{S}(h_{ion})$ (c) $\mathbb{E}(h_{ion})$
(d) $\mathbb{S}(h_{ion})$

(e) $\mathbb{E}(h_{ion})$ (f) $\mathbb{S}(h_{ion})$ (g) $\mathbb{E}(h_{ion})$
(h) $\mathbb{S}(h_{ion})$

Figure 13: Expectations and variances of reduced ion distribution $h_{ion}$,
in Brio-Wu shock tube with density uncertainty at $t=0.1$ (row 1:
$r_{g}=0.003$, row 2: $r_{g}=0.01$, row 3: $r_{g}=0.1$, row 4: $r_{g}=1$).

(a) $\mathbb{E}(N)$ (b) $\mathbb{E}(U)$ (c) $\mathbb{E}(B_{y})$ (d)
$\mathbb{E}(N)$

(e) $\mathbb{E}(U)$ (f) $\mathbb{E}(B_{y})$ (g) $\mathbb{E}(N)$ (h)
$\mathbb{E}(U)$

(i) $\mathbb{E}(B_{y})$ (j) $\mathbb{E}(N)$ (k) $\mathbb{E}(U)$ (l)
$\mathbb{E}(B_{y})$

(m) $\mathbb{E}(N)$ (n) $\mathbb{E}(U)$ (o) $\mathbb{E}(B_{y})$

Figure 14: Expectation values of $N$, $U$ and $B_{y}$ in Brio-Wu shock tube
with magnetic uncertainty at $t=0.1$ (row 1: $r_{g}=0.003$, row 2:
$r_{g}=0.01$, row 3: $r_{g}=0.1$, row 4: $r_{g}=1$, row 5: $r_{g}=100$).

(a) $\mathbb{S}(N)$ (b) $\mathbb{S}(U)$ (c) $\mathbb{S}(B_{y})$ (d)
$\mathbb{S}(N)$

(e) $\mathbb{S}(U)$ (f) $\mathbb{S}(B_{y})$ (g) $\mathbb{S}(N)$ (h)
$\mathbb{S}(U)$

(i) $\mathbb{S}(B_{y})$ (j) $\mathbb{S}(N)$ (k) $\mathbb{S}(U)$ (l)
$\mathbb{S}(B_{y})$

(m) $\mathbb{S}(N)$ (n) $\mathbb{S}(U)$ (o) $\mathbb{S}(B_{y})$

Figure 15: Standard deviations of $N$, $U$ and $B_{y}$ in Brio-Wu shock tube
with magnetic uncertainty at $t=0.1$ (row 1: $r_{g}=0.003$, row 2:
$r_{g}=0.01$, row 3: $r_{g}=0.1$, row 4: $r_{g}=1$, row 5: $r_{g}=100$).