**affiliationtext: Oden Institute for Computational Sciences and Engineering,
University of Texas at Austin${\ddagger}$${\ddagger}$affiliationtext:
Institute for Fusion Studies, University of Texas at
Austin${\dagger}$${\dagger}$affiliationtext: Department of Mechanical
Engineering, Eindhoven University of Technology

# A conservative Galerkin solver for the quasilinear diffusion model in
magnetized plasmas

Kun Huang  Michael Abdelmalik  Boris Breizman  Irene M. Gamba

###### Abstract

The quasilinear theory describes the resonant interaction between particles
and waves with two coupled equations: one for the evolution of the particle
probability density function(pdf), the other for the wave spectral energy
density(sed). In this paper, we propose a conservative Galerkin scheme for the
quasilinear model in three-dimensional momentum space and three-dimensional
spectral space, with cylindrical symmetry.

We construct an unconditionally conservative weak form, and propose a
discretization that preserves the unconditional conservation property, by
”unconditional” we mean that conservation is independent of the singular
transition probability. The discrete operators, combined with a consistent
quadrature rule, will preserve all the conservation laws rigorously. The
technique we propose is quite general: it works for both relativistic and non-
relativistic systems, for both magnetized and unmagnetized plasmas, and even
for problems with time-dependent dispersion relations.

We represent the particle pdf by continuous basis functions, and use
discontinuous basis functions for the wave sed, thus enabling the application
of a positivity-preserving technique. The marching simplex algorithm, which
was initially designed for computer graphics, is adopted for numerical
integration on the resonance manifold. We introduce a semi-implicit time
discretization, and discuss the stability condition. In addition, we present
numerical examples with a ”bump on tail” initial configuration, showing that
the particle-wave interaction results in a strong anisotropic diffusion effect
on the particle pdf.

Keywords: computational kinetic systems, quasilinear theory for magnetized
plasmas, mean-field effect, weak turbulence model.

##  1 Introduction

The Vlasov-Maxwell system and the Vlasov-Poisson system are widely used to
describe the collective(mean-field) effect of particles. Although a lot of
work has been done in numerical methods for these systems, in practice a
reduced model is often preferred when the problem is in high dimension and
some loss of details is justified from physics consideration.

For example, the electron runaway problem, which is the motivation of our
research, has attracted a lot of attention as a risk factor for magnetic
confinement fusion reactors like ITER [14, 4]. Runaway electrons are a group
of extremely fast electrons generated inside the tokamak, the release of which
can damage the wall. Therefore, it is important to have an answer to the
questions like how they are generated and how to mitigate them.

The dynamics of runaway electrons is determined by external electromagnetic
fields, collision, and particle-wave interaction. The quasilinear theory, as a
reduced mean-field model governing particle-wave interaction, arises from
averaging and linearizing over the original Vlasov-Maxwell system in weak
turbulence regime.

The quasilinear theory for unmagnetized plasmas was proposed by Vedenov et al.
[19] and Drummond et al. [6]. It was later generalized by Shapiro et al. [15]
to model the magnetized plasma. The same idea has been used extensively in the
following years, for example in the work of Kennel [10], Lerche [11], and
Kaufman [9], etc. The validity of such a model reduction was studied
numerically by Besse et al. [3], and analytically by Bardos and Besse [2] for
Vlasov-Poisson system. The existence of global weak solutions in one
dimensional electrostatic case has been proved in [8].

Since the quasilinear theory studies the spectrum of waves and the averaged
particle distribution function, it does not require a small time step to
characterize the high wave frequency numerically. However, the numerical
computation of the particle-wave resonance system is still challenging, due to
the resonance condition described with the Dirac delta function, the
complicated dispersion relation, high dimension, nonlinearity, and
conservation laws consisting of integrals in two different spaces. Therefore,
although the theory has been widely used in physics, there is no preceding
work focusing on the numerical method for quasilinear theory in magnetized
plasmas.

In this paper, we propose a conservative Galerkin solver for the homogeneous
quasilinear particle-wave interaction system.

Despite being a paradigm approach in the analysis and discretization of other
kinetic equations, the weak formulation of the quasilinear model has not
gained enough attention, partly because the equation for particles was usually
written in a nonlinear diffusion form, and the equation for waves was treated
as independent first-order ODEs with parameters. There are infinitely many
equivalent forms to the same equation because of the resonance condition.
Among all the equivalent forms, some are superior to the others, the reason is
as follows.

The quasilinear theory inherits the conservation laws from the original
Vlasov-Maxwell system. However, generally the conservation is conditional,
which means the gain and loss parts only offset each other on the resonance
manifold. When the resonance manifold is broadened or approximated,
conservation laws are no longer guaranteed. In this paper, we propose a novel
integro-differential form and the corresponding unconditionally conservative
weak form.

It is desired that the discrete weak form will preserve the unconditional
conservation property above, unfortunately, naive standard finite element
discretizations turn out to fail. We located the cause of discretization
errors by analyzing the weak form, and managed to construct a perfect
discretization by replacing some quantities with their projection in the
discrete finite element spaces.

Apart from that, for numerical integration on resonance manifold, we adopt the
marching simplex algorithm[5, 12], which enables us to deal with arbitrary
wave modes.

This paper is organized as follows. In section 2 we review the relativistic
quasilinear model for magnetized plasmas and introduce the integro-
differential system with its weak form. The conservative semi-discrete system,
as the main result of this paper, will be presented in section 3. In section
4, we derive the nonlinear ODE system associated with our conservative semi-
discrete form, and the relation between two interaction tensors is proved.
Stability and positivity will be discussed in section 5. The numerical results
are presented in section 6.

##  2 The Quasilinear Particle-Wave Interaction System

The quasilinear particle-wave interaction system consists of a diffusion
equation for electron pdf(probability distribution function) and a reaction
equation for wave sed(spectral energy density). They couple with each other
through the coefficients. As a reduced model for the Vlasov-Maxwell system,
the quasilinear theory inherits the conservation properties: mass, momentum,
and energy. Moreover, the entropy of particles dissipates as a result of
diffusion.

In this section, we are going to show that the system can be written in a
novel integro-differential form, which will finally lead to a conservative
discrete scheme.

###  2.1 The Integro-Differential Strong Form

There are two ways to interpret the quasilinear theory. In classical language,
it is a model reduction of the Vlasov-Maxwell system in a weak turbulence
regime. Meanwhile, in quantum mechanical language, the waves in a plasma can
also be regarded as a group of plasmons (wave packets, analogous to photons).
Hence the interaction between particles and waves can be interpreted as a
stochastic process of particles emitting/absorbing plasmons. For a derivation
in quantum mechanical language, we refer the readers to the review paper of
Vedenov[18] and the book of Thorne et al.[17]. In what follows, we will rely
more on the latter interpretation.

The particle pdf is a non-negative function defined on the particle momentum
space $\mathbb{R}_{p}^{3}$, i.e.

$f(\mathbf{p},t):\mathbb{R}_{p}^{3}\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}.$

Each particle carries mass $m$, momentum $\mathbf{p}$, and energy
$E(\mathbf{p})$. The particle velocity is
$\mathbf{v}(\mathbf{p})=\nabla_{\mathbf{p}}E(\mathbf{p})$.

At the same time, the wave sed is a non-negative function defined on the wave
momentum space $\mathbb{R}_{k}^{3}$ (the spectral space), i.e.

$W(\mathbf{k},t):\mathbb{R}_{k}^{3}\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}.$

Each plasmon carries momentum $\hbar\mathbf{k}$, energy $\hbar\omega$ and no
mass. The relation between wave frequency $\omega$ and wave-vector
$\mathbf{k}$ is called the dispersion relation, $\omega:=\omega(\mathbf{k})$.
The wave spectral energy density can be expressed in terms of plasmon number
density, $W(\mathbf{k})=N(\mathbf{k})\hbar\omega(\mathbf{k})$. Therefore it
makes no difference whether we use $W(\mathbf{k})$ or $N(\mathbf{k})$ as the
unknown.

The equations for quasilinear particle-wave interaction share the following
structure,

$\begin{split}\frac{\partial}{\partial
t}f(\mathbf{p},t)&=\nabla_{p}\cdot\left(D[W](\mathbf{p},t)\cdot\nabla_{p}f(\mathbf{p},t)\right),\\\
\frac{\partial}{\partial
t}W(\mathbf{k},t)&=\Gamma[f](\mathbf{k},t)W(\mathbf{k},t).\end{split}$ (2.1)

Both relations,
$D\\!:\\!L^{1}(\mathbb{R}_{k}^{3})\rightarrow(L^{\infty}(\mathbb{R}_{p}^{3}))^{3\times
3}$ and
$\Gamma\\!:\\!H^{1}(\mathbb{R}_{p}^{3})\\!\rightarrow\\!L^{\infty}(\mathbb{R}_{k}^{3})$,
are determined by transition probabilities of the stochastic
emission/absorption process. The transition probabilities per se, depend
solely on pre-interaction and post-interaction kinetic variables: particle
momentum $\mathbf{p}$, particle energy $E(\mathbf{p})$, plasmon momentum
$\hbar\mathbf{k}$ and plasmon energy $\hbar\omega(\mathbf{k})$. Hence the
particle energy relation $E(\mathbf{p})$ and plasmon dispersion relation
$\omega(\mathbf{k})$ must be specified before numerical simulation.

######  Remark 1.

The wave dispersion relation $\omega(\mathbf{k})$ depends on the medium, i.e.
the plasma itself, which is evolving. Since the computational cost for an
accurate dispersion relation $\omega:=\omega(\mathbf{k})$ can be quite high,
there is, in practice, a tendency to use low-order approximations based on
appropriate assumptions, for example, the cold plasma assumption (see
Appendix).

######  Remark 2.

In a plasma, there can be multiple wave modes, i.e. multiple ”species” of
plasmons, each having a distinct dispersion relation $\omega(\mathbf{k})$. In
our numerical experiment, we use the dispersion relation of whistler waves in
a cold magnetized plasma. Nevertheless, our numerical method is compatible
with any dispersion relation, and can be used to simulate multiple wave modes
at the same time.

In the previous paragraphs, we introduced the general governing equations for
particles absorbing/emitting plasmons. Next, let us focus on a specific
example that the proposed method is designed for.

#### The Particle-Wave Interaction System for Magnetized Plasmas with
Cylindrical Symmetry

We are interested in particle-wave interaction for plasmas embedded in the
magnetic field $\mathbf{B}(t)=\mathbf{B}_{0}$, as a consequence, we will focus
on gyro-averaged particle distribution functions $f(\mathbf{p},t)$ having
cylindrical symmetry. For simplicity, it is further assumed that
$W(\mathbf{k},t)$ is cylindrically symmetric.

The background magnetic field $\mathbf{B}_{0}$ induces an axis direction
$\mathbf{e}_{\parallel}=\mathbf{B}_{0}/|\mathbf{B}_{0}|$ and the associated
cylindrical coordinates, where
$u_{\parallel}=\mathbf{u}\cdot\mathbf{e}_{\parallel}$ and
$u_{\perp}=({u^{2}-u_{\parallel}^{2}})^{1/2}$ for any vector
$\mathbf{u}\in\mathbb{R}^{3}$. Then we have
$f(\mathbf{p},t)=f(p_{\parallel},p_{\perp},t)$ and
$W(\mathbf{k},t)=W(k_{\parallel},k_{\perp},t)$.

Particles with charge $q$ and mass $m$ has gyro-frequency
$\omega_{c}=qB_{0}/mc$. Note that for electrons, $q=-e$. Particles with
momentum $\mathbf{p}$ has relativistic energy
$E(\mathbf{p})=\gamma(\mathbf{p})mc^{2}$, where $\gamma$ is the Lorentz
factor.

To model the runaway electrons in a tokamak, Breizman et al.[4] used the
following equations written in spherical coordinates(i.e.
$\theta=\arccos{p_{\parallel}/p}$),

$\begin{split}\frac{\partial f(\mathbf{p},t)}{\partial
t}&=\frac{1}{p^{2}}\frac{\partial}{\partial p}p^{2}(D_{pp}\frac{\partial
f}{\partial p}+D_{p\theta}\frac{\partial
f}{p\partial\theta})+\frac{1}{p\sin\theta}\frac{\partial}{\partial\theta}\sin\theta(D_{\theta
p}\frac{\partial f}{\partial p}+D_{\theta\theta}\frac{\partial
f}{p\partial\theta}),\\\ \frac{\partial W(\mathbf{k},t)}{\partial
t}&=\Gamma(\mathbf{k},t)W,\end{split}$ (2.2)

where the diffusion coefficients are weighted integrals of wave sed
$W(\mathbf{k},t)$,

