# Anti-symmetric and Positivity Preserving Formulation of a Spectral Method
for Vlasov-Poisson Equations

Opal Issan<EMAIL_ADDRESS>Oleksandr Koshkarov  Federico D. Halpern  Boris
Kramer  Gian Luca Delzanno  Department of Mechanical and Aerospace
Engineering, University of California San Diego, La Jolla, CA, USA  T-5
Applied Mathematics and Plasma Physics Group, Los Alamos National Laboratory,
Los Alamos, NM, USA  General Atomics, P.O. Box 85608, San Diego, CA, USA

###### Abstract

We analyze the anti-symmetric properties of spectral discretization for the
one-dimensional Vlasov-Poisson equations. The discretization is based on a
spectral expansion in velocity with the symmetrically weighted Hermite basis
functions, central finite differencing in space, and an implicit Runge Kutta
integrator in time. The proposed discretization preserves the anti-symmetric
structure of the advection operator in the Vlasov equation, resulting in a
stable numerical method. We apply such discretization to two formulations: the
canonical Vlasov-Poisson equations and their continuously transformed square-
root representation. The latter preserves the positivity of the particle
distribution function. We derive analytically the conservation properties of
both formulations, including particle number, momentum, and energy, which are
verified numerically on the following benchmark problems: manufactured
solution, linear and nonlinear Landau damping, two-stream instability, and
bump-on-tail instability.

###### keywords:

Vlasov-Poisson equations , spectral methods , conservation laws , anti-
symmetry

##  1 Introduction

Many plasma physics phenomena require a kinetic description, rather than a
fluid description, to accurately represent the system’s dynamics. Important
examples include the Earth’s magnetosphere, the solar corona, galactic
dynamics, and laboratory fusion devices. The kinetic Vlasov equation describes
the evolution of the particle (e.g., ions, electrons, etc) distribution
function in six-dimensional phase space (three spatial and three velocity
coordinates) coupled to self-consistent electromagnetic fields.

Numerically solving the Vlasov equation for collisionless plasmas is
challenging since it is high dimensional, nonlinearly coupled to the Maxwell
or Poisson equations, has a rich geometric structure, and has large spatial
and temporal scale separation. The scale separation occurs in the microscopic
setting since ions are typically at least three orders of magnitude heavier
than electrons, and in the macroscopic setting since the system’s scales are
much larger than the local plasma scales. For instance, in the solar corona
the Debye length, a kinetic scale, is approximately $10\unit{}$, whereas the
solar corona’s system length scale is approximately $10^{11}\unit{}$ [39].
Moreover, the Vlasov equation has infinitely many invariants due to its
Hamiltonian structure. Retaining such conservation properties at the fully
discrete level is an important yet nontrivial task. Numerically handling the
high-dimensional phase space of the Vlasov equation is also an ongoing
challenge. Standard grid discretization, such as finite differencing, of the
six-dimensional phase space results in a prohibitively large number of degrees
of freedom. For example, discretizing each direction in phase space using
$10^{3}$ grid points results in $10^{18}$ degrees of freedom in each time
step, which is equivalent to 8 exabytes of memory in double precision. At the
time of this writing, the world’s fastest supercomputer (Frontier) has 9.2
petabytes (0.0092 exabytes) of memory [12].

Two approaches for solving the kinetic equations are prevalent: Lagrangian
(particle-based) and Eulerian (grid or spectral-based). The first approach
discretizes phase space via macroparticles and follows their characteristics.
Particle-based methods are the most commonly used methods in the community due
to their simplicity, robustness, and success during early adoption. However,
particle-based methods introduce statistical errors commonly referred to as
‘discrete particle noise’, which are due to random sampling of the particle
distribution function and its interpolation to the grid. The most
straightforward way to reduce the noise is to increase the number of
macroparticles, in which the noise decreases only as $1/\sqrt{N}$, where $N$
is the number of macroparticles [53]. Other denoising strategies include
filtering [53], variance reduction techniques [11], and phase space remapping
[55, 36]. This renders particle-based methods effective only for problems
where a low signal-to-noise ratio is acceptable.

The second approach is to solve the Vlasov equation in an Eulerian coordinate
frame using a grid or spectral-based discretization. Eulerian methods do not
have particle noise and therefore it is easier to resolve fine-scale
structures. Grid-based methods have been successful typically in reduced
dimensions, with a variety of discretization techniques such as finite
difference [8, 48], finite volume [13, 54], finite element [58, 57], and
discontinuous Galerkin [22, 9]. Alternatively, spectral methods expand the
particle distribution function in velocity using an orthogonal basis. This can
be achieved using different basis sets, including Fourier [25], Chebyshev
[49], and Legendre [32] expansions, but most of the spectral method literature
for the Vlasov equation uses a Hermite spectral expansion in velocity space
[44, 28, 6, 51, 52, 10, 24, 37, 38, 14, 41, 4, 27, 40].

Hermite-based velocity expansions originally proposed by [17] are particularly
advantageous since they can capture near-Maxwellian distribution functions
with a few basis functions, enabling 3D-3V simulations [10, 52, 44]. This is
because the Hermite basis is defined by the Hermite polynomials with a
Maxwellian weight. Therefore, an exact Maxwellian distribution can be
represented by the zeroth-order Hermite basis function. There are two types of
Hermite expansions: asymmetrically weighted (AW) and symmetrically weighted
(SW). The two expansions have different properties. Although the AW expansion
conserves particle number, momentum, and energy, it suffers from numerical
instability [6]. Conversely, the SW expansion is stable, yet conservation of
particle number, momentum, and energy is limited [33, 24, 27]. Additionally,
both expansions do not enforce the positivity of the particle distribution
function.

An important property of the Vlasov equation is the anti-symmetric (also known
as skew-adjoint, anti-self-adjoint, or skew-symmetry) structure of the phase
space advection operator. Since an anti-symmetric operator conserves square
norms, a discretization that preserves the anti-symmetric structure of the
advection operator in the Vlasov equation is guaranteed to be numerically
stable. The exploitation of this concept dates back to a classic paper by
Arakawa [2], which introduced the first numerically stable 2D finite-
difference method for non-linear incompressible fluid flows. Numerical
stability was achieved by introducing a linear combination of different
representations of the advection operator with different numerical properties.
This allows retaining some of the symmetries of the dynamics, such as the
anti-symmetry of the flow operator and two square quantities – kinetic energy
and enstrophy – to be conserved in discrete space, resulting in a numerically
stable method. Recently, Halpern et al. [21, 19, 20] formulated an anti-
symmetric representation of the plasma and fluid equations, ranging from the
two-fluid Braginskii model down to the Navier-Stokes model. This formulation
first transforms the continuous equations by introducing a novel set of
dependent variables in which the continuous force operator becomes anti-
symmetric. The new dependent variables include the square root of the density,
similar to Roe variables for the Euler equations [43], which by construction
is positivity preserving in the discrete setting. The spatial discretization
of the transformed equations is based on central finite differencing, which
results in a discrete anti-symmetric force operator, leading to a stable,
robust, and conservative method. The numerical stability and conservation
properties are guaranteed in the anti-symmetric formulation by drawing
connections between the continuous and discrete forms, instead of through
direct discretization of the Hamiltonian function.

This paper has four main contributions. First, we show that the SW expansion
of the particle distribution function preserves the anti-symmetric structure
of the advection operator, or equivalently the canonical Poisson bracket, in
the Vlasov equation. Second, we show that central finite differencing in space
can preserve the anti-symmetric structure and result in a discrete anti-
symmetric advection operator acting on the expansion coefficients. Preserving
the anti-symmetry of the advection operator results in a structure-preserving
and stable method. Additionally, central finite differencing scales well in
high-dimensional parallel code due to its sparsity and allows for more
flexibility with non-periodic boundary conditions, unlike the previously
proposed Fourier spectral expansion. Third, we propose a new SW square-root
formulation, which solves for the square root of the particle distribution
function, with the same SW velocity expansion and central finite differencing
spatial discretization. By construction, the SW square-root formulation is
positivity preserving and retains the numerical stability and anti-symmetry of
the SW formulation. Fourth, we derive the conservation properties of the SW
and SW square-root formulations for the one-dimensional Vlasov-Poisson
equations. We verify numerically such conservation properties on the following
benchmark problems: manufactured solution, linear and nonlinear Landau
damping, two-stream instability, and bump-on-tail instability.

This paper is organized as follows. Section 2 introduces the one-dimensional
Vlasov-Poisson equations and the velocity discretization for the SW and SW
square-root formulations. Section 3 presents the proposed anti-symmetric
discretization of the Vlasov equation, leading to a stable method. Section 4
derives the conservation properties of the SW and SW square-root formulations.
In section 5, we test the two formulations on a number of classical benchmark
problems. Lastly, section 6 offers conclusions and an outlook for future work.

##  2 Vlasov-Poisson Equations: Hermite Spectral Discretization in Velocity

We present the Vlasov-Poisson equations in section 2.1. Section 2.2 presents
the SW Hermite basis functions. We introduce the SW expansion in section 2.3
and the SW square-root expansion in section 2.4.

###  2.1 Vlasov-Poisson Equations

We study the one-dimensional Vlasov-Poisson equations which model the coupling
between collisionless plasmas and a self-consistent electric field. The plasma
is described by the non-negative and bounded particle distribution function
$f^{s}(x,v,t)$ in phase-space, which represents the number of particles
located at the spatial position between $x\in\Omega_{x}$ and
$x+\mathrm{d}x\in\Omega_{x}$ with velocity between $v\in\Omega_{v}$ and
$v+\mathrm{d}v\in\Omega_{v}$ at a given time
$t\in\Omega_{t}\subseteq\mathbb{R}_{\geq 0}$. The plasma species are denoted
by $s$, where $s=`e$’ denotes electrons and $s=`i$’ denotes singly-charged
ions. We assume the spatial coordinate $x$ is periodic in the spatial domain
$\Omega_{x}=[0,\ell]$, where $\ell$ is the length of the spatial domain,
$\Omega_{t}=[0,t_{f}]$, where $t_{f}\in\mathbb{R}_{\geq 0}$ is the final time,
and allow the velocity coordinate to be unbounded, i.e.
$\Omega_{v}=\mathbb{R}$. Let $L^{2}(\Omega_{x}\times\Omega_{v})$ denote the
real Hilbert space of square-integrable functions on
$\Omega_{x}\times\Omega_{v}$ with inner product $\langle
f,g\rangle=\int_{\Omega_{x}}\int_{\Omega_{v}}f(x,v,t)g(x,v,t)\mathrm{d}x\mathrm{d}v$.
We next introduce the spaces of functions that are periodic in space $x$ and
tend to zero as $|v|\to\infty$ as

$\displaystyle\mathcal{V}^{d}\coloneqq$
$\displaystyle\\{[f_{1},\ldots,f_{d}]^{\top}(\cdot,\cdot,t)\ \text{s.t.}\
f_{i}(\cdot,\cdot,t)\in L^{2}(\Omega_{x}\times\Omega_{v})\ \text{and}\
f_{i}\geq 0,$ $\displaystyle\hskip
5.0ptf_{i}(0,\cdot,t)=f_{i}(\ell,\cdot,t),\hskip 5.0ptf_{i}(\cdot,v,t)\to
0\hskip 5.0pt\mathrm{as}\hskip 5.0pt|v|\to\infty,\ i=1,\ldots,d,\hskip
5.0pt\forall t\in\Omega_{t}\\},$ $\displaystyle\mathcal{W}^{d}\coloneqq$
$\displaystyle\\{[g_{1},\ldots,g_{d}]^{\top}(\cdot,t)\ \text{s.t.}\
g_{i}(\cdot,t)\in L^{2}(\Omega_{x})\ \text{and}\ g_{i}(0,t)=g_{i}(\ell,t),\
i=1,\ldots,d,\hskip 5.0pt\forall t\in\Omega_{t}\\},$

where we abbreviate $\mathcal{V}=\mathcal{V}^{1}$ and
$\mathcal{W}=\mathcal{W}^{1}$. The one-dimensional (normalized) Vlasov-Poisson
equations are given by

$\displaystyle\left(\frac{\partial}{\partial t}+v\frac{\partial}{\partial
x}+\frac{q^{s}}{m^{s}}E(x,t)\frac{\partial}{\partial v}\right)f^{s}(x,v,t)$
$\displaystyle=0,\qquad$ $\displaystyle\text{in
}\Omega_{x}\times\Omega_{v}\times\Omega_{t},\hskip 5.0pt\forall s,$ (1)
$\displaystyle\frac{\partial E(x,t)}{\partial x}$
$\displaystyle=\rho(x,t)\qquad$ $\displaystyle\text{in
}\Omega_{x}\times\Omega_{t},$ (2)

where
$\rho(x,t)\coloneqq\sum_{s}q^{s}\int_{\Omega_{v}}f^{s}(x,v,t)\mathrm{d}v$ is
the charge density, $E\in\mathcal{W}$ is the normalized electric field, and
$m^{s}$ and $q^{s}$ are the normalized mass and charge of species $s$. The
electric current $J(x,t)$ is obtained from the particle distribution functions
as $J(x,t)\coloneqq\sum_{s}q^{s}\int_{\Omega_{v}}vf^{s}(x,v,t)\mathrm{d}v$.
The uniqueness of the solution is ensured by imposing that
$\int_{\Omega_{x}}E(x,t)\mathrm{d}x=0,\hskip 5.0pt\forall t\in\Omega_{t}$. To
complete the problem formulation, we set a suitable initial condition to the
particle distribution function $f^{s}(x,v,t=0)=f^{s}_{0}(x,v)$ for each of the
species $s$, and compute the initial electric field $E(x,t=0)$ by solving the
Poisson equation (2) at $t=0$.

##### Normalization

In the above, we normalized the physical quantities as follows. Time $t$ is
normalized to the electron plasma frequency
$\omega_{pe}=\sqrt{(q^{e})^{2}n_{0}^{e}/(\varepsilon_{0}m^{e})}$, where
$q^{e}$ is the elementary charge, $m^{e}$ is the electron mass,
$\varepsilon_{0}$ is the permittivity of vacuum, and $n_{0}^{e}$ is the
reference electron density. The velocity coordinate $v$ is normalized to the
electron thermal velocity $v_{t}=\sqrt{kT_{e}/m^{e}}$, where $k$ is the
Boltzmann constant and $T_{e}$ is a reference electron temperature; the
spatial coordinate $x$ is normalized to the electron Debye length
$\lambda_{D}=\sqrt{\varepsilon_{0}kT_{e}/((q^{e})^{2}n_{0}^{e})}$; the
electric field $E(x,t)$ is normalized to $m^{e}v_{t}\omega_{pe}/q^{e}$; the
particle distribution function $f^{s}(x,v,t)$ is normalized to
$n^{e}_{0}/v_{t}$; the mass and charge of species $s$, i.e. $m^{s}$ and
$q^{s}$, are normalized by electron mass $m^{e}$ and charge $q^{e}$.

###  2.2 Symmetrically-Weighted Hermite Basis Functions and Their Properties

We introduce the symmetrically-weighted (SW) Hermite basis functions
$\\{\psi_{n}(\xi^{s})\\}_{n\in\mathbb{N}_{0}}$, which are a function of the
scaled velocity coordinate $\xi^{s}$ and are defined as

$\psi_{n}(\xi^{s})\coloneqq(\sqrt{\pi}2^{n}n!)^{-\frac{1}{2}}\mathcal{H}_{n}(\xi^{s})\exp{\left(-\frac{(\xi^{s})^{2}}{2}\right)},\qquad\xi^{s}(v)\coloneqq\frac{v-u^{s}}{\alpha^{s}},$
(3)

where the velocity coordinate $v$ is shifted by $u^{s}\in\mathbb{R}$ and
scaled by $\alpha^{s}\in\mathbb{R}_{>0}$ [50]. Here, we treat parameters
$u^{s}$ and $\alpha^{s}$ as constants. Moreover, $\mathcal{H}_{n}(\xi^{s})$ is
the “physicist” Hermite polynomial [1, §22] of degree $n\in\mathbb{N}_{0}$,
i.e.

$\mathcal{H}_{n}(\xi^{s})\coloneqq(-1)^{n}\exp{\left((\xi^{s})^{2}\right)}\frac{\mathrm{d}^{n}}{\mathrm{d}(\xi^{s})^{n}}\exp{\left((-\xi^{s})^{2}\right)},$

with the following recursive properties

$\displaystyle\frac{\mathrm{d}\mathcal{H}_{n}(\xi^{s})}{\mathrm{d}\xi^{s}}$
$\displaystyle=2n\mathcal{H}_{n-1}(\xi^{s})\qquad$
$\displaystyle\mathrm{for}\qquad n\geq 1,$ (4)
$\displaystyle\xi^{s}\mathcal{H}_{n}(\xi^{s})$
$\displaystyle=\frac{1}{2}\mathcal{H}_{n+1}(\xi^{s})+n\mathcal{H}_{n-1}(\xi^{s})\qquad$
$\displaystyle\mathrm{for}\qquad n\geq 2.$ (5)

The SW-Hermite basis functions satisfy the orthogonality relation

$\displaystyle\int_{\mathbb{R}}\psi_{n}(\xi^{s})\psi_{m}(\xi^{s})\mathrm{d}\xi^{s}=\delta_{n,m},$
(6)

where $\delta_{n,m}$ is the Kronecker delta function. Additional properties of
the SW-Hermite basis functions that we leverage later include

