# Kinetic theory of particle-in-cell simulation plasma and the ensemble
averaging technique

Michaël<EMAIL_ADDRESS>Marvel Fusion Gmbh,
Blumenstrasse 28, Munich, Germany  Romain Codur  Department of Physics and
Astronomy, University of California Los Angeles, Los Angeles, USA  Frank Tsung
Department of Physics and Astronomy, University of California Los Angeles, Los
Angeles, USA  Viktor K Decyk  Department of Physics and Astronomy, University
of California Los Angeles, Los Angeles, USA  Warren B Mori  Department of
Physics and Astronomy, University of California Los Angeles, Los Angeles, USA
Luis O Silva  GOLP/Instituto de Plasma e Fusão Nuclear, Instituto Superior
Técnico, Universidade de Lisboa, Lisbon, Portugal

(August 2022)

###### Abstract

We derive the kinetic theory of fluctuations in physically and numerically
stable particle-in-cell (PIC) simulations of electrostatic plasmas. The
starting point is the single-time correlations at the simulation start between
the statistical fluctuations of weighted densities of macroparticle centers in
the plasma particle phase-space. The fluctuations are associated with
different initial conditions, typically due to the random initial conditions
(in velocity space) of the macroparticles/simulation plasma, assigned
according to their initial distribution of probability. The single-time
correlations at all time steps and in each spatial grid cell are then
determined from the Laplace-Fourier transforms of the discretized
Klimontovich-like equation for the macroparticles and Maxwell’s equations for
the fields as computed by modern PIC codes. We recover the expressions for the
electrostatic field and the plasma particle density fluctuation
autocorrelations spectra as well as the kinetic equations describing the
average evolution of PIC-simulated plasma particles, first derived by
[Langdon, 1970b] using a macroparticle test approach perturbing a discretized
Vlasovian plasma and then averaging the obtained physical quantity over the
initial macroparticle velocity distribution. We generalize and extend these
results to the modern algorithms in PIC codes and using arbitrary
macroparticle weights. Analytical estimates of statistical fluctuations
single-time correlation amplitudes are derived as a function of the plasma
simulation parameters, using the central limit theorem in the limit of a large
number of macroparticles per cell. The theory is then used to analyze the
ensemble averaging technique of PIC simulations where statistical averages are
performed over ensembles of PIC simulations, modeling the same plasma physics
problem but using different statistical realizations of the initial
distribution functions of the macroparticles. This method is illustrated with
linear Landau damping uncovering (out of what is usually considered numerical,
noise) the physical fluctuations driven by a single small amplitude
electrostatic wave perturbing a PIC simulation plasma in equilibrium.
keywords: kinetic theory, electrostatic plasma, particle-in-cell codes,
numerical noise, fluctuations correlations, ensemble averaging technique

###### Contents

  1. 1 Introduction
  2. 2 Electrostatic PIC simulations
    1. 2.1 Discretized particle and field equations in an electrostatic PIC code
    2. 2.2 Initialization of a fully ionized, infinite and homogeneous plasma in PIC simulations
    3. 2.3 Kinetic theory of PIC simulation electrostatic plasmas
  3. 3 Statistical fluctuations in electrostatic PIC simulations
    1. 3.1 Assumptions and notations
    2. 3.2 Estimate of statistical fluctuations in a PIC simulation electrostatic plasma
    3. 3.3 Kinetic equation of a PIC simulation electrostatic plasma
  4. 4 Discussion
    1. 4.1 Neglect of numerical consistency errors compared to statistical fluctuations
    2. 4.2 Stationarity of the plasma distribution function in PIC simulations
    3. 4.3 Evolution of the empirical entropy in electrostatic PIC simulation plasmas
  5. 5 The ensemble averaging technique
    1. 5.1 Reduction of statistical fluctuations via ensemble averaging
    2. 5.2 Kinetic equation for an ensemble average of electrostatic PIC simulations
    3. 5.3 Application to linear Landau damping
  6. 6 Conclusions

##  1 Introduction

A plasma is a set of charged particles consisting of electrons and ionized
atoms that behave collectively through their self-consistent long-distance
electromagnetic fields. In this paper, we consider (the, apparently, less
complex) non-degenerate, non-relativistic, and non-coupled collisionless
electrostatic plasmas without external forces. One possible approach to
describe the plasma considers all particle equations of motion under the
action of the plasma electrostatic fields, determined according to the self-
consistent Maxwell equations. This fully deterministic approach allows for
deducing the particle positions and velocities at every time (and thus the
temporal evolution of the plasma) if the full microscopic details of the
system are known at a given time. Differentiating the phase-space densities
associated with this detailed description with respect to time, the
Klimontovich equation for each plasma species [Klimontovich, 1958] is
obtained. The Klimontovich equation describes the evolution in time of all $N$
plasma particles in the 6-dimensional phase-space under the action of self-
consistent microscopic electrostatic fields. The detailed knowledge of the
positions and velocities of the particles is needed to evaluate the
microscopic electrostatic fields, which renders it considerably difficult due
to the challenge of computing simultaneously the particle dynamics and the
microscopic self-consistent fields, even resorting to today’s computer
technology. The PIC method relaxes this constraint by considering particles
with a finite spatial size, also called macroparticles [Hockney, 1965] [Yu et
al., 1965] [Hockney, 1966] [Burger et al., 1965] [Buneman, 1967] [Birdsall and
Fuss, 1969] [Dawson, 1983]. For dimensionless simulations mapped to a real
physical system (or for simulations that use real physical units), each
macroparticle usually represents several real particles – the relation between
real particles and their macroparticle representation depends on the
compromise between the plasma density of the system under study and the
available computer technology (for instance, restrictions on computer time and
memory / represent the plasma particles). In the PIC method, macroparticle
positions are interpolated onto a spatial grid used to compute the self-
consistent discretized Maxwell equations. Therefore, for simulations of a
plasma homogeneously distributed in the whole simulation box, the PIC
algorithm is much more efficient than the direct Klimontovich model previously
outlined, with a computational cost $C\sim\mathrm{O}\left(N\right)$ compared
to the computation of all binary Coulomb interactions between macroparticles
that would need a computational cost
$C^{\prime}\sim\mathrm{O}\left(N^{2}\right)$. In terms of the physics being
captured, the PIC method smooths the fields and thus underestimates close-
encounter Coulomb collisions [Okuda and Birdsall, 1970] [Langdon, 1970a].

Moreover, from classical thermodynamic theory, one can argue that the exact
positions and velocities of all particles are never known/measured
simultaneously. Instead, and for large N-body systems, knowledge of the system
can be obtained via quantities averaged over time, infinitesimally small
mesoscopic spatial volumes, or over different measurements performed under the
same conditions. This naturally motivates a probabilistic approach, based on
the temporal evolution of the probability distribution functions. Liouville
theorem describes the temporal evolution of the full $N$-body probability
distribution function of all particle position and velocity random variables
[Liouville, 1838]. As before, the information contained in the $N$-body
distribution function is too large to allow for the computation of the
Liouville equation coupled with the self-consistent microscopic Maxwell
equations even with today’s computer technology. Simplifications are also
possible in this approach: smaller order $s$-body distribution functions can
be obtained by integrating the Liouville equation over the position and
velocity random variables of all other particles $s+1,\,...,N$. In the
resulting equation, describing the temporal evolution of the $s$-body
distribution functions in the 6$s$-dimensional phase-space, the presence of
the ”collision” integral222We observe this term is used, even if, sometimes,
the ”collision” integral does not strictly describe close-encounter Coulomb
collisions, also poses challenges, and some additional assumptions are needed
to close the problem. The first order equation ($s=1$) in the BBGKY hierarchy
of equations [Bogoliubov, 1946, Born and Green, 1946, Kirkwood, 1946, Yvon,
1935], can be truncated to obtain the Vlasov equation by assuming an infinite
number of electrons in a Debye sphere, thus neglecting 2-body correlations
between 1-body distribution functions [Vlasov, 1938]. However, the number of
particles is necessarily finite in plasma, and particles are always correlated
due to their long-range electrostatic interactions. The Vlasov equation
remains valid only on time scales much smaller than this minimum interaction
time. On larger time scales, the Bogoliubov hypothesis [de Boer and Uhlenbeck,
1962] can be used to simplify the second order equation in the BBGKY hierarchy
and estimate the 2-body correlations between plasma particles. The Bogoliubov
hypothesis assumes the 2-body correlations relax on a time scale smaller than
the time scale on which the 1-body distribution function relaxes. Further
assumptions are still needed to get an analytical expression of correlations
e.g considering only binary electrostatic Coulomb interaction between plasma
particles and neglecting their 3-body correlations, it is possible to derive
the Vlasov-Boltzmann equation or the Vlasov-Lenard-Balescu equation depending
on the underlying assumptions [Lenard, 1960] [Balescu, 1960] [Guernsey, 1960].
Both the Vlasov-Boltzmann equation or the Vlasov-Lenard-Balescu equation
reduce to the Vlasov-Fokker-Planck-Landau equation when considering only non-
screened and small angle Coulomb collisions between plasma particles [Landau,
1937]. The numerical computation of these kinetic equations onto the resulting
7-dimensional phase-space-time grid is still challenging with modern computer
technology and further approximations are usually used to reduce the number of
variables [Tzoufras et al., 2011, Touati et al., 2014, Joglekar et al., 2014].

The theory of fluctuations in collisionless plasmas is related to (and
connects) both the deterministic (Klimontovich-based) and probabilistic
(Liouville-based) approaches. It is based on the study of single-time
correlations between exact microscopic plasma quantities fluctuations,
measuring the deviations of one statistical realization of the plasma
consisting of a discrete number of self-interacting particles around the
expected average, which is described by the continuous 1-body distribution
function. [Rostoker, 1961, Klimontovich and Silin, 1962, Dupree, 1963] have
shown that the study of single-time correlations between fluctuations also
leads to the Vlasov-Lenard-Balescu equation describing the 1-body distribution
function in the case of homogeneous and stationary collisionless plasmas. This
bears a striking resemblance with the fundamentals of PIC simulations: in PIC
simulations, the initial macroparticle velocities are determined by random
sampling of the (desired) probability distribution function. Therefore, this
random sampling (via a random number generator) naturally produces statistical
fluctuations around the initial (ensemble averaged) distribution function. The
main motivation for this paper lies in this observation: we apply the theory
of fluctuations to the discretized equations associated with electrostatic PIC
code (taking into account both the discretized equations of motion and the
discretized field equations), and we explore this theory to determine the
statistical properties (and possible advantages) of ensemble averages over
sets of PIC simulations. The paper is organized as follows. In section 2, we
recapitulate the discretized equations of an electrostatic PIC code and we
reformulate them within the mathematical framework of the kinetic theory of
plasmas. In section 3, we derive the theory of statistical fluctuations in
electrostatic PIC simulations, for stationary, homogeneous, infinite, and
fully ionized plasmas. We discuss the assumptions and the key results of the
theory in section 4. In section 5, we discuss the ”ensemble averaging
technique” [Shanny et al., 1967] within the theoretical framework developed in
the previous sections, and we determine the conditions under which this
technique relaxes the discrete particle effects inherent to the use of a
finite and relatively small number of macroparticles resulting from the limits
of computer technology. In section 6, the main conclusions of the paper are
presented.

##  2 Electrostatic PIC simulations

###  2.1 Discretized particle and field equations in an electrostatic PIC code

In a PIC simulation, the phase-space of each plasma species $a$ is sampled
according to $N_{a}$ macroparticles of mass $M_{a}$ and charge $Q_{a}$. When
the equations computed by a PIC code are renormalized to physical units, each
macroparticle represents $\delta N_{a}=M_{a}/m_{a}=Q_{a}/q_{a}$ real particles
of mass $m_{a}$ and electrical charge $q_{a}$. $\delta N_{a}$ is therefore a
number, a lineal density or an areal density of real particles from species
$a$ in 3D, 2D or 1D PIC simulations, respectively. Macroparticle trajectories
$\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,\mathbf{v}_{a,\ell}\left(t_{n}\right)\right)$
in phase-space are computed at each time step
$t_{n}=\left(n-1\right)\Delta_{t}$ according to their discretized equation of
motion

$\forall a,\,\forall\ell\in\left[1,\,N_{a}\right],\,\forall
n,\,\displaystyle\left\\{\begin{array}[]{lllr}{\displaystyle\left.\displaystyle\frac{d\mathbf{r}_{a,\ell}}{dt}\right|}^{n}&=&\mathbf{v}_{a,\ell}\left(t_{n}\right)&,\mathbf{r}_{a,\ell}\left(t_{1}\right)=\mathbf{R}_{a,\ell}\cr{\displaystyle\left.\displaystyle\frac{d\mathbf{v}_{a,\ell}}{dt}\right|}^{n}&=&\displaystyle\frac{Q_{a}}{M_{a}}\mathbf{E}_{s}\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)&,\mathbf{v}_{a,\ell}\left(t_{1}\right)=\mathbf{V}_{a,\ell}\end{array}\right.$
(1)

where

$\mathbf{E}_{s}\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)=\displaystyle\sum_{i,j,k}\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)S\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right)-\mathbf{r}_{i,j,k}\right)\Delta_{x}\Delta_{y}\Delta_{z}$
(2)

is the electric field, interpolated at macroparticle center locations. The
electrostatic field is deduced according to the discretized Maxwell equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}Q_{a}\displaystyle\sum_{\ell=1}^{N_{a}}S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\ell}\left(t_{n}\right)\right)\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right.$ (3)

where macroparticle centers $\mathbf{r}_{a,\ell}\left(t_{n}\right)$ are
interpolated at the grid points
$\forall\left\\{i,\,j,\,k\right\\}\in\left[1,N_{x}\right]\times\left[1,N_{y}\right]\times\left[1,N_{z}\right],\,\mathbf{r}_{i,\,j,\,k}=\displaystyle\left(\left(i-1\right)\Delta_{x},\,\displaystyle\left(j-1\right)\Delta_{y},\,\displaystyle\left(k-1\right)\Delta_{z}\right)$.
Here, $L_{x}=N_{x}\Delta_{x}$, $L_{y}=N_{y}\Delta_{y}$ and
$L_{z}=N_{z}\Delta_{z}$ are the simulation box dimensions and $\Delta_{x}$,
$\Delta_{y}$ and $\Delta_{z}$ are the spatial grid spacings that are usually
constant. The index $k$ or indices $j$ and $k$ for 2D or 3D simulations,
respectively, should be ommitted in the equations but we will keep writing
them in all the following to be the more general as possible. A smoothed
interpolating function