$\begin{split}D_{pp}(\mathbf{p},t)&=\sum_{l=-\infty}^{+\infty}\int
d^{3}k\left\\{W(\mathbf{k},t)U_{l}(\mathbf{p},\mathbf{k})\delta(\omega(\mathbf{k})-k_{\parallel}v(\mathbf{p})\cos\theta-l\omega_{c}/\gamma)\right\\},\\\
D_{p\theta}(\mathbf{p},t)=D_{p\theta}(\mathbf{p},t)&=\sum_{l=-\infty}^{+\infty}\int
d^{3}k\left\\{\left(\frac{l\omega_{c}/\gamma-\omega\sin^{2}\theta}{\omega\sin\theta\cos\theta}\right)W(\mathbf{k},t)U_{l}(\mathbf{p},\mathbf{k})\delta(\omega-
k_{\parallel}v\cos\theta-l\omega_{c}/\gamma)\right\\},\\\
D_{\theta\theta}(\mathbf{p},t)&=\sum_{l=-\infty}^{+\infty}\int
d^{3}k\left\\{\left(\frac{l\omega_{c}/\gamma-\omega\sin^{2}\theta}{\omega\sin\theta\cos\theta}\right)^{2}W(\mathbf{k},t)U_{l}(\mathbf{p},\mathbf{k})\delta(\omega-
k_{\parallel}v\cos\theta-l\omega_{c}/\gamma)\right\\},\end{split}$

and the growth rate is a weighted integral of particle pdf’s gradient,
$\nabla_{\mathbf{p}}f(\mathbf{p},t)$,

$\Gamma(\mathbf{k},t)=\sum_{l=-\infty}^{+\infty}\int
d^{3}p\left\\{v(\mathbf{p})(\frac{\partial f}{\partial
p}+\frac{l\omega_{c}/\gamma-\omega\sin^{2}\theta}{\omega(\bf
k)\sin\theta\cos\theta}\frac{\partial
f}{p\partial\theta})U_{l}(\mathbf{p},\mathbf{k})\delta(\omega-
k_{\parallel}v\cos\theta-l\omega_{c}/\gamma)\right\\}.$

The coefficients characterizing transition probability, is a function of
$\mathbf{p}$ and $\mathbf{k}$ that always takes finite non-negative values,

$U_{l}(\mathbf{p};\mathbf{k})=8\pi^{2}e^{2}\frac{\left\\{\frac{l\omega_{c}}{k_{\perp}p}J_{l}+E_{3}\cos\theta
J_{l}+iE_{2}\sin\theta
J_{l}^{\prime}\right\\}^{2}}{\left(1-E_{2}^{2}\right)\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}\varepsilon)+2iE_{2}\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}g)+E_{3}^{2}\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}\eta)},$
(2.3)

where the dielectric tensor components $\varepsilon,g$, $\eta$ and the wave
polarization vector components $E_{j}$ are functions of wave frequency
$\omega$ as defined in Equation(8.3) and Equation(8.4), argument of the Bessel
functions is $k_{\perp}p_{\perp}/m\omega_{c}$, for details, see Breizman et.
al.[4].

The above formulas are sufficient to perform a trivial numerical simulation:
treat Equation(2.1) as a normal diffusion equation and a normal reaction
equation with time-varying coefficients. The challenging part is to preserve
conservation, especially when there are integrals containing the Dirac delta
function. Numerical integrals are always performed by quadrature rules,
however, the quadrature points usually do not reside exactly on the resonance
manifold. In this paper, we propose an unconditionally conservative approach
by employing a novel equivalent integro-differential form instead of the
original equation. By ”unconditional” we mean that the scheme is conservative
no matter how the resonance manifold is discretized or broadened.

#### The Emission/Absorption Kernel and Directional Differential Operator

To rephrase Equation(2.2) in integral-differential form, let us introduce two
important concepts, the emission/absorption kernel and the directional
differential operator, along with some necessary notations.

Particles with momentum $\mathbf{p}$ do not emit or absorb plasmons with wave
vector $\mathbf{k}$ unless a certain resonance condition is satisfied. Define
the $l$-th resonance indicator function
$s_{l}(\mathbf{p},\mathbf{k})\coloneqq\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p})$,
$l\in\mathbb{Z}$, then the resonance condition reads
$s_{l}(\mathbf{p},\mathbf{k})=0$. Define the $l$-th resonance manifold as

$\mathcal{S}_{l}=\left\\{(\mathbf{p},\mathbf{k})\in\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}:s_{l}(\mathbf{p},\mathbf{k})=\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p})=0\right\\},$
(2.4)

then we can say that particles with momentum $\mathbf{p}$ emit or absorb
plasmons with wave vector $\mathbf{k}$ only when $(\mathbf{p},\mathbf{k})$
belongs to one of the resonance manifolds.

Analogous to the definition of collisional kernels in Boltzmann equations and
Fokker-Planck-Landau equations, we define the emission/absorption kernel which
characterizes the probability for a particle with momentum $\mathbf{p}$ to
absorb or emit a plasmon with wave vector $\mathbf{k}$,

$\mathcal{B}(\mathbf{p},\mathbf{k})=\sum_{l=-\infty}^{+\infty}U_{l}(\mathbf{p};\mathbf{k})\delta(s_{l}(\mathbf{p};\mathbf{k})).$
(2.5)

As we have mentioned above, interaction happens only if the resonance
condition is satisfied, so the emission/absorption kernel contains a Dirac
delta function. The coefficients $U_{l}(\mathbf{p},\mathbf{k})$ are given in
Equation(2.3). They take finite non-negative values for any coordinates
$(\mathbf{p},\mathbf{k})$.

Interaction with a plasmon results in diffusion of particle pdf
$f(\mathbf{p},t)$ along a particular direction
$\mathbf{\beta}(\mathbf{p},\mathbf{k})$, thus we define the directional
differential operator

$\mathcal{L}(g(\mathbf{p}))\coloneqq\mathbf{\beta}(\mathbf{p},\mathbf{k})\cdot\nabla_{p}g(\mathbf{p})=\frac{k_{\parallel}v_{\parallel}}{\omega}\frac{p}{p_{\parallel}}\frac{\partial
g}{\partial
p_{\parallel}}+(1-\frac{k_{\parallel}v_{\parallel}}{\omega})\frac{p}{p_{\perp}}\frac{\partial
g}{\partial p_{\perp}}$ (2.6)

Further, define the $L^{2}$ inner product in particle momentum space
$(u(\mathbf{p}),v(\mathbf{p}))_{p}\coloneqq\int_{\mathbb{R}_{p}^{3}}uvd^{3}p$,
and the $L^{2}$ inner product in wave spectral space
$(U(\mathbf{k}),V(\mathbf{k}))_{k}=\int_{\mathbb{R}_{k}^{3}}UVd^{3}k$.

Denote the adjoint operator of $\mathcal{L}$ by $\mathcal{L}^{*}$, then by
definition, we have

$(\mathcal{L}^{*}u,v)_{p}=(u,\mathcal{L}v)_{p}.$

#### The Bilinear Integro-Differential Operators

Now all the ingredients are prepared, we claim that the diffusion term and
reaction term can be rewritten as bilinear integro-differential operators.

######  Theorem 1.

The particle-wave interaction system in Equation(2.2) is equivalent to

$\begin{split}\frac{\partial f(\mathbf{p},t)}{\partial
t}&=\mathbb{B}(W,f)=-\int_{\mathbb{R}_{k}^{3}}\mathcal{L}^{*}\left(\mathcal{B}(\mathbf{p},\mathbf{k})W\mathcal{L}f\right)d\mathbf{k},\\\
\frac{\partial W(\mathbf{k},t)}{\partial
t}&=\mathbb{H}(f,W)=\int_{\mathbb{R}_{p}^{3}}(\mathcal{L}E(\mathbf{p}))W\mathcal{B}(\mathcal{L}f)d\mathbf{p}.\end{split}$
(2.7)

###### Proof.

By definition,

$\int\delta(s(\mathbf{x}))s(\mathbf{x})g(\mathbf{x})d\mathbf{x}=\frac{1}{|\nabla_{x}s|}\int_{\\{\mathbf{x}:s(\mathbf{x})=0\\}}s(\mathbf{x})g(\mathbf{x})d\pi_{s}=0$

Note that

$\frac{l\omega_{c}/\gamma-\omega\sin^{2}\theta}{\omega\sin\theta\cos\theta}=\frac{\omega\cos\theta-
k_{\parallel}v}{\omega\sin\theta},$

on the resonance manifold.

We obtain Equation(2.7) by substituting the above identity into Equation(2.2)
and rewriting everything in cylindrical coordinates(i.e. replacing $\theta$
with $\arccos{p_{\parallel}/p}$). ∎

######  Remark 3.

Both the particle diffusion operator $\mathbb{B}$ and the wave reaction
operator $\mathbb{H}$ mix particle momentum $\mathbf{p}$ and plasmon wave
vector $\mathbf{k}$ through the absorption/emission kernel
$\mathcal{B}(\mathbf{p},\mathbf{k})$.

######  Remark 4.

One might have noticed that $\mathcal{L}E(\mathbf{p})=v(\mathbf{p})$. The
reason we write $\mathcal{L}E$ rather than $v$ is to induce our conservative
semi-discrete form and to save preprocessing time. The details will be
addressed next.

###  2.2 The Unconditionally Conservative Weak Form and H-Theorem

For the purpose of either modeling or numerical implementation, the
absorption/emission kernel $\mathcal{B}(\mathbf{p},\mathbf{k})$ is usually
replaced with its approximation
$\mathcal{B}_{\varepsilon}(\mathbf{p},\mathbf{k})$. Here we present two
examples for such approximation.

  * •

Approximation to the identity.

The kernel is approximated with

$\mathcal{B}_{h}^{\text{ati}}(\mathbf{p},\mathbf{k})=\sum_{l}U_{l}(\mathbf{p},\mathbf{k})\frac{1}{\varepsilon}\varphi\left(\frac{s_{l}(\mathbf{p},\mathbf{k})}{\varepsilon}\right),$
(2.8)

where the compactly supported and positive function $\varphi(z)$ has unit
mass, i.e. $\int_{\mathbb{R}}\varphi(z)dz=1$.

Recall the definition of resonance manifold in Equation(2.4), it is a
hypersurface implicitly determined by the resonance condition. The above
approximation is equivalent to broadening of the resonance manifold, as the
approximated hypersurface has finite “width” proportional to $\varepsilon$.

  * •

Marching cube/simplex algorithm.

The kernel is approximated with

$\mathcal{B}_{\varepsilon}^{\text{msa}}(\mathbf{p},\mathbf{k})=\sum_{l}U_{l}(\mathbf{p},\mathbf{k})\delta\left(L_{\varepsilon}s_{l}(\mathbf{p},\mathbf{k})\right),$
(2.9)

where $L_{\varepsilon}s_{l}$ represents the piecewise linear interpolation of
$s_{l}$. As have been illustrated in [5, 12], such approximation discretizes
the resonance manifold, i.e. replaces the smooth hypersurface with a disjoint
union of simplices, thus enabling convenient numerical integration.

In what follows we derive the special weak form, and prove that even if the
emission/absorption kernel $\mathcal{B}$ is replaced/approximated, we can
still preserve mass, momentum and energy with the proposed form.

To obtain the weak formulation associated with the system(2.7), test it with
$\phi(\mathbf{p})$ and $\eta(\mathbf{k})$, we obtain that

$\begin{split}\int_{\mathbb{R}_{p}^{3}}\frac{\partial f}{\partial t}\phi
d\mathbf{p}&=-\int_{\mathbb{R}_{k}^{3}}d\mathbf{p}\int_{\mathbb{R}_{p}^{3}}d\mathbf{k}\\{\mathcal{L}\phi\mathcal{L}fW\mathcal{B}\\}\\\
\int_{\mathbb{R}_{k}^{3}}\frac{\partial W}{\partial t}\eta
d\mathbf{k}&=\int_{\mathbb{R}_{k}^{3}}d\mathbf{k}\int_{\mathbb{R}_{p}^{3}}d\mathbf{p}\\{\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}\\}\end{split}$

Note that the order of integration here is different on the right-hand side.
In what follows, assume that
$\int_{\mathbb{R}_{k}^{3}\times\mathbb{R}_{p}^{3}}d\mathbf{k}d\mathbf{p}|\mathcal{L}\phi\mathcal{L}fW\mathcal{B}|$
and
$\int_{\mathbb{R}_{k}^{3}\times\mathbb{R}_{p}^{3}}d\mathbf{k}d\mathbf{p}|\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}|$
are finite, therefore by Fubini’s theorem, the order of integration does not
matter,