$\displaystyle v\psi_{n}(\xi^{s})$
$\displaystyle=\alpha^{s}\sqrt{\frac{n+1}{{2}}}\psi_{n+1}(\xi^{s})+\alpha^{s}\sqrt{\frac{n}{2}}\psi_{n-1}(\xi^{s})+u^{s}\psi_{n}(\xi^{s}),$
(7) $\displaystyle v^{2}\psi_{n}(\xi^{s})$
$\displaystyle=(\alpha^{s})^{2}\frac{\sqrt{(n+1)(n+2)}}{2}\psi_{n+2}(\xi^{s})+2\alpha^{s}u^{s}\sqrt{\frac{n+1}{2}}\psi_{n+1}(\xi^{s})$
(8)
$\displaystyle+\left((\alpha^{s})^{2}\left(\frac{2n+1}{2}\right)+(u^{s})^{2}\right)\psi_{n}(\xi^{s})$
$\displaystyle+2\alpha^{s}u^{s}\sqrt{\frac{n}{2}}\psi_{n-1}(\xi^{s})+(\alpha^{s})^{2}\frac{\sqrt{n(n-1)}}{2}\psi_{n-2}(\xi^{s}),$
$\displaystyle\frac{\mathrm{d}\psi_{n}(\xi^{s})}{\mathrm{d}v}$
$\displaystyle=\frac{1}{\alpha^{s}}\sqrt{\frac{n}{2}}\psi_{n-1}(\xi^{s})-\frac{1}{\alpha^{s}}\sqrt{\frac{n+1}{2}}\psi_{n+1}(\xi^{s}).$
(9)

Identity (7) follows from the recursive relation in Eq. (5), identity (8)
follows from applying the recursive formula in Eq. (5) twice, and lastly,
identity (9) follows from the recursive relation in Eq. (4).

###  2.3 Symmetrically-Weighted Expansion

We approximate the particle distribution function $f^{s}(x,v,t)$ via a
truncated spectral decomposition in the scaled velocity coordinate $\xi^{s}$,
taking the following ansatz

$f^{s}(x,v,t)\approx
f^{s,N_{v}}(x,v,t)=\sum_{n=0}^{N_{v}-1}C_{n}^{s}(x,t)\psi_{n}(\xi^{s}),\qquad\xi^{s}(v)\coloneqq\frac{v-u^{s}}{\alpha^{s}},$
(10)

where $\psi_{n}(\xi^{s})$ is the SW-Hermite basis function defined in Eq. (3)
and $C_{n}^{s}\in\mathcal{W}$ is the expansion coefficient. The Hermite basis
is a natural choice for Maxwellian-type distributions because the Hermite
basis is defined by the Hermite polynomials with a Maxwellian weight, see Eq.
(3). By substituting the spectral approximation in Eq. (10) in the Vlasov
equation (1), multiplying against $\psi_{m}(\xi^{s})$, and integrating with
respect to $\xi^{s}$, we get

$\int_{\mathbb{R}}\psi_{m}(\xi^{s})\left[\sum_{n=0}^{N_{v}-1}\frac{\partial
C_{n}^{s}(x,t)}{\partial
t}\psi_{n}(\xi^{s})+v\sum_{n=0}^{N_{v}-1}\frac{\partial
C_{n}^{s}(x,t)}{\partial
x}\psi_{n}(\xi^{s})+\frac{q^{s}}{m^{s}}E(x,t)\sum_{n=0}^{N_{v}-1}C_{n}^{s}(x,t)\frac{\mathrm{d}\psi_{n}(\xi^{s})}{\mathrm{d}v}\right]\mathrm{d}\xi^{s}=0.$

By exploiting the orthogonality and recurrence relation of the Hermite basis
functions in Eqns. (6), (7) and (9), we derive a system of partial
differential equations (PDEs) for the expansion coefficients $C_{n}^{s}(x,t)$:

$\displaystyle\frac{\partial C_{n}^{s}(x,t)}{\partial t}$
$\displaystyle+\frac{\partial}{\partial
x}\underbrace{\left(\alpha^{s}\sqrt{\frac{n}{{2}}}C_{n-1}^{s}(x,t)+\alpha^{s}\sqrt{\frac{n+1}{2}}C_{n+1}^{s}(x,t)+u^{s}C_{n}^{s}(x,t)\right)}_{\coloneqq
Q^{s}_{n}(x,t)}$ (11)
$\displaystyle+\frac{q^{s}}{m^{s}\alpha^{s}}E(x,t)\left(-\sqrt{\frac{n}{2}}C_{n-1}^{s}(x,t)+\sqrt{\frac{n+1}{2}}C_{n+1}^{s}(x,t)\right)=0,$

for $n=0,1,\ldots,N_{v}-1$ with the convention that $C_{n<0}^{s}(x,t)=0$ and
imposing for the closure of the system to $C_{N_{v}}^{s}(x,t)=0$. Similarly,
after substituting the spectral approximation in Eq. (10) in the Poisson
equation (2), we get

$\frac{\partial E(x,t)}{\partial
x}=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}C_{n}^{s}(x,t)\qquad\mathrm{with}\qquad\mathcal{I}_{n}=\int_{\mathbb{R}}\psi_{n}(\xi^{s})\mathrm{d}\xi^{s}.$
(12)

The integral $\mathcal{I}_{n}$ has the following recursive property:

$\mathcal{I}_{0}=\sqrt{2}(\pi)^{\frac{1}{4}},\qquad\mathcal{I}_{1}=0,\qquad\mathcal{I}_{n}=\sqrt{\frac{n-1}{n}}\mathcal{I}_{n-2}\qquad\mathrm{for}\qquad
n\geq 2.$ (13)

###  2.4 Symmetrically-Weighted Square-Root Expansion

To preserve the non-negative property of the particle distribution function
$f^{s}(x,v,t)\geq 0$, we approximate the square root of the particle
distribution function, i.e. $\sqrt{f^{s}(x,v,t)}$, via an SW-Hermite spectral
expansion in velocity $v$, where

$\sqrt{f^{s}(x,v,t)}\approx\sqrt{f^{s,N_{v}}(x,v,t)}=\sum_{n=0}^{N_{v}-1}C_{n}^{s}(x,t)\psi_{n}(\xi^{s}).$
(14)

######  Proposition 1.

The Vlasov equation (1) in square-root form is

$\left(\frac{\partial}{\partial t}+v\frac{\partial}{\partial
x}+\frac{q^{s}}{m^{s}}E(x,t)\frac{\partial}{\partial
v}\right)\sqrt{f^{s}(x,v,t)}=0.$ (15)

A solution to Eq. (15) is a solution to the Vlasov equation (1) and vice
versa.

Proof: By the product rule and non-negative property of $f^{s}(x,v,t)$, we get

$\frac{\partial f^{s}(x,v,t)}{\partial
t}=\frac{\partial\left(\sqrt{f^{s}(x,v,t)}\right)^{2}}{\partial
t}=2\sqrt{f^{s}(x,v,t)}\frac{\partial\sqrt{f^{s}(x,v,t)}}{\partial t}.$ (16)

Equation (16) equivalently holds for the partial derivatives of $f^{s}(x,v,t)$
with respect to $x$ and $v$. We recast all partial derivatives in the Vlasov
equation (1) using Eq. (16), which gives

$2\sqrt{f^{s}(x,v,t)}\left(\frac{\partial}{\partial
t}+v\frac{\partial}{\partial
x}+\frac{q^{s}}{m^{s}}E(x,t)\frac{\partial}{\partial
v}\right)\sqrt{f^{s}(x,v,t)}=0.$

Thus, either $\sqrt{f^{s}(x,v,t)}=0$ or Eq. (15) holds. Since Eq. (15)
encompasses the trivial solution $\sqrt{f^{s}(x,v,t)}=0$, the result in
Proposition 1 is proven. $\Box$

By Proposition 1 we obtain the same system of PDEs for the expansion
coefficients $C_{n}^{s}(x,t)$ in Eq. (14) as the SW expansion coefficients in
Eq. (11). Moreover, the square-root spectral approximation (14) introduces
quadratic non-linearity in the Poisson equation (2), such that

$\frac{\partial E(x,t)}{\partial
x}=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}C_{n}^{s}(x,t)^{2}.$ (17)

The only differences between the SW and SW square-root formulations are the
velocity scaling coefficient $\alpha^{s}\in\mathbb{R}_{>0}$ in Eq. (3), the
right-hand side of the Poisson equation (12) and (17), and the expansion
coefficient’s initial condition.

##  3 Anti-symmetric Representation

The Vlasov-Poisson system has a rich geometric Hamiltonian structure.
Preserving such geometric structure on the discrete level results in accurate
long-term numerical simulations of plasma physics processes [29, 30, 56, 23].
Geometric structure includes conservation laws, such as the conservation of
mass, momentum, and energy, symmetries and constraints, such as the particle
distribution function being non-negative and bounded, and identities, for
example, those from vector calculus, $\text{curl grad}=0$ and $\text{div
curl}=0$. Recently, Halpern et al. [21, 19, 20] developed an anti-symmetric
structure-preserving formulation of the Navier-Stokes and magnetohydrodynamics
equations leading to simple, robust, and stable simulations. Here, we draw
insight from that work and focus on preserving the anti-symmetric property of
the advection operator in the Vlasov equation (1), thereby preventing
nonphysical dissipation. For ease of exposition, we present the anti-symmetric
formulation in the 1D-1V setting, however, the anti-symmetric formulation can
be extended to multidimensional problems with minor modifications.

Section 3.1 presents preliminary definitions and notation. In section 3.2, we
show that the semi-discrete Vlasov equation preserves the anti-symmetry of the
advection operator. We explain the proposed method’s stability properties in
section 3.3 and its time-reversibility in section 3.4. Lastly, section 3.5
relates the anti-symmetric properties to the Hamiltonian canonical and
noncanonical Poisson brackets.

###  3.1 Definitions and Notation

We start our discussion by recalling the definition of the adjoint of an
operator and an anti-symmetric operator along with simple examples.

######  Definition 1 (Adjoint of an operator [42, §8.5.8]).

Let $H$ be a complex Hilbert space, whose inner product is denoted by
$\langle\cdot,\cdot\rangle$. We consider a densely defined operator
$(\mathcal{D}(A),A)$ from a Hilbert space $H_{1}$ to a Hilbert space $H_{2}$.
The adjoint of the operator $(\mathcal{D}(A),A)$ is defined as the operator
$(\mathcal{D}(A^{*}),A^{*})$ from $H_{2}$ to $H_{1}$, such that $\langle
Ay,z\rangle=\langle y,A^{*}z\rangle,\forall y\in H_{1},\forall z\in H_{2}$.
The asterisk denotes the conjugate transpose of the operator; for real
operators, the conjugate transpose is the transpose, i.e. $A^{*}=A^{\top}$.

######  Example 1.

We give an example of the computation of an adjoint of a discrete operator.
Let $\mathcal{D}(A)=\mathbb{R}^{n}$, $A\in\mathbb{R}^{m\times n}$,
$y\in\mathbb{R}^{n}$, $z\in\mathbb{R}^{m}$, then the adjoint of the discrete
operator $(\mathcal{D}(A),A)$ is its transpose
$(\mathcal{D}(A^{\top}),A^{\top})$, where
$\mathcal{D}(A^{\top})=\mathbb{R}^{m}$ and $A^{\top}\in\mathbb{R}^{n\times
m}$.

######  Example 2.

We also give an example of the computation of an adjoint of a continuous
operator. Let $H=C^{1}(0,1)$, the space of continuous functions that have
continuous first derivatives on $[0,1]$, $\mathcal{D}(A)=\\{y\in H\hskip
5.0pt|\hskip 5.0pty(0)=0,\hskip 5.0pty(1)=0\\}$, $y,z\in H$, and
$A=\frac{\mathrm{d}}{\mathrm{d}x}$. We derive the adjoint of $A$ using
integration by parts

$\langle
Ay,z\rangle=\int_{0}^{1}zAy\mathrm{d}x=yz\Big{|}_{0}^{1}-\int_{0}^{1}yAz\mathrm{d}x=\langle
y,-Az\rangle$

such that $(\mathcal{D}(A),A^{\top})=(\mathcal{D}(A),-A)$.

Example 2 leads to the definition of an anti-symmetric operator.

######  Definition 2 (Anti-symmetric operator [20, §2]).

The densely defined operator $(\mathcal{D}(A),A)$ from $H$ to $H$ is anti-
symmetric111also referred to as “skew-adjoint”, “anti-self-adjoint”, or “skew-
symmetric” in literature. if $(\mathcal{D}(A),A)=(\mathcal{D}(A^{*}),-A^{*})$,
where $(\mathcal{D}(A^{*}),A^{*})$ is the adjoint of the operator
$(\mathcal{D}(A),A)$. An anti-symmetric operator $(\mathcal{D}(A),A)$ has the
following properties:

$\displaystyle\text{(anti-commutative) }\qquad\hskip 5.0pt\langle Ay,z\rangle$
$\displaystyle=-\langle y,Az\rangle,\qquad$ $\displaystyle\forall y,z\in H,$
(18) $\displaystyle\text{(preserves square-norms) }\langle Ay,y\rangle$
$\displaystyle=0,\qquad$ $\displaystyle\forall y\in H.$ (19)

###  3.2 Continuum and Discrete Anti-symmetric Representation of the Vlasov
Equation

The advection operator in the Vlasov equation (1) is anti-symmetric (in a
periodic domain), which we describe in detail in section 3.2.1. We prove in
section 3.2.2 that the anti-symmetric structure of the advection operator is
preserved after discretization in velocity using SW-Hermite expansion, and in
section 3.2.3, we show that anti-symmetry is preserved after discretization in
space using centered finite differencing.

####  3.2.1 Anti-symmetric Representation in the Continuum

We represent the Vlasov equation (1) and its square-root equivalent (15) in
anti-symmetric form, where the advection operator is anti-symmetric.

######  Theorem 1.

The Vlasov equation (1) can be written in an anti-symmetric form as

$\frac{\partial}{\partial
t}f^{s}(x,v,t)=A^{s}(x,v,t)f^{s}(x,v,t)\qquad\text{or}\qquad\frac{\partial}{\partial
t}\sqrt{f^{s}(x,v,t)}=A^{s}(x,v,t)\sqrt{f^{s}(x,v,t)},$ (20)

where

$A^{s}(x,v,t)=-\dot{\vec{x}}\cdot\nabla_{\vec{x}}$ (21)

and $\vec{x}=[x,v]^{\top}$,
$\dot{\vec{x}}=\frac{\mathrm{d}}{\mathrm{d}t}\vec{x}=[v,\frac{q^{s}}{m^{s}}E(x,t)]^{\top}$,
such that the advection operator $(\mathcal{V},A^{s}(x,v,t))$ is anti-
symmetric in a periodic domain.

Proof: The Vlasov equation (1) can be written as

$\frac{\partial}{\partial
t}f^{s}(x,v,t)=-\dot{\vec{x}}\cdot\nabla_{\vec{x}}f^{s}(x,v,t)=A^{s}(x,v,t)f^{s}(x,v,t),$
(22)

which proves the first equality in Eq. (20). By the non-negativity of
$f^{s}(x,v,t)$, the right hand side of Eq. (22) can be transformed to

$-\dot{\vec{x}}\cdot\nabla_{\vec{x}}f^{s}(x,v,t)=-2\sqrt{f^{s}(x,v,t)}\left(\dot{\vec{x}}\cdot\nabla_{\vec{x}}\sqrt{f^{s}(x,v,t)}\right)=2\sqrt{f^{s}(x,v,t)}A^{s}(x,v,t)\sqrt{f^{s}(x,v,t)}.$

The partial derivative of the distribution function in time, i.e.
$\frac{\partial f^{s}(x,v,t)}{\partial t}$, can be expressed in terms of its
square root, see Eq. (16), such that Eq. (22) becomes the second equality in
Eq. (20). The advection operator $\left(\mathcal{V},A^{s}(x,v,t)\right)$ in
Eq. (21) is anti-symmetric since $\langle A^{s}y,z\rangle=\langle
y,-A^{s}z\rangle$ holds for all $y,z\in\mathcal{V}$. The anti-symmetry is
obtained by the product rule, divergence theorem, and periodic boundary
conditions. Therefore, we get
$(\mathcal{V},A^{s}(x,v,t))=(\mathcal{V},-A^{s}(x,v,t)^{\top})$. $\Box$

####  3.2.2 Anti-symmetric Representation after Discretization in Velocity

We show that after discretization in velocity using the SW-Hermite spectral
expansion, we preserve the anti-symmetric structure of the advection operator
in the Vlasov equation (20) for the SW and SW square-root formulations. The
anti-symmetry of $A^{s}(x,v,t)$ is preserved since the SW-Hermite expansion
results in a self-adjoint Galerkin projection operator [16].

######  Theorem 2.

The SW-Hermite velocity spectral discretization in Eqns. (10) and (14)
preserves the anti-symmetric structure of the advection operator
$\left(\mathcal{V},A^{s}(x,v,t)\right)$ in Eq. (21).

Proof: The system of PDEs for the SW and SW square-root formulation expansion
coefficients $C_{n}^{s}(x,t)$ in Eq. (11) can be written in vector form as

$\frac{\partial}{\partial t}\Psi^{s}(x,t)=A_{v}^{s}(x,t)\Psi^{s}(x,t),$ (23)