$S\displaystyle\left(\mathbf{r}\right)=\displaystyle\left\\{\begin{array}[]{lll}S_{x}\displaystyle\left(x\right)&\,\mathrm{in}\,1\mathrm{D}\cr
S_{x}\displaystyle\left(x\right)S_{y}\displaystyle\left(y\right)&\,\mathrm{in}\,2\mathrm{D}\cr
S_{x}\displaystyle\left(x\right)S_{y}\displaystyle\left(y\right)S_{z}\displaystyle\left(z\right)&\,\mathrm{in}\,3\mathrm{D}\end{array}\right.\,\mathrm{with}\,\forall\xi\in[x,\,y,\,z],\,S_{\xi}(\xi)=\displaystyle\left[\Pi_{\Delta_{\xi}}^{(n)}*S_{0_{\xi}}\right]\displaystyle\left(\xi\right)$
(4)

is used where usually $\Pi_{\Delta}^{(n)}=\Pi_{\Delta}*...*\Pi_{\Delta}$ is
the $n$th-order B-spline, i.e. the centered square distribution of support
$\Delta$, $\Pi_{\Delta}\left(\xi\right)=1/\Delta$ if
$\left|\xi\right|\leq\Delta/2$ and $0$ else, convoluted $n$ times with itself.
$S_{0_{\xi}}$ are optional filters or smoothing functions that mitigate
spatial aliasing effects such as, for example, a Gaussian filter

$S_{0_{\xi}}\displaystyle\left(\xi\right)=\displaystyle\frac{1}{\sqrt{2\pi{a_{\xi}}^{2}}}\exp{\displaystyle\left[-\displaystyle\frac{\xi^{2}}{2{a_{\xi}}^{2}}\right]}$
(5)

for spectral Maxwell solvers or such as a compensated $N$th-order binomial
filter for finite-difference time-domain (FDTD) Maxwell solvers. Interpolating
functions and filters vary a lot across codes. For example, in PIC codes using
a FDTD scheme, the smoothing function is sometimes only applied to the
electrostatic field after computing it according to Maxwell equations (5). In
all the following, we will consider the same interpolating function in both
(2) and (3) so that the smoothed interpolating function (4) can be interpreted
physically as the spatial shape of macroparticles. The generalization of our
results to the case where different interpolating functions are used for (2)
and (3) is straightforward. In general, the interpolating and smoothing
functions have the following properties

$\displaystyle\left\\{\begin{array}[]{ll}\forall\mathbf{r},&\displaystyle\sum_{i,j,k}\Pi_{\mathbf{\boldsymbol{\Delta}}}^{(n)}\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}\right)\Delta_{x}\Delta_{y}\Delta_{z}=1,\\\
\cr\forall\mathbf{r},&\displaystyle\int_{\mathbb{R}^{3}}\Pi_{\mathbf{\boldsymbol{\Delta}}}^{(n)}\displaystyle\left(\mathbf{r}\right)d^{3}\mathbf{r}=1\\\
\cr\mathrm{and}&\displaystyle\left\\{\begin{array}[]{lll}\displaystyle\int_{\mathbb{R}^{3}}S_{0}\displaystyle\left(\mathbf{r}\right)d^{3}\mathbf{r}=1&\mathrm{for\,spectral\,Maxwell\,solver\,or}\\\
\cr\displaystyle\sum_{i,j,k}S_{0}\displaystyle\left(\mathbf{r}_{i,j,k}\right)\Delta_{x}\Delta_{y}\Delta_{z}=1&\mathrm{for\,FDTD\,Maxwell\,solvers}\\\
\end{array}\right.\end{array}\right.$ (6)

Here, we have noted
$\mathbf{\boldsymbol{\Delta}}={\displaystyle\left(\Delta_{x},\,\Delta_{y},\,\Delta_{z}\right)}^{t}$
with the uppercase t meaning the vector transposition,
$\Pi_{\mathbf{\boldsymbol{\Delta}}}^{(n)}\displaystyle\left(\mathbf{r}\right)=\Pi_{\Delta_{x}}^{(n)}\displaystyle\left(x\right)\Pi_{\Delta_{y}}^{(n)}\displaystyle\left(y\right)\Pi_{\Delta_{z}}^{(n)}\displaystyle\left(z\right)$
and
$S_{0}\displaystyle\left(\mathbf{r}\right)=S_{0_{x}}\displaystyle\left(x\right)S_{0_{y}}\displaystyle\left(y\right)S_{0_{z}}\displaystyle\left(z\right)$
to simplify the notations. Let us finally stress here that even in the case
where the PIC codes computes Maxwell equations by using the second order FDTD
scheme proposed by [Yee, 1966] coupled with the charge conserving scheme
proposed by [Villasenor and Buneman, 1992] for first order B-spline
macroparticle shapes or by [Esirkepov, 2001] for higher orders, the resulting
numerical scheme can be written with the form (3) if the discretized Maxwell-
Poisson equation is verified at the simulation start.

###  2.2 Initialization of a fully ionized, infinite and homogeneous plasma in
PIC simulations

Let us consider the electrostatic PIC simulation of a fully ionized, infinite
and homogeneous plasma of atomic number $Z$. We note $\bar{n}_{i}$ the ion
density and $\bar{n}_{e}=Z\bar{n}_{i}$ the electron density. We use $L_{x}$,
$L_{y}$ and $L_{z}$-periodic boundary conditions so that

$\forall a\in\\{e,\,i\\},\,\bar{n}_{a}=\displaystyle\frac{N_{a}\delta
N_{a}}{L_{x}L_{y}L_{z}}$ (7)

(In 1D simulations, $\bar{n}_{a}=N_{a}\delta N_{a}/L_{x}$ and in 2D
simulations, $\bar{n}_{a}=N_{a}\delta N_{a}/L_{x}L_{y}$). We choose a number
$N_{i}$ and a number $N_{e}=ZN_{i}$ of macroions and macroelectrons to sample
the plasma particle phase-spaces such that $\delta N_{e}=\delta N_{i}=\delta
N$. At the simulation start, we use a random number generator to initialize
macroparticle velocities $\mathbf{V}_{a,\ell}$ according to the desired
probability distribution $F_{a0}\left(\mathbf{V}_{a}\right)$. For example, we
have for the non-drifting Maxwell-Boltzmann equilibrium probability
distribution of velocities that maximizes entropy:

$F_{a0}\left(\mathbf{V}_{a}\right)=\displaystyle\frac{1}{{\left(2\pi{v_{T_{a}}}^{2}\right)}^{3/2}}\exp{\displaystyle\left(-\displaystyle\frac{{\mathbf{V}_{a}}^{2}}{2{v_{Ta}}^{2}}\right)}$
(8)

where $\mathbf{V}_{a}$ is the velocity random variable of macroparticles from
species $a$, $v_{T_{a}}=\sqrt{k_{B}T_{a}/m_{a}}$ is the thermal velocity of
particles, $k_{B}$ is the Boltzmann constant and $T_{a}$ are the species
temperatures. However, in order to ensure homogeneous electrical charge
neutrality in the whole PIC simulation plasma, macroparticle locations at the
simulation start $\mathbf{R}_{a,\ell}$ are not randomly chosen according to a
uniform probability law. Instead, we deposit each macroparticle
$\ell\in\left[1,\,N_{a}\right]$ homogeneously inside all spatial grid cells.
If we note $N_{a_{x}}$, $N_{a_{y}}$ and $N_{a_{z}}$ the number of
macroparticles distributed along the $x$, $y$ and $z$-axis respectively such
that $N_{a}=N_{a_{x}}N_{a_{y}}N_{a_{z}}$, the macroparticle locations read
$\forall\left(\ell_{x},\,\ell_{y},\,\ell_{z}\right)\in\left[1,\,N_{a_{x}}\right]\times\left[1,\,N_{a_{y}}\right]\times\left[1,\,N_{a_{z}}\right],$

$\mathbf{R}_{a,\ell}=\displaystyle\left(\displaystyle\begin{array}[]{c}X_{a,\ell_{x}}\cr
Y_{a,\ell_{y}}\cr
Z_{a,\ell_{z}}\end{array}\right)=\displaystyle\left(\displaystyle\begin{array}[]{c}\left(\ell_{x}-\displaystyle\frac{1}{2}\right)\displaystyle\frac{L_{x}}{N_{a_{x}}}\cr\left(\ell_{y}-\displaystyle\frac{1}{2}\right)\displaystyle\frac{L_{y}}{N_{a_{y}}}\cr\left(\ell_{z}-\displaystyle\frac{1}{2}\right)\displaystyle\frac{L_{z}}{N_{a_{z}}}\end{array}\right)$
(9)

where
$\ell=\ell_{x}+\left(\ell_{y}-1\right)N_{a_{x}}+\left(\ell_{z}-1\right)N_{a_{x}}N_{a_{y}}$
and where $-1/2$ is used to avoid having half the box have more charge than
the other half, which gives rise to non-physical larger amplitude oscillations
due to initial charge imbalances. The numbers of macroparticles are related to
the number of macroparticles per cell according to

$N_{a,\mathrm{mpc}}=\displaystyle\frac{N_{a}}{N_{x}N_{y}N_{z}}=\displaystyle\frac{N_{a_{x}}}{N_{x}}\displaystyle\frac{N_{a_{y}}}{N_{y}}\displaystyle\frac{N_{a_{z}}}{N_{z}}$
(10)

(In 1D simulations, $N_{a,\mathrm{mpc}}=N_{a}/N_{x}=N_{a_{x}}/N_{x}$ and in 2D
simulations,
$N_{a,\mathrm{mpc}}=N_{a}/N_{x}N_{y}=N_{a_{x}}N_{a_{y}}/N_{x}N_{y}$). The
latter is a key parameter of PIC simulations as we are going to see in the
next sections.

###  2.3 Kinetic theory of PIC simulation electrostatic plasmas

One can reformulate equations (1) and (3) computed by an electrostatic PIC
code by defining the phase-space densities of macroparticle centers multiplied
by their weights

$f_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)=\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}\delta\left(\mathbf{r}-\mathbf{r}_{a,\ell}\left(t\right)\right)\delta\left(\mathbf{v}-\mathbf{v}_{a,\ell}\left(t\right)\right)$
(11)

where we have noted
$\delta\left(\mathbf{r}\right)=\delta\left(x\right)\delta\left(y\right)\delta\left(z\right)$
($\delta\left(\mathbf{r}\right)=\delta\left(x\right)\delta\left(y\right)$ in
2D and $\delta\left(\mathbf{r}\right)=\delta\left(x\right)$ in 1D) the Dirac
distribution. Firstly, by differentiating the latter with time and then
discretizing it in space and time, we obtain

${\displaystyle\left.\displaystyle\frac{\partial f_{a_{c}}}{\partial
t}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}f_{a_{c}}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)=0$
(12)

at the grid points $\mathbf{r}_{i,j,k}$ and time steps $t_{n}$ according to
macroparticle discretized equations of motion (1). Secondly, one may express
Maxwell equations (3) as a function of weighted macroparticle center phase-
space densities (11) according to

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\int_{\mathbb{R}^{3}}f_{a_{c}}\displaystyle\left(\mathbf{r}^{\prime},\,\mathbf{v},\,t_{n}\right)S\left(\mathbf{r}_{i,j,k}-\mathbf{r}^{\prime}\right)d^{3}\mathbf{r}^{\prime}d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (13)

Here, we recognize the phase-space densities of real particles

$\begin{array}[]{lll}f_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)&=&\displaystyle\int_{\mathbb{R}^{3}}f_{a_{c}}\displaystyle\left(\mathbf{r}^{\prime},\,\mathbf{v},\,t\right)S\left(\mathbf{r}-\mathbf{r}^{\prime}\right)d^{3}\mathbf{r}^{\prime}\cr&=&\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}S\displaystyle\left(\mathbf{r}-\mathbf{r}_{a,\ell}\left(t\right)\right)\delta\left(\mathbf{v}-\mathbf{v}_{a,\ell}\left(t\right)\right)\end{array}$
(14)

as it is sampled according to the PIC method. By differentiating the latter
with time in order to get the deterministic description of electrostatic PIC
simulations, we obtain the discretized kinetic equation

${\begin{array}[]{l}{\displaystyle\left.\displaystyle\frac{\partial
f_{a}}{\partial
t}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}f_{a}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)=C_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right).\end{array}}$
(15)

It is self-consistently coupled with the discretized Maxwell equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}f_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (16)

The source term

$C_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)=-\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}\displaystyle\left[\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)-\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\right]S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\,\ell}\displaystyle\left(t_{n}\right)\right)\delta\displaystyle\left(\mathbf{v}-\mathbf{v}_{a,\,\ell}\displaystyle\left(t_{n}\right)\right)\right)$
(17)

is due to the sum of forces
$(q_{a}/m_{a})\left[\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)-\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\right]$
acting on each macroparticle $a,\ell$ because of the internal tension force
which is constraining the macroparticle charge and inertia as a unit, when
different parts of the macroparticle feel a different electrostatic field
[Decyk, 1982]. According to its mathematical expression (17), these non-
physical forces decreases with decreasing spatial grid spacings and/or with
increasing order of interpolation and smoothing functions. One must always
check

$\begin{array}[]{lll}&&\displaystyle\sum_{i,j,k}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{m_{a}\mathbf{v}^{2}}{2}C_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)d^{3}\mathbf{v}\Delta_{x}\Delta_{y}\Delta_{z}\cr&=&q_{a}\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}\displaystyle\left[\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)-\displaystyle\sum_{i,j,k}\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\,\ell}\displaystyle\left(t_{n}\right)\right)\Delta_{x}\Delta_{y}\Delta_{z}\right]\cdot\mathbf{v}_{a,\,\ell}\displaystyle\left(t_{n}\right)\end{array}$
(18)

remains sufficiently small in case of electrostatic oscillations on the
macroparticle size space scale by checking this energy flow associated with
this internal tension force at each time step $t_{n}$ during the whole
simulation.

The random number generator used to initialize macroparticle velocities
$\mathbf{V}_{a,\ell}$ introduces statistical fluctuations

$\delta
f_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)=f_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)-\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)$
(19)

between the weighted macroparticle center phase-space densities (11) from the
simulation, seen here as a statistical realization of the plasma that one
wants to simulate, and the corresponding weighted macroparticle center
distribution functions

$\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)=\mathbb{E}\displaystyle\left\\{f_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)\right\\}$
(20)

that ones would obtain on average over an infinite number of simulations and
where we have thus noted $\mathbb{E}\displaystyle\left\\{X\right\\}$ the
expectation of the random variable $X$ from the probability theory. According
to the previous section 2.2, the latter reads

$\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t_{1}\right)=n_{a_{c}}\displaystyle\left(\mathbf{r},\,t_{1}\right)F_{a_{0}}\displaystyle\left(\mathbf{v}\right)$
(21)

at the simulation start $t_{1}$ where

$n_{a_{c}}\displaystyle\left(\mathbf{r},\,t_{1}\right)=\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}\delta\displaystyle\left(\mathbf{r}-\mathbf{R}_{a,\ell}\right)$
(22)

is the weighted macroparticle center spatial density. They lead to
electrostatic field statistical fluctuations

$\delta\mathbf{E}\displaystyle\left(\mathbf{r},\,t\right)=\mathbf{E}\displaystyle\left(\mathbf{r},\,t\right)-\bar{\mathbf{E}}\displaystyle\left(\mathbf{r},\,t\right)$
(23)

between the electrostatic field from the simulation and its expected value

$\bar{\mathbf{E}}\displaystyle\left(\mathbf{r},\,t\right)=\mathbb{E}\displaystyle\left\\{\mathbf{E}\displaystyle\left(\mathbf{r},\,t\right)\right\\}$
(24)

corresponding to the plasma macroscopic state that one wants to simulate. We
remind the reader that we’re dealing here with the statistical fluctuations of
a statistical realization of the plasma that one wants to simulate around its
mean behaviour described by its distribution functions and not with spatial
nor temporal fluctuations. In this sense, a PIC simulation is seen here as a
numerical experiment. Taking the expected value of (12) and (13) at grid
points $\mathbf{r}_{i,j,k}$ and time steps $t_{n}$, we obtain without any
approximation

$\begin{array}[]{lll}&&{\displaystyle\left.\displaystyle\frac{\partial\bar{f}_{a_{c}}}{\partial
t}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}\bar{f}_{a_{c}}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\bar{\mathbf{E}}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\bar{f}_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)\cr&=&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}_{s}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}\end{array}$
(25)

self-consistently coupled with the discretized macroscopic Maxwell equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\bar{\mathbf{E}}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\int_{\mathbb{R}^{3}}\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r}^{\prime},\,\mathbf{v},\,t_{n}\right)S\left(\mathbf{r}_{i,j,k}-\mathbf{r}^{\prime}\right)d^{3}\mathbf{r}^{\prime}d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\bar{\mathbf{E}}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (26)

Here,

$\bar{\mathbf{E}}_{s}\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)=\displaystyle\sum_{i,j,k}\bar{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)S\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right)-\mathbf{r}_{i,j,k}\right)\Delta_{x}\Delta_{y}\Delta_{z}$
(27)

is the expected electrostatic field interpolated at macroparticle center
locations and

$\delta\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\displaystyle\sum_{i,j,k}\delta{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)S\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right)-\mathbf{r}_{i,j,k}\right)\Delta_{x}\Delta_{y}\Delta_{z}$
(28)

the interpolated electrostatic field fluctuating component from the
simulation. Subtracting (25) from (12) and (13) from (26), we obtain the self-
consistently coupled equations

$\begin{array}[]{lll}&&{\displaystyle\left.\displaystyle\frac{\partial\delta
f_{a_{c}}}{\partial
t}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}\delta
f_{a_{c}}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\bar{f}_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)\cr&=&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)-\mathbb{E}{\displaystyle\left\\{\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}}\right)\cr&&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\bar{\mathbf{E}}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)\end{array}$
(29)

and

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\delta\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\int_{\mathbb{R}^{3}}\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}^{\prime},\,\mathbf{v},\,t_{n}\right)S\left(\mathbf{r}_{i,j,k}-\mathbf{r}^{\prime}\right)d^{3}\mathbf{r}^{\prime}d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\delta\mathbf{E}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (30)

They describe the evolution with time of statistical fluctuations introduced
by different initial conditions implemented by random numbers. Knowing the
statistical fluctuations
$\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{1}\right)$ and $\delta
f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{1}\right)$ in all spatial
grid cells $\mathbf{r}_{i,j,k}$ at the simulation start $t_{1}$, they allow to
deduce their values at any time steps $t_{n}$.

Let us emphasize here that the weighted macroparticle center phase-space
density fluctuations (19) are not the statistical fluctuations

$\delta
f_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)=f_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)-\bar{f}_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)$
(31)

between the expected plasma distribution functions

$\bar{f}_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)=\mathbb{E}\displaystyle\left\\{f_{a}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t\right)\right\\}$
(32)