$\begin{split}\int_{\mathbb{R}_{p}^{3}}\left(\int_{\mathbb{R}_{k}^{3}}\mathcal{L}\phi\mathcal{L}fW\mathcal{B}d\mathbf{k}\right)d\mathbf{p}=\int_{\mathbb{R}_{k}^{3}}\left(\int_{\mathbb{R}_{p}^{3}}\mathcal{L}\phi\mathcal{L}fW\mathcal{B}d\mathbf{p}\right)d\mathbf{k}=\iint_{\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}}d\mathbf{k}d\mathbf{p}\left\\{\mathcal{L}\phi\mathcal{L}fW\mathcal{B}\right\\}\\\
\int_{\mathbb{R}_{k}^{3}}\left(\int_{\mathbb{R}_{p}^{3}}\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}d\mathbf{p}\right)d\mathbf{k}=\int_{\mathbb{R}_{p}^{3}}\left(\int_{\mathbb{R}_{k}^{3}}\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}d\mathbf{k}\right)d\mathbf{p}=\iint_{\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}}d\mathbf{k}d\mathbf{p}\left\\{\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}\right\\}.\end{split}$

On the right-hand side is inner products of bilinear integro-differential
operators with test functions. Therefore to simplify the notation, we can
define trilinear forms $B$ and $H$ as follows:

$\begin{split}B(f,W,\phi)&\coloneqq\iint_{\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}}d\mathbf{k}d\mathbf{p}\left\\{\mathcal{L}\phi\mathcal{L}fW\mathcal{B}\right\\},\\\
H(W,f,\eta)&\coloneqq\iint_{\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}}d\mathbf{k}d\mathbf{p}\left\\{\eta\mathcal{L}E\mathcal{L}fW\mathcal{B}\right\\}.\end{split}$
(2.10)

As a result, the weak form of system(2.7) can be written as,

$\begin{split}\left(\frac{\partial f}{\partial
t},\phi\right)_{p}&=-B(f,W,\phi),\\\ \left(\frac{\partial W}{\partial
t},\eta\right)_{k}&=H(W,f,\eta).\end{split}$ (2.11)

Due to resonance, there are infinitely many equivalent forms for the same
equation, for example,

$\sum_{l}U_{l}(\mathbf{p},\mathbf{k})\delta(\omega-
k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma)\left(\frac{k_{\parallel}v_{\parallel}}{\omega}\frac{p}{p_{\parallel}}\frac{\partial\phi}{\partial
p_{\parallel}}+(1-\frac{k_{\parallel}v_{\parallel}}{\omega})\frac{p}{p_{\perp}}\frac{\partial\phi}{\partial
p_{\perp}}\right)$

is always equal to

$\sum_{l}U_{l}(\mathbf{p},\mathbf{k})\delta(\omega-
k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma)\left(\left(\frac{k_{\parallel}v_{\parallel}+l\omega_{c}/\gamma}{\omega}\right)^{\alpha}\frac{k_{\parallel}v_{\parallel}}{\omega}\frac{p}{p_{\parallel}}\frac{\partial\phi}{\partial
p_{\parallel}}+(1-\frac{k_{\parallel}v_{\parallel}}{\omega})\frac{p}{p_{\perp}}\frac{\partial\phi}{\partial
p_{\perp}}\right),$

for any constant $\alpha>0$.

In the following theorem, we prove the superiority of the proposed form, i.e.
the unconditional conservation property.

######  Theorem 2 (unconditional conservation).

If $f(\mathbf{p},t)$ and $W(\mathbf{k},t)$ solve the system(2.11) with
emission/absorption kernel being replaced by $\mathcal{B}_{\varepsilon}$, then
for any $\mathcal{B}_{\varepsilon}$ we have the following conservation laws,

  * •

Mass Conservation

$\frac{\partial}{\partial t}\mathcal{M}_{tot}=\frac{\partial}{\partial
t}\left((f,1)_{p}+(W,0)_{k}\right)=0$

  * •

Momentum Conservation

$\frac{\partial}{\partial
t}\mathcal{P}_{\parallel}^{tot}=\frac{\partial}{\partial
t}\left((f,p_{\parallel})_{p}+(\frac{W}{\hbar\omega},\hbar
k_{\parallel})_{k}\right)=0,$

  * •

Energy Conservation

$\frac{\partial}{\partial t}\mathcal{E}_{tot}=\frac{\partial}{\partial
t}((f,E)_{p}+(\frac{W}{\hbar\omega},\hbar\omega)_{k})=0$

###### Proof.

Sum the two rows in Equation(2.11), we have

$\left(\frac{\partial f}{\partial t},\phi\right)_{p}+\left(\frac{\partial
W}{\partial
t},\eta\right)_{k}=\int_{\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}}d\mathbf{k}d\mathbf{p}\left\\{\left(\eta\mathcal{L}E-\mathcal{L}\phi\right)\mathcal{L}fW\mathcal{B}_{\varepsilon}\right\\}$
(2.12)

Substitute the test functions $\\{\phi,\eta\\}$ with conservation quantity
pairs $\\{\phi_{c},\eta_{c}\\}$, i.e. mass pair $\\{1,0\\}$, parallel momentum
pair $\\{p_{\parallel},\frac{k_{\parallel}}{\omega}\\}$ and energy pair
$\\{E,1\\}$. In principle the condition for conservation is
$\eta_{c}\mathcal{L}E-\mathcal{L}\phi_{c}=0$ on the approximate resonance
manifold
$\mathcal{S}=\left\\{(\mathbf{p},\mathbf{k})\in\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}:\mathcal{B}_{\varepsilon}(\mathbf{p},\mathbf{k})\neq
0\right\\}$. However, due to our particular definition of the directional
differential operator $\mathcal{L}$, we actually have
$\eta_{c}\mathcal{L}E-\mathcal{L}\phi_{c}=0$ on the whole domain
$\mathbb{R}_{p}^{3}\times\mathbb{R}_{k}^{3}$. Therefore the conservation laws
hold regardless of emission/absorption kernel $\mathcal{B}_{\varepsilon}$. ∎

######  Remark 5.

The unconditionally conservative form also exists for unmagnetized plasmas
with cylindrical symmetry, where we just replace the emission/absorption
kernel with
$\mathcal{B}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\omega^{2}}{k^{2}v^{2}}\delta(\omega-
k_{\parallel}v_{\parallel}-k_{\perp}v_{\perp}\cos\alpha)d\alpha$.

######  Remark 6.

Since the conservation laws solely depend on
$\eta\mathcal{L}E-\mathcal{L}\phi=0$, the unconditional conservative form and
the scheme we are going to propose can be generalized for time-dependent
dispersion relation $\omega(\mathbf{k};t)$ with no extra effort. An important
example is the self-consistent dispersion relation
$\omega=\omega(\mathbf{k};f_{1}(t),f_{2}(t),\cdots)$. The only obstacle is the
extra computational cost of updating the interaction tensors in each step. As
will be shown in section 4, that calculation can be expensive.

Recall the definition of emission/absorption kernel $\mathcal{B}$,

$\mathcal{B}(\mathbf{p},\mathbf{k})=\sum_{l=-\infty}^{+\infty}U_{l}(\mathbf{p};\mathbf{k})\delta(s_{l}(\mathbf{p};\mathbf{k})).$

Test the equation for particle pdf with $\phi(\mathbf{p})=\log f(\mathbf{p})$,
since $U_{l}$ and $W$ are non-negative, the right-hand side will be non-
positive,

$(\frac{\partial f}{\partial t},\log f)_{p}=-\iint
d\mathbf{k}d\mathbf{p}\frac{1}{f}(\mathcal{L}f)^{2}W\mathcal{B}\leq 0,$

thus we can prove the dissipation of entropy, i.e. H-theorem for the particle
pdf,

$\frac{\partial}{\partial t}(f,\log f)_{p}=(\frac{\partial f}{\partial t},\log
f)_{p}+(\frac{\partial f}{\partial t},1)_{p}\leq 0.$ (2.13)

##  3 The Conservative Discretization

This section aims to find a semi-discrete problem that consistently
approximates the original system, and at the same time preserves discrete
conservation laws. So in the following subsections, we will first introduce
our finite element discretization, the necessary projection operators, and
then elaborate on the conservation technique.

###  3.1 The Finite Element Discretization

#### The Cut-Off Domain and Boundary Conditions

Analogous to existing work on kinetic equations, for example, the papers of
Zhang et al.[20, 21], we assume that given any $0<\epsilon_{p}\ll 1$ and
$0<\epsilon_{k}\ll 1$, there exists finite cylindrical domains
$\Omega_{p}^{L}\subsetneq\mathbb{R}_{p}^{3}$ and
$\Omega_{k}^{L}\subsetneq\mathbb{R}_{k}^{3}$ such that for any $t\geq 0$,

$\left|1-\frac{\int_{\Omega^{L}_{p}}f(\mathbf{p},t)d^{3}p}{\int_{\mathbb{R}_{p}^{3}}f(\mathbf{p},t)d^{3}p}\right|\leq\epsilon_{p},$

and

$\left|1-\frac{\int_{\Omega^{L}_{k}}W(\mathbf{k},t)d^{3}k}{\int_{\mathbb{R}_{k}^{3}}W(\mathbf{k},t)d^{3}k}\right|\leq\epsilon_{k}.$

The particle momentum cut-off domain $\Omega_{p}^{L}$ is supposed to be
adaptive, while in our numerical experiments it turns out that, as a result of
anisotropic diffusion, there is no need to extend it.

Then it is reasonable to solve the equations in cut-off domains
$\Omega^{L}_{p}$ and $\Omega^{L}_{k}$. For the wave sed $W(\mathbf{k})$, there
is no need for a boundary condition since there is no flux in wave vector
space. For the particle pdf, we have the following choices, and when the
domain $\Omega^{L}_{p}$ is large enough, they are actually equivalent.

On the boundary $\partial\Omega^{L}_{p}$ of cut-off domain $\Omega^{L}_{p}$,
$|f|$ and $|\nabla_{\mathbf{p}}f|$ are nearly zero, two types of boundary
conditions can be applied,

  1. 1.

The zero-value boundary condition

$f=0,\ \forall\mathbf{p}\in\partial\Omega^{L}_{p}$

  2. 2.

The zero-flux boundary condition

$\left(D[W]\nabla_{\mathbf{p}}f\right)\cdot\mathbf{n}=0,\
\forall\mathbf{p}\in\partial\Omega^{L}_{p}$

Suppose we test the diffusion equation with $\phi_{h}\in V_{h}$. With
Neumann’s boundary condition, i.e. in the zero-flux case, the semi-discrete
weak form reads:

$(\frac{\partial f_{h}}{\partial
t},\phi_{h})+(D[W_{h}]\nabla_{\mathbf{p}}f_{h},\nabla_{\mathbf{p}}\phi_{h})=0$

For Dirichlet’s boundary conditions given by to zero-value on the discretized
boundary, i.e. $f_{\text{ini}}\mid_{\partial\Omega^{L}_{p}}\equiv 0$,
Nitsche’s method [13] applies, hence the weak the semi-discrete form reads

$(\frac{\partial f_{h}}{\partial
t},\phi_{h})+(D[W_{h}]\nabla_{\mathbf{p}}f_{h},\nabla_{\mathbf{p}}\phi_{h})-\langle\left(D[W_{h}]\nabla_{\mathbf{p}}f_{h}\right)\cdot\mathbf{n}_{p},\phi_{h}\rangle_{\partial\Omega^{L}_{p}}+\langle\left(D[W_{h}]\nabla_{\mathbf{p}}\phi_{h}\right)\cdot\mathbf{n}_{p},f_{h}\rangle_{\partial\Omega^{L}_{p}}=0$

The only difference between them is the boundary integral, which can be below
machine epsilon for large enough $\Omega^{L}_{p}$, because $D$ and $\phi_{h}$
are finite, while $|f_{h}|$ and $|\nabla_{\mathbf{p}}f_{h}|$ goes to zero as
we enlarge the domain. Stability can be proved for both formulations, in the
rest of the article, for simplicity, we will use the zero-flux boundary
condition.

#### The Finite Element Spaces

Since we have assumed cylindrical symmetry, the 3P-3K problem actually becomes
2P-2K.

$\begin{split}&\mathbf{p}=(p_{1},p_{2},p_{3})\in\Omega^{L}_{p}\Leftrightarrow(p_{\parallel},p_{\perp})\in\tilde{\Omega}^{L}_{p}\subset\mathbb{R}\times\mathbb{R}^{+}\\\
&\mathbf{k}=(k_{1},k_{2},k_{3})\in\Omega^{L}_{k}\Leftrightarrow(k_{\parallel},k_{\perp})\in\tilde{\Omega}^{L}_{k}\subset\mathbb{R}\times\mathbb{R}^{+}\end{split}$