where

$\Psi^{s}(x,t)\coloneqq\begin{bmatrix}C_{0}^{s}(x,t)\\\ C_{1}^{s}(x,t)\\\
\vdots\\\ C_{N_{v}-1}^{s}(x,t)\\\ \end{bmatrix}\in\mathcal{W}^{N_{v}},$
$A_{v}^{s}(x,t)=\begin{bmatrix}-u^{s}\frac{\partial}{\partial
x}&-\sqrt{\frac{1}{2}}G_{+}^{s}(x,t)&0&\ldots&0\\\
\sqrt{\frac{1}{2}}G_{-}^{s}(x,t)&-u^{s}\frac{\partial}{\partial
x}&-\sqrt{\frac{2}{2}}G_{+}^{s}(x,t)&\ldots&0\\\ &&&&\\\
&\ddots&\ddots&\ddots&\\\ &&&&\\\
0&0&\sqrt{\frac{N_{v}-2}{2}}G_{-}^{s}(x,t)&-u^{s}\frac{\partial}{\partial
x}&-\sqrt{\frac{N_{v}-1}{2}}G_{+}^{s}(x,t)\\\
0&0&0&\sqrt{\frac{N_{v}-1}{2}}G_{-}^{s}(x,t)&-u^{s}\frac{\partial}{\partial
x}\end{bmatrix}\in\mathcal{W}^{N_{v}\times N_{v}},$

such that

$G_{+}^{s}(x,t)=\frac{q^{s}}{m^{s}\alpha^{s}}E(x,t)+\alpha^{s}\frac{\partial}{\partial
x}\qquad\text{and}\qquad
G_{-}^{s}(x,t)=\frac{q^{s}}{m^{s}\alpha^{s}}E(x,t)-\alpha^{s}\frac{\partial}{\partial
x}.$

Since the operator $(\mathcal{W},\frac{\partial}{\partial x})$ is anti-
symmetric in a periodic domain, the operator
$\left(\mathcal{W}^{N_{v}},A_{v}^{s}(x,t)\right)$ is also anti-symmetric in a
periodic domain, i.e.
$(\mathcal{W}^{N_{v}},A_{v}^{s}(x,t))=(\mathcal{W}^{N_{v}},-A_{v}^{s}(x,t)^{\top})$.
$\Box$

####  3.2.3 Anti-symmetric Representation after Discretization in Space

After using the SW-Hermite spectral expansion in velocity, we discretize the
spatial domain using centered finite differencing. We show that the centered
finite differencing scheme preserves the anti-symmetric structure of the
advection operator in the Vlasov equation (20) for the SW and SW square-root
formulations. In the semi-discrete level, the dynamics are described as
infinitesimal rotations since the advection operator becomes an anti-symmetric
matrix whose nonzero eigenvalues are purely imaginary in the form of complex
conjugate pairs [20]. It is particularly advantageous to use centered finite
differencing since aside from retaining the anti-symmetric property of the
advection operator it also preserves important calculus identities such as
$\text{curl grad}=0$ and $\text{div curl}=0$ in multi-dimensional settings
[20].

######  Theorem 3.

Centered finite difference representation of spatial derivatives in Eq. (23)
preserves the anti-symmetric structure of the advection operator
$(\mathcal{W}^{N_{v}},A^{s}_{v}(x,t))$ in Eq. (23), and thus preserves the
anti-symmetric structure of the advection operator
$(\mathcal{V},A^{s}(x,v,t))$ in Eq. (20).

Proof: We discretize the spatial direction uniformly with $\Delta x$ mesh
spacing and denote the discretized expansion coefficients as
$\mathbf{C}^{s}_{n}(t)=[C^{s}_{n}(x^{(1)},t),C^{s}_{n}(x^{(2)},t),\ldots,C^{s}_{n}(x^{(N_{x})},t)]^{\top}\in\mathbb{R}^{N_{x}}$,
where $N_{x}$ is the number of mesh points in space and $x^{(j)}$ denotes the
spatial grid location with index $j=1,2,\ldots,N_{x}$. Discretizing Eq. (23)
using centered finite differencing results in the following semi-discrete
system of ordinary differential equations (ODEs):

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Psi}^{s}(t)=\mathbf{A}^{s}(t)\mathbf{\Psi}^{s}(t)$
(24)

where

$\mathbf{\Psi}^{s}(t)=\begin{bmatrix}\mathbf{C}^{s}_{0}(t)\\\
\mathbf{C}^{s}_{1}(t)\\\ \vdots\\\
\mathbf{C}^{s}_{N_{v}-1}(t)\end{bmatrix}\in\mathbb{R}^{N_{v}N_{x}}$ (25)
$\mathbf{A}^{s}(t)=\begin{bmatrix}-u^{s}\mathbf{D}^{p}&-\sqrt{\frac{1}{2}}\mathbf{G}_{+}^{s,p}(t)&0&\ldots&0\\\
\sqrt{\frac{1}{2}}\mathbf{G}_{-}^{s,p}(t)&-u^{s}\mathbf{D}^{p}&-\sqrt{\frac{2}{2}}\mathbf{G}_{+}^{s,p}(t)&\ldots&0\\\
&&&&\\\ &\ddots&\ddots&\ddots&\\\ &&&&\\\
0&0&\sqrt{\frac{N_{v}-2}{2}}\mathbf{G}_{-}^{s,p}(t)&-u^{s}\mathbf{D}^{p}&-\sqrt{\frac{N_{v}-1}{2}}\mathbf{G}_{+}^{s,p}(t)\\\
0&0&0&\sqrt{\frac{N_{v}-1}{2}}\mathbf{G}_{-}^{s,p}(t)&-u^{s}\mathbf{D}^{p}\end{bmatrix}\in\mathbb{R}^{N_{v}N_{x}\times
N_{v}N_{x}},$

such that
$\mathbf{G}_{+}^{s,p}(t)=\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)+\alpha^{s}\mathbf{D}^{p}$,
$\mathbf{G}_{-}^{s,p}(t)=\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)-\alpha^{s}\mathbf{D}^{p}$,
and the discrete form of the electric field is denoted by
$\mathbf{E}(t)=\left[E(x^{(1)},t),E(x^{(2)},t),\ldots,E(x^{(N_{x})},t)\right]^{\top}\in\mathbb{R}^{N_{x}}$.
The semi-discrete Vlasov equation (24) in the SW formulation is coupled to the
following semi-discrete Poisson equation

$\mathbf{D}^{p}\mathbf{E}(t)=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{C}_{n}^{s}(t),$
(26)

and similarly, the equivalent semi-discrete Vlasov equation (24) in the SW
square-root formulation is coupled to

$\mathbf{D}^{p}\mathbf{E}(t)=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)^{2},$
(27)

The centered finite difference operator of order $p$ is denoted by

$\mathbf{D}^{p}=\begin{bmatrix}d_{0}^{p}&\ldots&d_{p/2}^{p}&&-d_{p/2}^{p}\\\
\ldots&d_{0}^{p}&\ldots&d_{p/2}^{p}&\\\ &&\ddots&&\\\
&-d_{p/2}^{p}&\ldots&d_{0}^{p}&\ldots\\\
d_{p/2}^{p}&&-d_{p/2}^{p}&\ldots&d_{0}^{p}\end{bmatrix}\in\mathbb{R}^{N_{x}\times
N_{x}}.$ (28)

The stencil coefficients $-d_{p/2}^{p},\ldots,d_{0}^{p},\ldots,d_{p/2}^{p}$
are derived from the Taylor series approximation of order $p$. For example,
the stencil coefficients for $p=2,4,6$ are defined as

$\displaystyle d^{2}$ $\displaystyle\coloneqq\frac{1}{\Delta
x}\left[-\frac{1}{2}\hskip 5.0pt0\hskip 5.0pt\frac{1}{2}\right]\qquad\qquad$
(2nd order) $\displaystyle d^{4}$ $\displaystyle\coloneqq\frac{1}{\Delta
x}\left[\frac{1}{12}\hskip 5.0pt-\frac{2}{3}\hskip 5.0pt0\hskip
5.0pt\frac{2}{3}\hskip 5.0pt-\frac{1}{12}\right]\qquad\qquad$ (4th order)
$\displaystyle d^{6}$ $\displaystyle\coloneqq\frac{1}{\Delta
x}\left[-\frac{1}{60}\hskip 5.0pt\frac{3}{20}\hskip 5.0pt-\frac{3}{4}\hskip
5.0pt0\hskip 5.0pt\frac{3}{4}\hskip 5.0pt\hskip 5.0pt-\frac{3}{20}\hskip
5.0pt\frac{1}{60}\right]\qquad\qquad$ (6th order)

The finite difference operator $\mathbf{D}^{p}$, defined in Eq. (28), is a
singular matrix, thus it is not invertible. Upon reflection around
$d_{0}^{p}=0$, we can see that the centered finite difference operator
$\mathbf{D}^{p}$ is anti-symmetric since
$(\mathbb{R}^{N_{x}},\mathbf{D}^{p})=(\mathbb{R}^{N_{x}},-(\mathbf{D}^{p})^{\top})$.
Thus, the discretized advection operator is anti-symmetric as
$(\mathbb{R}^{N_{v}N_{x}},\mathbf{A}^{s}(t))=(\mathbb{R}^{N_{v}N_{x}},-\mathbf{A}^{s}(t)^{\top})$.
$\Box$

######  Remark 1.

A non-uniform spacing in the spatial direction $x$ breaks anti-symmetry, since
then
$(\mathbb{R}^{N_{x}},\mathbf{D}^{p})\neq-(\mathbb{R}^{N_{x}},(\mathbf{D}^{p})^{\top})$.

######  Remark 2.

Any spectral discretization in space, such as Fourier spectral expansion, also
preserves the anti-symmetric structure of the advection operator
$(\mathcal{W}^{N_{v}},A^{s}_{v}(x,t))$ in Eq. (23), and thus preserves the
anti-symmetric structure of the advection operator
$(\mathcal{V},A^{s}(x,v,t))$ in Eq. (20). A derives the semi-discrete Vlasov-
Poisson equations using a spatial Fourier spectral expansion and discusses in
detail its anti-symmetric property.

###  3.3 Numerical Stability

The SW and SW square-root formulations are algebraically stable as
$\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)\right)\leq
0$, see [16]. In particular,
$\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)$, is a quadratic invariant of
the SW and SW square root formulations. This is due to the anti-symmetric
structure of the semi-discrete Vlasov equation (24) and the following
proposition.

######  Proposition 2.

[18, §IV.1] Consider an ODE of the form

$\frac{\mathrm{d}}{\mathrm{d}t}y=A(y)y,$

where $y(t)\in\mathbb{R}^{n}$ and $(\mathbb{R}^{n},A(y))$ is an anti-symmetric
operator, then $y^{\top}y$ is an invariant of the ODE.

Proof: We recall the proof here to keep it self-consistent:
$\frac{\mathrm{d}}{\mathrm{d}t}(y^{\top}y)=(\frac{\mathrm{d}}{\mathrm{d}t}y)^{\top}y+y^{\top}\frac{\mathrm{d}}{\mathrm{d}t}y=y^{\top}A(y)^{\top}y+y^{\top}A(y)y=y^{\top}\left(A(y)^{\top}+A(y)\right)y=0.$
$\Box$
Thus, since the square of the $L^{2}$-norm of the state vector
$\mathbf{\Psi}^{s}(t)$ is conserved in the semi-discrete level, it guarantees
that the system is algebraically stable for both SW and SW square-root
formulations. Note that $\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)$ in
the SW formulation is equivalent to the $L^{2}$-norm of the distribution
function, i.e.
$\int_{\Omega_{x}}\int_{\Omega_{v}}f^{s}(x,v,t)^{2}\mathrm{d}v\mathrm{d}x$,
also known as the enstrophy, and in the SW square-root formulation it is the
discrete equivalent of the number of particles of species $s$, i.e.
$\int_{\Omega_{x}}\int_{\Omega_{v}}f^{s}(x,v,t)\mathrm{d}v\mathrm{d}x$. The SW
conservation of enstrophy is a well-known result shown by [24, 46, 27] with a
Fourier expansion in space, however, this work is the first to casts it under
the umbrella of anti-symmetry and central finite differencing.

###  3.4 Time Reversibility

The Vlasov-Poisson system of equations (1)-(2) are time-reversible, i.e. if
all velocities and time were reversed, the plasma’s prior motions would be
replicated. The anti-symmetric formulation is approximately time reversible
when using temporal integrators that are not time-reversal symmetric [19]. For
example, consider an explicit Euler scheme applied to Eq. (24), such that the
solution at each time iteration $k$ is

$\mathbf{\Psi}^{s}_{k+1}=\left(\mathbf{I}+\Delta
t\mathbf{A}^{s}_{k}\right)\mathbf{\Psi}^{s}_{k},$ (29)

where $\mathbf{I}$ is the identity matrix and $\Delta t$ is the time step.
When $\Delta t$ is sufficiently small, we have

$\left(\mathbf{I}+\Delta
t\mathbf{A}^{s}_{k}\right)^{\top}\left(\mathbf{I}+\Delta
t\mathbf{A}^{s}_{k}\right)=\left(\mathbf{I}-\Delta
t\mathbf{A}^{s}_{k}\right)\left(\mathbf{I}+\Delta
t\mathbf{A}^{s}_{k}\right)=\mathbf{I}+\mathcal{O}(\Delta t^{2}).$

Multiplying Eq. (29) by $(\mathbf{I}+\Delta t\mathbf{A}^{s}_{k})^{\top}$,
results in $\left(\mathbf{I}-\Delta
t\mathbf{A}^{s}_{k}\right)\mathbf{\Psi}^{s}_{k+1}=\mathbf{\Psi}^{s}_{k}+\mathcal{O}(\Delta
t^{2})$, which is approximately (up to second-order accuracy) the time-
reversed version of the explicit Euler method, see Eq. (29). Note that any
time-symmetric Runge-Kutta method is implicit [18, §V.2], and the most
straightforward among these methods is the implicit midpoint method. Although
explicit temporal integrators are not time-symmetric, the anti-symmetric
formulation renders such explicit integrators approximately time-reversible.
We demonstrate the (approximate) time-reversibility of the anti-symmetric
formulation on the linear Landau damping problem in section 5.3 with a third-
order explicit Runge-Kutta method.

###  3.5 Relation to the Anti-symmetric Poisson Bracket

The canonical Poisson bracket [34, 35] is defined as

$\left\\{f,g\right\\}\coloneqq\frac{\partial g}{\partial v}\frac{\partial
f}{\partial x}-\frac{\partial g}{\partial x}\frac{\partial f}{\partial v},$
(30)

for two arbitrary functions $f,g$, which is equivalent to the Jacobian
determinant $\partial(f,g)/\partial(x,v)$. The canonical Poisson bracket in
Eq. (30) is an anti-symmetric operator since
$\left\\{f,g\right\\}=-\left\\{g,f\right\\}$. The Vlasov equation can then be
written as

$\frac{\partial f^{s}}{\partial
t}=\left\\{f^{s},\frac{1}{m^{s}}\frac{\delta\mathcal{E}}{\delta
f^{s}}\right\\},$

where the functional

$\mathcal{E}[f^{s}]\coloneqq\frac{1}{2}\int_{\Omega_{x}}\int_{\Omega_{v}}\sum_{s}m^{s}v^{2}f^{s}(x,v,t)\mathrm{d}v\mathrm{d}x+\frac{1}{2}\int_{\Omega_{x}}E(x,t)^{2}\mathrm{d}x$
(31)

is the total energy of the system, i.e., the Hamiltonian. The variational
derivative222The variational derivative $\frac{\delta\mathcal{E}}{\delta
f^{s}}$ of a functional $\mathcal{E}[f^{s}]$ is defined as
$\int_{\Omega_{v}}\int_{\Omega_{x}}\frac{\delta\mathcal{E}}{\delta
f^{s}}\delta f^{s}\mathrm{d}x\mathrm{d}v\coloneqq\lim_{\varepsilon\to
0}\frac{\mathcal{E}[f^{s}+\varepsilon\delta
f^{s}]-\mathcal{E}[f^{s}]}{\varepsilon}$ . of the total energy is
$\frac{\delta\mathcal{E}}{\delta f^{s}}=\frac{1}{2}m^{s}v^{2}+q^{s}\phi(x,t)$,
where $\phi(x,t)$ is the electric potential, i.e.
$E(x,t)=-\frac{\partial\phi(x,t)}{\partial x}$. The variational derivative is
computed via integration by parts and the Poisson equation (2). The anti-
symmetric advection operator $A^{s}(x,v,t)$, defined in Eq. (21), can also be
written as the canonical Poisson bracket, that is

$A^{s}(x,v,t)=\left\\{\cdot,\frac{1}{m^{s}}\frac{\delta\mathcal{E}}{\delta
f^{s}}\right\\},$

Thus, by preserving the anti-symmetry of the advection operator in semi-
discrete form, the proposed method can also preserve the anti-symmetry of the
canonical Poisson bracket.

The non-canonical Poisson bracket of the Vlasov-Poisson system [34] is defined
as

$\\{\\{\mathcal{F},\mathcal{G}\\}\\}\coloneqq\int_{\Omega_{v}}\int_{\Omega_{x}}\frac{f^{s}}{m^{s}}\left\\{\frac{\delta\mathcal{F}}{\delta
f^{s}},\frac{\delta\mathcal{G}}{\delta f^{s}}\right\\}\mathrm{d}x\mathrm{d}v,$