that one wants to simulate and the particle phase-space densities (14) from
the PIC simulation. Therefore, one has to take the expected value of (15) and
(16) at grid points $\mathbf{r}_{i,j,k}$ and time steps $t_{n}$ in order to
get the probabilistic description of electrostatic PIC simulations. One thus
obtains without any approximation the discretized kinetic equation

$\begin{array}[]{lll}&&{\displaystyle\left.\displaystyle\frac{\partial\bar{f}_{a}}{\partial
t}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}\bar{f}_{a}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\bar{\mathbf{E}}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\bar{f}_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right)\cr&=&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}_{s}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}+\mathbb{E}{\displaystyle\left\\{C_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}}\end{array}$
(33)

that is self-consistently coupled with the discretized macroscopic Maxwell
equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\bar{\mathbf{E}}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\bar{f}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\bar{\mathbf{E}}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (34)

##  3 Statistical fluctuations in electrostatic PIC simulations

###  3.1 Assumptions and notations

We consider a stationary, homogeneous, infinite and fully ionized plasma
consisting of a stable PIC simulation plasma that is initialized as detailed
in the previous section 2.2 and that is reproduced periodically an infinite
number of times in all directions and that evolves perpetually from the
simulation start $t_{1}$ to infinity. Finite simulation box $L_{x}\times
L_{y}\times L_{z}$ and time duration $L_{t}$ would only modify the spectral
quantities by making them discrete with a resolution
$\mathbf{\Delta}_{k}=\left(2\pi/L_{x},\,2\pi/L_{y},\,2\pi/L_{z}\right)^{t}$
and $\Delta_{\omega}=2\pi/L_{t}$ due to the use of discrete inverse Laplace-
Fourier transforms instead of the continuous ones (42) and (51). In addition,
we assume

  1. 1.

Numerical consistency errors of numerical schemes used to compute the
equations are small compared to statistical fluctuations introduced by
different initial conditions implemented by random numbers.

  2. 2.

The weighted macroparticle centers spatial densities are homogeneous at the
simulation start

$\forall
i,\,j,\,k,\,n_{a_{c}}\left(\mathbf{r}_{i,j,k},\,t_{1}\right)=\bar{n}_{a}$ (35)

so that

$\forall
i,\,j,\,k,\,f_{a_{c}}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{1}\right)=\displaystyle\frac{\bar{n}_{a}}{N_{a,\mathrm{mpc}}}\displaystyle\sum_{\ell=1}^{N_{a,\mathrm{mpc}}}\delta\left(\mathbf{v}-\mathbf{V}_{a,\ell_{i,j,k}}\right)$
(36)

where
$\ell_{i,j,k}=\ell+\left[\left(i-1\right)+\left(j-1\right)N_{x}+\left(k-1\right)N_{x}N_{y}\right]N_{a,\mathrm{mpc}}$
and

$\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{1}\right)=\bar{n}_{a}F_{a_{0}}\left(\mathbf{v}\right).$
(37)

Therefore, we have also in this case

$\forall
i,\,j,\,k,\,\bar{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{1}\right)=\mathbf{0}$
(38)

according to macroscopic Maxwell equations (26), the macroparticle shape
properties (6) and Assumption 1.

  3. 3.

We neglect the quadratic terms

$\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\delta\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)-\mathbb{E}{\displaystyle\left\\{\delta\mathbf{E}_{s}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}}\right)$
(39)

describing the influence of macroparticle close-encounter Coulomb collisions
on the weighted macroparticle center discretized phase-space density
statistical fluctuations (29).

  4. 4.

The weighted macroparticle center distribution functions and macroscopic
electrostatic fields are stationary

$\forall
n,\,\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)=\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{1}\right)\,\mathrm{and}\,\bar{\mathbf{E}}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\bar{\mathbf{E}}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{1}\right).$
(40)

In order to check Assumption 2, let us introduce the Fourier transform in
space of discrete functions

$\widehat{\mathrm{F}}\left(\mathbf{k}\right)=\Delta_{x}\displaystyle\sum_{i=-\infty}^{\infty}\Delta_{y}\sum_{j=-\infty}^{\infty}\Delta_{z}\sum_{k=-\infty}^{\infty}F\displaystyle\left(\mathbf{r}_{i,j,k}\right)\exp{\displaystyle\left(-\iota\mathbf{k}\cdot\mathbf{r}_{i,j,k}\right)}.$
(41)

such that

$F\displaystyle\left(\mathbf{r}_{i,j,k}\right)=\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\widehat{\mathrm{F}}\left(\mathbf{k}\right)\exp{\displaystyle\left(\iota\mathbf{k}\cdot\mathbf{r}_{i,j,k}\right)}$
(42)

where the integration domain is the first Brillouin zone
${V_{\mathbf{k}_{g}}}=\left[-k_{g_{x}}/2,k_{g_{x}}/2\right]\times\left[-k_{g_{y}}/2,k_{g_{y}}/2\right]\times\left[-k_{g_{z}}/2,k_{g_{z}}/2\right]$.
The inverse Fourier transforms (42) can be seen as Fourier series components
of (41). Indeed, according to its mathematical expression, the latter is
necessarily $\mathbf{k}_{g}$-periodic where

$\mathbf{k}_{g}=\displaystyle\left(\begin{array}[]{c}2\pi/\Delta_{x}\cr
2\pi/\Delta_{y}\cr 2\pi/\Delta_{z}\end{array}\right)$ (43)

is the Nyquist frequency. It is defined as the highest spatial frequency that
one can resolve in a PIC simulation due to discretization of space. According
to (7), (9), (22), (41) and the Poisson summation formula applied first to the
infinite sums over macroparticle indices $\ell_{x}$, $\ell_{y}$ and $\ell_{z}$
and then to the infinite sums over spatial indices $i$, $j$ and $k$, one finds

$\widehat{n}_{a_{c}}\displaystyle\left(\mathbf{k},\,t_{1}\right)={\left(2\pi\right)}^{3}\bar{n}_{a}\displaystyle\sum_{p,q,r}{\left(-1\right)}^{p+q+r}\delta\left(\mathbf{k}-\mathbf{I}_{p,q,r}\mathbf{k}_{d},\,\mathbf{k}_{g}\right)$

where we have noted

$\delta\left(\mathbf{k},\,\mathbf{k}_{g}\right)=\displaystyle\sum_{p,q,r}\delta\displaystyle\left(\mathbf{k}-\mathbf{I}_{p,q,r}\cdot\mathbf{k}_{g}\right)=\delta\displaystyle\left(k_{x},k_{g_{x}}\right)\delta\displaystyle\left(k_{y},k_{g_{y}}\right)\delta\displaystyle\left(k_{z},k_{g_{z}}\right)$

the 3D Dirac comb distribution,

$\delta\displaystyle\left(k_{x},k_{g_{x}}\right)=\displaystyle\sum_{p=-\infty}^{\infty}\delta\displaystyle\left(k_{x}-pk_{g_{x}}\right),\,\mathbf{I}_{p,q,r}=\displaystyle\left(\begin{array}[]{ccc}p&0&0\cr
0&q&0\cr
0&0&r\end{array}\right)\,\mathrm{and}\,\mathbf{k}_{d}=\displaystyle\left(\begin{array}[]{c}2\pi
N_{a_{x}}/L_{x}\cr 2\pi N_{a_{y}}/L_{y}\cr 2\pi
N_{a_{z}}/L_{z}\end{array}\right)$

to simplify the notations. Considering only the first Brillouin zone due to
the discretization of space to compute Maxwell equations and a number of
macroparticles per cell distributed in all directions greater than one

$\forall\xi\in\displaystyle\left\\{x,\,y,\,z\right\\},\,\displaystyle\frac{N_{a_{\xi}}}{N_{\xi}}>1,$
(44)

the Nyquist frequencies $\mathbf{k}_{g_{\xi}}$ are necessarily lower than
$\mathbf{k}_{d_{\xi}}$ and we thus have necessarily

$\forall\mathbf{k}\in
V_{\mathbf{k}_{g}},\,\widehat{n}_{a_{c}}\displaystyle\left(\mathbf{k},\,t_{1}\right)={\left(2\pi\right)}^{3}\bar{n}_{a}\delta\displaystyle\left(\mathbf{k},\mathbf{k}_{g}\right).$
(45)

The inverse Fourier transform (42) of the latter validates consequently
Assumption 2 if the condition (44) is respected. A consequence of the
condition (44) is that the number of macroparticles per cell must necessarily
be $N_{a,\mathrm{mpc}}\geq 2$ in 1D, $\geq 4$ in 2D and $\geq 8$ in 3D
according to its definition (10) in such a way that, by depositing the same
amount strictly greater than 1 of macroparticles per cell in each directions,
the code necessarily sees homogeneous weighted macroparticle center spatial
densities.

Figure 1: Coulomb electrostatic field generated by one macroparticle $a,\ell$
compared to one ”point” particle of same electrical charge in 1D (sheet), 2D
(line) and 3D (point).

In Figure 1, we have plotted the Coulomb electrostatic field generated by one
macroparticle $a,\ell$ in the cases it has a Gaussian shape function (5) with
$a_{x}=a_{y}=a_{z}=\Delta$ or a cloud shape function (linear interpolating
function $\Pi_{\mathbf{\boldsymbol{\Delta}}}^{(1)}$) with
$\Delta=\Delta_{x}=\Delta_{y}=\Delta_{z}$. One can see that it is zero at the
macroparticle location. Then, it increases until the distance
$\left|\mathbf{r}-\mathbf{r}_{a,\ell}\right|\sim\Delta$ and it finally tends
to the classical Coulomb electrostatic field generated by one point particle
of equivalent electrical charge $Q_{a}$ when
$\left|\mathbf{r}-\mathbf{r}_{a,\ell}\right|\gtrapprox\Delta$. Therefore, long
range electrostatic interactions between macroparticles are correctly taken
into account in electrostatic PIC simulations while short range electrostatic
interaction between macroparticles when
$\left|\mathbf{r}-\mathbf{r}_{a,\ell}\right|<\Delta$ are underestimated [Okuda
and Birdsall, 1970] [Langdon and Birdsall, 1970]. We can thus neglect Coulomb
collisions between macroparticles and consider Assumption 3 to be valid. We
will check Assumptions 1 and 4 later in section 4.1 and section 4.2,
respectively. The following derivations are based on the theory of
fluctuations [Rostoker, 1961] [Klimontovich and Silin, 1962] [Dupree, 1963]
well summarized in the Landau and Lifschitz textbook [Landau and Lifshitz,
1981]. Notations used are recapitulated in Table 1.

$a$ | index of plasma species
---|---
$m_{a}$ | particle mass
$q_{a}$ | particle charge
$M_{a}$ | macroparticle mass
$Q_{a}$ | macroparticle charge
$\delta N_{a}$ | macroparticle weight
$S$ | macroparticle shape function
$N_{a}$ | total number of macroparticles
$\ell$ | index of macroparticles
$\mathbf{r}_{a,\ell}(t)$ | macroparticle center location
$\mathbf{R}_{a,\ell}$ | macroparticle center location at the simulation start
$\mathbf{v}_{a,\ell}(t)$ | macroparticle velocity
$\mathbf{V}_{a,\ell}$ | macroparticle velocity at the simulation start
$F_{a_{0}}(\mathbf{V}_{a})$ | probability distribution of macroparticle velocities at the simulation start
$\mathbf{r}$ | space variable
$\mathbf{v}$ | velocity variable
$t$ | time variable
$\mathbf{E}(\mathbf{r},\,t)$ | self-consistent electrostatic field
$\mathbf{E}_{s}(\mathbf{r}_{a,\ell}(t),\,t)$ | self-consistent electrostatic field interpolated at macroparticle center locations
$f_{a}(\mathbf{r},\,\mathbf{v}_{a},\,t)$ | plasma species phase-space density sampled according to the PIC method
$n_{a}(\mathbf{r},\,t)$ | plasma species density
$f_{a_{c}}(\mathbf{r},\,\mathbf{v}_{a},\,t)$ | weighted phase-space density of macroparticle centers
$n_{a_{c}}(\mathbf{r},\,t)$ | weighted density of macroparticle centers
$N_{x},\,N_{y},\,N_{z}$ | number of spatial grid points along $x$, $y$ and $z$-axis
$i,\,j,\,k$ | indices of spatial grid points along $x$, $y$ and $z$-axis
$\Delta_{x},\,\Delta_{y},\,\Delta_{z}$ | spatial cell size along $x$, $y$ and $z$-axis
$\mathbf{r}_{i,j,k}$ | spatial grid point location
$N_{t}$ | number of time iterations
$n$ | index of time iterations
$\Delta_{t}$ | time step
$t_{n}$ | iteration time
$N_{a,\mathrm{mpc}}$ | number of macroparticles per cell
$\widehat{\mathrm{F}}\left(\mathbf{k}\right)$ | Fourier transform of discrete function $F\left(\mathbf{r}_{i,j,k}\right)$
$\overset{\frown}{\mathrm{F}}\left(\mathbf{k}\right)$ | Fourier transform of continuous function $F\left(\mathbf{r}\right)$
$\widehat{\widehat{\mathrm{F}}}\left(\omega,\,\mathbf{k}\right)$ | Laplace-Fourier transform of discrete function $F\left(\mathbf{r}_{i,j,k},\,t_{n}\right)$
$k_{g,x},\,k_{g,y},\,k_{g,z}$ | Nyquist spatial angular frequency along $k_{x}$, $k_{y}$ and $k_{z}$-axis
$p,\,q,\,r$ | indices of spatial angular frequency aliases along $k_{x}$, $k_{y}$ and $k_{z}$-axis
$\mathbf{k}_{p,q,r}$ | spatial angular frequency vector alias
$\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}$ |  $\mathbf{k}$-bracket depending on the Maxwell solver
$\omega_{g}$ | Nyquist temporal angular frequency
$m$ | index of temporal angular frequency aliases
$\omega_{m}$ | temporal angular frequency aliases
$\displaystyle\left\\{\omega\right\\}_{t}$ |  $\omega$-bracket depending on the macroparticle pusher
$\bar{F}=\mathbb{E}\displaystyle\left\\{F\right\\}$ | expectation of stochastic process $F$ (sometimes called ensemble average)
$\delta{F}=F-\bar{F}$ | fluctuations of stochastic process $F$
$\mathbb{E}\displaystyle\left\\{\delta F(t)^{2}\right\\}$ | single-time autocorrelation of fluctuations $\delta F$
$\mathbb{E}\displaystyle\left\\{\delta F(t)\delta G(t)\right\\}$ | single-time correlation between fluctuations of stochastic processes $F$ and $G$
$N_{\text{\Large\textcalligra{s}}\,}$ | total number of statistical realizations
s | index of statistical realizations
$F^{(\text{\Large\textcalligra{s}}\,)}$ | statistical realization of $F$
$\langle F\rangle$ | finite ensemble average over $N_{\text{\Large\textcalligra{s}}\,}$ statistical realizations
Table 1: Recapitulation of main notations used in the text.

###  3.2 Estimate of statistical fluctuations in a PIC simulation
electrostatic plasma

According to their definitions (19), (20), (23) and (24), there is no
electrostatic field statistical fluctuations

$\mathbb{E}\displaystyle\left\\{\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\right\\}=0$
(46)

nor weighted macroparticle center phase-space density statistical fluctuations

$\mathbb{E}\displaystyle\left\\{\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}=0$
(47)

on average. Therefore, in order to get an estimate of the non-zero
electrostatic field -usually called ”numerical noise”- that one observes in
PIC simulations due to the macroparticle initial conditions implemented by
random numbers, we propose here to estimate it according to the square root of
its statistical fluctuation single-time autocorrelation

${\displaystyle\left|\mathbf{E}_{\mathrm{fluct}}\right|}=\displaystyle\sqrt{\mathbb{E}\displaystyle\left\\{\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)^{2}\right\\}}.$
(48)

Indeed, as it will be justified later in section 5.1, if the number of
macroparticles per cell is sufficiently large, the electrostatic field
statistical fluctuation
$\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,,t_{n}\right)$ at
grid points $\mathbf{r}_{i,j,k}$ and time steps $t_{n}$ follows a normal
distribution law with mean (46) and variance (48) according to the central
limit theorem. As a consequence,

$\displaystyle\left|\displaystyle\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\right|\underset{N_{e,\mathrm{mpc}}\gg
1}{\leq}\displaystyle\sqrt{2}\,\mathrm{erf}^{-1}\displaystyle\left(\displaystyle\frac{p}{100}\right)\displaystyle\left|\mathbf{E}_{\mathrm{fluct}}\right|$
(49)

with a $p\%$-confidence and the square root of electrostatic field statistical
fluctuation autocorrelation (48) thus provides a good estimate of statistical
fluctuations.

To carry out such a program, we use the Laplace-Fourier transforms of discrete
functions