Let $\mathcal{T}_{h}^{p}=\\{R_{p}\\}$, $\mathcal{T}_{h}^{k}=\\{R_{k}\\}$ be
rectangular partitions of $\tilde{\Omega}^{L}_{p}$ and
$\tilde{\Omega}^{L}_{k}$ respectively. We define the meshsize for momentum
space as
$h_{p}=\textnormal{max}_{R_{p}\in\mathcal{T}_{h}^{p}}\textnormal{diam}(R_{p})$
and the meshsize for wave vector space as
$h_{k}=\textnormal{max}_{R_{k}\in\mathcal{T}_{h}^{k}}\textnormal{diam}(R_{k})$.

The test space for particle pdf consists of continuous piecewise polynomials
with degree $\alpha_{1}$,

$\mathcal{G}_{h}^{\alpha_{1}}=\\{f(p_{\parallel},p_{\perp})\in
C^{0}(\Omega_{p}):f|_{R_{p}}\in
Q^{\alpha_{1}}(R_{p}),\forall\leavevmode\nobreak\
R_{p}\in\mathcal{T}_{h}^{p}\\}.$ (3.1)

The test space for wave sed consists of discontinuous piecewise polynomials
with degree $\alpha_{2}$,

$\mathcal{W}_{h}^{\alpha_{2}}=\\{W(k_{\parallel},k_{\perp}):W|_{R_{k}}\in
Q^{\alpha_{2}}(R_{k}),\forall\leavevmode\nobreak\
R_{k}\in\mathcal{T}_{h}^{k}\\}.$ (3.2)

To ensure positivity of $W_{h}$, it is required that $\alpha_{2}=0$ or
$\alpha_{2}=1$, the reason will be addressed later.

As will be shown in the next section, one of the key points to conservation is
replacing $v_{\parallel}$, $v_{\perp}$ and $k_{\parallel}/\omega$ with
$\frac{\partial E_{h}(\mathbf{p})}{\partial p_{\parallel}}$, $\frac{\partial
E_{h}(\mathbf{p})}{\partial p_{\perp}}$ and $N_{\parallel,h}$, where
$E_{h}=\Pi_{p,h}E(\mathbf{p})$ is the discrete particle kinetic energy, and
$N_{\parallel,h}=\Pi_{k,h}N_{\parallel}$ is the discrete refraction index. The
projection operators can be arbitrarily chosen as long as they satisfy the
following conditions:

  1. 1.

The projection $\Pi_{p,h}$ into test space $\mathcal{G}_{h}^{\alpha_{1}}$ must
satisfy that

$\lim_{h\rightarrow
0}\|\Pi_{p,h}g(\mathbf{p})-g(\mathbf{p})\|_{L^{2}(\Omega_{p}^{L})}=0,\ \forall
g\in L^{2}(\Omega_{p}^{L}),$

and

$\lim_{h\rightarrow
0}\|\Pi_{p,h}E(\mathbf{p})-E(\mathbf{p})\|_{H^{1}(\Omega_{p}^{L})}=0.$

  2. 2.

The projection $\Pi_{k,h}$ into test space $\mathcal{W}_{h}^{\alpha_{2}}$ must
satisfy that

$\lim_{h\rightarrow
0}\|\Pi_{k,h}\xi(\mathbf{k})-\xi(\mathbf{k})\|_{L^{2}(\Omega_{k}^{L})}=0,\
\forall\xi\in L^{2}(\Omega_{k}^{L}).$

There is no need to specify particular projections until we implement them in
the numerical examples, our method works with any of them.

###  3.2 The Conservative Semi-Discrete Form

Adopting the zero-flux boundary condition, testing the system on the cut-off
domain with $\phi_{h}\in\mathcal{G}_{h}^{\alpha_{1}}$ and
$\eta_{h}\in\mathcal{W}_{h}^{\alpha_{2}}$, we write the following semi-
discrete weak form,

$\begin{split}(\frac{\partial f_{h}}{\partial
t},\phi_{h})_{p}&=-B^{u}_{L}(f_{h},W_{h},\phi_{h})\coloneqq-\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\mathcal{L}\phi_{h}\mathcal{L}f_{h}W_{h}\mathcal{B}\\},\\\
(\frac{\partial W_{h}}{\partial
t},\eta_{h})_{k}&=H^{u}_{L}(W_{h},f_{h},\eta_{h})\coloneqq\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\eta_{h}\mathcal{L}E\mathcal{L}f_{h}W_{h}\mathcal{B}\\},\end{split}$

where the subscript $L$ means integral on cut-off domain, the superscript $u$
means unconservative. We will first analyze the source of conservation errors
and then present our conservative semi-discrete trilinear forms $B_{L}$ and
$H_{L}$.

#### The Source of Conservation Errors

Suppose different quadrature rules $R_{1}$ and $R_{2}$ are used for different
equations,

$\begin{split}(\frac{\partial f_{h}}{\partial
t},\phi_{h})_{p}&=R_{1}\left[-\int_{\Omega_{p}^{L}}d\mathbf{p}\int_{\Omega_{k}^{L}}d\mathbf{k}\\{\mathcal{L}\phi_{h}\mathcal{L}f_{h}W_{h}\mathcal{B}\\}\right]\\\
(\frac{\partial W_{h}}{\partial
t},\eta_{h})_{k}&=R_{2}\left[\int_{\Omega_{k}^{L}}d\mathbf{k}\int_{\Omega_{p}^{L}}d\mathbf{p}\\{\eta_{h}\mathcal{L}E\mathcal{L}f_{h}W_{h}\mathcal{B}\\}\right]\end{split}$

The error of conservation laws can be decomposed into three terms,

$\begin{split}&\frac{\partial}{\partial
t}((f_{h},\Pi_{p,h}\phi)_{p}+(W_{h},\Pi_{k,h}\eta)_{k})\\\
=&R_{1}\left[-\int_{\Omega_{p}^{L}}d\mathbf{p}\int_{\Omega_{k}^{L}}d\mathbf{k}\left\\{\left(\mathcal{L}\Pi_{p,h}\phi\right)\mathcal{L}f_{h}W_{h}\mathcal{B}\right\\}\right]+R_{2}\left[\int_{\Omega_{k}^{L}}d\mathbf{k}\int_{\Omega_{p}^{L}}d\mathbf{p}\left\\{\left(\Pi_{k,h}\eta\mathcal{L}E\right)\mathcal{L}f_{h}W_{h}\mathcal{B}\right\\}\right]\\\
=&A_{1}+A_{2}+A_{3},\end{split}$

where

$\begin{split}A_{1}&=\left(R_{1}-I\right)\left[\int_{\Omega_{k}^{L}}d\mathbf{k}\int_{\Omega_{p}^{L}}d\mathbf{p}\left\\{\left(-\mathcal{L}\Pi_{p,h}\phi\right)\mathcal{L}f_{h}W_{h}\mathcal{B}\right\\}\right],\\\
A_{2}&=\left(I-R_{2}\right)\left[\int_{\Omega_{k}^{L}}d\mathbf{k}\int_{\Omega_{p}^{L}}d\mathbf{p}\left\\{\left(-\mathcal{L}\Pi_{p,h}\phi\right)\mathcal{L}f_{h}W_{h}\mathcal{B}\right\\}\right],\\\
A_{3}&=R_{2}\left[\int_{\Omega_{k}^{L}}d\mathbf{k}\int_{\Omega_{p}^{L}}d\mathbf{p}\left\\{\left(\Pi_{k,h}\eta\mathcal{L}E-\mathcal{L}\Pi_{p,h}\phi\right)\mathcal{L}f_{h}W_{h}\mathcal{B}\right\\}\right].\end{split}$

The error terms $A_{1}$ and $A_{2}$ are caused by inconsistent numerical
integration on the resonance manifold. Suppose that $R_{1}-I$ is of the same
order as $O(h^{a})$, and quadrature rule $R_{2}$ has error $O(h^{b})$, then
the sum will be roughly $O(h^{\min\\{a,b\\}})$. The last error term $A_{3}$ is
a result of projection error, whose order depends on the degree of test
spaces, $\alpha_{1}$ and $\alpha_{2}$.

Note that $A_{1}$ and $A_{2}$ cancel out when we use the same quadrature
rules, i.e. $R_{1}=R_{2}$. In what follows, we will introduce a conservative
semi-discrete form such that $A_{3}$ disappears.

#### The Conservative Semi-Discrete Form

Recall the definition of directional differential operator $\mathcal{L}$,

$\mathcal{L}g=\frac{k_{\parallel}v_{\parallel}}{\omega}\frac{p}{p_{\parallel}}\frac{\partial
g}{\partial
p_{\parallel}}+(1-\frac{k_{\parallel}v_{\parallel}}{\omega})\frac{p}{p_{\perp}}\frac{\partial
g}{\partial p_{\perp}}=N_{\parallel}\frac{\partial E}{\partial
p_{\perp}}\frac{p}{p_{\perp}}\frac{\partial g}{\partial
p_{\parallel}}+(1-N_{\parallel}\frac{\partial E}{\partial
p_{\parallel}})\frac{p}{p_{\perp}}\frac{\partial g}{\partial p_{\perp}}.$

We propose a discretized operator $\mathcal{L}_{h}$ defined as follows,

$\mathcal{L}_{h}g\coloneqq N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\perp}}\frac{p}{p_{\perp}}\frac{\partial g}{\partial
p_{\parallel}}+(1-N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\parallel}})\frac{p}{p_{\perp}}\frac{\partial g}{\partial p_{\perp}},$
(3.3)

where the discretized kinetic energy is defined as
$E_{h}=\Pi_{p,h}E(\mathbf{p})$, and the discretized wave refractive index is
defined as
$N_{\parallel,h}=\Pi_{k,h}N_{\parallel}=\Pi_{k,h}\frac{k_{\parallel}}{\omega(\mathbf{k})}$.

The main result of this paper is stated in the following theorem.

######  Theorem 3.

If $f_{h}(\mathbf{p},t)$ and $W_{h}(\mathbf{k},t)$ are solutions of the
following semi-discrete weak form,

$\begin{split}(\frac{\partial f_{h}}{\partial
t},\phi_{h})_{p}&=-B_{L}(f_{h},W_{h},\phi_{h})\coloneqq-\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\mathcal{L}_{h}\phi_{h}\mathcal{L}_{h}f_{h}W_{h}\mathcal{B}_{h}\\},\\\
(\frac{\partial W_{h}}{\partial
t},\eta_{h})_{k}&=H_{L}(W_{h},f_{h},\eta_{h})\coloneqq\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\eta_{h}\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}W_{h}\mathcal{B}_{h}\\},\end{split}$
(3.4)

then the following discrete conservation laws hold,

$\begin{split}\frac{\partial}{\partial
t}\mathcal{M}_{tot,h}&=\frac{\partial}{\partial
t}((f_{h},\Pi_{p,h}1)_{p}+(W_{h},0)_{k})=0,\\\ \frac{\partial}{\partial
t}\mathcal{P}^{\parallel}_{tot,h}&=\frac{\partial}{\partial
t}\left((f_{h},\Pi_{p,h}p_{\parallel})_{p}+(W_{h},\Pi_{k,h}N_{\parallel})_{k}\right)=0,\\\
\frac{\partial}{\partial t}\mathcal{E}_{tot,h}&=\frac{\partial}{\partial
t}((f_{h},\Pi_{p,h}E(\mathbf{p}))_{p}+(W_{h},\Pi_{k,h}1)_{k})=0.\end{split}$

###### Proof.

Substitute the discrete conservation pairs $\\{\Pi_{p,h}1,0\\}$,
$\\{\Pi_{p,h}p_{\parallel},\Pi_{k,h}N_{\parallel}\\}$ and
$\\{\Pi_{p,h}E(\mathbf{p}),\Pi_{k,h}1\\}$ into semi-discrete form (3.4) and
use the definition of $\mathcal{L}_{h}$. ∎

######  Corollary 1.

If in addition to the assumptions of Theorem(3), the projections are $L^{2}$
orthogonal projections, i.e.

$\left(u-\Pi_{p,h}u,v\right)_{p}=0,\forall v\in\mathcal{G}_{h}^{\alpha_{1}}$

and

$\left(U-\Pi_{k,h}U,V\right)_{k}=0,\forall V\in\mathcal{W}_{h}^{\alpha_{2}},$

then the exact conservation laws are preserved, i.e.