where $\mathcal{F}[f^{s}]$ and $\mathcal{G}[f^{s}]$ are arbitrary functionals
acting on $f^{s}$. The Vlasov equation (1) can then be represented as

$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{F}=\\{\\{\mathcal{F},\mathcal{E}\\}\\},$

where $\mathcal{F}[f^{s}]$ is any functional of $f^{s}$ and $\mathcal{E}$ is
the total energy from Eq. (31). The noncanonical Poisson bracket is also anti-
symmetric, such that
$\\{\\{\mathcal{F},\mathcal{E}\\}\\}=-\\{\\{\mathcal{E},\mathcal{F}\\}\\}$.
Therefore, by the anti-symmetry of the noncanonical Poisson bracket we know
that
$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{E}=\\{\\{\mathcal{E},\mathcal{E}\\}\\}=-\\{\\{\mathcal{E},\mathcal{E}\\}\\}=0$.
Any discretization that preserves the anti-symmetry of the noncanonical
Poisson bracket will automatically conserve the Hamiltonian of the system,
$\mathcal{E}[f^{s}]$, and other functional invariants, such as Casimirs
$\mathcal{C}[f^{s}]$ (where $\\{\\{\mathcal{C},\mathcal{E}\\}\\}=0$). As we
show in section 4, the SW and SW square-root formulations do not conserve the
total energy $\mathcal{E}[f^{s}]$ of the system and therefore do not preserve
the anti-symmetry of the noncanonical Poisson bracket. The work by [47]
suggests that developing a discretization of the Vlasov-Poisson system with a
discrete (Jacobi) noncanonical Poisson bracket can not be achieved via an
Eulerian solver.

##  4 Conservation Properties of the Symmetrically Weighted and Symmetrically
Weighted Square-Root Formulations

An ideal discretization method should conserve all the invariants of the
Vlasov-Poisson system. This is a challenging problem since the system has an
infinite number of invariants. The Vlasov-Poisson system invariants include
the number of particles, total momentum, total energy, and infinitely many
Casimirs [35], i.e.
$\mathcal{C}[f^{s}]=\int_{\Omega_{v}}\int_{\Omega_{x}}h(f^{s}(x,v,t))\mathrm{d}x\mathrm{d}v$,
where $h(\cdot)$ is any function of the particle distribution function
$f^{s}$. For example, by choosing $h(f^{s})=f^{s}$, the Casimir invariant is
the number of particles of species $s$, and by choosing
$h(f^{s})=f^{s}\log(f^{s})$, the Casimir invariant is the kinetic entropy, and
lastly, by choosing $h(f^{s})=(f^{s})^{p}$ with $p\in\mathbb{N}_{>0}$, the
Casimir invariant is the $L^{p}$-norm of the distribution function [24]. We
focus on analyzing the classical invariants, namely, the number of particles,
total momentum, and total energy.

The conservation constraints for the SW and SW square-root formulations are
due to the SW expansion in velocity, i.e. analogous conservation properties
hold prior to finite difference discretization in space. For each formulation,
we derive the drift rate as the rate at which conservation no longer holds,
which is a quadratic expression of the expansion coefficients and electric
field. The drift rates (aside from total energy in the SW square-root
formulation) depend only on the last expansion coefficient, which we will see
in Theorem 4–8 below, such that in the limit of an infinite number of velocity
spectral terms $N_{v}\to\infty$ the drift rate approaches zero.

Table 1 summarizes the conservation and stability properties of the SW and SW
square-root formulations from section 3 and section 4. We discuss in greater
detail the conservation properties of the SW formulation in section 4.1 and
the SW square-root formulation in section 4.2. In all derivations we use the
trapezoidal rule333The trapezoid rule for approximating the spectral expansion
coefficients $C_{n}^{s}(x,t)$ spatial integral in a periodic domain is given
by $\int_{0}^{\ell}C_{n}^{s}(x,t)\mathrm{d}x\approx\Delta
x\mathbf{1}^{\top}\mathbf{C}_{n}^{s}(t)$, where
$\mathbf{1}\in\mathbb{R}^{N_{x}}$ is a column vector of all ones. to
approximate all integrals in space.

Table 1: A comparison between the SW and SW square-root formulation properties. Properties | SW | SW square-root
---|---|---
conservation of the number of particles | if $N_{v}$ is odd | ✓
conservation of total momentum | if $N_{v}$ is even & $u^{s}=0,\hskip 5.0pt\forall s$ | ✗
conservation of total energy | if $N_{v}$ is odd & $u^{s}=0,\hskip 5.0pt\forall s$ | ✗
positivity preserving | ✗ | ✓
algebraically stable | ✓ | ✓

###  4.1 Symmetrically Weighted Formulation Conservation Properties

We show below that the SW formulation conserves the number of particles if
$N_{v}$ is odd, total momentum if $N_{v}$ is even and $u^{s}=0$ for all $s$,
and total energy if $N_{v}$ is odd and $u^{s}=0$ for all $s$. Thus, the SW
formulation conserves either total momentum or particle number and total
energy, requiring $u^{s}=0$ for all $s$.

######  Theorem 4 (SW formulation particle number conservation).

Let $f^{s}_{\mathrm{sw}}(x,v,t)$ be the solution of the semi-discrete Vlasov-
Poisson equations (24)-(26) solved via the SW formulation. Then, the number of
particles of each plasma species $s$, i.e.

$\mathcal{N}^{s}(t)\coloneqq\int_{0}^{\ell}\int_{\mathbb{R}}f^{s}_{\mathrm{sw}}(x,v,t)\mathrm{d}v\mathrm{d}x,$

is conserved in time if $N_{v}$ is odd, where the drift rate is

$\displaystyle\frac{\mathrm{d}\mathcal{N}^{s}}{\mathrm{d}t}=\begin{cases}0\qquad&\text{if
}N_{v}\text{ is odd},\\\ -\Delta
x\frac{q^{s}}{m^{s}}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)^{\top}\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is even.}\\\ \end{cases}$ (32)

Proof: The number of particles of each plasma species $s$,
$\mathcal{N}^{s}(t)$, in the SW formulation is given by

$\mathcal{N}^{s}(t)\coloneqq\int_{0}^{\ell}\int_{\mathbb{R}}f^{s}_{\mathrm{sw}}(x,v,t)\mathrm{d}v\mathrm{d}x=\Delta
x\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{1}^{\top}\mathbf{C}_{n}^{s}(t),$
(33)

where $\mathbf{1}\in\mathbb{R}^{N_{x}}$ is a column vector of all ones and
$\mathcal{I}_{n}$ is defined in Eq. (13). The derivative in time of the
particle number is

$\displaystyle\frac{\mathrm{d}\mathcal{N}^{s}}{\mathrm{d}t}$
$\displaystyle=\Delta
x\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{1}^{\top}\frac{\mathrm{d}\mathbf{C}_{n}^{s}(t)}{\mathrm{d}t}$
$\displaystyle=\Delta
x\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\left[-\alpha^{s}\underbrace{\mathbf{1}^{\top}\mathbf{D}^{p}\mathbf{Q}^{s}_{n}(t)}_{=0}+\frac{q^{s}}{m^{s}}\mathbf{E}(t)^{\top}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\right]\hskip
5.0pt$ $\displaystyle\text{from Eq. }\eqref{vlasov-ode}$ $\displaystyle=\Delta
x\frac{q^{s}}{m^{s}}\mathbf{E}(t)^{\top}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\hskip
5.0pt$ $\displaystyle\text{from Eqns. }\eqref{anti-
commutative}\leavevmode\nobreak\ \&\leavevmode\nobreak\ \eqref{central-finite-
differencing}$

where
$\mathbf{Q}_{n}^{s}(t)=\alpha^{s}\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)+\alpha^{s}\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)+u^{s}\mathbf{C}_{n}^{s}(t)$.
Since the summation over $n$ is telescopic and from the recurrence relation of
$\mathcal{I}_{n}$ in Eq. (13), we obtain Eq. (32). $\Box$

######  Theorem 5 (SW formulation total momentum conservation).

Let $f^{s}_{\mathrm{sw}}(x,v,t)$ be the solution of the semi-discrete Vlasov-
Poisson equations (24)-(26) solved via the SW formulation. Then, the total
momentum, $\mathcal{P}(t)$, defined as

$\mathcal{P}(t)\coloneqq\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}vf^{s}_{\mathrm{sw}}(x,v,t)\mathrm{d}v\mathrm{d}x,$

is conserved in time if $N_{v}$ is even and $u^{s}=0$ for all $s$, where the
drift rate is

$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}=\begin{cases}-\Delta
x(N_{v}-1)\mathcal{I}_{N_{v}-1}\mathbf{E}(t)^{\top}\sum_{s}\alpha^{s}q^{s}\mathbf{C}_{N_{v}-1}^{s}(t)\qquad&\text{if
}N_{v}\text{ is odd},\\\ -\Delta
x\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)^{\top}\sum_{s}u^{s}q^{s}\mathbf{C}_{N_{v}-1}^{s}(t)\qquad&\text{if
}N_{v}\text{ is even}.\\\ \end{cases}$ (34)

Proof: The total momentum $\mathcal{P}(t)$ in the SW formulation is given by

$\mathcal{P}(t)\coloneqq\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}vf^{s}_{\mathrm{sw}}(x,v,t)\mathrm{d}v\mathrm{d}x=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{1}\mathbf{1}^{\top}\mathbf{C}_{n}^{s}(t),$
(35)

where the integral $\mathcal{I}^{1}_{n}$ is defined as

$\mathcal{I}^{1}_{n}\coloneqq\int_{\mathbb{R}}v\psi_{n}(\xi^{s})\mathrm{d}\xi^{s}=\begin{cases}\alpha^{s}\sqrt{\frac{n+1}{2}}\mathcal{I}_{n+1}+\alpha^{s}\sqrt{\frac{n-1}{2}}\mathcal{I}_{n-1}\qquad&\text{if
}n\geq 1\text{ and odd,}\\\ u^{s}\mathcal{I}_{n}\qquad&\text{if }n\geq 0\text{
and even,}\\\ 0&\text{otherwise.}\end{cases}$

The derivative in time of the total momentum is

$\displaystyle\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}$ $\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{1}\mathbf{1}^{\top}\frac{\mathrm{d}\mathbf{C}_{n}^{s}(t)}{\mathrm{d}t}$
(36) $\displaystyle=\Delta
x\sum_{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{1}\left[-m^{s}\alpha^{s}\underbrace{\mathbf{1}^{\top}\mathbf{D}^{p}\mathbf{Q}^{s}_{n}(t)}_{=0}+q^{s}\mathbf{E}(t)^{\top}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\right]\hskip
5.0pt$ $\displaystyle\text{from Eq. }\eqref{vlasov-ode}$ $\displaystyle=\Delta
x\mathbf{E}(t)^{\top}\sum_{s}q^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{1}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\hskip
5.0pt$ $\displaystyle\text{from Eqns. }\eqref{anti-
commutative}\leavevmode\nobreak\ \&\leavevmode\nobreak\ \eqref{central-finite-
differencing}.$

We divide the summation into even and odd indices and use the recursive
relation of $\mathcal{I}_{n}$ from Eq. (13), such that the even indices result
in

$\sum_{\begin{subarray}{c}n=0\\\ n\text{ is
even}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n}^{1}\left(\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)=\begin{cases}0\qquad&\text{if
}N_{v}\text{ is odd,}\\\
-\sqrt{\frac{N_{v}-1}{2}}u^{s}\mathcal{I}_{N_{v}-2}\mathbf{C}_{N_{v}-1}(t)&\text{if
}N_{v}\text{ is even,}\end{cases}$ (37)

and the odd indices result in

$\displaystyle\sum_{\begin{subarray}{c}n=1\\\ n\text{ is
odd}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n}^{1}\left(\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)$
$\displaystyle=\sum_{\begin{subarray}{c}n=1\\\ \text{n is
odd}\end{subarray}}^{N_{v}-1}\alpha^{s}\mathcal{I}_{n-1}\mathbf{C}_{n-1}^{s}(t)$
(38)
$\displaystyle+\begin{cases}-(N_{v}-1)\mathcal{I}_{N_{v}-1}\alpha^{s}\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is odd,}\\\ 0&\text{if }N_{v}\text{ is even.}\end{cases}$

Thus, by inserting Eqns. (37) and (38) into Eq. (36), we get

$\displaystyle\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}$ $\displaystyle=\Delta
x\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\sum_{\begin{subarray}{c}n=1\\\
n\text{ is
odd}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n-1}\mathbf{C}_{n-1}^{s}(t)$
$\displaystyle+\begin{cases}-\Delta
x(N_{v}-1)\mathcal{I}_{N_{v}-1}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is odd,}\\\ -\Delta
x\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}u^{s}\mathbf{C}^{s}_{N_{v}-1}(t)&\text{if
}N_{v}\text{ is even.}\end{cases}$

The first term in the expression above is equal to zero from the Poisson
equation (26), the recursive relation of $\mathcal{I}_{n}$ in Eq. (13), and
the anti-symmetry properties (19) of the centered finite differencing
derivative operator $(\mathbb{R}^{N_{x}},\mathbf{D}^{p})$ defined in Eq. (28),
i.e.

$\Delta
x\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\sum_{\begin{subarray}{c}n=1\\\
n\text{ is
odd}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n-1}\mathbf{C}_{n-1}^{s}(t)=\Delta
x\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{C}_{n}^{s}(t)=\Delta
x\mathbf{E}(t)^{\top}\mathbf{D}^{p}\mathbf{E}(t)=0,$

which implies that Eq. (34) holds. $\Box$

######  Theorem 6 (SW formulation total energy conservation).

Let $f^{s}_{\mathrm{sw}}(x,v,t)$ be the solution of the semi-discrete Vlasov-
Poisson equations (24)-(26) solved via the SW formulation. Consider the total
energy
$\mathcal{E}(t)\coloneqq\mathcal{E}_{\mathrm{kin}}(t)+\mathcal{E}_{\mathrm{pot}}(t)$,
where

$\mathcal{E}_{\mathrm{kin}}(t)\coloneqq\frac{1}{2}\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}v^{2}f_{\mathrm{sw}}^{s}(x,v,t)\mathrm{d}v\mathrm{d}x\qquad\mathrm{and}\qquad\mathcal{E}_{\mathrm{pot}}(t)\coloneqq\frac{1}{2}\int_{0}^{\ell}E(x,t)^{2}\mathrm{d}x.$

The total energy $\mathcal{E}(t)$ is conserved in time if $N_{v}$ is odd and
$u^{s}=0$ for all $s$, such that the drift rate is

$\frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}=\begin{cases}-\Delta
x(N_{v}-1)\mathcal{I}_{N_{v}-1}\mathbf{E}(t)^{\top}\sum_{s}u^{s}q^{s}\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is odd},\\\ -\Delta
x\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\left((2N_{v}-1)(\alpha^{s})^{2}+(u^{s})^{2}+\frac{q^{s}}{m^{s}}(\mathbf{D}^{p})^{\dagger}\mathbf{E}(t)\odot\right)\mathbf{C}^{s}_{N_{v}-1}(t)&\text{if
}N_{v}\text{ is even,}\end{cases}$ (39)

where $(\mathbf{D}^{p})^{\dagger}\in\mathbb{R}^{N_{x}\times N_{x}}$ is the
pseudoinverse of the finite difference operator $\mathbf{D}^{p}$ defined in
Eq. (28).

Proof: The kinetic energy $\mathcal{E}_{\mathrm{kin}}(t)$ in the SW
formulation is given by

$\mathcal{E}_{\mathrm{kin}}(t)\coloneqq\frac{1}{2}\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}v^{2}f_{\mathrm{sw}}^{s}(x,v,t)\mathrm{d}v\mathrm{d}x=\frac{\Delta
x}{2}\sum_{s}\alpha^{s}m^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{2}\mathbf{1}^{\top}\mathbf{C}_{n}^{s}(t),$
(40)

where the integral $\mathcal{I}^{2}_{n}$ is defined as

$\displaystyle\mathcal{I}_{n}^{2}\coloneqq\begin{cases}2\alpha^{s}u^{s}\left(\sqrt{\frac{n+1}{2}}\mathcal{I}_{n+1}+\sqrt{\frac{n}{2}}\mathcal{I}_{n-1}\right)=2u^{s}\mathcal{I}_{n}^{1}\qquad&\text{if
}n\geq 1\text{ and odd,}\\\
(\alpha^{s})^{2}\left(\sqrt{\frac{(n+1)(n+2)}{4}}\mathcal{I}_{n+2}+\left(\frac{2n+1}{2}+\left(\frac{u^{s}}{\alpha^{s}}\right)^{2}\right)\mathcal{I}_{n}+\sqrt{\frac{n(n-1)}{4}}\mathcal{I}_{n-2}\right)\qquad&\text{if
}n\geq 0\text{ and even,}\\\ 0&\text{otherwise.}\end{cases}$

The derivative in time of the kinetic energy is