$\begin{array}[]{lll}\widehat{\widehat{\mathrm{F}}}\left(\omega,\,\mathbf{k}\right)&=&\Delta_{t}\displaystyle\sum_{n=1}^{\infty}\Delta_{x}\displaystyle\sum_{i=-\infty}^{\infty}\Delta_{y}\sum_{j=-\infty}^{\infty}\Delta_{z}\sum_{k=-\infty}^{\infty}\displaystyle
F\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\exp{\displaystyle\left[\iota\displaystyle\left(\omega
t_{n}-\mathbf{k}\cdot\mathbf{r}_{i,j,k}\right)\right]}\cr&=&\Delta_{t}\displaystyle\sum_{n=1}^{\infty}\,\,\widehat{\mathrm{F}}\displaystyle\left(\mathbf{k},\,t_{n}\right)\exp{\displaystyle\left(\iota\omega
t_{n}\right)}\end{array}$ (50)

such that

$F\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\,\,\widehat{\widehat{\mathrm{F}}}\left(\omega,\,\mathbf{k}\right)\exp{\displaystyle\left[-\iota\displaystyle\left(\omega
t_{n}-\mathbf{k}\cdot\mathbf{r}_{i,j,k}\right)\right]}$ (51)

where $\nu\in\mathbb{R}$ such that the contour path of integration is in the
region of convergence. Again, according to their definitions (50), one can see
that they are necessarily $\left(\omega_{g},\,\mathbf{k}_{g}\right)$-periodic
where the Nyquist frequencies are (43) and

$\omega_{g}=2\pi/\Delta_{t}.$ (52)

This is not the case for the macroparticle shape function (4) whose support is
necessarily continuous since macroparticle positions
$\mathbf{r}_{a,\ell}\left(t_{n}\right)$ may have spatial variations less than
the grid spacing. We will therefore note

$\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)=\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{r}\displaystyle
S\displaystyle\left(\mathbf{r}\right)\exp{\displaystyle\left(-\iota\mathbf{k}\cdot\mathbf{r}\right)}$
(53)

its continuous Fourier transform such that

$\displaystyle
S\displaystyle\left(\mathbf{r}\right)=\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)\exp{\displaystyle\left(\iota\mathbf{k}\cdot\mathbf{r}\right)}$
(54)

to avoid misunderstandings. For example, in a spectral simulation using a
$n$th-order B-spline interpolating function and a Gaussian filter, it simply
reads

$\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)=\prod_{\xi=x,y,z}{\displaystyle\left[\displaystyle\frac{\sin{\displaystyle\left(k_{\xi}\Delta_{\xi}/2\right)}}{k_{\xi}\Delta_{\xi}/2}\right]}^{n+1}\exp{\displaystyle\left(-\displaystyle\frac{{a_{\xi}}^{2}{k_{\xi}}^{2}}{2}\right)}.$

By linearizing the leap-frog macroparticle pusher algorithm

$\displaystyle\left\\{\begin{array}[]{clllr}\displaystyle\frac{\mathbf{r}_{a,\ell}\displaystyle\left(t_{n+1}\right)-\mathbf{r}_{a,\ell}\displaystyle\left(t_{n}\right)}{\Delta_{t}}&=&\mathbf{v}_{a,\ell}\displaystyle\left(t_{n+1/2}\right),&\mathbf{r}_{a,\ell}\displaystyle\left(t_{1}\right)&=\mathbf{R}_{a,\ell}\cr\displaystyle\frac{\mathbf{v}_{a,\ell}\displaystyle\left(t_{n+1/2}\right)-\mathbf{v}_{a,\ell}\displaystyle\left(t_{n-1/2}\right)}{\Delta_{t}}&=&\displaystyle\frac{Q_{a}}{M_{a}}\mathbf{E}_{s}\left(\mathbf{r}_{a,\ell}^{n},\,t_{n}\right),&\mathbf{v}_{a,\ell}\displaystyle\left(t_{1/2}\right)&=\mathbf{V}_{a,\ell}\end{array}\right.,$
(55)

that is commonly used to compute the discretized macroparticle equations of
motions (1), and by associating the non-perturbed straight line orbits
$(\mathbf{r}_{a,\ell}^{(0)}\left(t_{n}\right),\,\mathbf{v}_{a,\ell}^{(0)}\left(t_{n}\right))=\left(\mathbf{R}_{a,\ell}+\mathbf{V}_{a,\ell}t,\,\mathbf{V}_{a,\ell}\right)$
with $\bar{f}_{a_{c}}$ and the first order perturbed orbits
$\left(\delta\mathbf{r}_{a,\ell}\left(t_{n}\right),\,\delta\mathbf{v}_{a,\ell}\left(t_{n}\right)\right)$
with $\delta f_{a_{c}}$ according to Assumptions 1, 2, 3 and 4 [Lindman,
1970], the discrete Laplace-Fourier transform (50) of the resulting
discretized equation (29) gives

$\iota\displaystyle\left\\{\mathbf{k}\cdot\mathbf{v}-\omega\right\\}_{t}\widehat{\widehat{\delta
f}}_{a_{c}}\displaystyle\left(\omega,\,\mathbf{k},\,\mathbf{v}\right)+\bar{n}_{a}\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\sum_{p,q,r}\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)\widehat{\widehat{\delta\mathbf{E}}}\displaystyle\left(\omega,\,\mathbf{k}_{p,q,r}\right)\cdot\displaystyle\frac{dF_{a_{0}}}{d\mathbf{v}}=\widehat{\delta
f}_{a_{c}}\displaystyle\left(\mathbf{k},\,\mathbf{v},\,t_{1}\right)$ (56)

where

$\displaystyle\frac{1}{\displaystyle\left\\{\mathbf{k}\cdot\mathbf{v}-\omega\right\\}_{t}}=\displaystyle\frac{\Delta_{t}}{2}\cot{\displaystyle\left[\displaystyle\left(\mathbf{k}\cdot\mathbf{v}-\omega\right)\displaystyle\frac{\Delta_{t}}{2}\right]}=\displaystyle\sum_{m}\displaystyle\frac{1}{\mathbf{k}\cdot\mathbf{v}-\omega_{m}}$
(57)

[Langdon, 1970b] [Birdsall and Langdon, 1991] [Abramowitz and Stegun, 1965].
Here, we have noted

$\omega_{m}=\omega-m\omega_{g}$ (58)

and

$\mathbf{k}_{p,q,r}=\mathbf{k}-\mathbf{I}_{p,q,r}\cdot\mathbf{k}_{g}$ (59)

the different temporal and spatial angular frequency aliases, following the
same notation as in [Langdon, 1970a] and [Birdsall and Langdon, 1991]. Maxwell
equations (13) are usually computed by using a spectral solver or the second
order FDTD method proposed by [Yee, 1966] coupled with the charge conserving
scheme proposed by [Villasenor and Buneman, 1992] or [Esirkepov, 2001] rather
than the Poisson equation solver considered in [Langdon, 1970b]. In all cases,
the Laplace-Fourier transform (50) of the resulting Maxwell equations (30)
linearized according to Assumptions 1, 2 and 4 gives

$\displaystyle\left\\{\begin{array}[]{lllll}\iota\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}\,\,\cdot\,\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\,\mathbf{k}\right)&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\sum_{p^{\prime},q^{\prime},r^{\prime}}\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\right)\displaystyle\int_{\mathbb{R}^{3}}\widehat{\widehat{\delta
f}}_{a_{c}}\displaystyle\left(\omega,\,\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}},\,\mathbf{v}\right)d^{3}\mathbf{v}\cr\iota\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}\times\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\,\mathbf{k}\right)&=&\mathbf{0}\end{array}\right.$
(60)

where we have noted $\forall\mathbf{k}\in V_{\mathbf{k}_{g}},$

${\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}=\displaystyle\left\\{\begin{array}[]{lcl}&{\displaystyle\left(k_{x},\,k_{y},\,k_{z}\right)}^{t}&\mathrm{\,for\,the\,spectral\,solver}\cr&{\displaystyle\left(\displaystyle\frac{\sin{\left(k_{x}\Delta_{x}/2\right)}}{\Delta_{x}/2},\,\displaystyle\frac{\sin{\left(k_{y}\Delta_{y}/2\right)}}{\Delta_{y}/2},\,\displaystyle\frac{\sin{\left(k_{z}\Delta_{z}/2\right)}}{\Delta_{z}/2}\right)}^{t}&\mathrm{\,for\,the\,2nd\,order\,FDTD\,scheme}\end{array}\right.}$
(61)

and $\forall p,\,q,\,r$,
$\displaystyle\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}=\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}$.
By substituting $\widehat{\widehat{\delta
f}}_{a_{c}}\displaystyle\left(\omega,\,\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}},\,\mathbf{v}\right)$
from (56) into the expression of
$\delta\widehat{\widehat{\mathbf{E}}}\left(\omega,\,\mathbf{k}\right)$ from
(60), replacing $\mathbf{k}$ by $\mathbf{k}_{p,q,r}$, changing the indices
$\left(p^{\prime},q^{\prime},r^{\prime}\right)$ into
$\left(p^{\prime}-p,q^{\prime}-q,r^{\prime}-r\right)$ in the resulting
equation, one can impose then
$\left(p,q,r\right)=\left(-p^{\prime},-q^{\prime},-r^{\prime}\right)$ and use
the $\mathbf{k}_{g}$-periodicity of
$\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}$,
$\widehat{\widehat{\delta\mathbf{E}}}$ and $\delta\widehat{f}_{a_{c}}$ to
obtain the discrete Laplace-Fourier transform (50) of electrostatic field
statistical fluctuations [Langdon, 1970a]. It reads

$\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\,\mathbf{k}\right)=\displaystyle\frac{-4\pi\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right){\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}^{2}}\displaystyle\sum_{p,q,r}\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)\displaystyle\sum_{a}q_{a}\sum_{m}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{\delta\widehat{f}_{a_{c}}\left(\mathbf{k},\,\mathbf{v},\,t_{1}\right)}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}}d^{3}\mathbf{v}$
(62)

where the longitudinal permittivity of the PIC simulation plasma reads

$\begin{array}[]{lll}\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)&=&1-\displaystyle\sum_{a}\displaystyle\frac{{\omega_{p_{a}}}^{2}}{{\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}^{2}}\displaystyle\sum_{p,q,r}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}\displaystyle\sum_{m}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}}\cdot\displaystyle\frac{dF_{a0}}{d\mathbf{v}}d^{3}\mathbf{v}\end{array}$
(63)

with

${\omega_{p_{a}}}=\displaystyle\sqrt{\displaystyle\frac{4\pi\bar{n}_{a}{q_{a}}^{2}}{m_{a}}}$
(64)

the Langmuir angular frequencies of each plasma species. We can check that it
tends to the longitudinal permittivity of collisionless plasmas in the limit
$\forall\xi\in\left\\{x,\,y,\,z,\,t\right\\},\,\Delta_{\xi}\rightarrow 0$. The
study of poles $\omega=\omega\left(\mathbf{k}\right)$ given by the dispersion
relation $\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)=0$
provides the stability criterion of the simulated plasma. We have assumed the
simulated plasma remains stable, both physically and numerically. It means
that all these poles have a negative imaginary part
$\mathrm{Im}\left\\{\omega\left(\mathbf{k}\right)\right\\}<0$. In general, the
stability criterion depends on the macroparticle shape and initial
macroparticle velocity probability distributions $F_{a_{0}}$. For example, it
has been shown that electrostatic PIC simulations of plasmas at equilibrium
using the first or greater orders B-spline interpolating function (and same
spatial spacings in all directions for simplicity) are stable if the condition

$\Delta_{t}\leq\Delta_{x}/v_{T_{e}}<1/\omega_{p}$ (65)

is respected [Birdsall and Langdon, 1991] where

${\omega_{p}}^{2}={\omega_{p_{e}}}^{2}+{\omega_{p_{i}}}^{2}\approx{\omega_{p_{e}}}^{2}$
(66)

is the Langmuir plasma angular frequency in the limit $Zm_{e}/m_{i}\ll 1$.

One can now deduce the expression of electrostatic field statistical
fluctuations autocorrelation from the discrete Laplace-Fourier transform of
electrostatic field statistical fluctuations (62). It reads

$\begin{array}[]{lll}\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}^{2}\right\\}&=&\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu^{\prime}-\omega_{g}/2}^{\iota\nu^{\prime}+\omega_{g}/2}\displaystyle\frac{d\omega^{\prime}}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}^{\prime}}{{\left(2\pi\right)}^{3}}\cr&&\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta\mathbf{E}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}\exp{\displaystyle\left\\{-\iota\displaystyle\left[\left(\omega+\omega^{\prime}\right)t_{n}-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{r}_{i,j,k}\right]\right\\}}\end{array}$
(67)

where

$\begin{array}[]{ll}\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta\mathbf{E}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}&=\displaystyle\frac{{\left(4\pi\right)}^{2}}{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)\varepsilon_{L}\displaystyle\left(\omega^{\prime},\mathbf{k}^{\prime}\right)}\displaystyle\frac{{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}\cdot{\left\\{\mathbf{k}^{\prime}\right\\}_{\mathbf{r}}}}{{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}^{2}{\left\\{\mathbf{k}^{\prime}\right\\}_{\mathbf{r}}}^{2}}\displaystyle\sum_{p,q,r}\displaystyle\sum_{p^{\prime},q^{\prime},r^{\prime}}\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)\overset{\frown}{\mathrm{S}}\left(\mathbf{k}^{\prime}_{p^{\prime},q^{\prime},r^{\prime}}\right)\cr&\displaystyle\sum_{a,b}{q_{a}}{q_{b}}\displaystyle\sum_{m,m^{\prime}}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{\mathbb{E}\displaystyle\left\\{\delta\widehat{f}_{a_{c}}\left(\mathbf{k},\,\mathbf{v},\,t_{1}\right)\delta\widehat{f}_{b_{c}}\left(\mathbf{k}^{\prime},\,\mathbf{v}^{\prime},\,t_{1}\right)\right\\}}{\displaystyle\left(\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}\right)\displaystyle\left({\mathbf{k}^{\prime}}_{p^{\prime},q^{\prime},r^{\prime}}\cdot\mathbf{v}^{\prime}-{\omega^{\prime}}_{m^{\prime}}\right)}d^{3}\mathbf{v}d^{3}\mathbf{v}^{\prime}.\end{array}$
(68)

It depends on the weighted macroparticle center phase-space density
fluctuations Fourier transforms correlation at the simulation start $t_{1}$.
One can derive the exact analytical expression of the latter. One finds
according to section 2.2

$\begin{array}[]{l}\mathbb{E}\displaystyle\left\\{\delta
f_{a_{c}}\left(\mathbf{r},\,\mathbf{v},\,t_{1}\right)\delta
f_{b,c}\left(\mathbf{r}^{\prime},\,\mathbf{v}^{\prime},\,t_{1}\right)\right\\}=\delta_{ab}\delta
N_{a}\delta\left(\mathbf{r}-\mathbf{r}^{\prime}\right)n_{a_{c}}\displaystyle\left(\mathbf{r},\,t_{1}\right)F_{a0}\displaystyle\left(\mathbf{v}\right)\displaystyle\left[\delta\displaystyle\left(\mathbf{v}-\mathbf{v}^{\prime}\right)-F_{a0}\displaystyle\left(\mathbf{v}^{\prime}\right)\right].\end{array}$
(69)

The first term in the square bracket corresponds to the case where there is
only one macroparticle in the infinitezimal phase-space volume located between
$\left(\mathbf{r},\,\mathbf{v}\right)$ and
$\left(\mathbf{r}+d^{3}\mathbf{r},\,\mathbf{v}+d^{3}\mathbf{v}\right)$.
Replacing consequently $\delta\left(\mathbf{r}-\mathbf{r}^{\prime}\right)$ by
its equivalent expression
$\delta_{i,i^{\prime}}\delta_{j,j^{\prime}}\delta_{k,k^{\prime}}/\Delta_{x}\Delta_{y}\Delta_{z}$
in a discretized space, one can compute the Fourier transform (41) of the
initial weighted macroparticle center phase-space densities correlation (69)
expressed at the grid points $\mathbf{r}_{i,j,k}$ and
$\mathbf{r}_{i^{\prime},j^{\prime},k^{\prime}}$. It reads

$\begin{array}[]{l}\mathbb{E}\displaystyle\left\\{\delta\widehat{f}_{a_{c}}\left(\mathbf{k},\,\mathbf{v},\,t_{1}\right)\delta\widehat{f}_{b,c}\left(\mathbf{k}^{\prime},\,\mathbf{v}^{\prime},\,t_{1}\right)\right\\}=\delta_{ab}\delta
N_{a}{\displaystyle\left(2\pi\right)}^{3}\bar{n}_{a}\delta\left(\mathbf{k}+\mathbf{k}^{\prime},\,\mathbf{k}_{g}\right)F_{a0}\displaystyle\left(\mathbf{v}\right)\displaystyle\left[\delta\displaystyle\left(\mathbf{v}-\mathbf{v}^{\prime}\right)-F_{a0}\displaystyle\left(\mathbf{v}^{\prime}\right)\right]\end{array}$
(70)