$\begin{split}\frac{\partial}{\partial
t}\mathcal{M}_{tot,h}&=\frac{\partial}{\partial
t}((f_{h},1)_{p}+(W_{h},0)_{k})=0,\\\ \frac{\partial}{\partial
t}\mathcal{P}^{\parallel}_{tot,h}&=\frac{\partial}{\partial
t}\left((f_{h},p_{\parallel})_{p}+(W_{h},N_{\parallel})_{k}\right)=0,\\\
\frac{\partial}{\partial t}\mathcal{E}_{tot,h}&=\frac{\partial}{\partial
t}((f_{h},E(\mathbf{p}))_{p}+(W_{h},1)_{k})=0.\end{split}$

###### Proof.

Use the fact that $f_{h}\in\mathcal{G}_{h}^{\alpha_{1}}$ and
$W_{h}\in\mathcal{W}_{h}^{\alpha_{2}}$. ∎

######  Remark 7.

Same as stated in Theorem(2), our semi-discrete weak form is also
unconditionally conservative, i.e. the conservation does not depend on a
particular discrete emission/absorption kernel $\mathcal{B}_{h}$.

##  4 The Sparse Interaction Tensors

Suppose that the test spaces are spanned by basis functions, i.e.
$\mathcal{G}_{h}^{\alpha_{1}}=\text{span}\\{\phi_{i}\\}$ and
$\mathcal{W}_{h}^{\alpha_{2}}=\text{span}\\{\eta_{j}\\}$. Then we can express
the discrete particle pdf $f_{h}$ and wave sed $W_{h}$ as a linear combination
of basis functions.

$\begin{split}f_{h}(\mathbf{p},t)&=\sum_{i=1}^{N_{f}}a_{i}(t)\phi_{i}(\mathbf{p})\\\
W_{h}(\mathbf{k},t)&=\sum_{j=1}^{N_{w}}w_{j}(t)\eta_{j}(\mathbf{k})\end{split}$

By definition, $E_{h}=\Pi_{p,h}E\in\mathcal{G}_{h}^{\alpha_{1}}$, therefore it
is also a linear combination of basis functions,
$E_{h}=\sum_{q=1}^{N_{f}}E_{q}\phi_{q}$.

Substitute the above expressions into Equation(3.4), then the semi-discrete
system becomes a first-order finite dimension ODE system:

$\begin{split}\sum_{i=1}^{N_{f}}\frac{\partial a_{i}}{\partial
t}\int\phi_{i}\phi_{m}d^{3}p&=-\sum_{n=1}^{N_{f}}\sum_{k=1}^{N_{w}}a_{n}w_{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\int_{\Omega^{L}_{k}}d\mathbf{k}\\{\mathcal{L}_{h}\phi_{m}\mathcal{L}_{h}\phi_{n}\eta_{k}\mathcal{B}(\mathbf{p};\mathbf{k})\\}\\\
\sum_{j=1}^{N_{w}}\frac{\partial w_{j}}{\partial
t}\int\eta_{j}\eta_{q}d^{3}k&=\sum_{k=1}^{N_{w}}\sum_{n=1}^{N_{f}}w_{k}a_{n}\int_{\Omega^{L}_{k}}d\mathbf{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\eta_{q}\mathcal{L}_{h}E_{h}\mathcal{L}_{h}\phi_{n}\eta_{k}\mathcal{B}(\mathbf{p};\mathbf{k})\\}\end{split}$

Denote the mass matrix for particle pdf as $A_{im}=(\phi_{i},\phi_{m})_{p}$,
and denote the mass matrix for wave sed as $G_{jq}=(\eta_{j},\eta_{q})_{k}$.

Analogously, define the interaction tensors $B$ and $H$ corresponding to the
trilinear forms.

$\begin{split}B_{nkm}&=B(\phi_{n},\eta_{k},\phi_{m})\\\
H_{knq}&=H(\eta_{k},\phi_{n},\eta_{q})\end{split}$

As a result, we obtain the nonlinear ODE system corresponding to semi-discrete
weak form(3.4):

$\begin{split}\frac{\partial a_{i}}{\partial t}A_{im}&=-a_{n}w_{k}B_{nkm}\\\
\frac{\partial w_{j}}{\partial t}G_{jq}&=w_{k}a_{n}H_{knq}\end{split}$ (4.1)

The interaction tensors $B$ and $H$ are both sparse tensors for two reasons:
compactly supported basis and the resonant feature of trilinear forms. Taking
particle interaction tensor $B$ as an example, $B_{nkm}=0$ when

  1. 1.

$\phi_{m}$ and $\phi_{n}$ are not in neighboring elements.

  2. 2.

$\phi_{m}$ and $\eta_{k}$ do not ”resonate”, i.e.
$\text{supp}(\phi_{m})\times\text{supp}(\eta_{k})$ does not intersect with the
resonant manifold.

Suppose in each dimension we have $O(n)$ meshes, then the shape of particle
interaction tensor $B$ is roughly $O(n^{2})\times O(n^{2})\times O(n^{2})$,
while the number of nonzero elements will be only $O(n^{3})$, i.e. the
sparsity of tensor $B$ is about $1-\frac{1}{O(n^{3})}$. A similar analysis can
also be applied to the wave interaction tensor $H$.

We observed that the trilinear forms $B$ and $H$ defined in Equation(2.10)
have similar structures. Therefore one might wonder if there is any relation
between the interaction tensors $B$ and $H$. It turns out that when
$\alpha_{2}=0$, i.e piecewise constant basis functions are used for wave sed
$W_{h}$, we can infer any nonzero element of wave interaction tensor $H$ from
particle interaction tensor $B$. In practice, the interaction tensors are
precomputed and saved for later use. Taking advantage of this relation, we can
save half the time of preprocessing. The derivation is as follows.

When $\alpha_{2}=0$, $\mathcal{W}_{h}^{\alpha_{2}}=\text{span}\\{\eta_{j}\\}$
are piecewise constant functions, we have

$\eta_{i}(\mathbf{k})\eta_{j}(\mathbf{k})=\delta_{ij}\eta_{i}(\mathbf{k}).$

Then the mass matrix for wave sed is diagonal,

$G_{jq}=(\eta_{j},\eta_{q})_{k}=\int_{R_{k}^{j}}\delta_{jq}d^{3}k=\text{diag}(\mu(R_{k}^{j})),$

where $\mu(R_{k}^{j})=\int_{R_{k}^{j}}1d^{3}k$ is the measure of $j$-th
element in $\Omega^{L}_{k}$.

Moreover, note that if we define a $4$-th order tensor

$\tilde{H}_{mknq}\coloneqq\int_{\Omega^{L}_{k}}d\mathbf{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\mathcal{L}_{h}\phi_{m}\mathcal{L}_{h}\phi_{n}\eta_{k}\eta_{q}\mathcal{B}(\mathbf{p};\mathbf{k})\\}$

Recall the expansion $E_{h}=\sum_{q=1}^{N_{f}}E_{q}\phi_{q}$, and substitute
it into the definition of wave interaction tensor $H$, we obtain the relation
between $H_{knq}$ and $\tilde{H}_{mknq}$,

$H_{knq}=H(\eta_{k},\phi_{n},\eta_{q})=\sum_{m=1}^{N_{f}}E_{m}\tilde{H}_{mknq}$

It can be observed that the form of $\tilde{H}_{mknq}$ is almost identical to
the definition of particle tensor $B_{mnk}$, except for the extra $\eta_{q}$.
Replace $\eta_{k}(\mathbf{k})\eta_{q}(\mathbf{k})$ with
$\delta_{kq}\eta_{k}(\mathbf{k})$, we obtain the relation between
$\tilde{H}_{mknq}$ and $B_{mnk}$,

$\tilde{H}_{mknq}=\delta_{kq}\int_{\Omega^{L}_{k}}d\mathbf{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\mathcal{L}_{h}\phi_{m}\mathcal{L}_{h}\phi_{n}\eta_{k}\mathcal{B}(\mathbf{p};\mathbf{k})\\}=\delta_{kq}B_{mnk}.$

Therefore $B_{mnk}$ and $E_{m}$ is all we need to calculate $H_{knq}$,

$H_{knq}=\sum_{m=1}^{N_{f}}E_{m}\tilde{H}_{mknq}=\sum_{m=1}^{N_{f}}E_{m}\delta_{kq}B_{mnk}=\delta_{kq}\sum_{m=1}^{N_{f}}E_{m}B_{mnk}=\begin{cases}0,&k\neq
q\\\ \sum_{m=1}^{N_{f}}E_{m}B_{mnk}&k=q\end{cases}$

##  5 Stability and Positivity

In this section, we investigate the stability of the fully discretized
nonlinear system. With semi-implicit time discretization, there is no
constraint on time step size from the CFL condition. However, the stability
will rely on the positivity of $W_{h}$, which results in a condition for the
time step size, relevant to the gradient of particle pdf $f_{h}$. The
condition will not cause any trouble for implementation, because we can always
adapt the step size a posteriori.

###  5.1 Stability of the Semi-Discrete Form

Consider the equation for particle pdf only, it has the form of a diffusion
equation, thus its stability relies on the fact that the diffusion coefficient
is positive semi-definite, which further relies on the positivity of wave sed
$W_{h}$.

######  Lemma 1 ($L^{2}$ stability of $f_{h}(\mathbf{p})$ and $L^{1}$ bound of
$W_{h}(\mathbf{k})$).

Suppose $f_{h}(\mathbf{p},t)$ and $W_{h}(\mathbf{k},t)$ are the solution of
equation(3.4) with the following initial condition:

$\begin{split}f_{h}(\mathbf{p},0)&=f_{h}^{0}(\mathbf{p})\\\
W_{h}(\mathbf{k},0)&=W_{h}^{0}(\mathbf{k})\end{split}$

If $W_{h}$ always takes non-negative values, i.e. $W_{h}(\mathbf{k},t)\geq
0,\forall\leavevmode\nobreak\ \mathbf{k}\in\Omega^{L}_{k},\forall t\geq 0$,
then $f_{h}$ has $L^{2}$ stability

$\|f_{h}\|_{L^{2}(\Omega^{L}_{p})}\leq\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}$

and $W_{h}$ has bounded $L^{1}$ norm.

$\|W_{h}\|_{L^{1}(\Omega^{L}_{k})}\leq\mathcal{E}_{tot,h}^{0}+\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}\cdot\|\Pi_{p,h}E\|_{L^{2}(\Omega^{L}_{p})}$

###### Proof.

Since $f_{h}$ belongs to the test space $\mathcal{G}_{h}^{\alpha_{1}}$, we
test the equation for particles with $f_{h}$, we obtain that

$(\frac{\partial f_{h}}{\partial
t},f_{h})_{p}=-\int_{\Omega^{L}_{k}}d\mathbf{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\\{(\mathcal{L}_{h}f_{h})^{2}W_{h}\mathcal{B}\\}.$

The right hand side is non-positive as long as $W_{h}$ always take non-
negative values, therefore the $L^{2}$ norm of $f_{h}$ always decreases,

$\frac{1}{2}\frac{\partial}{\partial
t}\|f_{h}\|_{L^{2}(\Omega^{L}_{p})}^{2}\leq
0\Rightarrow\|f_{h}\|_{L^{2}(\Omega^{L}_{p})}\leq\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}$

Now consider $W_{h}$, by definition,

$\|W_{h}\|_{L^{1}(\Omega^{L}_{k})}=(W_{h},\text{sgn}(W_{h}))_{k}=(W_{h},1)_{k}$

Recall the energy conservation property in Theorem 3:

$\|W_{h}\|_{L^{1}(\Omega^{L}_{k})}+(f_{h},\Pi_{p,h}E)_{p}=\mathcal{E}_{tot,h}^{0}$

Use Holder’s inequality

$\begin{split}\|W_{h}\|_{L^{1}(\Omega^{L}_{k})}&=\mathcal{E}_{tot,h}^{0}-(f_{h},\Pi_{p,h}E)_{p}\\\
&\leq\mathcal{E}_{tot,h}^{0}+|(f_{h},\Pi_{p,h}E)_{p}|\\\
&\leq\mathcal{E}_{tot,h}^{0}+\|f_{h}\|_{L^{2}(\Omega^{L}_{p})}\cdot\|\Pi_{p,h}E\|_{L^{2}(\Omega^{L}_{p})}\end{split}$

By the $L^{2}$-stability of $f_{h}$, we obtain the upper bound of $W_{h}$’s
$L^{1}$ norm,

$\|W_{h}\|_{L^{1}(\Omega^{L}_{k})}\leq\mathcal{E}_{tot,h}^{0}+\|f_{h}\|_{L^{2}(\Omega^{L}_{p})}\cdot\|E\|_{L^{2}(\Omega^{L}_{p})}\leq\mathcal{E}_{tot,h}^{0}+\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}\cdot\|\Pi_{p,h}E\|_{L^{2}(\Omega^{L}_{p})}$