$\displaystyle\frac{\mathrm{d}\mathcal{E}_{\mathrm{kin}}}{\mathrm{d}t}$
$\displaystyle=\frac{\Delta
x}{2}\sum_{s}\alpha^{s}m^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{2}\mathbf{1}^{\top}\frac{\mathrm{d}\mathbf{C}_{n}^{s}(t)}{\mathrm{d}t}$
$\displaystyle=\frac{\Delta
x}{2}\sum_{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{2}\left[-m^{s}\alpha^{s}\underbrace{\mathbf{1}^{\top}\mathbf{D}^{p}\mathbf{Q}^{s}_{n}(t)}_{=0}+q^{s}\mathbf{E}(t)^{\top}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\right]\hskip
5.0pt$ $\displaystyle\text{from Eq. }\eqref{vlasov-ode}$
$\displaystyle=\frac{\Delta
x}{2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}^{2}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\hskip
5.0pt$ $\displaystyle\text{from Eqns. }\eqref{anti-
commutative}\leavevmode\nobreak\ \&\leavevmode\nobreak\ \eqref{central-finite-
differencing}$

We first examine the odd indices, where from Eq. (38) we get that

$\displaystyle\frac{\Delta
x}{2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\sum_{\begin{subarray}{c}n=1\\\ n\text{
is odd}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n}^{2}$
$\displaystyle\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)=$
(41) $\displaystyle\begin{cases}-\Delta
x(N_{v}-1)\mathcal{I}_{N_{v}-1}\mathbf{E}(t)^{\top}\sum_{s}u^{s}q^{s}\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is odd,}\\\ 0&\text{if }N_{v}\text{ is even.}\end{cases}$

Similarly, the even indices, result in

$\displaystyle\frac{\Delta x}{2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}$
$\displaystyle\sum_{\begin{subarray}{c}n=1\\\ n\text{ is
even}\end{subarray}}^{N_{v}-1}\mathcal{I}_{n}^{2}\left(\sqrt{\frac{n}{2}}\mathbf{C}^{s}_{n-1}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)=\Delta
x\mathbf{E}(t)^{\top}\mathbf{J}(t)$ (42)
$\displaystyle+\begin{cases}0\qquad&\text{if }N_{v}\text{ is odd,}\\\
-\frac{\Delta
x}{2}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\left((2N_{v}-1)(\alpha^{s})^{2}+(u^{s})^{2}\right)\mathbf{C}_{N_{v}-1}^{s}(t)&\text{if
}N_{v}\text{ is even,}\end{cases}$

where
$\mathbf{J}(t)=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{Q}_{n}^{s}(t)$
is the discretized electric current. The potential energy
$\mathcal{E}_{\mathrm{pot}}(t)$ in the SW formulation is given by

$\mathcal{E}_{\mathrm{pot}}(t)=\frac{1}{2}\int_{0}^{\ell}E(x,t)^{2}\mathrm{d}x=\frac{\Delta
x}{2}\mathbf{E}(t)^{\top}\mathbf{E}(t),$ (43)

and its time derivative is

$\frac{\mathrm{d}\mathcal{E}_{\mathrm{pot}}}{\mathrm{d}t}=\Delta
x\mathbf{E}(t)^{\top}\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}.$ (44)

By taking the time derivative of the semi-discrete Poisson equation (12), we
get

$\displaystyle\mathbf{D}^{p}\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}$
$\displaystyle=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\frac{\mathrm{d}\mathbf{C}_{n}^{s}(t)}{\mathrm{d}t}$
$\displaystyle=-\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\left(\mathbf{D}^{p}\mathbf{Q}_{n}^{s}(t)+\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)\odot\left(\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\right)\hskip
5.0pt$ $\displaystyle\text{from Eq. }\eqref{vlasov-ode}$
$\displaystyle=-\mathbf{D}^{p}\left(\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}\mathbf{Q}_{n}^{s}(t)\right)+\begin{cases}0\qquad&\text{if
}N_{v}\text{ is odd},\\\
-\frac{(q^{s})^{2}}{m^{s}}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)\odot\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is even.}\\\ \end{cases}$
$\displaystyle=-\mathbf{D}^{p}\mathbf{J}(t)+\begin{cases}0\qquad&\text{if
}N_{v}\text{ is odd},\\\
-\frac{(q^{s})^{2}}{m^{s}}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)\odot\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is even,}\\\ \end{cases}$

where $\odot$ denotes the Hadamard (element-wise) product. Thus,

$\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}=-\mathbf{J}(t)+c\mathbf{1}+\begin{cases}0\qquad&\text{if
}N_{v}\text{ is odd},\\\
-(\mathbf{D}^{p})^{\dagger}\frac{(q^{s})^{2}}{m^{s}}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)\odot\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is even.}\\\ \end{cases}$

where $c\in\mathbb{R}$ is an arbitrary constant and
$(\mathbf{D}^{p})^{\dagger}\in\mathbb{R}^{N_{x}\times N_{x}}$ is the
pseudoinverse of the finite difference operator $\mathbf{D}^{p}$ defined in
Eq. (28), and Eq. (44) becomes

$\frac{\mathrm{d}\mathcal{E}_{\mathrm{pot}}}{\mathrm{d}t}=\Delta
x\mathbf{E}(t)^{\top}\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}=-\Delta
x\mathbf{E}(t)^{\top}\mathbf{J}(t)+\begin{cases}0\qquad&\text{if }N_{v}\text{
is odd},\\\ -\Delta
x\mathbf{E}(t)^{\top}(\mathbf{D}^{p})^{\dagger}\frac{(q^{s})^{2}}{m^{s}}\sqrt{\frac{N_{v}-1}{2}}\mathcal{I}_{N_{v}-2}\mathbf{E}(t)\odot\mathbf{C}^{s}_{N_{v}-1}(t)\qquad&\text{if
}N_{v}\text{ is even.}\\\ \end{cases}$ (45)

Summing up Eqns. (41), (42), and (45) results in Eq. (39). $\Box$

###  4.2 Symmetrically Weighted Square-Root Formulation Conservation
Properties

We show below that the SW square-root formulation conserves the number of
particles without restrictions on the number of velocity spectral terms
$N_{v}$ or velocity coordinate scaling and shifting parameters. However, the
total momentum and energy are not conserved in time.

######  Theorem 7 (SW square-root formulation particle number conservation).

Let $f^{s}_{\mathrm{swsr}}(x,v,t)$ be the solution of the semi-discrete
Vlasov-Poisson equations (24) and (27) solved via the SW square-root
formulation. Then, the number of particles of each plasma species $s$, i.e.

$\mathcal{N}^{s}(t)\coloneqq\int_{0}^{\ell}\int_{\mathbb{R}}f^{s}_{\mathrm{swsr}}(x,v,t)\mathrm{d}v\mathrm{d}x,$

is conserved in time, such that
$\frac{\mathrm{d}\mathcal{N}^{s}}{\mathrm{d}t}=0$.

Proof: The number of particles of each plasma species $s$,
$\mathcal{N}^{s}(t)$, in the SW square-root formulation is given by

$\mathcal{N}^{s}(t)\coloneqq\int_{0}^{\ell}\int_{\mathbb{R}}f^{s}_{\mathrm{swsr}}(x,v,t)\mathrm{d}v\mathrm{d}x=\Delta
x\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)=\Delta
x\alpha^{s}\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t).$ (46)

From the anti-symmetric property of the advection operator
$\mathbf{A}^{s}(t)\in\mathbb{R}^{N_{v}N_{x}\times N_{v}N_{x}}$ in semi-
discrete form, see Eq. (24) and Proposition 2, the quadratic term
$\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)$ is an invariant of the semi-
discrete equations. Thus, the particle number $\mathcal{N}^{s}(t)$ of species
$s$ is conserved in the SW square-root formulation. $\Box$

######  Theorem 8 (SW square-root formulation total momentum conservation).

Let $f^{s}_{\mathrm{swsr}}(x,v,t)$ be the solution of the semi-discrete
Vlasov-Poisson equations (24) and (27) solved via the SW square-root
formulation. Then, the total momentum, $\mathcal{P}(t)$, which is defined as

$\mathcal{P}(t)\coloneqq\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}vf^{s}_{\mathrm{swsr}}(x,v,t)\mathrm{d}v\mathrm{d}x,$

is not conserved in time, where the drift rate is

$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}=-\Delta
xN_{v}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\mathbf{C}_{N_{v}-1}^{s}(t)^{2}$
(47)

Proof: The total momentum $\mathcal{P}(t)$ in the SW square-root formulation
is given by

$\displaystyle\mathcal{P}(t)$
$\displaystyle\coloneqq\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}vf^{s}_{\mathrm{swsr}}(x,v,t)\mathrm{d}v\mathrm{d}x$
(48) $\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\left(u^{s}\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)+2\alpha^{s}\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right).$

The time derivative of the total momentum is

$\displaystyle\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}$ $\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\frac{\mathrm{d}}{\mathrm{d}t}\left(u^{s}\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)+2\alpha^{s}\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\left(2u^{s}\underbrace{\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)}_{=0\text{
from Proposition}\leavevmode\nobreak\ \ref{thm:anti-symmetric-
invariants}}+2\alpha^{s}\sum_{n=0}^{N_{v}-1}\sqrt{\frac{n}{2}}\left(\mathbf{C}_{n}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n-1}^{s}(t)+\mathbf{C}_{n-1}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n}^{s}(t)\right)\right).$

From the semi-discrete form of the Vlasov equation (24), we get that

$\displaystyle\sum_{n=0}^{N_{v}-1}\sqrt{\frac{n}{2}}\mathbf{C}_{n}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n-1}^{s}(t)$
$\displaystyle=\sum_{n=0}^{N_{v}-1}-\sqrt{\frac{n}{2}}\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{D}^{p}\left(\alpha^{s}\sqrt{\frac{n-1}{2}}\mathbf{C}_{n-2}^{s}(t)+\alpha^{s}\sqrt{\frac{n}{2}}\mathbf{C}_{n}^{s}(t)+u^{s}\mathbf{C}_{n-1}^{s}(t)\right)$
$\displaystyle+\sum_{n=0}^{N_{v}-1}\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)^{\top}\left(\frac{\sqrt{n(n-1)}}{2}\mathbf{C}_{n-2}^{s}(t)\odot\mathbf{C}_{n}^{s}(t)-\frac{n}{2}\mathbf{C}_{n}^{s}(t)^{2}\right),$
$\displaystyle\sum_{n=0}^{N_{v}-1}\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n}^{s}(t)$
$\displaystyle=\sum_{n=0}^{N_{v}-1}-\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)^{\top}\mathbf{D}^{p}\left(\alpha^{s}\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)+\alpha^{s}\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)+u^{s}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle+\sum_{n=0}^{N_{v}-1}\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)^{\top}\left(\frac{n}{2}\mathbf{C}_{n-1}^{s}(t)^{2}-\frac{\sqrt{n(n+1)}}{2}\mathbf{C}_{n-1}^{s}(t)\odot\mathbf{C}_{n+1}^{s}(t)\right).$

The telescoping sum results in

$\displaystyle\sum_{n=0}^{N_{v}-1}\sqrt{\frac{n}{2}}\left(\mathbf{C}_{n}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n-1}^{s}(t)+\mathbf{C}_{n-1}^{s}(t)^{\top}\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle=\frac{q^{s}}{2m^{s}\alpha^{s}}\mathbf{E}(t)^{\top}\left(\left(\sum_{n=0}^{{N_{v}-1}}\mathbf{C}_{n}^{s}(t)^{2}\right)-N_{v}\mathbf{C}_{N_{v}-1}^{s}(t)^{2}\right).$

Thus,

$\displaystyle\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}t}$ $\displaystyle=\Delta
x\mathbf{E}(t)^{\top}\underbrace{\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)^{2}}_{=\mathbf{D}^{p}\mathbf{E}(t)\text{
from Eq.}\leavevmode\nobreak\ \eqref{poisson-ode-sw-sqrt}}-\Delta
xN_{v}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\mathbf{C}_{N_{v}-1}^{s}(t)^{2}$
$\displaystyle=-\Delta
xN_{v}\mathbf{E}(t)^{\top}\sum_{s}q^{s}\alpha^{s}\mathbf{C}_{N_{v}-1}^{s}(t)^{2}\qquad$
$\displaystyle\text{from Eqns.}\leavevmode\nobreak\ \eqref{preserves-square-
norms}\leavevmode\nobreak\ \&\leavevmode\nobreak\ \eqref{central-finite-
differencing}.$

$\Box$

######  Theorem 9 (SW square-root formulation total energy conservation).

Let $f^{s}_{\mathrm{swsr}}(x,v,t)$ be the solution of the semi-discrete
Vlasov-Poisson equations (24) and (27) solved via the SW square-root
formulation. Consider the total energy
$\mathcal{E}(t)\coloneqq\mathcal{E}_{\mathrm{kin}}(t)+\mathcal{E}_{\mathrm{pot}}(t)$,
where

$\mathcal{E}_{\mathrm{kin}}(t)\coloneqq\frac{1}{2}\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}v^{2}f_{\mathrm{swsr}}^{s}(x,v,t)\mathrm{d}v\mathrm{d}x\qquad\mathrm{and}\qquad\mathcal{E}_{\mathrm{pot}}(t)\coloneqq\frac{1}{2}\int_{0}^{\ell}E(x,t)^{2}\mathrm{d}x.$

The total energy $\mathcal{E}(t)$ is not conserved in time, where the drift
rate is

$\displaystyle\frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}=$
$\displaystyle-\Delta
x\frac{N_{v}}{2}\sqrt{\frac{N_{v}-1}{2}}\sum_{s}(\alpha^{s})^{2}\left(m^{s}(\alpha^{s})^{2}\mathbf{C}_{N_{v}-1}^{s}(t)^{\top}\mathbf{D}^{p}\mathbf{C}_{N_{v}-2}^{s}(t)+q^{s}\mathbf{E}(t)^{\top}\left(\mathbf{C}_{N_{v}-2}^{s}(t)\odot\mathbf{C}_{N_{v}-1}^{s}(t)\right)\right)$
$\displaystyle-\Delta
xN_{v}\mathbf{E}(t)^{\top}\sum_{s}q^{s}u^{s}\alpha^{s}\mathbf{C}^{s}_{N_{v}-1}(t)^{2}+\Delta
x\mathbf{E}(t)^{\top}\sum_{s}\sum_{n=0}^{N_{v}-1}q^{s}\alpha^{s}\mathbf{C}_{n}^{s}(t)\odot\mathbf{Q}_{n}^{s}(t)$
(49) $\displaystyle-2\Delta
x\mathbf{E}(t)^{\top}(\mathbf{D}^{p})^{\dagger}\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)\odot\mathbf{D}^{p}\mathbf{Q}_{n}^{s}(t).$

Proof: The kinetic energy $\mathcal{E}_{\mathrm{kin}}(t)$ in the SW square-
root formulation is given by

$\displaystyle\mathcal{E}_{\mathrm{kin}}(t)$
$\displaystyle\coloneqq\frac{1}{2}\sum_{s}m^{s}\int_{0}^{\ell}\int_{\mathbb{R}}v^{2}\frac{\partial
f^{s}_{\mathrm{swsr}}(x,v,t)}{\partial t}\mathrm{d}v\mathrm{d}x$ (50)
$\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\frac{1}{2}\left((u^{s})^{2}+(\alpha^{s})^{2}n+\frac{(\alpha^{s})^{2}}{2}\right)\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle+\sum_{n=0}^{N_{v}-1}2u^{s}\alpha^{s}\sqrt{\frac{n}{2}}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-1}^{s}(t)\right)+\frac{(\alpha^{s})^{2}}{2}\sum_{n=0}^{N_{v}-1}\sqrt{n(n-1)}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-2}^{s}(t)\right).$

The derivative in time of the kinetic energy is

$\displaystyle\frac{\mathrm{d}\mathcal{E}_{\mathrm{kin}}}{\mathrm{d}t}$
$\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\frac{1}{2}\left((u^{s})^{2}+(\alpha^{s})^{2}n+\frac{(\alpha^{s})^{2}}{2}\right)\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle+\sum_{n=0}^{N_{v}-1}2u^{s}\alpha^{s}\sqrt{\frac{n}{2}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-1}^{s}(t)\right)+\frac{(\alpha^{s})^{2}}{2}\sum_{n=0}^{N_{v}-1}\sqrt{n(n-1)}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-2}^{s}(t)\right)\hskip
5.0pt$ $\displaystyle\text{from Eq.}\leavevmode\nobreak\ \eqref{identity-2}$
$\displaystyle=\Delta
x\sum_{s}m^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\frac{n(\alpha^{s})^{2}}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right)+\underbrace{\left(\frac{(u^{s})^{2}}{2}+\frac{(\alpha^{s})^{2}}{4}\right)\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Psi}^{s}(t)^{\top}\mathbf{\Psi}^{s}(t)}_{=0}$
from Proposition2
$\displaystyle+\sum_{n=0}^{N_{v}-1}2u^{s}\alpha^{s}\sqrt{\frac{n}{2}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-1}^{s}(t)\right)+\frac{(\alpha^{s})^{2}}{2}\sum_{n=0}^{N_{v}-1}\sqrt{n(n-1)}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-2}^{s}(t)\right).$

We expand each term in the expression above, such that