according to Assumption 2. By performing successively the integration over
$\mathbf{k}^{\prime}$, $\mathbf{v}^{\prime}$ and $\omega^{\prime}$ in (67)
(see Appendix A), one finds

$\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta\mathbf{E}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}\underset{\omega_{p}t_{n}\gg
1}{=}{\displaystyle\left(2\pi\right)}^{4}\delta\displaystyle\left(\omega+\omega^{\prime},\,\omega_{g}\right)\delta\displaystyle\left(\mathbf{k}+\mathbf{k}^{\prime},\,\mathbf{k}_{g}\right)\widehat{\widehat{{\delta\mathbf{E}}^{2}}}\left(\omega,\mathbf{k}\right)$
(71)

with

$\begin{array}[]{lcl}\widehat{\widehat{{\delta\mathbf{E}}^{2}}}\left(\omega,\mathbf{k}\right)&=&\displaystyle\int_{\iota\nu^{\prime}-\omega_{g}/2}^{\iota\nu^{\prime}+\omega_{g}/2}\displaystyle\frac{d\omega^{\prime}}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}^{\prime}}{{\left(2\pi\right)}^{3}}\cr&&\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta\mathbf{E}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}\exp{\displaystyle\left\\{-\iota\displaystyle\left[\left(\omega+\omega^{\prime}\right)t_{n}-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{r}_{i,j,k}\right]\right\\}}\cr&\underset{\omega_{p}t_{n}\gg
1}{=}&\displaystyle\frac{32\pi^{3}}{{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)}^{2}}\displaystyle\frac{1}{{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}^{2}}\displaystyle\sum_{p,q,r}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}\displaystyle\sum_{a}\delta
N_{a}\bar{n}_{a}{q_{a}}^{2}\displaystyle\sum_{m}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}F_{a0}\displaystyle\left(\mathbf{v}\right)\delta\displaystyle\left(\omega_{m}-\mathbf{k}_{p,q,r}\cdot\mathbf{v}\right).\end{array}$
(72)

in the expression of single-time electrostatic field statistical fluctuations
autocorrelation (67). For example, we find for PIC simulation plasmas at
equilibrium the finite estimate (see Appendix B)

$\displaystyle\left|\mathbf{E}_{\mathrm{fluct}}\right|\underset{\omega_{p}t_{n}\gg
1}{\sim}\displaystyle\sqrt{\delta
Nk_{B}T\eta\left(k_{g}\lambda_{D}\right)}\,\mathrm{with}\,\eta\left(k_{g}\lambda_{D}\right)=\displaystyle\left\\{\begin{array}[]{cll}\displaystyle\frac{4\arctan{\displaystyle\left(k_{g}\lambda_{D}/2\right)}}{{\lambda_{D}}}&\,\mathrm{in}\,1\mathrm{D}\cr\displaystyle\frac{1}{2}\displaystyle\frac{\ln{\left[1+{(k_{g}\lambda_{D}/2)}^{2}\right]}}{{\lambda_{D}}^{2}}&\,\mathrm{in}\,2\mathrm{D}\cr\displaystyle\frac{2}{3}\displaystyle\frac{k_{g}\lambda_{D}-2\arctan{\displaystyle\left(k_{g}\lambda_{D}/2\right)}}{\pi{\lambda_{D}}^{3}}&\,\mathrm{in}\,3\mathrm{D}\end{array}\right..$
(73)

These estimates are compared with the amplitude of electrostatic field
statistical fluctuations in electrostatic PIC plasma simulations in Figure 2
considering $\Delta_{x}=\Delta_{y}=\Delta_{z}$ for simplicity and using

$\delta
N=\displaystyle\frac{\bar{n}_{e}{\Delta_{x}}^{\text{dim}}}{N_{e,\text{mpc}}}$

according to (7) and (10).

Figure 2: Average over all spatial cells of electrostatic field fluctuations
amplitudes at the discrete time $t_{n}=10/\omega_{p}$ (Left) and their
analytical estimates (73) (Right) in 1D (Top), 2D (Middle) and 3D (Bottom)
periodic electrostatic PIC plasmas at Maxwell-Boltzmann equilibrium simulated
with the spectral electrostatic PIC codes from the UPIC framework [Decyk,
1995] [Decyk, 2007] [Decyk and Singh, 2014]; each simulation out of the 30,000
uses a total of 4096 spatial cells and interpolation functions
$\Pi^{(1)}_{\Delta_{x}}$ without filtering.

By performing the similar derivation for the single-time PIC simulation plasma
spatial density statistical fluctuations autocorrelation for each species, one
also finds

$\begin{array}[]{lll}\mathbb{E}\displaystyle\left\\{{\delta
n_{a}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}^{2}\right\\}&=&\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu^{\prime}-\omega_{g}/2}^{\iota\nu^{\prime}+\omega_{g}/2}\displaystyle\frac{d\omega^{\prime}}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}^{\prime}}{{\left(2\pi\right)}^{3}}\cr&&\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta
n_{a}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta
n_{a}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}\exp{\displaystyle\left\\{-\iota\displaystyle\left[\left(\omega+\omega^{\prime}\right)t_{n}-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{r}_{i,j,k}\right]\right\\}}\end{array}$
(74)

with

$\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta
n_{a}}}\left(\omega,\mathbf{k}\right)\widehat{\widehat{\delta
n_{a}}}\left(\omega^{\prime},\,\mathbf{k}^{\prime}\right)\right\\}\underset{\omega_{p}t_{n}\gg
1}{=}{\displaystyle\left(2\pi\right)}^{4}\delta\displaystyle\left(\omega+\omega^{\prime},\,\omega_{g}\right)\delta\displaystyle\left(\mathbf{k}+\mathbf{k}^{\prime},\,\mathbf{k}_{g}\right)\widehat{\widehat{{\delta
n_{a}}^{2}}}\left(\omega,\mathbf{k}\right)$ (75)

and

$\widehat{\widehat{{\delta
n_{a}}^{2}}}\left(\omega,\mathbf{k}\right)\underset{\omega_{p}t_{n}\gg
1}{=}\displaystyle\frac{2\pi\delta
N_{a}\bar{n}_{a}}{{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)}^{2}}\displaystyle\sum_{p,q,r}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}F_{a0}\displaystyle\left(\mathbf{v}\right)\delta\displaystyle\left(\omega-\mathbf{k}_{p,q,r}\cdot\mathbf{v},\omega_{g}\right).$
(76)

Note that these lasts are the single-time autocorrelation of plasma particle
densities statistical fluctuations

$\delta
n_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\displaystyle\int_{\mathbb{R}^{3}}\delta
f_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)d^{3}\mathbf{v}$
(77)

and not the one of weighted macroparticle center spatial densities statistical
fluctuations

$\delta
n_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\displaystyle\int_{\mathbb{R}^{3}}\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)d^{3}\mathbf{v}.$
(78)

The single-time plasma electrical charge statistical fluctuations
autocorrelation can be obtained starting from (72) or from (76). One can note
that initial single-time statistical fluctuations correlation of differing
plasma species quantities necessarily cancels according to (69). Fluctuations
spectra (72) and (76) are the same as the ones obtained by [Langdon, 1970b]
[Langdon, 1979], however generalizing them to simulations of multiple species
plasmas using today’s algorithms and arbitrary macroparticle weights factor
$\delta N_{a}$. To derive these spectra, the authors used a test macroparticle
formalism similar to the method used by [Thompson and Hubbard, 1960] [Hubbard,
1961] and later justified by [Rostoker, 1964] from the classical kinetic
theory of plasmas. It consists in perturbing a discretized Vlasovian plasma by
one macroparticle propagating along a straight-line orbit and averaging the
resulting linear PIC simulated plasma quantity over all possible macroparticle
initial conditions. The validity of their method is ensured by the first term,
corresponding to the case where there is only one macroparticle in the
infinitesimal volume of phase-space, in our results (69) of initial weighted
macroparticle center phase-space density correlations. While providing the
physical interpretation of the results up to the first order, it may lead to
confusion between Vlasov simulations and PIC simulations. Therefore, it seems
to us that our method relates with fewer ambiguity equations computed by an
electrostatic PIC code with the kinetic theory of plasmas. In addition, our
approach adds the information that the statistical fluctuation spectra are
only valid on a time scale $\omega_{p}t_{n}\gg 1$. It means that one must let
the simulation run during a time greater than $\sim 10/\omega_{p}$ before
obtaining comparable spectra by directly computing it from the PIC simulation
results.

###  3.3 Kinetic equation of a PIC simulation electrostatic plasma

Performing the same reasoning as in the previous section, one can also find
the single-time correlation between interpolated electrostatic field
statistical fluctuations and weighted macroparticle center phase-space density
statistical fluctuations at the grid points $\mathbf{r}_{i,j,k}$ and time
steps $t_{n}$ for each species. We find starting from (56), (62) and (70),
assuming again Assumptions 1, 2, 3, 4 and considering again only non-damped
terms on time scales $\omega_{p}t_{n}\gg 1$

$\begin{array}[]{l}\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}_{s}\left(\omega^{\prime},\mathbf{k}^{\prime}\right)\widehat{\widehat{\delta
f}}_{a_{c}}\left(\omega,\,\mathbf{k},\,\mathbf{v}_{a}\right)\right\\}\underset{\omega_{p}t_{n}\gg
1}{=}{\displaystyle\left(2\pi\right)}^{4}\delta\displaystyle\left(\omega+\omega^{\prime},\,\omega_{g}\right)\delta\displaystyle\left(\mathbf{k}+\mathbf{k}^{\prime},\,\mathbf{k}_{g}\right)\widehat{\widehat{\delta\mathbf{E}_{s}\delta
f}}_{a_{c}}\left(\omega,\mathbf{k},\,\mathbf{v}_{a}\right)\end{array}$ (79)

where

$\begin{array}[]{lcl}\widehat{\widehat{\delta\mathbf{E}_{s}\delta
f}}_{a_{c}}\left(\omega,\mathbf{k},\,\mathbf{v}_{a}\right)&\underset{\omega_{p}t_{n}\gg
1}{=}&-\iota\displaystyle\sum_{p,q,r}\displaystyle\frac{8\pi^{2}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}}{{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)}}\displaystyle\frac{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}{{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{2}}\delta
N_{a}\bar{n}_{a}{q_{a}}\displaystyle\sum_{m}F_{a0}\displaystyle\left(\mathbf{v}\right)\delta\displaystyle\left(\omega_{m}-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a}\right)\cr&+&\hskip
6.49994pt\iota\displaystyle\sum_{p,q,r}\bar{n}_{a}\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\sum_{m}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\displaystyle\left(\mathbf{k}_{p,q,r}\right)}^{2}}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a}-\omega_{m}}\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}\cdot\displaystyle\frac{dF_{a0}}{d\mathbf{v}_{a}}\widehat{\widehat{{\delta\mathbf{E}}^{2}}}\left(\omega,\mathbf{k}\right)\displaystyle\frac{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}{{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{2}}.\end{array}$
(80)

It allows for estimating the discretized right-hand side of kinetic equations
(25) depending on

$\begin{array}[]{lcl}&&\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{a},\,t_{n}\right)\right\\}\cr&=&\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu^{\prime}-\omega_{g}/2}^{\iota\nu^{\prime}+\omega_{g}/2}\displaystyle\frac{d\omega^{\prime}}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}^{\prime}}{{\left(2\pi\right)}^{3}}\cr&&\mathbb{E}\displaystyle\left\\{\widehat{\widehat{\delta\mathbf{E}}}_{s}\left(\omega^{\prime},\mathbf{k}^{\prime}\right)\widehat{\widehat{\delta
f}}_{a_{c}}\left(\omega,\,\mathbf{k},\,\mathbf{v}_{a}\right)\right\\}\exp{\displaystyle\left\\{-\iota\displaystyle\left[\left(\omega+\omega^{\prime}\right)t_{n}-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{r}_{i,j,k}\right]\right\\}}\cr&\underset{\omega_{p}t_{n}\gg
1}{=}&\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\widehat{\widehat{\delta\mathbf{E}_{s}\delta
f}}_{a_{c}}\left(\omega,\mathbf{k},\,\mathbf{v}\right).\end{array}$ (81)

One finds after integration over $\omega$ (see Appendix C)

$\begin{array}[]{lll}&&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}_{a}}\cdot\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}\delta
f_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}\cr&\underset{\omega_{p}t_{n}\gg
1}{=}&-\displaystyle\frac{1}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}_{a}}\cdot\displaystyle\sum_{b}\bar{n}_{a}\bar{n}_{b}\displaystyle\int_{\mathbb{R}^{3}}\mathbf{Q}\displaystyle\left(\mathbf{v}_{a},\,\mathbf{v}_{b}\right)\cdot\displaystyle\left[\displaystyle\frac{\delta
N_{b}}{m_{a}}F_{b0}\left(\mathbf{v}_{b}\right)\displaystyle\frac{dF_{a0}}{d\mathbf{v}_{a}}-\displaystyle\frac{\delta
N_{a}}{m_{b}}F_{a0}\left(\mathbf{v}_{a}\right)\displaystyle\frac{dF_{b0}}{d\mathbf{v}_{b}}\right]d^{3}\mathbf{v}_{b}\end{array}$
(82)

depending on the tensor

$\mathbf{Q}\displaystyle\left(\mathbf{v}_{a},\,\mathbf{v}_{b}\right)=-2\,{q_{a}}^{2}{q_{b}}^{2}\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}\otimes{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}}{{\varepsilon_{L}\left(\mathbf{k}\cdot\mathbf{v}_{a},\,\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}\right\\}_{\mathbf{r}}}^{4}}\displaystyle\sum_{p,q,r}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}\displaystyle\sum_{m}\delta\displaystyle\left(\mathbf{k}\cdot\mathbf{v}_{a}-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{b}-m\omega_{g}\right)d^{3}\mathbf{k}$
(83)

for 3D simulations ($d^{3}\mathbf{k}$ and the entire 3D space of integration
$\mathbb{R}^{3}$ should be replaced by ${(2\pi)}^{2}dk_{x}$ and $\mathbb{R}$
or by $2\pi\,dk_{x}dk_{y}$ and $\mathbb{R}^{2}$ for 1D or 2D simulations
respectively). It tends to the Lenard-Balescu tensor for collisionless plasmas
in the limit
$\forall\xi\in\left\\{x,\,y,\,z,\,t\right\\},\,\Delta_{\xi}\rightarrow 0$.
According to Assumptions 2, 4 and the macroparticle shape function properties
(6), we have necessarily

$\bar{n}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\delta
N_{a}\displaystyle\sum_{i,j,k}S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\ell}\left(t_{1}\right)\right)=\bar{n}_{a}.$
(84)

and therefore

$\bar{f}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)=\bar{n}_{a}F_{a_{0}}\displaystyle\left(\mathbf{v}\right)=\bar{f}_{a_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right).$
(85)

Neglecting the expectation of the term depending on the internal tension
forces of macroparticles (17), kinetic equations corresponding to an
electrostatic PIC simulation follow immediately. They read

$\begin{array}[]{l}\hskip
20.00003pt{\displaystyle\left.\displaystyle\frac{\partial{\bar{f}}_{a}}{\partial{t}}\right|}^{i,j,k,n}+{\displaystyle\left.\displaystyle\frac{\partial}{\partial{\mathbf{r}}}\cdot\displaystyle\left({\mathbf{v}}_{a}{\bar{f}}_{a}\right)\right|}^{i,j,k,n}+\displaystyle\frac{\partial}{\partial{\mathbf{v}}_{a}}\cdot\displaystyle\left(\displaystyle\frac{{q}_{a}}{{m}_{a}}{\bar{\mathbf{E}}}_{s}\left({\mathbf{r}}_{i,j,k},\,{t}_{n}\right){\bar{f}}_{a}\left({\mathbf{r}}_{i,j,k},\,{\mathbf{v}}_{a},\,{t}_{n}\right)\right)\cr\underset{\omega_{p}{t}_{n}\gg
1}{=}-\displaystyle\frac{1}{{m}_{a}}\displaystyle\frac{\partial}{\partial{\mathbf{v}}_{a}}\cdot\displaystyle\sum_{b}\displaystyle\int_{\mathbb{R}^{3}}{\mathbf{Q}}\displaystyle\left({\mathbf{v}}_{a},\,{\mathbf{v}}_{b}\right)\cdot{\displaystyle\left.\displaystyle\left(\displaystyle\frac{{\delta
N}_{b}}{{m}_{a}}{\bar{f}}_{b}\left({\mathbf{v}}_{b}\right)\displaystyle\frac{\partial{\bar{f}}_{a}}{\partial{\mathbf{v}}_{a}}-\displaystyle\frac{{\delta
N}_{a}}{{m}_{b}}{\bar{f}}_{a}\left({\mathbf{v}}_{a}\right)\displaystyle\frac{\partial{\bar{f}}_{b}}{\partial{\mathbf{v}}_{b}}\right)\right|}^{i,j,k,n}d^{3}{\mathbf{v}}_{b}\end{array}$
(86)