∎

###  5.2 Time Discretization

Recall our conservative semi-discrete weak form,

$\begin{split}(\frac{\partial f_{h}}{\partial
t},\phi_{h})_{p}&=-B_{L}(f_{h},W_{h},\phi_{h})\coloneqq-\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\mathcal{L}_{h}\phi_{h}\mathcal{L}_{h}f_{h}W_{h}\mathcal{B}_{h}\\},\\\
(\frac{\partial W_{h}}{\partial
t},\eta_{h})_{k}&=H_{L}(W_{h},f_{h},\eta_{h})\coloneqq\iint_{\Omega^{L}_{k}\times\Omega^{L}_{p}}d\mathbf{k}d\mathbf{p}\\{\eta_{h}\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}W_{h}\mathcal{B}_{h}\\},\end{split}$

The time step size of the explicit scheme for diffusion equations is
restricted by the CFL condition. Two reasons urge us to avoid explicit
schemes,

  1. 1.

The CFL bound of step size may be too restrictive, and we might lose
efficiency.

  2. 2.

The upper bound depends on the eigenvalues of time-varying diffusion
coefficients. However, in the proposed scheme, we never calculate the
diffusion coefficient explicitly, instead, we compute the interaction tensor
$B$ associated with the trilinear form $B$.

On the other hand, due to nonlinearity, a fully implicit scheme requires
fixed-point iteration involving both particle pdf $f_{h}$ and wave sed
$W_{h}$, which can be time-consuming. Therefore, the objective is to find a
scheme that is only implicit for $f_{h}$, and at the same time preserves
discrete conservation laws.

We propose the following semi-implicit scheme,

$\begin{split}(\frac{f_{h}^{s+1}-f_{h}^{s}}{\Delta
t},\phi_{h})_{p}&+B_{L}(f_{h}^{s+1},W_{h}^{s},\phi_{h})=0\\\
(\frac{W_{h}^{s+1}-W_{h}^{s}}{\Delta
t},\eta_{h})_{k}&-H_{L}(W_{h}^{s},f_{h}^{s+1},\eta_{h})=0\end{split}$ (5.1)

The scheme is implicit for particle pdf $f_{h}$ if we focus on the first line,
meanwhile it is explicit for wave sed $W_{h}$, considering the second line.
For implementation, we solve the first row and then substitute the next step
particle pdf $f^{s+1}_{h}$ into the second row. It can be easily verified that
the discrete conservation laws still hold, i.e. we have

$(f_{h}^{s+1},\phi_{c,h})_{p}+(W_{h}^{s+1},\eta_{c,h})_{k}=(f_{h}^{s},\phi_{c,h})_{p}+(W_{h}^{s},\eta_{c,h})_{k}$

The following theorem is the fully discrete version of Lemma 1, giving the
unconditional $L^{2}$-stability of $f^{s}_{h}$ when $W^{s}_{h}$ is non-
negative.

######  Theorem 4.

Suppose $f^{s}_{h}(\mathbf{p})$ and $W^{s}_{h}(\mathbf{k})$ are the solution
of Equation(5.1).

If $W^{s}_{h}$ always takes non-negative values, i.e.
$W^{s}_{h}(\mathbf{k})\geq 0,\forall\leavevmode\nobreak\
\mathbf{k}\in\Omega^{L}_{k},\forall s\geq 0$, then $f^{s}_{h}$ has $L^{2}$
stability

$\|f^{s}_{h}\|_{L^{2}(\Omega^{L}_{p})}\leq\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}$

and $W^{s}_{h}$ has bounded $L^{1}$ norm.

$\|W^{s}_{h}\|_{L^{1}(\Omega^{L}_{k})}\leq\mathcal{E}_{tot,h}^{0}+\|f_{h}^{0}\|_{L^{2}(\Omega^{L}_{p})}\cdot\|\Pi_{p,h}E\|_{L^{2}(\Omega^{L}_{p})}$

###### Proof.

Given that $W_{h}^{s}(\mathbf{k})\geq 0,\forall\mathbf{k}\in\Omega^{L}_{k}$,
we have $B_{L}(f_{h}^{s+1},W_{h}^{s},f_{h}^{s+1})\geq 0$. Therefore,
$f^{s}_{h}$ has unconditional $L^{2}$-stability.

Since the scheme(5.1) preserves energy conservation, the $L^{1}$ bound of
$W_{h}$ can be proved in the same approach as we have done in Lemma 1. ∎

Note that the stability depends on our assumption that $W^{s}_{h}$ is non-
negative. Therefore, in what follows, we will discuss the positivity-
preserving technique of $W^{s}_{h}$.

###  5.3 Positivity-Preserving Technique for the Wave SED

To ensure positivity of wave sed $W^{s}_{h}$, we draw the strategy from Zhang
et al.[22]:

  1. 1.

Use a small enough time step to ensure positive cell-average of a temporary
wave sed $W^{s+1,*}_{h}$, given that we have pointwise positivity of last step
wave sed $W^{s}_{h}$.

  2. 2.

Apply a slope limiter on $W^{s+1,*}_{h}$ which preserves cell-average at the
same time, then we obtain a pointwise positive $W^{s+1}_{h}$ as our solution
of the next step wave sed. (Obviously, if we use piecewise constant basis
functions, this step is not necessary).

Firstly we will derive the constraint on time step size. After that, we
explain why the slope limiter will not break discrete conservation laws.

Suppose $\eta_{j,0}$ is the characteristic function of the $j$-th element
$R_{k}^{j}\subset\Omega^{L}_{k}$, i.e. $\eta_{j,0}=1_{\mathbf{k}\in
R_{k}^{j}}$, which belongs to the test space $\mathcal{W}_{h}^{\alpha_{2}}$.
According to the time discretization in Equation(5.1),

$\begin{split}(\frac{W_{h}^{s+1,*}-W_{h}^{s}}{\Delta
t},\eta_{j,0})_{k}=\int_{\Omega^{L}_{k}}\frac{W_{h}^{s+1,*}-W_{h}^{s}}{\Delta
t}\eta_{j,0}d\mathbf{k}=\int_{\Omega^{L}_{k}}d\mathbf{k}\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}W_{h}^{s}\eta_{j,0}\mathcal{B}_{h}\\}\end{split},$

which is equivalent to

$\int_{R_{k}^{j}}W_{h}^{s+1,*}d^{3}k=\int_{R_{k}^{j}}W_{h}^{s}(1+\Delta
t\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}\mathcal{B}_{h}\\})d\mathbf{k}$

To ensure positive cell-average, i.e.
$\int_{R_{k}^{j}}W_{h}^{s+1,*}d\mathbf{k}\geq 0$, we require that there exists
a constant $\epsilon>0$ such that

$1+\Delta
t\int_{\Omega^{L}_{p}}d\mathbf{p}\\{\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}\mathcal{B}_{h}\\}\geq\epsilon,\
\forall\mathbf{k}\in R_{k}^{j}$ (5.2)

As long as the time step size $\Delta t$ satisfy condition(5.2), we have
$\int_{R_{k}^{j}}W_{h}^{s+1,*}d\mathbf{k}\geq\epsilon\int_{R_{k}^{j}}W_{h}^{s}d\mathbf{k}\geq
0$.

The following theorem guarantees that our bound for $\Delta t$ will not shrink
over time.

######  Theorem 5.

For any $\epsilon>0$, given a regular enough discrete emission/absorption
kernel $\mathcal{B}_{h}$, there exists a constant $\Delta t_{M}$ determined by
$\epsilon$, $f^{0}$, $\Omega^{L}_{p}$, $\Omega^{L}_{k}$ and $h$, such that any
$\Delta t<\Delta t_{M}$ satisfies condition (5.2).

###### Proof.

By Hölder’s inequality, the ”growth rate” is bounded as follows,

$\begin{split}\left|\int_{\Omega_{p}^{L}}d\mathbf{p}\\{\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}\mathcal{B}_{h}\\}\right|\leq&\int_{\Omega_{p}^{L}}d\mathbf{p}\left|\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}\mathcal{B}_{h}\right|\\\
\leq&\left\|\frac{p_{\perp}}{p}\mathcal{L}_{h}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}\cdot\left\|\frac{p}{p_{\perp}}\mathcal{B}_{h}\mathcal{L}_{h}E_{h}\right\|_{L^{1}(\Omega_{p}^{L})}\\\
=&2\pi\left\|\frac{p_{\perp}}{p}\mathcal{L}_{h}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}\cdot\left(\int
p\mathcal{B}_{h}\mathcal{L}_{h}E_{h}dp_{\perp}dp_{\parallel}\right)\end{split}$
(5.3)

Firstly, consider the $L^{\infty}$-norm factor from inequality (5.3). Recall
the definition of $\mathcal{L}_{h}$,

$\frac{p_{\perp}}{p}\mathcal{L}_{h}f_{h}^{s+1}\coloneqq
N_{\parallel,h}\frac{\partial E_{h}}{\partial p_{\perp}}\frac{\partial
f_{h}^{s+1}}{\partial p_{\parallel}}+(1-N_{\parallel,h}\frac{\partial
E_{h}}{\partial p_{\parallel}})\frac{\partial f_{h}^{s+1}}{\partial
p_{\perp}}.$

Both of the coefficients $N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\perp}}$ and $(1-N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\parallel}})$ are bounded by some constant $C$ dependent on
$\Omega_{p}^{L}$ and $\Omega_{k}^{L}$, hence it follows that,

$\begin{split}\left\|\frac{p_{\perp}}{p}\mathcal{L}_{h}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}&=\left\|N_{\parallel,h}\frac{\partial
E_{h}}{\partial p_{\perp}}\frac{\partial f_{h}^{s+1}}{\partial
p_{\parallel}}+(1-N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\parallel}})\frac{\partial f_{h}^{s+1}}{\partial
p_{\perp}}\right\|_{L^{\infty}(\Omega_{p}^{L})}\\\ &\leq
C(\Omega_{p}^{L},\Omega_{k}^{L})\left\|\frac{\partial f_{h}^{s+1}}{\partial
p_{\parallel}}\right\|_{L^{\infty}(\Omega_{p}^{L})}+C(\Omega_{p}^{L},\Omega_{k}^{L})\left\|\frac{\partial
f_{h}^{s+1}}{\partial p_{\perp}}\right\|_{L^{\infty}(\Omega_{p}^{L})}\\\ &\leq
C_{1}(\Omega_{p}^{L},\Omega_{k}^{L})\cdot\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})},\end{split}$
(5.4)

We claim that
$\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}$ is bounded
uniformly in time. Indeed, since the domain $\Omega_{p}^{L}$ is finite, all
$L^{r}$ norms are equivalent, therefore,

$\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}\leq
C_{2}(\Omega_{p}^{L})\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{2}(\Omega_{p}^{L})}.$

Moreover, the inverse inequality for finite element spaces,

$\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{2}(\Omega_{p}^{L})}\leq\frac{C_{3}}{h_{p}}\left\|f_{h}^{s+1}\right\|_{L^{2}(\Omega_{p}^{L})},$

and the $L^{2}$ stability estimate from Theorem 4

$\left\|f_{h}^{s+1}\right\|_{L^{2}(\Omega_{p}^{L})}\leq\left\|f_{h}^{0}\right\|_{L^{2}(\Omega_{p}^{L})},$

leads to the following estimate for $\nabla_{p}f_{h}^{s+1}$,

$\left\|\nabla_{p}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}\leq
C_{2}(\Omega_{p}^{L})\frac{C_{3}}{h_{p}}\left\|f_{h}^{0}\right\|_{L^{2}(\Omega_{p}^{L})}.$
(5.5)

Therefore, the $L^{\infty}$-norm factor from inequality (5.3) is bounded as
follows,

$\left\|\frac{p_{\perp}}{p}\mathcal{L}_{h}f_{h}^{s+1}\right\|_{L^{\infty}(\Omega_{p}^{L})}\leq
C_{1}(\Omega_{p}^{L},\Omega_{k}^{L})C_{2}(\Omega_{p}^{L})\frac{C_{3}}{h_{p}}\left\|f_{h}^{0}\right\|_{L^{2}(\Omega_{p}^{L})}.$

Next, consider the $L^{1}$-norm factor from inequality (5.3), and write it as
follows,

$\int
p\mathcal{B}_{h}\mathcal{L}_{h}E_{h}dp_{\perp}dp_{\parallel}=\int\left[pU_{l}(\mathbf{p};\mathbf{k})\delta_{h}(\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p}))\frac{p}{p_{\perp}}\frac{\partial
E_{h}}{\partial p_{\perp}}\right]dp_{\perp}dp_{\parallel},$

where $\delta_{h}$ represents an approximation of Dirac delta, see Equation
(2.8) and 2.9.