$\displaystyle\sum_{n=0}^{N_{v}-1}\frac{n(\alpha^{s})^{2}}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n}^{s}(t)\right)$
$\displaystyle=\sum_{n=0}^{N_{v}-1}-(\alpha^{s})^{3}\sqrt{\frac{n}{2}}\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{D}^{p}\mathbf{C}^{s}_{n-1}+\frac{q^{s}\alpha^{s}}{m^{s}}\sqrt{\frac{n}{2}}\mathbf{E}^{\top}\left(\mathbf{C}_{n-1}^{s}(t)\odot\mathbf{C}_{n}^{s}(t)\right),$
$\displaystyle\sum_{n=0}^{N_{v}-1}2u^{s}\alpha^{s}\sqrt{\frac{n}{2}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-1}^{s}(t)\right)$
$\displaystyle=\frac{q^{s}u^{s}}{m^{s}}\mathbf{E}(t)^{\top}\left(\left(\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)^{2}\right)-N_{v}\mathbf{C}_{N_{v}-1}^{s}(t)^{2}\right),$
$\displaystyle\sum_{n=0}^{N_{v}-1}\frac{(\alpha^{s})^{2}}{2}\sqrt{n(n-1)}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathbf{C}_{n}^{s}(t)^{\top}\mathbf{C}_{n-2}^{s}(t)\right)$
$\displaystyle=\sum_{n=0}^{N_{v}-1}-(\alpha^{s})^{3}\sqrt{\frac{n-1}{2}}\mathbf{C}^{s}_{n-2}(t)^{\top}\mathbf{D}^{p}\mathbf{C}_{n-1}^{s}(t)$
$\displaystyle+\sum_{n=0}^{N_{v}-1}\frac{q^{s}\alpha^{s}}{m^{s}}\sqrt{\frac{n-1}{2}}\mathbf{E}(t)^{\top}\left(\mathbf{C}_{n-1}^{s}(t)\odot\mathbf{C}_{n-2}^{s}(t)\right)$
$\displaystyle-\frac{q^{s}\alpha^{s}}{2m^{s}}(N_{v}-2)\sqrt{\frac{N_{v}-1}{2}}\mathbf{E}(t)^{\top}\left(\mathbf{C}_{N_{v}-2}^{s}(t)\odot\mathbf{C}_{N_{v}-1}^{s}(t)\right)$
$\displaystyle-\frac{(\alpha^{s})^{3}}{2}(N_{v}-2)\sqrt{\frac{N_{v}-1}{2}}\mathbf{C}_{N_{v}-1}^{s}(t)^{\top}\mathbf{D}^{p}\mathbf{C}_{N_{v}-2}^{s}(t).$

From the anti-commutative property of the anti-symmetric operator
$(\mathbb{R}^{N_{x}},\mathbf{D}^{p})$, see Eqns. (18) and (28), we end up with

$\displaystyle\frac{\mathrm{d}\mathcal{E}_{\mathrm{kin}}}{\mathrm{d}t}=$
$\displaystyle-\Delta
x\frac{N_{v}}{2}\sqrt{\frac{N_{v}-1}{2}}\sum_{s}(\alpha^{s})^{2}\left(m^{s}(\alpha^{s})^{2}\mathbf{C}_{N_{v}-1}^{s}(t)^{\top}\mathbf{D}^{p}\mathbf{C}_{N_{v}-2}^{s}(t)+q^{s}\mathbf{E}(t)^{\top}\left(\mathbf{C}_{N_{v}-2}^{s}(t)\odot\mathbf{C}_{N_{v}-1}^{s}(t)\right)\right)$
$\displaystyle-\Delta
xN_{v}\mathbf{E}(t)^{\top}\sum_{s}q^{s}u^{s}\alpha^{s}\mathbf{C}^{s}_{N_{v}-1}(t)^{2}+\Delta
x\mathbf{E}(t)^{\top}\sum_{s}\sum_{n=0}^{N_{v}-1}q^{s}\alpha^{s}\mathbf{C}_{n}^{s}(t)\odot\mathbf{Q}_{n}^{s}(t).$

The potential energy $\mathcal{E}_{\mathrm{pot}}(t)$ is given by

$\frac{\mathrm{d}\mathcal{E}_{\mathrm{pot}}}{\mathrm{d}t}\coloneqq\frac{1}{2}\int_{0}^{\ell}E(x,t)^{2}\mathrm{d}x=\frac{\Delta
x}{2}\mathbf{E}(t)^{\top}\mathbf{E}(t).$ (51)

The derivative in time of the potential energy is

$\frac{\mathrm{d}\mathcal{E}_{\mathrm{pot}}}{\mathrm{d}t}=\Delta
x\mathbf{E}(t)^{\top}\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}.$

By taking the time derivative of the semi-discrete Poisson equation (27), we
get

$\displaystyle\mathbf{D}^{p}\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}$
$\displaystyle=\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\frac{\mathrm{d}\mathbf{C}_{n}^{s}(t)^{2}}{\mathrm{d}t}$
$\displaystyle=-2\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)\odot\left(\mathbf{D}^{p}\mathbf{Q}_{n}^{s}(t)+\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}(t)^{\top}\left(\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)\
-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right)\right)$
$\displaystyle=-2\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)\odot\mathbf{D}^{p}\mathbf{Q}_{n}^{s}(t).$

Thus,

$\frac{\mathrm{d}\mathbf{E}(t)}{\mathrm{d}t}=-2(\mathbf{D}^{p})^{\dagger}\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathbf{C}_{n}^{s}(t)\odot\mathbf{D}^{p}\mathbf{Q}_{n}^{s}(t),$
(52)

resulting in Eq. (49). $\Box$

######  Remark 3.

The finite difference operator $\mathbf{D}^{p}\in\mathbb{R}^{N_{x}\times
N_{x}}$ defined in Eq. (28) does not preserve the product rule, i.e.
$\mathbf{D}^{p}(\mathbf{y}\odot\mathbf{z})\neq(\mathbf{D}^{p}\mathbf{y})\odot\mathbf{z}+\mathbf{y}\odot(\mathbf{D}^{p}\mathbf{z})$.
Otherwise, the last two terms in the energy drift in Eq. (49) cancel out and
the energy drift is a function of the electric field and the last two spectral
coefficient terms. On the other hand, any spectral discretization preserves
the product rule, such as the Fourier spectral discretization, see A.

##  5 Numerical Results

We numerically investigate the accuracy, convergence, and conservation
properties of the SW and SW square-root formulations. The numerical setup for
the test cases is described in section 5.1. We test the SW and SW square-root
formulations on classical benchmark problems: (1) manufactured solution in
section 5.2, (2) linear and nonlinear Landau damping in section 5.3, (3) two-
stream instability in section 5.4, and (4) bump-on-tail instability in section
5.5.

###  5.1 Numerical Setup

The number of particles, total momentum, and total energy are either linear or
quadratic invariants of the system, see Eqns. (33)(35)(40)(43) for the SW
formulation and Eqns. (46)(48)(50)(51) for the SW square-root formulation.
These can be conserved at the fully discrete level using implicit Gauss-
Legendre temporal integrators. We use the second-order implicit midpoint
temporal integrator, which is part of the Gauss-Legendre family of
integrators. The nonlinear system at each time step is solved using an
unpreconditioned Jacobian-Free-Newton-Krylov (JFNK) method [26] implemented in
Python’s scipy package with relative tolerance set to $10^{-8}$ and absolute
tolerance set to $10^{-14}$. The Krylov linear solver that is embedded in the
JFNK method is the default linear solver restarted generalized minimum
residual method (LGMRES) [3] with relative and absolute tolerance set to
$10^{-5}$. The semi-discrete Poisson equations (26)-(27) are solved using the
generalized minimal residual method (GMRES) [45] with a relative tolerance and
absolute tolerance of $10^{-15}$. We set the ions as a static neutralizing
background since the ions are approximately motionless on the electron time
scales. The number of grid points in space is set to $N_{x}=100$ and the time
step to $\Delta t=10^{-2}$, unless specified otherwise. Additionally, we use a
second-order central finite differencing derivative approximation in space. In
contrast to asymmetrically-weighted solvers [6], we do not consider an
artificial collisional operator in this work. The initial electron
distribution for all test cases, aside from the manufactured solution, is
represented by

$f^{e}(x,v,t=0)=\frac{n_{0}(1+\epsilon\cos(kx))}{\sqrt{\pi}}\begin{cases}\frac{1}{\sqrt{2}\alpha^{e}_{\mathrm{sw}}}\exp\left(-\frac{1}{2}\left(\frac{v-u^{e}}{\alpha^{e}_{\mathrm{sw}}}\right)^{2}\right),\qquad&\text{SW
formulation}\\\
\frac{1}{\alpha^{e}_{\mathrm{swsr}}}\exp\left(-\left(\frac{v-u^{e}}{\alpha^{e}_{\mathrm{swsr}}}\right)^{2}\right),\qquad&\text{SW
square-root formulation}\end{cases}$

where $\epsilon$ is the amplitude of the initial perturbation, $n_{0}$ is the
average density, $\alpha^{e}_{\mathrm{sw}}$ is the electron velocity scaling
for the SW formulation, $\alpha^{e}_{\mathrm{swsr}}$ is the electron velocity
scaling for the SW square-root formulation, $u^{e}$ is the electron velocity
shifting, and $k$ is the perturbation’s wavenumber. The parameters for each
numerical test (aside from the manufactured solution test) are shown in Table
2. Moreover, we set a realistic ion-to-electron mass ratio $m_{i}/m_{e}=1836$.

Table 2: Parameter setup for the linear Landau damping, nonlinear Landau damping, two-stream instability, and bump-on-tail instability. The subscripts ‘sw’ and ‘swsr’ correspond to the SW and SW square-root formulations, respectively. Parameter | Linear Landau | Nonlinear Landau | Two-stream | Bump-on-tail
---|---|---|---|---
$\ell$ | $2\pi$ | $4\pi$ | $2\pi$ | $20\pi$
$k$ | $1$ | $0.5$ | $1$ | $0.3$
$\epsilon$ | $0.01$ | $0.5$ | $10^{-3}$ | $0.03$
$t_{f}$ | $10$ | $10$ | $45$ | $20$
$n_{0}$ | $1$ | $1$ |  $0.5$ & $0.5$ |  $0.9$ & $0.1$
$u^{e}$ | 0 | 0 |  $-1$ & $1$ |  $0$ & $4.5$
$\alpha^{e}_{\mathrm{sw}}$ | $1$ | $1$ |  $0.5/\sqrt{2}$ & $0.5/\sqrt{2}$ |  $1$ & $0.5$
$\alpha^{e}_{\mathrm{swsr}}$ | $\sqrt{2}$ | $\sqrt{2}$ |  $0.5$ & $0.5$ |  $\sqrt{2}$ & $1/\sqrt{2}$

###  5.2 Manufactured Solution

We test the accuracy of the SW and SW square-root formulations using the
method of manufactured solutions. We consider the transport of uncharged
particles corresponding to

$\left(\frac{\partial}{\partial t}+v\frac{\partial}{\partial
x}\right)f^{m}(x,v,t)=0\qquad\text{or}\qquad\left(\frac{\partial}{\partial
t}+v\frac{\partial}{\partial x}\right)\sqrt{f^{m}(x,v,t)}=0\qquad\text{in
}\Omega_{x}\times\Omega_{v}\times\Omega_{t}.$

We set $\ell=2\pi$, $t_{f}=1$, and the initial condition to

$f^{m}(x,v,t=0)=f_{0}^{m}(x,v)=\pi^{-\frac{1}{4}}(2-\cos(x))^{2}\exp{\left(-\left(\frac{v-1}{2}\right)^{2}\right)}.$

From the method of characteristics, the analytic solution to the linear
advection equation is given by $f^{m}(x,v,t)=f_{0}^{m}(x-vt,v)$. We set
$\Delta t=10^{-3}$ for both formulations. The initial condition for the SW
formulation is parameterized with $\alpha_{\mathrm{sw}}^{m}=\sqrt{2}$,
$u_{\mathrm{sw}}^{m}=1$,
$\mathbf{C}^{m}_{0,\mathrm{sw}}(t=0)=\left(2-\cos(\mathbf{x})\right)^{2}$, and
$\mathbf{C}^{m}_{n,\mathrm{sw}}(t=0)=0$ for $n=1,2,\ldots,N_{v}-1$, where
$\mathbf{x}=[x^{(1)},x^{(2)},\ldots,x^{(N_{x})}]^{\top}\in\mathbb{R}^{N_{x}}$.
Similarly, the SW square-root formulation is parameterized with
$\alpha_{\mathrm{swsr}}^{m}=2$, $u_{\mathrm{swsr}}^{m}=1$,
$\mathbf{C}^{m}_{0,\mathrm{swsr}}(t=0)=\pi^{\frac{1}{8}}(2-\cos(\mathbf{x}))$,
and $\mathbf{C}^{m}_{n,\mathrm{swsr}}(t=0)=0$ for $n=1,2,\ldots N_{v}-1$. The
governing equations for the expansion coefficients in the SW and SW square-
root formulations are of the same form, i.e.

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}^{m}_{n}(t)=-\mathbf{D}^{p=2}\left(\alpha^{m}\sqrt{\frac{n}{2}}\mathbf{C}^{m}_{n-1}(t)+\alpha^{m}\sqrt{\frac{n+1}{2}}\mathbf{C}^{m}_{n+1}(t)+u^{m}\mathbf{C}^{m}_{n}(t)\right).$

Figure 1 shows the SW and SW square-root numerical results at $t=1$, along
with the relative error with $N_{v}=N_{x}=100$. The numerical results show
that the two formulations are in good agreement with the analytic solution and
have comparable accuracy. The relative error in the SW formulation indicates
the interaction of higher-order modes in comparison to the SW square-root
formulation. Figure 2(a) shows that indeed the SW formulation requires $\sim
44$ Hermite modes to approximate the manufactured solution, whereas the SW
square-root formulation only requires $\sim 36$ modes to approximate the
manufactured solution. Figure 2(b) presents the $L^{2}$-norm absolute error of
the two formulations as a function of the number of spatial grid points
$N_{x}$. The $L^{2}$-norm absolute error is computed using $10^{3}$ grid
points uniformly spaced in the velocity direction from $v\in[-5,5]$ with
$N_{x}=50,100,200,400$. The spatial resolution convergence rate is $2$, which
agrees with the analytic second-order central finite difference convergence.
Figure 2(b) also shows that the SW square-root formulation is slightly more
accurate than the SW formulation for the manufactured test case.

(a) SW formulation

(b) SW square-root formulation

(c) SW formulation relative error

(d) SW square-root formulation relative error

Figure 1: The (a) SW and (b) SW square-root formulations manufactured solution
distribution function $f^{m}(x,v,t=1)$ results with $N_{x}=N_{v}=100$.
Subfigures (c) and (d) show the relative error between the analytic solution
and the SW and SW square-root numerical results, respectively. The numerical
results show that as expected the SW and SW square-root formulation accuracy
is comparable, yet the SW formulation has more interaction between higher-
order modes.

(a) Velocity spectral convergence

(b) Spatial finite difference convergence

Figure 2: Subfigure (a) shows the $L^{2}$-norm of the expansion coefficients
with $N_{v}=N_{x}=100$. The spectral convergence illustrates that the SW
square-root formulation requires fewer spectral coefficients in comparison to
the SW formulation to approximate the manufactured solution. Subfigure (b)
shows the $L^{2}$-norm absolute error of the SW and SW square-root
manufactured solutions. In subfigure (b), we set $N_{v}=100$ and
$N_{x}=50,100,200,400$. The numerical convergence rate is $2$ for both
formulations, which agrees with the second-order central finite difference
convergence rate.

###  5.3 Landau Damping

Landau damping is a classic benchmark problem for kinetic plasma codes caused
by wave-particle resonance [31], where particles gain energy from or lose
energy to the wave. In a near-Maxwellian distribution, more particles have
velocities slightly slower than the wave phase velocity so there is a net
transfer of energy from the wave to the particles, and as a result, the
electric field damps and particles gain energy. We examine the linear and
nonlinear Landau damping test cases in section 5.3.1 and section 5.3.2,
respectively.

####  5.3.1 Linear Landau Damping

For a small perturbation in the initial electron distribution, Landau damping
can be described using linear theory. We set the initial perturbation
amplitude to $\epsilon=0.01$. Figure 3 shows the conservation of particle
number, momentum, and energy of the SW and SW square-root formulations with
$N_{v}=100$ and $N_{v}=101$. Since we set the electron velocity shifting
parameter to $u^{e}=0$, the particle number and energy are conserved in the SW
formulation when $N_{v}$ is odd, as shown in Figure 3(a), and momentum is
conserved when $N_{v}$ is even, as shown in Figure 3(b). On the contrary, the
particle number is conserved in the SW square-root formulation for any
$N_{v}$. Figure 3(c) and Figure 3(d) show the SW square-root formulation
conservation properties for $N_{v}=101$ and $N_{v}=100$, respectively, which
are almost identical, unlike the SW formulation conservation properties which
depend on $N_{v}$ being even or odd. Figure 3 also shows the drift rate in the
conservation for each setup. The numerical results show that the analytic
conservation properties derived in section 4 match the numerical drift rate in
the conservation of particle number, momentum, and energy. Despite the
momentum drift rates given by Eqns. (34) and (47), momentum is conserved close
to machine precision in the SW square-root formulation (for both $N_{v}=100$
and $N_{v}=101$) and the SW formulation (with $N_{v}=101$). This is because
the momentum drift rate is a function of the product of the electric field and
the last spectral coefficient, which are both small due to wave damping and
spectral convergence.