according to (25) and (82), in agreement with [Langdon, 1970b] [Birdsall and
Langdon, 1991]. The component
$\left(m,\,p,\,q,\,r\right)=\left(0,\,0,\,0,\,0\right)$ in (82) is physical.
As shown by [Dawson, 1983], it describes the average friction and diffusion of
macroparticles in the fluctuating component of the electrostatic force which
are direct consequences from the emission [Decyk, 1987] and absorption
[Landau, 1946] of plasma waves by the macroparticles in the simulation.
However, if each macroparticle represent more than one real particle such that
$\delta N_{a}=\delta N_{b}=\delta N\gg 1$, these physical processes are
greatly overestimated by the weight factor $\delta N$ compared to reality
according to the expression of the ”collision operator” (82). While the
component $\left(m,\,p,\,q,\,r\right)=\left(0,\,0,\,0,\,0\right)$ of the
latter cancels for single species 1D plasmas [Dawson, 1964] and for plasmas at
equilibrium for which $\delta N_{e}=\delta N_{i}$ such that the diffusion of
macroparticles compensates exactly their friction (in a plasma in thermal
equilibrium, the rate of emission of the plasma waves is exactly balanced by
their energy loss by Landau damping), this is never the case for its time
aliases components $m\neq 0$ and its spatial aliases $p\neq 0$, $q\neq 0$ or
$r\neq 0$. Worst, they may even be at the origin of numerical instabilities
[Langdon, 1970b] [Birdsall and Langdon, 1991]. [Langdon, 1970b] [Birdsall and
Langdon, 1991] have shown that, in general, the ”collision tensor” (83) has
the following properties

$\displaystyle\left\\{\begin{array}[]{l}\forall\mathbf{X},\,\mathbf{X}\cdot\mathbf{Q}\left(\mathbf{v}_{a},\,\mathbf{v}_{b}\right)\cdot\mathbf{X}<0\cr\mathbf{Q}\left(\mathbf{v}_{a},\,\mathbf{v}_{b}\right)=\mathbf{Q}\left(\mathbf{v}_{b},\,\mathbf{v}_{a}\right)\end{array}\right..$
(87)

According to the kinetic equation (86), a consequence is that the average PIC
simulation plasma entropy

$\bar{\mathrm{H}}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=-\displaystyle\sum_{a}\displaystyle\int_{\mathbb{R}^{3}}\bar{f}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{a},\,t_{n}\right)\ln{\bar{f}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{a},\,t_{n}\right)}d^{3}\mathbf{v}_{a},$
(88)

following average macroparticle straight line trajectories in space,
necessarily increases

${\displaystyle\left.\displaystyle\frac{d\bar{\mathrm{H}}}{dt}\displaystyle\right|}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)>0.$
(89)

We conclude that whatever the velocity distributions
$F_{a0}\left(\mathbf{v}\right)$ that are used to initialize macroparticle
velocities using a random number generator, if one lets the simulation running
during a sufficiently long time duration, the PIC simulation plasma will
necessarily tend on average to a quasi equilibrium state that maximizes
entropy.

##  4 Discussion

###  4.1 Neglect of numerical consistency errors compared to statistical
fluctuations

In this subsection, we discuss the simulation conditions for which Assumption
1 is valid. First, (73) has to be compared with the numerical consistency
errors associated with the leap-frog scheme (55) and with the Maxwell solver.
For Maxwell spectral solvers, there is no numerical consistency error in
space. The counterpart is the compulsory use of periodic boundary conditions.
While it is not problematic here because we are considering an infinite and
homogeneous plasma, it may become constraining under other circumstances. For
the 2nd order FDTD scheme, the numerical consistency error in space is given
by

$\mathbf{\boldsymbol{\varepsilon}}_{\mathbf{r}}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\,\mathrm{with}\,\varepsilon_{\mathbf{r},\xi}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=\displaystyle\frac{{\Delta_{\xi}}^{3}}{24}{\displaystyle\left.\displaystyle\frac{\partial^{3}E_{\xi}}{\partial\xi^{3}}\right|}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)+O\left({\Delta_{\xi}}^{5}\right).$
(90)

The numerical consistency error on macroparticle locations due to the leap-
frog scheme is

$\displaystyle\frac{{\Delta_{t}}^{3}}{24}{\displaystyle\left.\displaystyle\frac{d^{3}\mathbf{r}_{a,\ell}}{dt^{3}}\right|}\displaystyle\left(t_{n}\right)+O\left({\Delta_{t}}^{5}\right).$
(91)

Therefore, injecting it in Maxwell equations, we can estimate the consequent
numerical consistency error on the electrostatic field as

$\mathbf{\boldsymbol{\varepsilon}}_{t}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=4\pi\displaystyle\frac{{\Delta_{t}}^{3}}{24}\displaystyle\sum_{a}Q_{a}\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\sum_{\ell=1}^{N_{a}}{\displaystyle\left.\displaystyle\frac{\partial\mathbf{E}_{s}}{\partial
t}\right|}\displaystyle\left(\mathbf{r}_{a,\ell}\displaystyle\left(t_{n}\right),\,t_{n}\right)S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\ell}\displaystyle\left(t_{n}\right)\right)+O\left({\Delta_{t}}^{5}\right).$
(92)

By normalizing time, space, electrical charge and mass according to
$\underline{t}=\omega_{p_{e}}t$,
$\underline{\mathbf{r}}=\mathbf{r}/\lambda_{D}$ where
$\lambda_{D}=v_{T_{e}}/\omega_{p_{e}}$ is the electron Debye screening length,
$\underline{q_{a}}=q_{a}/e$ and $\underline{m_{a}}=m_{a}/m_{e}$, respectively,
and by assuming reasonably $Zm_{e}/m_{i}\ll 1$, we find the scalings

$\displaystyle\left|\underline{\mathbf{E}}_{\mathrm{fluct}}\right|\underset{\underline{t_{n}}\gg
1}{\sim}\displaystyle\sqrt{\displaystyle\frac{\delta
N}{4\pi\bar{n}_{e}{\lambda_{D}}^{3}}\eta\left(2\pi/\underline{\Delta_{x}}\right)},\,\displaystyle\left|\underline{\mathbf{\boldsymbol{\varepsilon}}_{\mathbf{r}}}\right|\leq\alpha\displaystyle\frac{{\underline{\Delta_{x}}}^{3}}{24}\mathrm{and}\displaystyle\left|\underline{\mathbf{\boldsymbol{\varepsilon}}_{t}}\right|\leq\beta\displaystyle\frac{{\underline{\Delta_{t}}}^{3}}{24}$
(93)

where

$\alpha\left(t_{n}\right)=\displaystyle\max_{i,j,k}{\displaystyle\left\\{\displaystyle\sqrt{\displaystyle\sum_{\underline{\xi}=\underline{x},\underline{y},\underline{z}}{\displaystyle\left({\displaystyle\left.\displaystyle\frac{\partial^{3}\underline{E}_{\xi}}{\partial\underline{\xi}^{3}}\right|}\displaystyle\left(\underline{\mathbf{r}}_{i,j,k},\,\underline{t}_{n}\right)\right)}^{2}}\right\\}}$

and

$\beta\left(t_{n}\right)=\displaystyle\max_{i,j,k}{\displaystyle\left\\{\displaystyle\sqrt{\displaystyle\sum_{\underline{\xi}=\underline{x},\underline{y},\underline{z}}{\displaystyle\left({\displaystyle\left.\displaystyle\frac{\partial\underline{E}_{s,\xi}}{\partial\underline{t}}\right|}\displaystyle\left(\underline{\mathbf{r}}_{i,j,k},\,\underline{t}_{n}\right)\right)}^{2}}\right\\}}\displaystyle\max_{i,j,k}{\displaystyle\left\\{\displaystyle\frac{n_{e}\displaystyle\left(\underline{\mathbf{r}}_{i,j,k},\,\underline{t}_{n}\right)}{\bar{n}_{e}}\right\\}}.$

Assumption 1 is therefore valid if
$\displaystyle\left|\underline{\mathbf{E}}_{\mathrm{fluct}}\right|\gg\displaystyle\left|\underline{\mathbf{\boldsymbol{\varepsilon}}_{\mathbf{r}}}\right|$
and
$\displaystyle\left|\underline{\mathbf{E}}_{\mathrm{fluct}}\right|\gg\displaystyle\left|\underline{\mathbf{\boldsymbol{\varepsilon}}_{t}}\right|$.
Considering roughly $\alpha\sim\beta\sim 1$ and
$\underline{\Delta_{t}}\sim\underline{\Delta_{x}}$, this is the case whenever

$\displaystyle\frac{{\underline{\Delta_{x}}}^{6}}{\eta\left(2\pi/\underline{\Delta_{x}}\right)24^{2}}\ll\displaystyle\frac{\delta
N}{4\pi\bar{n}_{e}{\lambda_{D}}^{3}}$ (94)

which is not constraining in many cases. It means that statistical
fluctuations plays a crucial role in PIC simulations using a random number
generator.

###  4.2 Stationarity of the plasma distribution function in PIC simulations

In this subsection, we discuss the simulation conditions for which Assumption
4 is valid. It seems paradoxical that the right-hand side of kinetic equations
(86) that we have obtained according to Assumptions 1, 2, 3 and 4 never
cancels while their left-hand side cancels according to Assumptions 1, 2 and
4. The paradox is resolved if Assumption 4 remains valid only on a time scale
$\omega_{p}t_{n}$ during which the right-hand side of kinetic equations (86)
is sufficiently small thus ensuring the distribution functions (32) remain
homogeneous and stationary. Above a certain simulation time duration much
greater than $\omega_{p}t_{n}$, the distribution functions begin to evolve
with time. To estimate this simulation time threshold, let us consider the 3D
electrostatic spectral PIC simulation of a plasma consisting of immobile ions
and thermal electrons that uses equal weights $\delta N_{e}=\delta
N_{i}=\delta N$ and a smoothing function filtering all spatial frequency
aliases. We also neglect time aliasing effects assuming a sufficiently small
time step $\Delta_{t}$ and we assume $\Delta_{x}=\Delta_{y}=\Delta_{z}=\Delta$
and $a_{x}=a_{y}=a_{z}=a\leq\Delta$ such that
$\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)\approx\exp{\left(-a^{2}\mathbf{k}^{2}/2\right)}$.
In this case, the PIC Coulomb logarithm, defined according to

$\ln{\Lambda\left(\mathbf{v}_{e}\right)}=\displaystyle\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{4}dk}{k\,{\varepsilon_{L}\left(\mathbf{k}\cdot\mathbf{v}_{e},\,\mathbf{k}\right)}^{2}},$
(95)

only depends on macroelectron velocities $\mathbf{v}_{e}$ if they have a
velocity much greater than the thermal velocity
$\left|\mathbf{v}_{e}\right|\gg v_{T_{e}}$, implying the PIC simulation plasma
permittivity is close to unity
$\varepsilon_{L}\left(\mathbf{k}\cdot\mathbf{v}_{e},\,\mathbf{k}\right)\approx
1$. Here, we have noted $k=\left|\mathbf{k}\right|$ and
$\mathbf{\boldsymbol{\kappa}}=\mathbf{k}/k={\left(\cos{\theta},\,\sin{\theta}\cos{\varphi},\,\sin{\theta}\sin{\varphi}\right)}^{t}$.
We have also introduced arbitrary cutoffs $k_{\mathrm{min}}$ and
$k_{\mathrm{max}}$ that are not necessary for the convergence of the Coulomb
logarithm contrary to the one obtained for point particles. In this particular
case, the Lenard-Balescu-like tensor (83) reduces to the Landau collision
tensor

$\begin{array}[]{lll}\mathbf{Q}\left(\mathbf{v}_{e},\,\mathbf{v}_{b}\right)&=&-2\displaystyle\frac{{e}^{2}{q_{b}}^{2}\ln{\Lambda\left(\mathbf{v}_{e}\right)}}{\left|\mathbf{v}_{e}-\mathbf{v}_{b}\right|}\displaystyle\int_{0}^{2\pi}d\varphi\displaystyle\int_{0}^{\pi}\sin{\theta}d\theta\,\mathbf{\boldsymbol{\kappa}}\otimes\mathbf{\boldsymbol{\kappa}}\,\delta\displaystyle\left(\mathbf{\boldsymbol{\kappa}}\cdot\displaystyle\frac{\mathbf{v}_{e}-\mathbf{v}_{b}}{\left|\mathbf{v}_{e}-\mathbf{v}_{b}\right|}\right)\cr&=&-2\displaystyle\frac{{e}^{2}{q_{b}}^{2}\ln{\Lambda\left(\mathbf{v}_{e}\right)}}{\left|\mathbf{v}_{e}-\mathbf{v}_{b}\right|}\pi\displaystyle\left[\mathbf{I}-\displaystyle\frac{\left(\mathbf{v}_{e}-\mathbf{v}_{b}\right)\otimes\left(\mathbf{v}_{e}-\mathbf{v}_{b}\right)}{{\left|\mathbf{v}_{e}-\mathbf{v}_{b}\right|}^{2}}\right]\end{array}$
(96)

that describes small angle binary Coulomb collisions between point particles.
It seems very strange to obtain such a collision integral since we have
neglected close-encounter Coulomb collisions between macroparticles and
consequently between real particles. In the textbook [Landau and Lifshitz,
1981], they explain ”it is sufficient to consider fluctuations in a
collisionless plasma to calculate the collision integral because the important
Fourier components of the electric field in collisions in plasmas are those
with $\left|\mathbf{k}\right|\gtrapprox 1/\lambda_{D}\gg
1/\lambda_{\mathrm{mfp}}$, so that collisions may be neglected” where
$\lambda_{\mathrm{mfp}}$ denotes real particle mean free paths due to close-
encounter binary Coulomb collisions. As a consequence, ”collisions” are
naturally taken into account in electrostatic PIC simulations of collisionless
plasmas using a random number generator to initialize macroparticle
velocities. Knowing from the kinetic theory of plasmas that the smallest
relaxation rate is given by the electron-ion collision frequency, let us
neglect in a first attempt the electron-electron fluctuations correlation term
in the kinetic equation (86). Considering then the limit $Zm_{e}/m_{i}\ll 1$,
one can neglect the friction experienced by such electrons compared to their
diffusion so that the right hand-side of the kinetic equation (86) reduces to
a Lorentz-like collision operator whose time scale is governed by the
overestimated electron-ion collision frequency

$\nu_{ei}\left(\mathbf{v}_{e}\right)=\delta
N\displaystyle\frac{4\pi\bar{n}_{i}Z^{2}e^{4}}{{m_{e}}^{2}{\left|\mathbf{v}_{e}\right|}^{3}}\ln{\Lambda}\left(\mathbf{v}_{e}\right)\lesssim\bar{\nu}_{ei}=\delta
N\displaystyle\frac{4\pi\bar{n}_{e}Ze^{4}}{{m_{e}}^{2}{v_{T_{e}}}^{3}}\ln{\Lambda}\left(\mathbf{v}_{T_{e}}\right).$
(97)

We conclude that Assumption 4 remains valid on a time scale $t_{n}$ such that

$1\ll\omega_{p}t_{n}\ll\displaystyle\frac{\omega_{p}}{\bar{\nu}_{ei}}=\displaystyle\frac{4\pi\bar{n}_{e}{\lambda_{D}}^{3}}{Z\delta
N\ln{\Lambda}\left(\mathbf{v}_{T_{e}}\right)}\approx\displaystyle\frac{\lambda_{\mathrm{mfp}}}{\delta
N\lambda_{D}}$ (98)

that is to say for not too long electrostatic PIC simulations of ”not too
dense” collisionless plasmas; ”not too dense” referring to the value of
macroparticle weight factors $\delta N$ constrained by today’s computer
technology. Although overestimating them by the weight factor $\delta N$, the
term $\left(m,\,p,\,q,\,r\right)=\left(0,\,0,\,0,\,0\right)$ of (82) describes
single-time correlations between the positions of real plasma particles due to
their long-range Coulomb interactions. It is consequently not a numerical
noise contrary to what is usually considered. Contrary to ”not too dense”
collisionless plasmas, the distribution functions evolve faster than the
statistical fluctuations correlation in electrostatic PIC simulations of
”dense” collisionless plasmas or collisional plasmas for which
$\omega_{p}/\bar{\nu}_{ei}\lessapprox 1$. Assumption 4. is therefore not valid
in this case. In addition, the physics of such fast evolution is not taken
into account by electrostatic PIC codes because they underestimate close-
encounter Coulomb collisions in this regime where $1/\lambda_{D}\lessapprox
1/\delta N\lambda_{\mathrm{mfp}}<1/\lambda_{\mathrm{mfp}}$. To model correctly
this fast evolution of distribution functions, a Langevin force based on a
Monte-Carlo model of multiple binary Coulomb collisions between real particles
per time step $\Delta_{t}$ should be added into (1) such as the one proposed
by [Nanbu, 1997]. Such a Monte-Carlo algorithm emulates statistical
fluctuations due to binary Coulomb collisions between particles and ensures
Bogoliubov hypothesis. However, these PIC simulations are out of the scope of
this paper.