We discuss the following two cases,

  * •

When $l\neq 0$, since

$\lim_{p_{\perp}\rightarrow
0^{+}}U_{l}(\mathbf{p};\mathbf{k})\frac{p^{2}}{p_{\perp}}=\lim_{p_{\perp}\rightarrow
0^{+}}8\pi^{2}e^{2}\frac{\left(iE_{2}J_{l}^{\prime}\right)^{2}\left(\frac{p_{\perp}}{p}\right)^{2}}{\left(1-E_{2}^{2}\right)\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}\varepsilon)+2iE_{2}\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}g)+E_{3}^{2}\frac{1}{\omega}\frac{\partial}{\partial\omega}(\omega^{2}\eta)}\frac{p^{2}}{p_{\perp}}=0,$

the integral is bounded as follows,

$\begin{split}\int
p\mathcal{B}_{h}\mathcal{L}_{h}E_{h}dp_{\perp}dp_{\parallel}&=\int\left[\left(U_{l}(\mathbf{p};\mathbf{k})\frac{p^{2}}{p_{\perp}}\right)\frac{\partial
E_{h}}{\partial
p_{\perp}}\delta_{h}(\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p}))\right]dp_{\perp}dp_{\parallel}\\\
&\leq\sup_{\Omega_{p}^{L}}\left(U_{l}(\mathbf{p};\mathbf{k})\frac{p^{2}}{p_{\perp}}\right)\frac{\partial
E_{h}}{\partial
p_{\perp}}\int\left[\delta_{h}(\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p}))\right]dp_{\perp}dp_{\parallel}\\\
&\leq C_{4}(\Omega_{p}^{L},\Omega_{k}^{L}).\end{split}$ (5.6)

  * •

When $l=0$, the above trick does not work, because $\lim_{p_{\perp}\rightarrow
0^{+}}U_{0}(\mathbf{p};\mathbf{k})>0$. For this special case, as an
alternative to the original operator

$\mathcal{L}_{h}g\coloneqq N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\perp}}\frac{p}{p_{\perp}}\frac{\partial g}{\partial
p_{\parallel}}+(1-N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\parallel}})\frac{p}{p_{\perp}}\frac{\partial g}{\partial p_{\perp}},$

we adopt a new discrete operator,

$\mathcal{L}^{0}_{h}g\coloneqq N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\perp}}\frac{p}{E_{h}\frac{\partial E_{h}}{\partial
p_{\perp}}}\frac{\partial g}{\partial
p_{\parallel}}+(1-N_{\parallel,h}\frac{\partial E_{h}}{\partial
p_{\parallel}})\frac{p}{E_{h}\frac{\partial E_{h}}{\partial
p_{\perp}}}\frac{\partial g}{\partial p_{\perp}}.$

The operator is still a consistent discretization since
$E=\sqrt{1+p_{\parallel}^{2}+p_{\perp}^{2}}$. It can also be easily verified
that the $L^{\infty}$ bound in inequality (5.4) is still true with this new
operator $\mathcal{L}^{0}_{h}$.

In addition, the $L^{1}$-norm factor becomes,

$\int
p\mathcal{B}_{h}\mathcal{L}_{h}^{0}E_{h}dp_{\perp}dp_{\parallel}=\int\left[pU_{l}(\mathbf{p};\mathbf{k})\delta_{h}(\omega(\mathbf{k})-k_{\parallel}v_{\parallel}-l\omega_{c}/\gamma(\mathbf{p}))\frac{p}{E_{h}}\right]dp_{\perp}dp_{\parallel}.$

Therefore inequality (5.6) still holds.

Combine inequalities (5.4), (5.5), and (5.6) to obtain

$\left|\int_{\Omega_{p}^{L}}d\mathbf{p}\\{\mathcal{L}_{h}E_{h}\mathcal{L}_{h}f_{h}^{s+1}\mathcal{B}_{h}\\}\right|\leq
2\pi\frac{C_{1}\cdot C_{2}\cdot C_{3}\cdot
C_{4}}{h_{p}}\left\|f_{h}^{0}\right\|_{L^{2}(\Omega_{p}^{L})},$

which enables us to define the uniform-in-time upper bound,

$\Delta t_{M}\coloneqq\frac{(1-\epsilon)h_{p}}{2\pi\cdot C_{1}\cdot C_{2}\cdot
C_{3}\cdot C_{4}\left\|f_{h}^{0}\right\|_{L^{2}(\Omega_{p}^{L})}}.$

It can be easily verified that any $\Delta t<\Delta t_{M}$ satisfies condition
(5.2). ∎

The condition does not need to be calculated explicitly, because we can adapt
time step size a posteriori in the code: monitor the cell averages, if any
cell average of the temporary solution $W_{h}^{s+1,*}$ is non-positive,
replace $\Delta t$ with $0.5\Delta t$ and calculate $W_{h}^{s+1,*}$ again.

Now let us discuss the effect of slope limiters on conservation laws. If
$\alpha_{2}=0$, there is no need for any slope limiter. If $\alpha_{2}=1$, we
apply the slope limiter $\theta$ and obtain
$W_{h}^{s+1}=\theta(W^{s+1,*}_{h})$. According to Zhang et al.[22], the cell
average is preserved, i.e.
$\int_{R_{k}^{j}}W_{h}^{s+1,*}d\mathbf{k}=\int_{R_{k}^{j}}W_{h}^{s+1}d\mathbf{k}$.
In other words,
$\left(W_{h}^{s+1,*},\eta\right)_{k}=\left(W_{h}^{s+1},\eta\right)_{k}$, for
any piecewise constant test function, i.e.
$\forall\eta\in\mathcal{W}_{h}^{0}$. Therefore, to preserve discrete
conservation laws, in the definition of the discrete directional differential
operator $\mathcal{L}_{h}$, we need to pick a projection $\Pi_{k,h}$ such that
$\Pi_{k,h}U$ belongs to
$\mathcal{W}_{h}^{0}\subset\mathcal{W}_{h}^{\alpha_{2}}$ for any function $U$.

##  6 Numerical Results

###  6.1 Problem Setting

Although the emission/absorption kernel contains a summation from $l=-\infty$
to $l=+\infty$, it is not practical to perform that numerically. In practice,
we keep the dominant part of those terms. In the following example, we will
only consider one term with $l=1$, associated with the anomalous Doppler
resonance. We used the dispersion relation $\omega(\mathbf{k})$ of the
whistler wave in cold magnetized plasma(see Appendix), with electron gyro-
frequency $\omega_{c}=-2\omega_{p}$.

Set the cut-off computational domain as follows,

$\begin{split}\Omega^{L}_{p}&=\\{(p_{\parallel},p_{\perp}):p_{\parallel}\in(-5mc,25mc),\
p_{\perp}\in(0,15mc)\\}\\\
\Omega^{L}_{k}&=\\{(k_{\parallel},k_{\perp}):k_{\parallel}\in(0.05\frac{\omega_{p}}{c},0.65\frac{\omega_{p}}{c}),\
k_{\perp}\in(0,0.6\frac{\omega_{p}}{c})\\}\end{split}$

Take piecewise linear quadrilateral basis
$\mathcal{G}_{h}^{1}=\\{f(p_{\parallel},p_{\perp})\in
C^{0}(\Omega^{L}_{p}):f|_{R_{p}}\in Q^{1}(R_{p}),\forall\leavevmode\nobreak\
R_{p}\in\mathcal{T}_{h}^{p}\\}$ and piecewise constant basis
$\mathcal{W}_{h}^{0}=\\{W(k_{\parallel},k_{\perp}):W|_{R_{k}}\in
Q^{0}(R_{k}),\forall\leavevmode\nobreak\ R_{k}\in\mathcal{T}_{h}^{k}\\}$ as
our test spaces. Choose the $L^{2}$ orthogonal projections $\Pi_{p,h}$ and
$\Pi_{k,h}$ as stated in Corollary 1.

The numerical experiment is performed with $75\times 75$ elements in
$\Omega_{p}^{L}$, and $40\times 40$ elements in $\Omega_{k}^{L}$. The initial
time step size is set as $\Delta t=1.0\times 10^{4}\frac{1}{2\pi\omega_{p}}$.

The integration on resonance manifold is performed with Gauss-Legendre
quadrature on $\Omega^{L}_{p}$ and the marching simplex method[7, 12] on
$\Omega^{L}_{k}$.

Consider the following initial conditions, which is the so-called ’bump on
tail instability’ configuration.

$\begin{cases}\begin{aligned}
f(p_{\parallel},p_{\perp})|_{t=0}&=\left[10^{-5}\frac{1}{\sqrt{\pi}}\exp\left(-\left(\frac{p_{\parallel}}{mc}-20\right)^{2}-\left(\frac{p_{\perp}}{mc}\right)^{2}\right)\right]\frac{n_{0}}{m^{3}c^{3}}\\\
W(k_{\parallel},k_{\perp})|_{t=0}&=10^{-5}\frac{n_{0}mc^{2}}{(\omega_{p}/c)^{3}}\end{aligned}\end{cases}$

######  Remark 8.

The bump on tail configuration actually refers to the sum of a bulk and a
bump, i.e. $f(\mathbf{p},t)=f_{c}(\mathbf{p})+f_{b}(\mathbf{p},t)$, where the
cold bulk
$f_{c}(\mathbf{p})\approx\frac{n_{0}}{m^{3}c^{3}}\delta(\mathbf{p})$, and the
bump $f_{b}(\mathbf{p},t)$ is a peak with a much smaller population, centered
far from the origin. However, as shown in the following equation,

$\begin{split}&\frac{\partial(f-f_{c})}{\partial
t}=\mathbb{B}(W,f)=\mathbb{B}(W,f-f_{c})+\mathbb{B}(W,f_{c})=\mathbb{B}(W,f-f_{c}),\\\
&\frac{\partial W}{\partial
t}=\mathbb{H}(f,W)=\mathbb{H}(f-f_{c},W)+\mathbb{H}(f_{c},W)=\mathbb{H}(f-f_{c},W),\end{split}$

we do not have to really compute the contribution from $f_{c}$.

###  6.2 Temporal Evolution

In analogous to the analysis done by Kennel and Engelmann [10], for a given
wave vector $\mathbf{k}$, the characteristics associated with directional
differential operator $\mathcal{L}$ is

$z(p_{\parallel},p_{\perp})\coloneqq\frac{\omega}{k_{\parallel}}p_{\parallel}-E(p_{\parallel},p_{\perp})=\frac{\omega}{k_{\parallel}}p_{\parallel}-\sqrt{m^{2}c^{4}+p_{\parallel}^{2}c^{2}+p_{\perp}^{2}c^{2}}=\text{const},$
(6.1)

which is the isoenergy contour in the reference frame moving at the wave’s
phase velocity.

When the wave sed $W(\mathbf{k},t)$ is concentrated around the given
$\mathbf{k}$, the contours as illustrated in Figure(2) indicates the principal
diffusion direction. For the specific problem setting,
$\frac{\omega}{k_{\parallel}}$ is small, hence the contour lines are almost
concentric circles.

In Figure(1) we show the evolution of electron pdf
$f(p_{\parallel},p_{\perp},t)$ and wave sed $W(k_{\parallel},k_{\perp},t)$. It
can be observed that the bump on tail results in the excitation of the
approximate waves in a narrow region of spectral space $\Omega^{L}_{k}$, and
as predicted by Equation(6.1), the whistler waves in turn cause anisotropic
diffusion of electron pdf almost along the contour lines in Figure(2).

Figure 1: Temporal evolution of the electron pdf and wave sed. Figure 2: The
characteristics of directional differential operator $\mathcal{L}$ given
$(k_{\parallel},k_{\perp})=(0.2,0.3)$.

###  6.3 Verification of Conservation

To verify the discrete conservation property of the proposed scheme, we define
the relative error for conserved quantity as follows,

$e_{rel}\left(\mathcal{Q}_{tot,h}\right)\coloneqq\frac{\|\mathcal{Q}_{tot,h}-\mathcal{Q}_{tot,h}^{0}\|_{L^{\infty}(0,T_{max})}}{\mathcal{Q}_{tot,h}^{0}},$

where $\mathcal{Q}$ is the conserved quantities defined in Theorem 2.

Then with $T_{max}=8.0\times 10^{7}\frac{1}{2\pi\omega_{p}}$, we have

$\begin{split}e_{rel}(\mathcal{M}_{tot,h})&=4.96\times 10^{-14},\\\
e_{rel}(\mathcal{P}_{\parallel,tot,h})&=4.58\times 10^{-14},\\\
e_{rel}(\mathcal{E}_{tot,h})&=4.69\times 10^{-14}.\end{split}$