As mentioned in section 3.4, the anti-symmetry of the semi-discrete Eq. (24)
results in explicit Runge-Kutta temporal integrators being approximately time-
reversal symmetric. We numerically demonstrate this with the SW square-root
formulation. Figure 4 shows the electric field amplitude as a function of time
for the forward and backward in-time simulations using the non-adaptive 3rd-
order explicit Runge-Kutta integrator of Bogacki-Shampine [5] with $\Delta
t=10^{-3}$ and $N_{v}=100$. Both the forward and backward in-time simulations
agree with the linear theory electric field damping rate444The linear theory
damping rate $\gamma=-0.851$ for the linear Landau damping test is obtained by
solving the dispersion equation numerically with wavenumber $k=1$ [7, Table
1], i.e. $2+\zeta Z(\zeta)=0$, where $\zeta=(w_{r}+\gamma i)/\sqrt{2}$ and
$Z(\zeta)=\pi^{-1/2}\int_{-\infty}^{\infty}(x-\zeta)^{-1}\exp(-x^{2})\mathrm{d}x$
is the plasma dispersion function, see [15, §2] for a detailed derivation. .
The numerical results show that although explicit Runge-Kutta temporal
integrators are not time-reversal symmetric, the anti-symmetric structure of
the semi-discrete equations (24) makes such temporal integrators approximately
time-reversible.

(a) SW formulation $N_{v}=101$ (odd)

(b) SW formulation $N_{v}=100$ (even)

(c) SW square-root formulation $N_{v}=101$ (odd)

(d) SW square-root formulation $N_{v}=100$ (even)

Figure 3: Conservation of particle number, momentum, and energy via the SW
formulation (a/b) and the SW square-root formulation (c/d) for linear Landau
damping. Subfigures (a/c) have an odd number of spectral terms ($N_{v}=101$)
and subfigures (b/d) have an even number of spectral terms ($N_{v}=100$). The
numerical drift rate (the rate at which conservation no longer holds) in each
setup matches the analytic drift rate, which we derive in Eqns. (32), (34),
(39), (47), and (49). Figure 4: The linear Landau damping electric field
amplitude as a function of time computed via the SW square-root formulation.
The numerical damping rate agrees with the theoretical damping rate from [7,
Table 1]. The forward and backward simulations are evolved via the non-
adaptive 3rd-order explicit Runge-Kutta integrator of Bogacki-Shampine [5]
with $\Delta t=10^{-3}$ and $N_{v}=100$. Note that such an explicit Runge-
Kutta temporal integrator is not time-reversal symmetric, yet the anti-
symmetric structure of the equations makes this method approximately time-
reversible.

####  5.3.2 Nonlinear Landau Damping

The linear theory fails to describe the damping rate if the magnitude of the
initial perturbation is large as nonlinear effects become important. We set
the initial perturbation amplitude for the nonlinear Landau damping test case
to $\epsilon=0.5$. The electron distribution function $f^{e}(x,v,t)$ at
various time instances is shown in Figure 5. Figure 5(a) illustrates that the
SW formulation results become negative in phase space starting from $t=5$. By
construction, the SW square-root formulation preserves the distribution
function’s non-negative property, as shown in Figure 5(b). At $t=10$, it
appears that the SW formulation results display a slightly higher degree of
filamentation than the SW formulation results, as shown in Figure 6.

Figure 7 shows the conservation of particle number, momentum, and energy for
the nonlinear Landau damping test case. Similar to the linear Landau damping
test case, we set the velocity shifting parameter $u^{e}=0$, thus the SW
formulation conserves particle number and energy when $N_{v}$ is odd as shown
in Figure 7(a) and conserves momentum when $N_{v}$ is even as shown in Figure
7(b). As predicted, the SW square-root formulation conserves particle number,
independently of $N_{v}$, as shown in Figure 7(a) (with $N_{v}=101$) and
Figure 7(d) (with $N_{v}=100$). Interestingly, although energy and momentum
are not conserved in the SW square-root formulation, their drift rate seems to
grow very slowly in time.

(a) SW formulation

(b) SW square-root formulation

Figure 5: The evolution of the electron distribution function $f^{e}(x,v,t)$
for nonlinear Landau damping with $N_{v}=100$. The distribution function
becomes negative in the (a) SW formulation whereas the (b) SW square-root
formulation is positivity preserving by construction.  Figure 6: A cross-
section of the nonlinear Landau damping electron particle distribution
function $f^{e}(x=2\pi,v,t=10)$ via the two formulations: SW and SW square-
root (with $N_{v}=101$). The SW formulation appears to have more filamentation
near regions close to zero.

(a) SW formulation $N_{v}=101$ (odd)

(b) SW formulation $N_{v}=100$ (even)

(c) SW square-root formulation $N_{v}=101$ (odd)

(d) SW square-root formulation $N_{v}=100$ (even)

Figure 7: Same as Figure 3 for nonlinear Landau damping. The numerical drift
rate for the SW square-root formulation is independent of $N_{v}$, whereas the
SW formulation drift rate depends on whether $N_{v}$ is either odd or even.

###  5.4 Two-stream Instability

The two-stream instability is a classic velocity-space microinstability
corresponding to two counter-streaming electron beams. Figure 8 shows the
electron distribution function in phase space at $t=45$ with $N_{v}=101$. The
SW formulation results become negative in certain regions of phase space,
which is relatively benign as such regions are far from the phase space
vortex. Figure 9(a) shows a cross-section of the electron distribution
function for the two formulations, which indicates along with Figure 8, that
the SW square-root formulation exhibits more filamentation and the negativity
of the distribution function in SW is insignificant. In contrast to the
nonlinear Landau damping test case, the SW square-root formulation results
display a higher degree of filamentation in comparison to the SW formulation
results. This suggests that the filamentation effects are problem-dependent.
The numerical electric field growth rate for the SW and SW square-root
formulations with $N_{v}=101$ is shown in Figure 9(b). The numerical growth
rate for both formulations agrees with the linear theory growth rate555The
linear theory growth rate $\gamma=0.185$ for the two-stream instability test
is obtained by solving numerically the dispersion equation for $k=1$
wavenumber, i.e. $9+4\zeta_{1}Z(\zeta_{1})+4\zeta_{2}Z(\zeta_{2})=0$, where
$\zeta_{1}=2(w_{r}-1+\gamma i)$, $\zeta_{2}=2(w_{r}+1+\gamma i)$, and
$Z(\zeta)=\pi^{-1/2}\int_{-\infty}^{\infty}(x-\zeta)^{-1}\exp(-x^{2})\mathrm{d}x$,
see [15, §3] for a detailed derivation..

Figure 10 compares the SW and SW square-root conservation laws for the two-
stream instability. Since the velocity shifting parameter $u^{s}\neq 0$, the
SW formulation only conserves particle number when $N_{v}$ is odd, see Figure
10(a), and does not conserve particle number, momentum, or energy when $N_{v}$
is even, see Figure 10(a). However, the momentum drift rate is very small in
both cases. The SW square-root formulation conserves particle number
regardless of whether $N_{v}$ is even or odd and momentum drift is negligible
as shown in Figure 10(c) and Figure 10(d). Here, the momentum drifts for both
the SW and SW square-root formulations, see Eqn (34) and (47), are negligible
due to the phase space symmetry.

(a) SW formulation

(b) SW square-root formulation

Figure 8: The two-stream instability electron distribution function
$f^{e}(x,v,t=45)$ with $N_{v}=101$ obtained by the (a) SW and the (b) SW
square-root formulations. The SW formulation becomes negative in certain
regions in phase space, yet such regions are far from the region of interest,
i.e. the phase space vortex.

(a) Two-stream instability filamentation

(b) Two-stream instability growth rate

Figure 9: Subfigure (b) shows a cross-section of the two-stream instability
electron particle distribution function $f^{e}(x=\frac{\pi}{4},v,t=45)$ via
the two formulations: SW and SW square-root (with $N_{v}=101$). The SW square-
root formulation appears to have more filamentation near the peak regions.
Subfigure (b) presents the two-stream instability electric field amplitude as
a function of time computed via the SW and SW square-root formulation with
$N_{v}=101$. The electric field growth rate of the SW and SW square-root
simulations agree with the theoretical growth rate [15, §3].

(a) SW formulation $N_{v}=101$ (odd)

(b) SW formulation $N_{v}=100$ (even)

(c) SW square-root formulation $N_{v}=101$ (odd)

(d) SW square-root formulation $N_{v}=100$ (even)

Figure 10: Same as Figure 3 for the two-stream instability. Due to the phase
space symmetry in the two-stream instability, momentum is approximately
conserved in both formulations.

###  5.5 Bump-on-tail Instability

The bump-on-tail instability is an asymmetric velocity-space instability.
Figure 11 shows the SW and SW square-root formulations distribution function
evolution. The SW formulation results become negative at $t=10$ and at $t=20$
the distribution function becomes negative in the the majority of phase space
along with the phase space vortex between the two beams. Figure 12(a) shows a
cross-section of the electron distribution function at $t=20$, indicating that
the SW and SW square-root formulations suffer from comparable filamentation
effects, manifesting as artificial oscillations of the distribution function.
Figure 12(b) shows the kinetic and potential energy exchange from the two
formulations. Such macroscopic quantities are nearly identical for the two
formulations.

The conservation of particle number, momentum, and energy is shown in Figure
13. The numerical drift rates match the analytically derived drift rates for
both formulations. Similar to the two-stream instability, the velocity
shifting parameter $u_{e}\neq 0$, thus only particle number is conserved in
the SW formulation if the number of spectral terms is odd. The numerical drift
rates for both formulations are comparable in magnitude. The momentum drift
rate for the SW formulation (with $N_{v}=100$) and SW square-root formulations
is larger than the two-stream instability drift rate due to the asymmetry in
the velocity space of the bump-on-tail test case.

(a) SW formulation

(b) SW square-root formulation

Figure 11: The evolution of the electron distribution function $f^{e}(x,v,t)$
for the bump-on-tail instability with $N_{v}=100$. The distribution function
becomes negative in the (a) SW formulation in the majority of phase space
whereas the (b) SW square-root formulation is positivity preserving by
construction.

(a) Bump-on-tail instability filamentation

(b) Bump-on-tail instability energy exchange

Figure 12: Subfigure (a) shows a cross-section of the bump-on-tail electron
particle distribution function $f^{e}(x=11\pi,v,t=20)$ via the two
formulations: SW and SW square-root (with $N_{v}=100$). The filamentation is
relatively comparable in the SW and SW square-root formulations. Subfigure (b)
shows the evolution of the transferred kinetic and potential energy for the
bump-on-tail instability with $N_{v}=100$. The SW and SW square-root
macroscopic quantities–kinetic and potential energy– are nearly identical.

(a) SW formulation $N_{v}=101$ (odd)

(b) SW formulation $N_{v}=100$ (even)

(c) SW square-root formulation $N_{v}=101$ (odd)

(d) SW square-root formulation $N_{v}=100$ (even)

Figure 13: Same as Figure 3 for the bump-on-tail instability. Due to the
asymmetry in the velocity space in the bump-on-tail instability, the momentum
drift in both formulations is more pronounced than in the previous symmetric
test cases.

##  6 Conclusions

We proposed the anti-symmetric discretization of the one-dimensional Vlasov-
Poisson equations, which is based on the symmetrically-weighted Hermite
spectral expansion in velocity, centered finite differencing in space, and an
implicit Runge-Kutta temporal integrator. We showed that the method preserves
the anti-symmetric structure of the advection operator (or equivalently the
canonical Poisson bracket) of the Vlasov equation, which results in a stable
method. The anti-symmetric discretization is applied to two formulations: SW
and SW square-root formulation. The former discretizes the standard Vlasov-
Poisson equations and the latter discretizes a continuous square-root variable
transformation of the governing equations. The SW square-root formulation is
positivity preserving by construction. We derived the conservation properties
of the two formulations. The SW formulation conserves particle number and
energy if the number of velocity spectral terms is odd and momentum if the
number of spectral terms is even. Moreover, the conservation of momentum and
energy only holds in the SW expansion if the velocity shifting parameter is
set to zero. Alternatively, the SW square-root formulation conserves particle
number independently of the number of spectral terms and velocity
parameterization. Overall, this study introduces the construction of
structure-preserving and stable Vlasov-Poisson solvers.

A future goal of ours is to study techniques to mitigate filamentation.
Filamentation is a common issue in Eulerian solvers due to their limited
velocity resolution. In particular, in the case of a spectral expansion in
velocity space, filamentation can cause discontinuities in the higher-order
expansion coefficients, resulting in numerical instabilities. Filamentation
can also cause recurrence effects, where the solution is artificially periodic
in time [41]. Typical filamentation strategies include (1) artificial
collisions, (2) filtering, and (3) flux-limiter techniques. However, applying
such standard strategies to the SW and SW square-root formulations will break
conservation laws and the method’s anti-symmetric structure since its
invariants depend on all the coefficients in the expansion. We plan to study
such existing strategies and their influence on the SW and SW square-root
formulations in the future.

##  Appendix A Anti-symmetric Representation via a Fourier Spectral Expansion
in Space

We derive the Hermite-Fourier spectral discretization of the SW and SW square-
root formulations and show that a Fourier discretization in space preserves
the anti-symmetric structure of the advection operator in the Vlasov equation
(21). We begin by introducing the Fourier basis functions, such that

$\eta_{k}(x)\coloneqq\exp\left(\frac{2\pi ikx}{\ell}\right),\qquad\forall
k\in\mathbb{N}_{0}.$ (53)

Since the spatial domain $\Omega_{x}\coloneqq[0,\ell]$ is bounded, the Fourier
basis functions $\\{\eta_{k}(x)\\}_{k\in\mathbb{N}_{0}}$ form an orthogonal
basis with the following orthogonality relation

$\frac{1}{\ell}\int_{0}^{\ell}\eta_{k}(x)\eta_{m}(x)\mathrm{d}x=\delta_{k+m,0},\qquad\forall
k,m\in\mathbb{Z}.$ (54)

We approximate the Hermite expansion coefficient $C_{n}^{s}\in\mathcal{W}$ in
Eqns. (10) and (14) and the electric field $E\in\mathcal{W}$ in Eqns. (1)-(2)
via a spatial Fourier expansion, i.e.

$C_{n}^{s}(x,t)\approx
C_{n}^{s,N_{x}}(x,t)=\sum_{k=-N_{x}}^{N_{x}}C_{n,k}^{s}(t)\eta_{k}(x)\qquad\text{and}\qquad
E(x,t)\approx E^{N_{x}}(x,t)=\sum_{k=-N_{x}}^{N_{x}}E_{k}(t)\eta_{k}(x),$

where the Fourier basis is defined in Eq. (53). The product
$E^{N_{x}}(x,t)C_{n}^{s,N_{x}}(x,t)$ can be described as a Fourier
convolution, such that

$E^{N_{x}}(x,t)C_{n}^{s,N_{x}}(x,t)=\sum_{k=-2N_{x}}^{2N_{x}}\eta_{k}(x)[\mathbf{E}(t)\ast\mathbf{C}_{n}^{s}(t)][k],$

where

$\displaystyle\mathbf{E}(t)$
$\displaystyle=[E_{-N_{x}}(t),E_{-N_{x}+1}(t),\ldots,E_{0}(t),\ldots,E_{N_{x}-1}(t),E_{N_{x}}(t)]^{\top}\in\mathbb{C}^{2N_{x}+1},$
$\displaystyle\mathbf{C}_{n}(t)$
$\displaystyle=[C_{n,-N_{x}}^{s}(t),C_{n,-N_{x}+1}^{s}(t),\ldots,C_{n,0}^{s}(t),\ldots,C_{n,N_{x}-1}^{s}(t),C^{s}_{n,N_{x}}(t)]^{\top}\in\mathbb{C}^{2N_{x}+1},$

and $\ast$ denotes the following convolution

$[\mathbf{E}(t)\ast\mathbf{C}_{n}^{s}(t)][k]=\sum_{j=-N_{x}}^{N_{x}}E_{k-j}(t)C_{n,j}^{s}(t)$

Multiplying the Vlasov equation (11) by $\eta_{-m}(x)$, integrating over the
spatial coordinate $x$, and using the orthogonality relation of the Fourier
basis functions in Eq. (54), we get

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{C}_{n}^{s}(t)$
$\displaystyle=-\mathbf{D}\left(\alpha^{s}\sqrt{\frac{n}{{2}}}\mathbf{C}_{n-1}^{s}(t)+\alpha^{s}\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)+u^{s}\mathbf{C}_{n}^{s}(t)\right)$
(55)
$\displaystyle+\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}_{\mathrm{conv}}(t)\left(\sqrt{\frac{n}{2}}\mathbf{C}_{n-1}^{s}(t)-\sqrt{\frac{n+1}{2}}\mathbf{C}_{n+1}^{s}(t)\right),$

where