In order to give orders of magnitude, let us consider a laser-generated Helium
plasma ($Z=2$), heated to temperatures up to $T_{e}\approx 1\mathrm{keV}$ and
where the electron density is about $\bar{n}_{e}\approx 5\times
10^{19}\mathrm{cm}^{-3}$ ($\omega_{p}\approx 4\times
10^{14}\mathrm{rad.s}^{-1}$), the Debye length is about $\lambda_{D}\approx
0.03\,\mu\mathrm{m}$. The resulting number of electrons in a Debye cube is
about $\bar{n}_{e}{\lambda_{D}}^{3}\approx 1776$. If we choose reasonably
$\Delta=0.01\,\mu\mathrm{m}\approx\lambda_{D}/3$, we have necessarily
$N_{e,\mathrm{mpc}}{\left(\lambda_{D}/\Delta\right)}^{\mathrm{3}}\gtrapprox
266$ for the assumption 2. to be valid. It leads to a maximum possible value
of macroelectron weights $\delta N_{e}\lessapprox 7$ in 3D simulations for
matching the correct plasma electron density. We deduce Assumption 4. remains
valid on a time scale $\omega_{p}t_{n}\ll 1670$. Under such simulation
conditions, assumption 1. is valid if one chooses
$\omega_{p}\Delta_{t}\lessapprox\Delta/\lambda_{D}$ to ensure the numerical
stability of the simulation. Indeed,
$\underline{\Delta}^{6}/\eta\left(2\pi/\underline{\Delta}\right)24^{2}\approx
1\times 10^{-15}\eta\left(2\pi/\underline{\Delta}\right)\ll\delta
N_{e}/4\pi\bar{n}_{e}{\lambda_{D}}^{3}\approx 3\times 10^{-5}$ in this case.
Considering a reasonable number of macroparticles simulated with today’s
computer technology $N_{e}\approx 1\times 10^{11}$ using a highly parallelized
PIC code and several $10^{5}$ CPU $\times$ hours with a computational cost
typical of $\approx 100\,\mathrm{ns}/\mathrm{timestep}/\mathrm{macroparticle}$
and a simulation time duration $L_{t}=1000/\omega_{p}$, we conclude that we
are limited to plasma simulation boxes of $L_{x}=L_{y}=L_{z}\approx
40\,\mu\mathrm{m}$ by using such a macroelectron weight $\delta N_{e}\approx
7$.

Performing the same reasoning on the main term considering the same
assumptions for 2D electrostatic PIC simulations but defining a Coulomb
logarithm

$\displaystyle\frac{\ln{\Lambda^{\mathrm{2D}}\left(\mathbf{v}_{e}\right)}}{\lambda_{D}}=\displaystyle\int_{k_{\mathrm{min}}}^{k_{\mathrm{max}}}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{4}dk}{{\varepsilon_{L}\left(\mathbf{k}\cdot\mathbf{v}_{e},\,\mathbf{k}\right)}^{2}}$
(99)

instead, we find Assumption 4 remains valid on a time scale $t_{n}$ such that

$1\ll\omega_{p}t_{n}\ll\displaystyle\frac{\omega_{p}}{\bar{\nu}_{ei}^{\mathrm{2D}}}=\displaystyle\frac{2\pi\bar{n}_{e}{\lambda_{D}}^{2}}{Z\delta
N\ln{\Lambda^{\mathrm{2D}}\left(\mathbf{v}_{T_{e}}\right)}}$ (100)

in agreement with [Montgomery and Nielson, 1970], [Hockney, 1971] and [Hockney
and Eastwood, 1981]. For 1D electrostatic PIC simulations, one finds the main
term $\left(m,\,p\right)=\left(0,\,0\right)$ is identically zero owing to the
Dirac delta function similarly as for the 1D Lenard-Balescu or the 1D Landau
collision terms in agreement with [Dawson, 1964]. Indeed, according to the
latter, the observed relaxation time

$\displaystyle\frac{\omega_{p}}{\bar{\nu}^{\mathrm{1D}}_{\mathrm{Dawson}}}=\displaystyle\frac{1}{10}\,{\displaystyle\left(\displaystyle\frac{\bar{n}_{e}{\lambda_{D}}}{\delta
N}\right)}^{2}$ (101)

of 1D electrostatic PIC simulation plasmas is due to the 3-body collisions
between macroparticles and we should have considered the fluctuations
equations (29) and (30) without neglecting the quadratic terms (Assumption 3)
to obtain such a result. However, the aliasing terms do not cancel in
agreement with [Birdsall and Langdon, 1991] and [Montgomery and Nielson,
1970]. One finds for example a relaxation time

$\displaystyle\frac{\omega_{p}}{\bar{\nu}_{(0,1)}^{\mathrm{1D}}}=\displaystyle\frac{2{(\mathbf{k}_{g}\lambda_{D})}^{2}V}{v_{T_{e}}}\displaystyle\frac{\bar{n}_{e}{\lambda_{D}}}{\delta
N}$ (102)

where

$\displaystyle\frac{1}{V}=\displaystyle\frac{1}{V_{e}}+\displaystyle\frac{Zm_{e}/m_{i}}{V_{i}}\,\mathrm{and}\,\displaystyle\frac{1}{V_{b}}=\displaystyle\int_{-\infty}^{\infty}{\displaystyle\left(\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\displaystyle\frac{-k_{g}v^{\prime}_{b,x}}{v_{e,x}-v^{\prime}_{b,x}}\right)}{\overset{\frown}{\mathrm{S}}\left(\displaystyle\frac{-k_{g}v_{e,x}}{v_{e,x}-v^{\prime}_{b,x}}\right)}}{\varepsilon_{L}\displaystyle\left(\displaystyle\frac{-k_{g}v_{e,x}v^{\prime}_{b,x}}{v_{e,x}-v^{\prime}_{b,x}},\,\displaystyle\frac{-k_{g}v^{\prime}_{b,x}}{v_{e,x}-v^{\prime}_{b,x}}\right)}\right)}^{2}{\displaystyle\left(1-\displaystyle\frac{v_{e,x}}{v^{\prime}_{b,x}}\right)}^{2}\displaystyle\frac{F_{b_{0}}(v^{\prime}_{b,x})}{v^{\prime}_{b,x}}dv^{\prime}_{b,x}$

considering only the spatial aliasing term
$\left(m,\,p\right)=\left(0,\,1\right)$ when deriving the collision operator
(82) that reads

$-\displaystyle\frac{q_{e}}{m_{e}}\displaystyle\frac{\partial}{\partial\mathbf{v}_{e}}\cdot\mathbb{E}\displaystyle\left\\{{\delta\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}\delta
f_{e_{c}}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\right\\}\underset{\omega_{p}t_{n}\gg
1}{\approx}\displaystyle\frac{\partial}{\partial
v_{e,x}}\displaystyle\left({\nu}_{(0,1)}^{\mathrm{1D}}{v_{T_{e}}}^{2}\displaystyle\frac{\partial\bar{f_{e_{c}}}}{\partial
v_{e,x}}\right)$ (103)

in agreement with [Birdsall and Langdon, 1991] and the empirical diffusion
term proposed by [Virtamo and Tuomisto, 1979]. As a conclusion for 1D
electrostatic PIC simulations, we conclude Assumption 4 remains valid on a
time scale $t_{n}$ such that

$1\ll\omega_{p}t_{n}\ll\displaystyle\frac{\omega_{p}}{\bar{\nu}_{ei}^{\mathrm{1D}}}$
(104)

with $\bar{\nu}^{\mathrm{1D}}$ scaling as (102) for a large number of
macroparticle centers per Debye length and as (101) for a small number of
macroparticle centers per Debye length as already suggested by [Birdsall and
Langdon, 1991].

###  4.3 Evolution of the empirical entropy in electrostatic PIC simulation
plasmas

Landau damping of electrostatic waves is a time-reversible process [Landau,
1946]. [Decyk, 1987] have shown that the emission of electrostatic wakefield
is also a time-reversible process. Therefore, it may seem paradoxical that we
concluded that whatever velocity distributions $F_{a0}\left(\mathbf{v}\right)$
that are used to initialize macroparticle velocities using a random number
generator, if one lets the simulation running during a sufficiently long time
duration, the PIC-simulated plasma will necessarily tend on average to a
quasi-equilibrium state that maximizes the entropy. The physical reason for
this entropy increase is the average distribution function entropy increase:
any deterministic and time-reversible statistical fluctuations around the
average solution makes tend the average solution to be the one that maximizes
the entropy. Consequently, even if they are time-reversible, the statistical
realizations also tend to the solution that maximizes entropy since they are
supposed to be close to the average result in the limit $N_{a,\mathrm{mpc}}\gg
1$. Let us emphasize here that the macroscopic PIC-simulated plasma entropy
(88) is not the empirical PIC-simulated plasma entropy

$\mathrm{H}_{\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\,t_{n}\right)=-\displaystyle\sum_{a}\displaystyle\sum_{\mathrm{l,m,n}}f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{\mathrm{l,m,n}},\,t_{n}\right)\ln{f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{\mathrm{l,m,n}},\,t_{n}\right)}{\Delta_{v}}^{3}$
(105)

that one can compute directly from the PIC simulation by discretizing the
velocity space according to bins of volume ${\Delta_{v}}^{3}$ located between
velocity grid points
$\mathbf{v}_{l,m,n}=\mathbf{I}_{l,m,n}\cdot\mathbf{\boldsymbol{\Delta}}_{v}$
where
$\mathbf{\boldsymbol{\Delta}}_{v}={\displaystyle\left(\Delta_{v},\,\Delta_{v},\,\Delta_{v}\right)}^{t}$
and by counting the number of macroparticles contained inside each bins

$\begin{array}[]{lll}f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{\mathrm{l,m,n}},t_{n}\right)&=&\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{r}f_{a,c}\displaystyle\left(\mathbf{r},\,\mathbf{v},\,t_{n}\right)S\left(\mathbf{r}-\mathbf{r}_{i,j,k}\right)\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}-\mathbf{v}_{\mathrm{l,m,n}}\right)\cr&=&\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}\hskip
5.0ptf_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}-\mathbf{v}_{\mathrm{l,m,n}}\right)\cr&=&\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}S\displaystyle\left(\mathbf{r}_{i,j,k}-\mathbf{r}_{a,\,\ell}\displaystyle\left(t_{n}\right)\right)\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}_{\mathrm{l,m,n}}-\mathbf{v}_{a,\ell}\left(t_{n}\right)\right).\end{array}$
(106)

Indeed, by differentiating the latter empirical distribution functions
$f_{a,\Pi}\displaystyle\left(\mathbf{r},\mathbf{v},t\right)$ with time, one
obtains the kinetic equation

$\displaystyle\frac{\partial f_{a,\Pi}}{\partial
t}+\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}f_{a,\Pi}\right)+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\mathbf{E}_{s}\left(\mathbf{r},\,t\right)f_{a,\Pi}\right)=C_{a}^{\prime}\left(\mathbf{r},\,\mathbf{v},\,t\right)$
(107)

where

$C_{a}^{\prime}\left(\mathbf{r},\,\mathbf{v},\,t\right)=C_{a}\left(\mathbf{r},\,\mathbf{v},\,t\right)-\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\delta
N_{a}\displaystyle\sum_{\ell=1}^{N_{a}}\displaystyle\left[\mathbf{v}_{a,\ell}\displaystyle\left(t\right)-\mathbf{v}\right]S\displaystyle\left(\mathbf{r}-\mathbf{r}_{a,\,\ell}\displaystyle\left(t\right)\right)\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}-\mathbf{v}_{a,\ell}\left(t\right)\right)\right)$
(108)

according to the continuous macroparticle equations of motion. If
$\Delta_{v_{x}}=\Delta_{v_{y}}=\Delta_{v_{z}}=\Delta_{v}$ is chosen
sufficiently small such that there is a maximum of one macroparticle per bin
volume, meaning
$\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}-\mathbf{v}_{l,m,n}\right)\underset{\Delta_{v}\rightarrow
0}{\rightarrow}\delta\left(\mathbf{v}-\mathbf{v}_{l,m,n}\right)$ and
$f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{l,m,n},t_{n}\right)\underset{\Delta_{v}\rightarrow
0}{\rightarrow}f_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{l,m,n},t_{n}\right)$,
the discretized empirical entropy (105) is only determined by the right-hand
side (17) of Klimontovich-like equations (15) describing the effects of
macroparticles internal tension forces. However, in practice, we do not have
sufficient computer memory and there is always more than one macroparticle per
velocity bin volume. As a consequence, we lose necessarily information about
the simulation results and the second term in the right-hand side (108) of
equation (107) cannot be neglected. Finally, let us stress here that, if one
does not use a random number generator to initialize macroparticle velocities
$\mathbf{V}_{a,\ell}$ at the simulation start $t_{1}$ by initializing, for
example, a cold plasma where $\forall
a,\,\ell,\,\mathbf{V}_{a,\ell}=\mathbf{0}$ instead, all equations are fully
deterministic and there is no statistical fluctuations. For all PIC-simulated
plasma quantities $F$, $\mathbb{E}{\left\\{F\right\\}}=F$ and the kinetic
equations (33) become meaningless. Only the deterministic description (15)
applies in this case.

##  5 The ensemble averaging technique

###  5.1 Reduction of statistical fluctuations via ensemble averaging

Let us consider $N_{\text{\Large\textcalligra{s}}}\,$ PIC simulations of the
same stationary, homogeneous, infinite and fully ionized plasma that we have
studied in the previous section 3, the simulations differing only by the
randomly generated numbers used to initialize macroparticle velocities. We
note for all PIC simulation plasma quantity
$F\in\\{\delta\left(\mathbf{v}-\mathbf{v}_{a,\ell}\left(t_{n}\right)\right)$,
$\bar{f}_{a_{c}}$, $f_{a_{c}}$, $\delta f_{a_{c}}$, $\bar{n}_{a_{c}}$,
$n_{a_{c}}$, $\delta n_{a_{c}}$, $\bar{f}_{a}$, $f_{a}$, $\delta f_{a}$,
$\bar{n}_{a}$, $n_{a}$, $\delta n_{a}$, $\bar{\mathbf{E}}$, $\mathbf{E}$,
$\delta\mathbf{E}$, $\bar{\mathbf{E}}_{s}$, $\mathbf{E}_{s}$ ,
$\delta\mathbf{E}_{s}\\}$,

$\displaystyle\langle
F\rangle=\displaystyle\frac{1}{N_{\text{\Large\textcalligra{s}}}\,}\displaystyle\sum_{\text{\Large\textcalligra{s}}\,=1}^{N_{\text{\Large\textcalligra{s}}}\,}F^{(\text{\Large\textcalligra{s}}\,)}\underset{N_{\text{\Large\textcalligra{s}}}\,\rightarrow\infty}{\rightarrow}\mathbb{E}\displaystyle\left\\{F\right\\}$
(109)

its finite ensemble average obtained from all PIC simulation results, each one
denoted by the superscript
$(\text{\Large\textcalligra{s}}\,),\,\text{\Large\textcalligra{s}}\,\in\left[1,\,N_{\text{\Large\textcalligra{s}}}\,\right]$.
By linearity, the expected value of the ensemble average of PIC simulation
plasma quantities is also their expected value

$\mathbb{E}\displaystyle\left\\{\displaystyle\langle
F\rangle\right\\}=\mathbb{E}\displaystyle\left\\{F\right\\}$ (110)

during one single PIC simulation. We thus have in particular

$\mathbb{E}\displaystyle\left\\{\displaystyle\langle
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle\right\\}=\bar{f}_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)$
(111)

and

$\mathbb{E}\displaystyle\left\\{\displaystyle\langle\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle\right\\}=\bar{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right).$
(112)

As a consequence, the deterministic description of the finite ensemble average
of PIC-simulations consists in the kinetic equation