For the evolution of the electron-plasmon system momentum and energy, see
solid lines in Figure(3).

###  6.4 Comparison of Different Dispersion Relations

The above results were obtained with the exact whistler wave dispersion
relation $\omega_{\text{imp}}(\mathbf{k})$ for cold magnetized plasma, given
implicitly by Equation(8.2). One might wonder what if we replace it with a
simpler explicit approximate relation, for instance,

$\omega_{\text{exp}}(\mathbf{k})=|\omega_{c}|\frac{|k_{\parallel}|kc^{2}}{\omega^{2}_{p}\sqrt{1+k^{2}c^{2}\omega_{c}^{2}/\omega_{p}^{4}}}$

which is asymptotic to the implicit relation $\omega_{\text{imp}}(\mathbf{k})$
when $k$ is small, i.e.

$\lim_{k\rightarrow
0}\left(\omega_{\text{exp}}(\mathbf{k})-\omega_{\text{imp}}(\mathbf{k})\right)=0.$

(a) Momentum evolution.

(b) Energy evolution.

Figure 3: Comparison between $\omega_{\text{exact}}$ and
$\omega_{\text{approx}}$

As shown in Figure(3), for both cases, energy and momentum are transferred
from particles to waves. Meanwhile, we do observe a different transfer rate
for the approximate whistler dispersion relation when compared to the exact
implicit dispersion relation derived from Equation(8.2) in the Appendix.

##  7 Summary

We studied the numerical method for the initial value problem associated with
the relativistic quasilinear diffusion model in magnetized plasma. We showed
that a conservative semi-discrete form can be derived by adopting a novel
integro-differential form of this wave-particle interaction system. We
incorporated the marching simplex algorithm in numerical integration on the
resonance manifold. A semi-implicit time discretization was introduced to
ensure the stability of the particle pdf and positivity of wave sed, which
also preserves conservation in the fully discrete form. In the end, we
presented our numerical results for the bump-on-tail instability, and the
conservation properties are verified.

In the future, we will consider the problem with the spatial non-uniform
setting, and other factors such as Landau collision operator and external
electric field will be included. Error estimates for the Galerkin scheme will
also be investigated.

##  8 Appendix : Waves in Cold Magnetized Plasma

The directional differential operator and the emission/absorption kernel both
vary for different wave modes. Since our numerical experiment is based on the
whistler mode, we introduce the wave modes in cold magnetized plasma here to
make the paper self-contained.

Inside the medium with conductivity tensor $\sigma$, we have the linear
relation between current and field $\mathbf{J}=\sigma\cdot\mathbf{E}$. Apply
the Fourier transform to the Maxwell equations, we obtain that

$\mathbf{k}\times\mathbf{k}\times\hat{\mathbf{E}}=\mu_{0}\varepsilon_{0}\omega^{2}(\frac{i\mathbf{\sigma}}{\omega\varepsilon_{0}}+I)\cdot\hat{\mathbf{E}}$

Write it with Einstein’s summation notation, we have the wave equation in
spectral form,

$(k_{\alpha}k_{\beta}-\delta_{\alpha\beta}k^{2}+\frac{\omega^{2}}{c^{2}}\varepsilon_{\alpha\beta})\hat{E}_{\beta}=0,\beta=1,2,3$
(8.1)

where the dielectric tensor $\varepsilon_{\alpha\beta}$ is dependent on
$\sigma$.

Let
$M(\omega,\mathbf{k})=k_{\alpha}k_{\beta}-\delta_{\alpha\beta}k^{2}+\frac{\omega^{2}}{c^{2}}\varepsilon_{\alpha\beta}$,
then the wave equation(8.1) has nontrivial solution if and only if

$Det[M(\omega,\mathbf{k})]=0$ (8.2)

Equation(8.2) gives the graph of implicit function $\omega(\mathbf{k})$, which
is known as the dispersion relation.

The above discussion works for any medium. Now we focus on the plasma.
Consider cold magnetized plasma with background field
$\mathbf{B}_{0}(\mathbf{x})=B_{0}\mathbf{b}$, where $B_{0}$ is constant and
$\mathbf{b}$ is a fixed unit vector. A plasma is ”cold” when the waves
propagate faster than its thermal speed.

The following dielectric tensor for cold magnetized plasmas can be found in
textbooks[16, 17]. (For simplicity, ion motion is neglected here.)

$\varepsilon_{\alpha\beta}(\omega)\equiv\varepsilon\delta_{\alpha\beta}+ige_{\alpha\beta\gamma}b_{\gamma}+(\eta-\varepsilon)b_{\alpha}b_{\beta},$

where

$\begin{split}\varepsilon&=\varepsilon^{H}\equiv
1-\frac{\omega_{pe}^{2}}{\omega^{2}-\omega_{ce}^{2}},\\\
g&=g^{H}\equiv-\frac{\omega_{ce}}{\omega}\frac{\omega_{pe}^{2}}{\omega^{2}-\omega_{ce}^{2}},\\\
\eta&=\eta^{H}\equiv 1-\frac{\omega_{pe}^{2}}{\omega^{2}}.\end{split}$ (8.3)

There are two parameters in the above formulas. The electron gyro-frequency
$\omega_{ce}=-|e|B/m_{e}c$ is proportional to the background magnetic field.
The plasma frequency $\omega_{pe}=\sqrt{4\pi n_{e}e^{2}/m_{e}}$ is
proportional to the square root of particle density.

Recall that the dispersion relation $\omega(\mathbf{k})$ is given implicitly
in Equation(8.2), one might wonder whether $\omega(\mathbf{k})$ is multi-
valued, and for a specific branch, whether it is well-defined for any
$\mathbf{k}\in\Omega_{k}$. Textbooks never elaborate on this issue, therefore
we provide an answer here.

Define the parallel component and perpendicular component of wave vector
$\mathbf{k}$, $k_{\parallel}\coloneqq\mathbf{k}\cdot\mathbf{B}/|\mathbf{B}|$,
$k_{\perp}\coloneqq\sqrt{\mathbf{k}^{2}-k_{\parallel}^{2}}$. Denote the
magnitude $\mathbf{k}$ as $k$, and define the cosine of polar angle as
$\xi\coloneqq{k_{\parallel}/k}$. If ion motion is neglected, then we have what
follows.

######  Proposition 1.

$\forall\leavevmode\nobreak\
(k_{\parallel},k_{\perp})\in(0,+\infty)\times(0,+\infty)$,
$\exists\leavevmode\nobreak\
0<\omega_{1}<\omega_{2}<\omega_{3}<\omega_{4}<\infty$ s.t.

$Det[M(\omega_{j},k_{\parallel},k_{\perp})]=0,\leavevmode\nobreak\ j=1,2,3,4$

i.e. the equation admits exactly $4$ positive single-value implicit functions
$\omega_{j}(k_{\parallel},k_{\perp})$, $j=1,2,3,4$, on domain
$(0,+\infty)\times(0,+\infty)$, moreover,
$\omega_{j}(k,\xi)\coloneqq\omega_{j}(k_{\parallel}(k,\xi),k_{\perp}(k,\xi))$
satisfy that $\frac{\partial}{\partial k}\omega_{j}(k,\xi)\geq
0,\leavevmode\nobreak\ \forall\leavevmode\nobreak\ \xi\in(0,1)$.

The first branch $\omega_{1}(\mathbf{k})$ is defined on the whole spectral
space $\mathbb{R}^{3}_{k}$. In a relatively strong magnetic
field($\omega_{p}^{2}/\omega_{c}^{2}\ll 1$), the whistler wave actually refers
to waves with wave vector
$\mathbf{k}\in\Omega_{w}\subsetneq\mathbb{R}^{3}_{k}$, and meanwhile has
frequency $\omega_{1}(\mathbf{k})$. Outside the region $\Omega_{w}$, the first
branch has another name. For details, see Aleynikov et al.[1]. In our
numerical experiment, the cut-off domain $\Omega^{L}_{k}\subset\Omega_{w}$,
therefore we say that we consider the whistler wave.

The polarization vector components in the emission/absorption kernel are given
below

$\begin{split}E_{1}(\mathbf{k})&=1,\\\
E_{2}(\mathbf{k})&=i\frac{g}{\varepsilon-N^{2}},\\\
E_{3}(\mathbf{k})&=-\frac{N_{\parallel}N_{\perp}}{\eta-
N_{\perp}^{2}},\end{split}$ (8.4)

where $\mathbf{N}=\frac{\mathbf{k}}{\omega}$ is the refractive index.

Acknowledgements. The authors thank and gratefully acknowledge the support
from the Oden Institute of Computational Engineering and Sciences and the
University of Texas Austin. This project was supported by funding from NSF
DMS: 2009736 and DOE DE-SC0016283 project Simulation Center for Runaway
Electron Avoidance and Mitigation.

## References

  * [1] Pavel Aleynikov and Boris Breizman.  Stability analysis of runaway-driven waves in a tokamak.  Nuclear Fusion, 55(4):043014, April 2015.
  * [2] Claude Bardos and Nicolas Besse.  Diffusion limit of the vlasov equation in the weak turbulent regime.  Journal of Mathematical Physics, 62(10):101505, 2021.
  * [3] Nicolas Besse, Yves Elskens, DF Escande, and Pierre Bertrand.  Validity of quasilinear theory: refutations and new numerical confirmation.  Plasma Physics and Controlled Fusion, 53(2):025012, 2011.
  * [4] Boris N. Breizman, Pavel Aleynikov, Eric M. Hollmann, and Michael Lehnen.  Physics of runaway electrons in tokamaks.  Nuclear Fusion, 59(8):083001, August 2019.
  * [5] Akio Doi and Akio Koide.  An efficient method of triangulating equi-valued surfaces by using tetrahedral cells.  IEICE TRANSACTIONS on Information and Systems, 74(1):214–224, 1991.
  * [6] WE Drummond and D Pines.  Non-linear stability of plasma oscillations.  1962\.
  * [7] André Guéziec and Robert Hummel.  Exploiting triangulated surface extraction using tetrahedral decomposition.  IEEE Transactions on visualization and computer graphics, 1(4):328–342, 1995.
  * [8] Kun Huang and Irene M. Gamba.  Existence of global weak solutions to quasilinear theory for electrostatic plasmas, 2023.
  * [9] Allan N Kaufman.  Resonant interactions between particles and normal modes in a cylindrical plasma.  The Physics of Fluids, 14(2):387–397, 1971.
  * [10] C. F. Kennel.  Velocity Space Diffusion from Weak Plasma Turbulence in a Magnetic Field.  Physics of Fluids, 9(12):2377, 1966.
  * [11] I Lerche.  Quasilinear theory of resonant diffusion in a magneto-active, relativistic plasma.  The Physics of Fluids, 11(8):1720–1727, 1968.
  * [12] Chohong Min and Frédéric Gibou.  Geometric integration over irregular domains with application to level-set methods.  Journal of Computational Physics, 226(2):1432–1443, 2007.
  * [13] Joachim Nitsche.  Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind.  In Abhandlungen aus dem mathematischen Seminar der Universität Hamburg, volume 36, pages 9–15. Springer, 1971.
  * [14] MN Rosenbluth and SV Putvinski.  Theory for avalanche of runaway electrons in tokamaks.  Nuclear fusion, 37(10):1355, 1997.
  * [15] VD Shapiro and VI Shevchenko.  On the nonlinear theory of interaction between charged particle beams and a plasma in a magnetic field.  Zhur. Eksptl’. i Teoret. Fiz., 42, 1962.
  * [16] Thomas H Stix.  Waves in plasmas.  Springer Science & Business Media, 1992.
  * [17] Kip S Thorne and Roger D Blandford.  Modern classical physics: optics, fluids, plasmas, elasticity, relativity, and statistical physics.  Princeton University Press, 2017.
  * [18] AA Vedenov.  Theory of a weakly turbulent plasma.  Reviews of plasma physics, pages 229–276, 1967.
  * [19] AA Vedenov, EP Velikhov, and RZ Sagdeev.  Nonlinear oscillations of rarified plasma.  Nuclear Fusion, 1(2):82, 1961.
  * [20] Chenglong Zhang and Irene M Gamba.  A conservative scheme for vlasov poisson landau modeling collisional plasmas.  Journal of Computational Physics, 340:470–497, 2017.
  * [21] Chenglong Zhang and Irene M Gamba.  A conservative discontinuous galerkin solver for the space homogeneous boltzmann equation for binary interactions.  SIAM Journal on Numerical Analysis, 56(5):3040–3070, 2018.
  * [22] Xiangxiong Zhang and Chi-Wang Shu.  On maximum-principle-satisfying high order schemes for scalar conservation laws.  Journal of Computational Physics, 229(9):3091–3120, 2010.