$\mathbf{D}=\frac{2\pi
i}{\ell}\text{diag}([-k,-k+1,\ldots,0,\ldots,k-1,k]\in\mathbb{C}^{(2N_{x}+1)\times(2N_{x}+1)}$
(56)

and

$\mathbf{E}_{\mathrm{conv}}(t)=\left[\begin{array}[]{cccccc}E_{0}(t)&E_{-1}(t)&\ldots&E_{-N_{x}}(t)&0&\ldots\\\
E_{1}(t)&E_{0}(t)&E_{-1}(t)&\ldots&E_{-N_{x}}(t)&0\\\
&\ddots&\ddots&\ddots&\ddots&\\\
E_{N_{x}}(t)&E_{N_{x}-1}(t)&\ldots&E_{0}(t)&\ldots&E_{-N_{x}}(t)\\\
&\ddots&\ddots&\ddots&\ddots&\\\
0&E_{N_{x}}(t)&E_{N_{x}-1}(t)&\ldots&E_{0}(t)&E_{-1}(t)\\\
\ldots&0&E_{N_{x}}(t)&E_{N_{x}-1}(t)&\ldots&E_{0}(t)\\\
\end{array}\right]\in\mathbb{C}^{(2N_{x}+1)\times(2N_{x}+1)}.$ (57)

Since the electric field $E\in\mathcal{W}$ is a real-valued function, then
$(\mathbb{C}^{2N_{x}+1},\mathbf{E}_{\mathrm{conv}}(t))=(\mathbb{C}^{2N_{x}+1},\mathbf{E}_{\mathrm{conv}}(t)^{*})$,
defined in Eq. (57), is a complex symmetric matrix. Similarly, the derivative
diagonal matrix
$(\mathbb{C}^{2N_{x}+1},\mathbf{D})=(\mathbb{C}^{2N_{x}+1},-\mathbf{D}^{*})$,
defined in Eq. (56), is a complex anti-symmetric matrix (also known as skew-
Hermitian).

The Fourier spectral discretization of the SW Poisson equation (12) is given
by

$E_{k}(t)=\frac{\ell}{2\pi
ik}\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}\mathcal{I}_{n}C_{n,k}^{s}(t),$

and the SW square-root Poisson equation (17) reads as

$E_{k}(t)=\frac{\ell}{2\pi
ik}\sum_{s}q^{s}\alpha^{s}\sum_{n=0}^{N_{v}-1}[\mathbf{C}_{n}^{s}(t)\ast\mathbf{C}_{n}^{s}(t)][k],$

with $E_{0}(t)=0$ from the uniqueness condition, i.e.
$\int_{0}^{\ell}E(x,t)\mathrm{d}x=0$.

######  Theorem 10.

A Fourier spatial spectral expansion in Eq. (23) preserves the anti-symmetric
structure of the advection operator $(\mathcal{W}^{N_{v}},A^{s}_{v}(x,t))$ in
Eq. (23), and thus preserves the anti-symmetric structure of the advection
operator $(\mathcal{V},A^{s}(x,v,t))$ in Eq. (20).

Proof: The semi-discrete Vlasov equation (55) can be written in vector form as

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Psi}^{s}(t)=\mathbf{A}^{s}(t)\mathbf{\Psi}^{s}(t)$

where

$\mathbf{\Psi}^{s}(t)=\begin{bmatrix}\mathbf{C}^{s}_{0}(t)\\\
\mathbf{C}^{s}_{1}(t)\\\ \vdots\\\
\mathbf{C}^{s}_{N_{v}-1}(t)\end{bmatrix}\in\mathbb{C}^{N_{v}(2N_{x}+1)}$

and

$\mathbf{A}^{s}(t)=\begin{bmatrix}-u^{s}\mathbf{D}&-\sqrt{\frac{1}{2}}\mathbf{G}_{+}^{s}(t)&0&\ldots&0\\\
\sqrt{\frac{1}{2}}\mathbf{G}_{-}^{s}(t)&-u^{s}\mathbf{D}&-\sqrt{\frac{2}{2}}\mathbf{G}_{+}^{s}(t)&\ldots&0\\\
&&&&\\\ &\ddots&\ddots&\ddots&\\\ &&&&\\\
0&0&\sqrt{\frac{N_{v}-2}{2}}\mathbf{G}_{-}^{s}(t)&-u^{s}\mathbf{D}&-\sqrt{\frac{N_{v}-1}{2}}\mathbf{G}_{+}^{s}(t)\\\
0&0&0&\sqrt{\frac{N_{v}-1}{2}}\mathbf{G}_{-}^{s}(t)&-u^{s}\mathbf{D}\end{bmatrix}\in\mathbb{C}^{N_{v}(2N_{x}+1)\times
N_{v}(2N_{x}+1)},$

such that

$\mathbf{G}_{+}^{s}(t)=\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}_{\mathrm{conv}}(t)+\alpha^{s}\mathbf{D}\qquad\mathrm{and}\qquad\mathbf{G}_{-}^{s}(t)=\frac{q^{s}}{m^{s}\alpha^{s}}\mathbf{E}_{\mathrm{conv}}(t)-\alpha^{s}\mathbf{D}.$

Since
$(\mathbb{C}^{2N_{x}+1},\mathbf{D})=(\mathbb{C}^{2N_{x}+1},-\mathbf{D}^{*})$
and
$(\mathbb{C}^{2N_{x}+1},\mathbf{E}_{\mathrm{conv}}(t))=(\mathbb{C}^{2N_{x}+1},\mathbf{E}_{\mathrm{conv}}(t)^{*})$,
the discretized advection operator is anti-symmetric where
$(\mathbb{C}^{N_{v}(2N_{x}+1)},\mathbf{A}^{s}(t))=(\mathbb{C}^{N_{v}(2N_{x}+1)},-\mathbf{A}^{s}(t)^{*})$.
$\Box$

## Code Availability

The public repository https://github.com/opaliss/Antisymmetric-Vlasov-Poisson
contains a collection of Jupyter notebooks in Python 3.9 with the code used in
this study.

## Acknowledgement

O.K. and G.L.D. are grateful for the useful and informative discussions with
Cecilia Pagliantini, Gianmarco Manzini, and Daniel Livescu. O.I. was partially
supported by the Los Alamos National Laboratory Vela Fellowship. O.I. and B.K.
were partially supported by the National Science Foundation under Award
2028125 for “SWQU: Composable Next Generation Software Framework for Space
Weather Data Assimilation and Uncertainty Quantification”. O.K. and G.L.D.
were supported by the Laboratory Directed Research and Development Program of
Los Alamos National Laboratory under project number 20220104DR. Los Alamos
National Laboratory is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. Department of Energy
(Contract No. 89233218CNA000001). F.D.H. was supported by the U.S. Department
of Energy, Office of Science, Office of Fusion Energy Sciences, Theory
Program, under Award No. DE-FG02-95ER54309.

## References

  * [1] M. Abramowitz and I. A. Stegun.  Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables.  Applied mathematics series. Dover Publications, 1965.
  * [2] A. Arakawa.  Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I.  Journal of Computational Physics, 1(1):119–143, 1966.
  * [3] A. H. Baker, E. R. Jessup, and T. Manteuffel.  A Technique for Accelerating the Convergence of Restarted GMRES.  SIAM Journal on Matrix Analysis and Applications, 26(4):962–984, 2005.
  * [4] M. Bessemoulin-Chatard and F. Filbet.  On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system.  Journal of Computational Physics, 451:110881, Feb 2022.
  * [5] P. Bogacki and L. Shampine.  A 3(2) pair of Runge-Kutta formulas.  Applied Mathematics Letters, 2(4):321–325, 1989.
  * [6] E. Camporeale, G. Delzanno, B. Bergen, and J. Moulton.  On the velocity space discretization for the Vlasov-Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cell methods.  Computer Physics Communications, 198:47–58, 2016.
  * [7] J. Canosa.  Numerical solution of Landau’s dispersion equation.  Journal of Computational Physics, 13(1):158–160, 1973.
  * [8] M. Carrié and B. Shadwick.  An unconditionally stable, time-implicit algorithm for solving the one-dimensional Vlasov-Poisson system.  Journal of Plasma Physics, 88(2):905880201, 2022.
  * [9] Y. Cheng, A. J. Christlieb, and X. Zhong.  Energy-conserving discontinuous Galerkin methods for the Vlasov-Ampére system.  Journal of Computational Physics, 256:630–655, 2014.
  * [10] G. Delzanno.  Multi-dimensional, fully-implicit, spectral method for the Vlasov-Maxwell equations with exact conservation laws in discrete form.  Journal of Computational Physics, 301:338–356, 2015.
  * [11] R. E. Denton and M. Kotschenreuther.  $\delta f$ Algorithm.  Journal of Computational Physics, 119(2):283–294, 1995.
  * [12] J. Dongarra and A. Geist.  Report on the Oak Ridge National Laboratory’s Frontier System.  Technical Report ICL-UT-22-05, Oak Ridge National Laboratory, 2022-05 2022.
  * [13] F. Filbet.  Convergence of a Finite Volume Scheme for the Vlasov–Poisson System.  SIAM Journal on Numerical Analysis, 39(4):1146–1169, 2001.
  * [14] F. Filbet and T. Xiong.  Conservative Discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System.  Communications on Applied Mathematics and Computation, 4:1–26, 09 2020.
  * [15] S. Gary.  Theory of Space Plasma Microinstabilities.  Cambridge Atmospheric and Space Science Series. Cambridge University Press, 1993.
  * [16] D. Gottlieb and S. Orszag.  Numerical Analysis of Spectral Methods: Theory and Applications.  Society for Industrial and Applied Mathematics, 1977.
  * [17] H. Grad.  On the kinetic theory of rarefied gases.  Communications on Pure and Applied Mathematics, 2(4):331–407, 1949.
  * [18] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, second edition, 2006.
  * [19] F. D. Halpern.  Anti-symmetric representation of the extended magnetohydrodynamic equations.  Physics of Plasmas, 27(4):042303, 04 2020.
  * [20] F. D. Halpern, I. Sfiligoi, M. Kostuk, R. Stefan, and R. E. Waltz.  Simulations of plasmas and fluids using anti-symmetric models.  Journal of Computational Physics, 445:110631, 2021.
  * [21] F. D. Halpern and R. E. Waltz.  Anti-symmetric plasma moment equations with conservative discrete counterparts.  Physics of Plasmas, 25(6), 06 2018.
  * [22] R. Heath, I. Gamba, P. Morrison, and C. Michler.  A discontinuous Galerkin method for the Vlasov-Poisson system.  Journal of Computational Physics, 231(4):1140–1174, 2012.
  * [23] J. Holloway.  On Numerical Methods for Hamiltonian PDEs and a Collocation Method for the Vlasov-Maxwell Equations.  Journal of Computational Physics, 129(1):121–133, 1996.
  * [24] J. Holloway.  Spectral velocity discretizations for the Vlasov-Maxwell equations.  Transport Theory and Statistical Physics, 25:1–32, 01 1996.
  * [25] A. Klimas and W. Farrell.  A Splitting Algorithm for Vlasov Simulation with Filamentation Filtration.  Journal of Computational Physics, 110(1):150–163, 1994.
  * [26] D. Knoll and D. Keyes.  Jacobian-free Newton-Krylov methods: a survey of approaches and applications.  Journal of Computational Physics, 193(2):357–397, 2004.
  * [27] K. Kormann and A. Yurova.  A generalized Fourier-Hermite method for the Vlasov-Poisson system.  BIT Numerical Mathematics, 61, 04 2021.
  * [28] O. Koshkarov, G. Manzini, G. Delzanno, C. Pagliantini, and V. Roytershteyn.  The multi-dimensional Hermite-discontinuous Galerkin method for the Vlasov-Maxwell equations.  Computer Physics Communications, 264:107866, 2021.
  * [29] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker.  GEMPIC: geometric electromagnetic particle-in-cell methods.  Journal of Plasma Physics, 83(4), July 2017.
  * [30] M. Kraus, E. Tassi, and D. Grasso.  Variational integrators for reduced magnetohydrodynamics.  Journal of Computational Physics, 321:435–458, 2016.
  * [31] L. D. Landau.  On the vibrations of the electronic plasma.  Journal of Physics, 10(1):25–34, 1946.
  * [32] G. Manzini, G. Delzanno, J. Vencels, and S. Markidis.  A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system.  Journal of Computational Physics, 317:82–107, July 2016.
  * [33] G. Manzini, D. Funaro, and G. L. Delzanno.  Convergence of Spectral Discretizations of the Vlasov–Poisson System.  SIAM Journal on Numerical Analysis, 55(5):2312–2335, Jan 2017.
  * [34] P. Morrison.  Hamiltonian field description of the one-dimensional Poisson-Vlasov equations.  Princeton Plasma Physics Laboratory, Technical Report PPPL-1788, 07 1981.
  * [35] P. J. Morrison.  Hamiltonian description of the ideal fluid.  Rev. Mod. Phys., 70:467–521, Apr 1998.
  * [36] A. Myers, P. Colella, and B. V. Straalen.  A 4th-Order Particle-in-Cell Method with Phase-Space Remapping for the Vlasov–Poisson Equation.  SIAM Journal on Scientific Computing, 39(3):B467–B485, 2017.
  * [37] C. Pagliantini, G. L. Delzanno, and S. Markidis.  Physics-based adaptivity of a spectral method for the Vlasov-Poisson equations based on the asymmetrically-weighted Hermite expansion in velocity space.  Journal of Computational Physics, 488:112252, Sep 2023.
  * [38] C. Pagliantini, G. Manzini, O. Koshkarov, G. L. Delzanno, and V. Roytershteyn.  Energy-conserving explicit and implicit time integration methods for the multi-dimensional Hermite-DG discretization of the Vlasov-Maxwell equations.  Computer Physics Communications, 284:108604, 2023.
  * [39] M. Palmroth, U. Ganse, Y. Pfau-Kempf, M. Battarbee, L. Turc, T. Brito, M. Grandin, S. Hoilijoki, A. Sandroos, and S. Alfthan.  Vlasov methods in space physics and astrophysics.  Living Reviews in Computational Astrophysics, 4, 08 2018.
  * [40] J. T. Parker and P. J. Dellar.  Fourier-Hermite spectral representation for the Vlasov-Poisson system in the weakly collisional limit.  Journal of Plasma Physics, 81(2):305810203, 2015.
  * [41] O. Pezzi, E. Camporeale, and F. Valentini.  Collisional effects on the numerical recurrence in Vlasov-Poisson simulations.  Physics of Plasmas, 23, 01 2016.
  * [42] M. Renardy and R. Rogers.  An Introduction to Partial Differential Equations.  Texts in Applied Mathematics. Springer New York, 2004.
  * [43] P. Roe.  Approximate Riemann solvers, parameter vectors, and difference schemes.  Journal of Computational Physics, 43(2):357–372, 1981.
  * [44] V. Roytershteyn and G. L. Delzanno.  Spectral Approach to Plasma Kinetic Simulations Based on Hermite Decomposition in the Velocity Space.  Frontiers in Astronomy and Space Sciences, 5, 2018.
  * [45] Y. Saad and M. H. Schultz.  GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.  SIAM Journal on Scientific and Statistical Computing, 7(3):856–869, 1986.
  * [46] J. W. Schumer and J. P. Holloway.  Vlasov Simulations Using Velocity-Scaled Hermite Representations.  Journal of Computational Physics, 144(2):626–661, 1998.
  * [47] C. Scovel and A. Weinstein.  Finite dimensional Lie-Poisson approximations to Vlasov-Poisson equations.  Communications on Pure and Applied Mathematics, 47(5):683–709, 1994.
  * [48] T. Shiroto, N. Ohnishi, and Y. Sentoku.  Quadratic conservative scheme for relativistic Vlasov-Maxwell system.  Journal of Computational Physics, 379:32–50, 2019.
  * [49] M. Shoucri and G. Knorr.  Numerical integration of the Vlasov equation.  Journal of Computational Physics, 14(1):84–92, 1974.
  * [50] T. Tang.  The Hermite Spectral Method for Gaussian-Type Functions.  SIAM Journal on Scientific Computing, 14(3):594–606, 1993.
  * [51] J. Vencels, G. Delzanno, A. Johnson, I. Peng, E. Laure, and S. Markidis.  Spectral Solver for Multi-Scale Plasma Physics Simulations with Dynamically Adaptive Number of Moments.  Procedia Computer Science, 51:1148–1157, 12 2015.
  * [52] J. Vencels, G. Delzanno, G. Manzini, S. Markidis, I. Peng, and V. Roytershteyn.  SpectralPlasmaSolver: a Spectral Code for Multiscale Simulations of Collisionless, Magnetized Plasmas.  Journal of Physics: Conference Series, 719:012022, 05 2016.
  * [53] J. P. Verboncoeur.  Particle simulation of plasmas: review and advances.  Plasma Physics and Controlled Fusion, 47(5A):A231, Apr 2005.
  * [54] G. Vogman, U. Shumlak, and P. Colella.  Conservative fourth-order finite-volume Vlasov-Poisson solver for axisymmetric plasmas in cylindrical $(r,v_{r},v_{\theta})$ phase space coordinates.  Journal of Computational Physics, 373:877–899, 2018.
  * [55] B. Wang, G. H. Miller, and P. Colella.  A Particle-in-cell Method with Adaptive Phase-space Remapping for Kinetic Plasmas.  SIAM Journal on Scientific Computing, 33(6):3509–3537, 2011.
  * [56] J. Xiao, H. Qin, and J. Liu.  Structure-preserving geometric particle-in-cell methods for Vlasov-Maxwell systems.  Plasma Science and Technology, 20(11):110501, Sep 2018.
  * [57] S. Zaki, T. Boyd, and L. Gardner.  A finite element code for the simulation of one-dimensional Vlasov plasmas. II. Applications.  Journal of Computational Physics, 79(1):200–208, 1988.
  * [58] S. Zaki, L. Gardner, and T. Boyd.  A finite element code for the simulation of one-dimensional Vlasov plasmas. I. Theory.  Journal of Computational Physics, 79(1):184–199, 1988.