$\begin{array}[]{lll}&&{\displaystyle\left.\displaystyle\frac{\partial\langle
f_{a}\rangle}{\partial
t}\right|}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)+{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\cdot\displaystyle\left(\mathbf{v}\langle
f_{a}\rangle\right)\right|}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)+\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left(\displaystyle\frac{q_{a}}{m_{a}}\langle\mathbf{E}_{s}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle\langle
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle\right)\cr&=&-\displaystyle\frac{q_{a}}{m_{a}}\displaystyle\frac{\partial}{\partial\mathbf{v}}\cdot\displaystyle\left[\displaystyle\langle{\delta\mathbf{E}_{s}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\delta
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle-\langle{\delta\mathbf{E}_{s}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle\langle\delta
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle\right]\cr&&+\displaystyle\langle
C\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle\end{array}$
(113)

self-consistently coupled with the finite ensemble average of discretized
Maxwell equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\langle{\mathbf{E}\rangle}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\langle{f_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle}d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\langle{\mathbf{E}\rangle}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (114)

according to (15), (16) and (109). Because all PIC simulations use different
randomly generated numbers to initialize macroparticle velocities, they are
all statistically independent. We thus deduce the relationship between the
finite ensemble average of PIC simulation plasma quantity variances with their
own variances from one single PIC simulation. For all PIC simulation plasma
quantity $G\in\left\\{\bar{n}_{a_{c}},\,n_{a_{c}},\,\delta
n_{a_{c}},\,\bar{n}_{a},\,\,n_{a},\,\delta
n_{a},\,\bar{\mathbf{E}},\,\mathbf{E},\,\delta\mathbf{E},\,\bar{\mathbf{E}}_{s},\,\mathbf{E}_{s},\,\delta\mathbf{E}_{s}\right\\}$,
we have

$\mathbb{E}\Big{\\{}{\Big{(}\displaystyle\langle
G\rangle-\mathbb{E}\displaystyle\left\\{\displaystyle\langle
G\rangle\right\\}\Big{)}}^{2}\Big{\\}}=\displaystyle\frac{1}{N_{\text{\Large\textcalligra{s}}}\,}\mathbb{E}\Big{\\{}{\Big{(}G-\mathbb{E}\displaystyle\left\\{G\right\\}\Big{)}}^{2}\Big{\\}}.$
(115)

We thus obtain in particular

$\mathbb{E}\Big{\\{}{\displaystyle\langle\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle}^{2}\Big{\\}}=\displaystyle\frac{\mathbb{E}\Big{\\{}{\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}^{2}\Big{\\}}}{N_{\text{\Large\textcalligra{s}}}\,}.$
(116)

If $N_{\text{\Large\textcalligra{s}}}\,$ is sufficiently large, the finite
ensemble average of electrostatic field statistical fluctuation
$\langle\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,,t_{n}\right)\rangle$
at grid points $\mathbf{r}_{i,j,k}$ and time steps $t_{n}$ follows a normal
distribution law with mean (112) and variance (116) according to the central
limit theorem. One deduces

$\displaystyle\left|\displaystyle\langle\delta\mathbf{E}\displaystyle\left(\mathbf{r}_{i,j,k},\,,t_{n}\right)\rangle\rangle\right|\underset{N_{\text{\Large\textcalligra{s}}}\,\gg
1}{\leq}\displaystyle\sqrt{2}\,\text{erf}^{-1}\displaystyle\left(\displaystyle\frac{p}{100}\right)\displaystyle\left|\displaystyle\langle\mathbf{E}_{\text{fluct}}\rangle\right|$
(117)

with a $p\%$-confidence where we have noted

$\displaystyle\left|\displaystyle\langle\mathbf{E}_{\text{fluct}}\rangle\right|=\displaystyle\frac{\displaystyle\left|\displaystyle\mathbf{E}_{\text{fluct}}\right|}{\sqrt{N_{\text{\Large\textcalligra{s}}}\,}}$
(118)

assuming again Assumptions 1, 2, 3 and 4 as in section 3.2. Considering again
the particular case of PIC simulations of a collisionless plasma at
equilibrium using equal weights $\delta N_{e}=\delta N_{i}=\delta N$ and
$\Delta_{t}<\Delta_{x}/v_{T}=\Delta_{y}/v_{T}=\Delta_{z}/v_{T}<1/\omega_{p}$,
it reads

$\displaystyle\left|\displaystyle\langle\mathbf{E}_{\text{fluct}}\rangle\right|\underset{\omega_{p}t_{n}\gg
1}{\sim}\displaystyle\sqrt{\displaystyle\frac{\delta
N}{N_{\text{\Large\textcalligra{s}}}\,}\displaystyle\frac{k_{B}T}{{\lambda_{D}}^{3}}\eta\left(k_{g}\lambda_{D}\right)}.$
(119)

A finite ensemble average of $N_{\text{\Large\textcalligra{s}}}\,$ PIC
simulations thus reduces statistical fluctuation amplitudes of electrostatic
fields and densities by a factor $\sqrt{N_{\text{\Large\textcalligra{s}}}\,}$
compared to one single electrostatic PIC simulation. Also, since $\delta
N\propto 1/N_{e,\text{mpc}}$ according to its definition (7), we deduce there
is an equivalence between running one PIC simulation using $N_{e,\text{mpc}}$
macroparticles and a finite ensemble average of
$N_{e,\text{mpc}}/N_{e,\text{mpc}}^{\prime}$ PIC simulations using each one
$N_{e,\text{mpc}}^{\prime}$ macroelectrons per cell. We can thus apply the
central limit theorem to estimate the statistical fluctuations in one single
PIC simulations using a sufficiently large number of macroparticle per cell
$N_{e,\text{mpc}}\gg N_{e,\text{mpc}}^{\prime}\geq 2$ in 1D, 4 in 2D and 8 in
3D as we did in section 3.2.

###  5.2 Kinetic equation for an ensemble average of electrostatic PIC
simulations

Similarly, if we note $F^{\prime}$ another PIC simulation plasma quantity
distinct from $F$, the statistical independence of PIC simulations makes the
finite ensemble average of PIC simulation plasma quantity covariance between
$F$ and $F^{\prime}$ decrease by a factor
$\sqrt{N_{\text{\Large\textcalligra{s}}}\,}$ compared to their covariance from
one single PIC simulation

$\mathbb{E}\Big{\\{}{\Big{(}\displaystyle\langle
F\rangle-\mathbb{E}\displaystyle\left\\{\displaystyle\langle
F\rangle\right\\}\Big{)}}{\Big{(}\displaystyle\langle
F^{\prime}\rangle-\mathbb{E}\displaystyle\left\\{\displaystyle\langle
F^{\prime}\rangle\right\\}\Big{)}}\Big{\\}}=\displaystyle\frac{1}{N_{\text{\Large\textcalligra{s}}}\,}\mathbb{E}\Big{\\{}{\Big{(}F-\mathbb{E}\displaystyle\left\\{F\right\\}\Big{)}}{\Big{(}F^{\prime}-\mathbb{E}\displaystyle\left\\{F^{\prime}\right\\}\Big{)}}\Big{\\}}.$
(120)

We thus obtain in particular

$\mathbb{E}\displaystyle\left\\{\displaystyle\langle\delta\widehat{f}_{a_{c}}\left(\mathbf{k},\,\mathbf{v},\,t_{1}\right)\rangle\displaystyle\langle\delta\widehat{f}_{b,c}\left(\mathbf{k}^{\prime},\,\mathbf{v}^{\prime},\,t_{1}\right)\rangle\right\\}=\delta_{ab}\displaystyle\frac{\delta
N_{a}}{N_{\text{\Large\textcalligra{s}}}\,}{\displaystyle\left(2\pi\right)}^{3}\bar{n}_{a}\delta\left(\mathbf{k}+\mathbf{k}^{\prime},\,\mathbf{k}_{g}\right)F_{a0}\displaystyle\left(\mathbf{v}\right)\displaystyle\left[\delta\displaystyle\left(\mathbf{v}-\mathbf{v}^{\prime}\right)-F_{a0}\displaystyle\left(\mathbf{v}^{\prime}\right)\right].$
(121)

Therefore, by defining statistical fluctuations

$\displaystyle\langle\delta
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle=\displaystyle\langle
f_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle-\displaystyle\langle\bar{f}_{a}\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle$
(122)

and

$\displaystyle\langle\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle=\displaystyle\langle\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle-\displaystyle\langle\bar{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle$
(123)

between the finite ensemble average of PIC simulation results and the expected
ones, by performing the same derivation as in section 3.3 using (121) instead
of (70) and assuming again Assumptions 1, 2, 3 and 4, we obtain the
probabilistic description of the finite ensemble average of PIC simulations.
The kinetic equation reads

$\begin{array}[]{l}\hskip
15.00002pt{\displaystyle\left.\displaystyle\frac{\partial\langle{\bar{f}}_{a}\rangle}{\partial{t}}\right|}\left({\mathbf{r}}_{i,j,k},\,{\mathbf{v}}_{a},\,{t}_{n}\right)+{\displaystyle\left.\displaystyle\frac{\partial}{\partial{\mathbf{r}}}\cdot\displaystyle\left({\mathbf{v}}_{a}\langle{\bar{f}}_{a}\rangle\right)\right|}\left({\mathbf{r}}_{i,j,k},\,{\mathbf{v}}_{a},\,{t}_{n}\right)+\displaystyle\frac{\partial}{\partial{\mathbf{v}}_{a}}\cdot\displaystyle\left(\displaystyle\frac{{q}_{a}}{{m}_{a}}\langle{\bar{\mathbf{E}}}_{s}\left({\mathbf{r}}_{i,j,k},\,{t}_{n}\right)\rangle\langle{\bar{f}}_{a}\left({\mathbf{r}}_{i,j,k},\,{\mathbf{v}}_{a},\,{t}_{n}\right)\rangle\right)\cr\underset{\omega_{p}{t}_{n}\gg
1}{=}-\displaystyle\frac{1}{{m}_{a}}\displaystyle\frac{\partial}{\partial{\mathbf{v}}_{a}}\cdot\displaystyle\sum_{b}\displaystyle\int_{\mathbb{R}^{3}}{\mathbf{Q}}\displaystyle\left({\mathbf{v}}_{a},\,{\mathbf{v}}_{b}\right)\cdot{\displaystyle\left.\displaystyle\left(\displaystyle\frac{{\delta
N}_{a}}{N_{\text{\Large\textcalligra{s}}}\,}\displaystyle\frac{\langle{\bar{f}}_{a}\left({\mathbf{v}}_{a}\right)\rangle}{m_{b}}\displaystyle\frac{\partial\langle{\bar{f}}_{b}\rangle}{\partial{\mathbf{v}}_{b}}-\displaystyle\frac{{\delta
N}_{b}}{N_{\text{\Large\textcalligra{s}}}\,}\displaystyle\frac{\langle{\bar{f}}_{b}\left({\mathbf{v}}_{b}\right)\rangle}{m_{a}}\displaystyle\frac{\partial\langle{\bar{f}}_{a}\rangle}{\partial{\mathbf{v}}_{a}}\right)\right|}\left({\mathbf{r}}_{i,j,k},\,{t}_{n}\right)d^{3}{\mathbf{v}}_{b}.\end{array}$
(124)

It is self-consistently coupled with the ensemble average of expected
discretized Maxwell equations

$\displaystyle\left\\{\begin{array}[]{lllll}{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\,\,\cdot\,\langle{\bar{\mathbf{E}}\rangle}\hskip
5.0pt\right|}^{i,j,k,n}&=&4\pi\displaystyle\sum_{a}q_{a}\displaystyle\int_{\mathbb{R}^{3}}\langle{\bar{f}_{a}\displaystyle\left(\mathbf{r}_{i,j,k},\,\mathbf{v},\,t_{n}\right)\rangle}d^{3}\mathbf{v}\cr{\displaystyle\left.\displaystyle\frac{\partial}{\partial\mathbf{r}}\times\langle{\bar{\mathbf{E}}\rangle}\hskip
5.0pt\right|}^{i,j,k,n}&=&\mathbf{0}\end{array}\right..$ (125)

via the ensemble average of expected electrostatic field interpolated at
macroparticle center locations

$\langle\bar{\mathbf{E}}_{s}\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right),\,t_{n}\right)\rangle=\displaystyle\sum_{i,j,k}\langle\bar{\mathbf{E}}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)\rangle
S\displaystyle\left(\mathbf{r}_{a,\ell}\left(t_{n}\right)-\mathbf{r}_{i,j,k}\right)\Delta_{x}\Delta_{y}\Delta_{z}.$
(126)

Consequently, the kinetic equations (124) of the ensemble average of
distribution functions tend to the Vlasov equation in the limit of an infinite
number of simulations $N_{\text{\Large\textcalligra{s}}}\,/\delta
N\rightarrow\infty$. More realistically, the ensemble average of
$N_{\text{\Large\textcalligra{s}}}\,=\delta N$ independent electrostatic PIC
simulations modeling the same macroscopic plasma with macroparticle weights
$\delta N=\delta N_{e}=\delta N_{i}$ is strictly equivalent to the numerical
resolution of the Vlasov-Lenard-Balescu/Maxwell set of equations discretized
in space and time only. Let us mitigate here these two important results (118)
and (124) by reminding the reader that they are only valid assuming
Assumptions 1, 2, 3 and 4, that is, only for simulation time durations $L_{t}$
of collisionless plasmas using sufficiently small cell sizes and time steps
such that

$\displaystyle\frac{{\left(\omega_{p_{e}}\Delta_{t}\right)}^{6}}{{\eta\left(k_{g}\lambda_{D}\right)}24^{2}}\ll\displaystyle\frac{\delta
N/N_{\text{\Large\textcalligra{s}}}\,}{4\pi\bar{n}_{e}{\lambda_{D}}^{3}}\ll\displaystyle\frac{1}{Z\omega_{p_{e}}L_{t}}$
(127)

similarly as we obtained in the previous section 4.

Precautions have been made in the previous paragraphs to distinguish between
PIC simulation plasma quantities $F$, $F^{\prime}$ and $G$. Indeed, the
variance of plasma particle phase-space densities at the grid point
$\mathbf{r}_{i,j,k}$ and velocity $\mathbf{v}$ naturally diverges due to the
contribution of the term expressing the case where there is only one
macroparticle in the phase-space volume located between
$\left(\mathbf{r}_{i,j,k},\,\mathbf{v}\right)$ and
$\left(\mathbf{r}_{i,j,k}+\mathbf{\boldsymbol{\Delta}},\,\mathbf{v}+d^{3}\mathbf{v}\right)$.
To the best of our knowledge, there is no unbiased estimator of distribution
function statistical realizations. However, to troubleshoot this issue, one
may use the biased statistical estimator introduced in section 3.3 for the
empirical distribution function (106) to express

$f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{l,m,n},t_{1}\right)=\displaystyle\frac{\bar{n}_{a}}{N_{a,mpc}}\displaystyle\sum_{\ell_{i,j,k}=1}^{N_{a,mpc}}\Pi_{\mathbf{\boldsymbol{\Delta}}_{v}}^{(0)}\left(\mathbf{v}_{l,m,n}-\mathbf{V}_{a,\ell_{i,j,k}}\right)$
(128)

according to Assumption 2. Assuming
$\Delta_{v_{x}}=\Delta_{v_{y}}=\Delta_{v_{z}}=\Delta_{v}\ll v_{T_{a}}$ for
simplicity, one obtains

$\begin{array}[]{lcl}\bar{f}_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{l,m,n},t_{1}\right)&=&\mathbb{E}{\displaystyle\left\\{f_{a,\Pi}\displaystyle\left(\mathbf{r}_{i,j,k},\mathbf{v}_{l,m,n},t_{1}\right)\right\\}}\cr&\underset{\Delta_{v}\ll
v_{T_{a}}}{=}&\bar{n}_{a}F_{a_{0}}\left(\mathbf{v}_{l,m,n}\right)+\bar{n}_{a}\displaystyle\frac{{\Delta_{v}}^{2}}{24}\displaystyle\sum_{\xi}\displaystyle\frac{d^{2}F_{a_{0}}}{dv_{\xi}^{2}}\left(\mathbf{v}_{l,m,n}\right)+O\left[{\displaystyle\left(\displaystyle\frac{\Delta_{v}}{v_{T_{a}}}\right)}^{4}\right]\end{array}$
(129)

and

$\mathbb{E}\displaystyle\left\\{\delta
f_{a,\Pi}\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{l,m,n},\,t_{1}\right)^{2}\right\\}\underset{\Delta_{v}\ll
v_{T_{a}}}{=}\displaystyle\frac{{\left.\bar{n}_{a}\right.}^{2}}{N_{a,\text{mpc}}}F_{a0}\displaystyle\left(\mathbf{v}_{l,m,n}\right)\displaystyle\left[\displaystyle\frac{1}{{\Delta_{v}}^{3}}-F_{a0}\displaystyle\left(\mathbf{v}_{l,m,n}\right)\right].$
(130)

Therefore, in the limit of large number of macroparticles per cell, the