statistical estimator of particle phase-space densities (128) follows a normal
distribution law with mean (129) and variance (130) according to the central
limit theorem. Assuming the PIC simulation plasma is stable, we deduce in the
limits $\omega_{p}t_{n}\gg 1$, $\Delta_{v}\ll v_{T_{a}}$, $N_{a,\text{mpc}}\gg
1$ and $N_{\text{\Large\textcalligra{s}}}\,\gg 1$

$\displaystyle\left|\displaystyle\delta
f_{a,\Pi}\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{l,m,n},\,t_{n}\right)\right|\leq\bar{n}_{a}\displaystyle\frac{{\Delta_{v}}^{2}}{24}\displaystyle\sum_{\xi}\displaystyle\left|\displaystyle\frac{d^{2}F_{a_{0}}}{dv_{\xi}^{2}}\left(\mathbf{v}_{l,m,n}\right)\right|+\displaystyle\sqrt{\displaystyle\frac{2}{N_{a,\text{mpc}}}}\,\text{erf}^{-1}\displaystyle\left(\displaystyle\frac{p}{100}\right)\bar{n}_{a}\displaystyle\sqrt{\displaystyle\frac{F_{a0}\displaystyle\left(\mathbf{v}_{l,m,n}\right)}{{\Delta_{v}}^{3}}}$
(131)

and

$\displaystyle\left|\displaystyle\langle\delta
f_{a,\Pi}\left(\mathbf{r}_{i,j,k},\,\mathbf{v}_{l,m,n},\,t_{n}\right)\rangle\right|\leq\bar{n}_{a}\displaystyle\frac{{\Delta_{v}}^{2}}{24}\displaystyle\sum_{\xi}\displaystyle\left|\displaystyle\frac{d^{2}F_{a_{0}}}{dv_{\xi}^{2}}\left(\mathbf{v}_{l,m,n}\right)\right|+\displaystyle\sqrt{\displaystyle\frac{2}{N_{a,\text{mpc}}N_{\text{\Large\textcalligra{s}}}\,}}\,\text{erf}^{-1}\displaystyle\left(\displaystyle\frac{p}{100}\right)\bar{n}_{a}\displaystyle\sqrt{\displaystyle\frac{F_{a0}\displaystyle\left(\mathbf{v}_{l,m,n}\right)}{{\Delta_{v}}^{3}}}$
(132)

with a $p\%$-confidence. Statistical fluctuations of phase-space densities
depends consequently on a scaling different from $1/\sqrt{N_{a,\text{mpc}}}$
in single PIC simulations. Indeed, by using velocity bin sizes
$\Delta_{v}\propto v_{T_{a}}/{N_{a,\text{mpc}}}^{n}$ that decrease with
decreasing number of macroparticles per cell, the optimal choice minimizing
the bias is $n=1/\left(4+\text{dim}\right)$ where ”dim” denotes here the PIC
simulation velocity space dimension. It provides instead a scaling
$1/\sqrt{{N_{a,\text{mpc}}}^{4/\left(4+\text{dim}\right)}}$ for single PIC
simulations and $1/\sqrt{{\left(N_{a,\text{mpc}}\times
N_{\text{\Large\textcalligra{s}}}\,\right)}^{4/\left(4+\text{dim}\right)}}$
for a finite ensemble average of $N_{\text{\Large\textcalligra{s}}}\,$ PIC
simulations.

###  5.3 Application to linear Landau damping

To illustrate how to decrease the amplitude of statistical fluctuations in PIC
simulations by ensemble averaging, we consider here the linear Landau damping
of a small amplitude electrostatic wave propagating in the 1D-3V PIC
simulation plasma described in section 2. In order to generate the
electrostatic wave, we perturb the PIC simulation plasma by driving an
electrostatic field $\mathbf{E}_{\text{d}}\left(x_{i},\,t_{n}\leq
0\right)=E_{\text{d}}\cos{\displaystyle\left[k_{0}x_{i}-\omega_{0}(t_{n}+\tau)\right]}\mathbf{e}_{x}$
from $\omega_{p_{e}}t_{n}=-\omega_{p_{e}}\tau=-3.6$ to $\omega_{p_{e}}t_{n}=0$
at a spatial angular frequency $k_{0}\lambda_{D}=0.319$. [Landau, 1946]
predicted that if particles move in a collisionless plasma at a velocity
$v_{x}$ close to the wave phase velocity $v_{\varphi}=\omega/k$ of an
electrostatic wave, they see an almost constant electrostatic field and can
thus interact with the wave. If $v_{x}\lessapprox v_{\varphi}$, trapped
particles earn energy from the wave and if $v_{x}\gtrapprox v_{\varphi}$, they
lose energy to it. Consequently, if one consider a PIC simulation plasma
consisting of immobile ions and thermal electrons with an initial Maxwellian
distribution of velocities (8), more macroelectrons will earn energy from the
wave than the electrostatic wave from macroelectrons since $\partial
F_{e_{0}}/\partial v_{x}\left(v_{x}>0\right)<0$ in this case. As a result, the
driven electrostatic wave is linearly damped at the Landau damping rate
$\gamma_{L}$ in a first-order approximation. If one uses a Gaussian smoothing
function (5) and chooses a sufficiently small time step
$\omega_{p_{e}}\Delta_{t}=0.2$, aliasing effects are small in the expression
of the PIC simulation plasma permittivity (63) and one may deduce from the
resulting dispersion relation
$\varepsilon_{L}\left(\omega-\iota\gamma_{L},\,k\right)=0$

$\displaystyle\frac{\omega\left(k\right)}{\omega_{p_{e}}}\underset{k\lambda_{D}\ll
1}{=}\overset{\frown}{\text{S}}\left(k\right)\displaystyle\left[1+\displaystyle\frac{3}{2}{\left(k\lambda_{D}\right)}^{2}\right]$
(133)

and

$\displaystyle\frac{\gamma_{L}\left(k\right)}{\omega_{p_{e}}}\underset{k\lambda_{D}\ll
1}{=}\displaystyle\sqrt{\displaystyle\frac{\pi}{8}}\displaystyle\frac{\overset{\frown}{\text{S}}\left(k\right)}{{\left(k\lambda_{D}\right)}^{3}}\exp{\displaystyle\left[-\displaystyle\frac{1}{2{\left(k\lambda_{D}\right)}^{2}}-\displaystyle\frac{3}{2}\right]}$
(134)

in the limit $k\lambda_{D}\ll 1$. Actually,
$\omega_{0}=\omega\left(k_{0}\right)\approx
1.18\omega_{p_{e}}\overset{\frown}{\text{S}}\left(k_{0}\right)$ and
$\gamma_{L_{0}}=\gamma_{L}\left(k_{0}\right)\approx 1.96\times
10^{-2}\omega_{p_{e}}\overset{\frown}{\text{S}}\left(k_{0}\right)$ at
$k_{0}\lambda_{D}=0.319$ if one solves numerically the dispersion relation;
cf. the left panel of Figure 3. By using spatial grid spacings
$\Delta_{x}=\lambda_{D}$, a quadratic interpolating function
$\Pi_{\Delta_{x}}^{(2)}$ and a smoothing parameter $a_{x}=0.866667\lambda_{D}$
such that we lose the less physical information from spatial angular frequency
spectra as possible, the macroparticle shape (4) is
$\overset{\frown}{\text{S}}\left(k_{0}\right)\approx 0.903$ at the considered
spatial angular frequency $k_{0}\lambda_{D}=0.319$ and it doesn’t affect too
much the theoretical results.

As macroelectrons damp the wave, they also earn energy. As a result, at the
end of this linear stage more macroelectrons finish with a velocity
$v_{x}\gtrapprox v_{\phi}$ compared to those who finish with a velocity
$v_{x}\lessapprox v_{\phi}$ so that the sign of the resulting distribution of
electron velocities changes and macroelectrons give back the received energy
to the wave. By studying the energy conservation between the electrostatic
wave and one trapped macroelectron, one deduces the non-linear Landau damping
occurs at the bounce angular frequency

$\omega_{b}\left(k\right)=\displaystyle\sqrt{\displaystyle\frac{eE_{0}k}{m_{e}}}$
(135)

where $E_{0}$ is the electrostatic wave amplitude. Consequently, the condition
for observing the linear Landau damping of an electrostatic wave is given by
the scaling $\omega\gg\gamma_{L}\gg\omega_{b}$ and the electrostatic wave
amplitude is consequently constrained according to

$\displaystyle\frac{eE_{0}}{m_{e}\omega_{p_{e}}v_{T_{e}}}\ll\displaystyle\frac{{\left(\gamma_{L}/\omega_{p_{e}}\right)}^{2}}{k\lambda_{D}}\ll\displaystyle\frac{{\left(\omega/\omega_{p_{e}}\right)}^{2}}{k\lambda_{D}}.$
(136)

Plots in Figure 3 stop at the spatial angular frequency of $k\lambda_{D}=1$.
Indeed, above this value, the linear Landau damping rate has the same order of
magnitude compared to the angular frequency of electrostatic wave
oscillations. Therefore, the concept of electrostatic waves is meaningless in
this regime. In the opposite regime where angular spatial frequencies
$k\lambda_{D}\lessapprox 0.2$, the linear Landau damping rate decreases
drastically with decreasing spatial angular frequencies by many order of
magnitudes according to equation (134). In this regime, Landau damping is
negligible and electrostatic waves may consequently propagate almost freely in
the PIC simulation plasma. In the right panel of Figure 3, we compare the
Landau damping rate $\gamma_{L}(k)$ and bounce angular frequencies (135) for
different values of initial perturbing electrostatic field amplitude $E_{0}$
in the interesting regime of Landau damping $0.2\lessapprox k\lambda_{D}\leq
1$. While the observation of non-linear Landau damping of a perturbing
electrostatic wave is relatively easy, the observation of linear Landau
damping is much more challenging due to the smaller initial amplitude of the
perturbing electrostatic field that is needed and that is usually below the
PIC simulation plasma electrostatic field statistical fluctuation amplitudes
(49). In our case, $k_{0}\lambda_{D}=0.319$ leads to a value of
${\left(\gamma_{L}/\omega_{p_{e}}\right)}^{2}/k_{0}\lambda_{D}\approx
1.2\times 10^{-3}$.

Figure 3: (Left) Comparison between electrostatic wave angular frequency
$\omega(k)$ and its Landau damping rate $\gamma_{L}(k)$ versus its spatial
angular frequency $k\lambda_{D}$ obtained using Taylor expansions in the limit
$k\lambda_{D}\ll 1$ (dashed curves) and exact theoretical ones (full curves).
(Right) Comparisons between Landau damping rates $\gamma_{L}$ and bounce
frequencies $\omega_{b}$ versus spatial angular frequency of electrostatic
waves $k\lambda_{D}$ in the interesting regime for different electrostatic
wave amplitudes $E_{0}$.

In a first approximation and neglecting statistical fluctuations, the self-
consistent PIC simulation driven plasma electrostatic wave increases linearly
with time $\mathbf{E}_{x}\left(x_{i},\,t_{n}\leq 0\right)\approx
E_{\text{d}}\alpha(t_{n}+\tau)\sin{\displaystyle\left[k_{0}x_{i}-\omega_{0}(t_{n}+\tau)\right]}\mathbf{e}_{x}$
under the action of the driver where $\alpha$ is the growing rate of the
field. In order to give orders of magnitude, let us consider a laser-generated
Hydrogen plasma ($Z=1$) at equilibrium with an electron density
$\bar{n}_{e}\approx 5.4\times 10^{19}\text{ cm}^{-3}$ consisting of immobile
ions ($m_{e}/m_{i}\ll 1\Rightarrow v_{T_{i}}\ll v_{T_{e}}$) and thermal
electrons at a temperature $T_{e}\approx 7$ keV. For such a 1D collisionless
plasma, the number of electrons in a Debye cube is about
$\bar{n}_{e}{\lambda_{D}}^{3}\approx 32,768$ and the Debye length about
$\lambda_{D}\approx 0.085\,\mu$m. We perform the simulations by using the open
source spectral electrostatic PIC code BEPS developed from the UPIC framework
[Decyk, 2007] available at https://github.com/UCLA-Plasma-Simulation-Group. We
use a macroelectron weight $\delta N_{e}=1$, a number of macroelectrons
$N_{e}=8,388,608$ and a simulation box $L_{x}=256\lambda_{D}$ leading to a
number $N_{e,\text{mpc}}=32,768$ of macroelectrons per cell and
$k_{0}\lambda_{D}=2\pi p_{0}/L_{x}$ with $p_{0}=13$. The choice of a
macroelectron weight $\delta N_{e}=1$ allows us to emphasize again that
physical statistical fluctuations are also present in dedicated experiments
and they are not a numerical noise. The input deck for single simulations can
be found in Appendix D. Given all these simulation parameters, we have

$\displaystyle\frac{4\pi{\left(\omega_{p_{e}}\Delta_{t}\right)}^{6}}{\arctan{\left(k_{g}\lambda_{D}\right)}24^{2}}\approx
9.8\times
10^{-7}\leq\displaystyle\frac{1}{4\pi\bar{n}_{e}{\lambda_{D}}^{3}}\approx
2.4\times 10^{-6}\ll\displaystyle\frac{1}{\omega_{p_{e}}L_{t}}=5\times
10^{-3}$ (137)

with $L_{t}=200/\omega_{p_{e}}$ ($N_{t}=1000$ time iterations). It means
according to the previous section 2 that all our estimates obtained in
sections 3, 4 and 5 are valid for these simulations. In particular,
consistency errors due to the leap-frog scheme are smaller than statistical
fluctuation amplitudes

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

with
$e\displaystyle\left|\mathbf{E}_{\text{fluct}}\right|/m_{e}\omega_{p_{e}}v_{T_{e}}\approx
5\times 10^{-4}$ and a $p\%$-confidence. We choose
$eE_{\text{d}}/m_{e}\omega_{p_{e}}v_{T_{e}}=10^{-3}$. According to
simulations, we measure a growing rate of the field $\alpha\approx
0.243\,\omega_{p_{e}}$; cf. Figure 4.

Figure 4: (Left) Comparison of linearly Landau-damped electrostatic wave
spatial Fourier mode $p_{0}$ from one single electrostatic 1D-3V PIC
simulation (black curve) and an average of 30 simulations (red curve). (Right)
Comparison of linearly Landau-damped electrostatic wave spatial Fourier mode
$p_{0}$ from Landau theory (green curve), from an average of 30 electrostatic
1D-3V PIC simulations (red curve) and from the subtraction of one nominal
simulation without perturbation from one PIC simulation with perturbation
(blue curve).

At $t_{n}=0$, the driver is stopped and we expect the driven traveling
electrostatic wave is Landau-damped according to

$\mathbf{E}_{x}\left(x_{i},\,t_{n}\geq 0\right)\underset{\omega_{b}t\ll
1}{\approx}E_{0}\sin{\displaystyle\left[k_{0}x_{i}-\omega_{0}(t_{n}+\tau)\right]}\exp{\displaystyle\left(-\gamma_{L_{0}}t_{n}\right)}\mathbf{e}_{x}$
(139)

with
$eE_{0}/m_{e}\omega_{p_{e}}v_{T_{e}}=eE_{\text{d}}\alpha\tau/m_{e}\omega_{p_{e}}v_{T_{e}}\approx
8.75\times 10^{-4}$. It is just slightly above the statistical fluctuations
level and leads to a bounce frequency
$\omega_{b_{0}}=\omega_{b}(k_{0})=1.67\times 10^{-2}\omega_{p_{e}}$ slightly
lower than the linear Landau damping rate $\gamma_{L_{0}}$ in agreement with
the simulations; cf Figure 3. The observation of linear Landau damping is
therefore difficult to observe in such electrostatic PIC simulations or
dedicated experiments. A method proposed by [Decyk, 1987] has already been
applied by [Grismayer et al., 2011] for resolving this problem. The
subtraction technique consists in running two identical PIC simulations that
use the same randomly generated numbers to initialize macroparticle velocities
at the simulation start. Then, the nominal simulation without perturbation is
subtracted from the one initialized with the small amplitude perturbing
electrostatic field. Here, we use the ensemble average of
$N_{\text{\Large\textcalligra{s}}}\,=30$ simulations and we compare the
different results. In the left panel of Figure 4, we have plotted the spatial
mode $p_{0}$ of electrostatic field discrete Fourier transform

$\widehat{E}_{x}\left(k_{0},\,t_{n}\right)=\displaystyle\frac{1}{N_{x}}\displaystyle\sum_{i=1}^{N_{x}}E_{x}\left(x_{i},\,t_{n}\right)\exp{\displaystyle\left(-\iota
2\pi p_{0}x_{i}/L_{x}\right)}$ (140)

from one single simulation and from their ensemble average. We can see
statistical fluctuation amplitudes (138) are effectively reduced by a factor
$\sqrt{30}\approx 5.5$ in agreement with our estimate (118). After
$\omega_{p_{e}}t_{n}\approx 20$, one single PIC simulation doesn’t allow for
the observation of the tiny electrostatic field amplitude (136) below
statistical fluctuations (138) while the electrostatic wave amplitude continue
to decrease until $\omega_{p_{e}}t_{n}\approx 100$ where it reaches the
reduced statistical fluctuation amplitudes
$\displaystyle\left|\displaystyle\langle\mathbf{E}_{\text{fluct}}\rangle\right|\approx
9\times 10^{-5}$. An advantage of the ensemble averaging technique is that
simulations do not need to be parallelized compared to an equivalent but more-
computationally-extensive simulation using
$N_{\text{\Large\textcalligra{s}}}\,$ times more macroparticles and
macroparticle weights $\delta N_{e}/N_{\text{\Large\textcalligra{s}}}\,$. In
the right panel of Figure 4, the ensemble average results is compared with the
substraction technique and the theoretical result

$\widehat{E}_{x}\left(k_{0},\,t_{n}>0\right)\underset{\omega_{b}t\ll
1}{\approx}\displaystyle\frac{E_{0}}{2}\exp{\left(-\gamma_{L_{0}}t_{n}\right)}.$
(141)

We find good agreements between the two different methods and the
Vlasov/Maxwell theory.

##  6 Conclusions

PIC simulations are one of the most widely used and most powerful tools in
plasma physics, being a core component and methodology of plasma physics as a
scientific discipline. However, the statistical physics principles of a PIC
code (and what is actually computed) are still a source of confusion among
many users and developers. Conflicting views assume that a PIC code evolves
the macroscopic Vlasov-Maxwell set of equations, others consider PIC
simulations describe the microscopic Klimontovich/Maxwell set of equations,
while others claim that PIC simulations cannot be compared to (or used with)
the kinetic theory of plasmas. In this paper, we clarify the kinetic
properties of PIC simulations. We have derived the deterministic description
of a PIC simulation electrostatic plasma, or PIC plasma. This description
consists of a Klimontovich-like equation discretized in space and time, but
not in velocity space, with a source term (in the right-hand side of the
Klimontovich-like equation) due to the internal tension force of
macroparticles. In the left-hand side of the Klimontovich-like equation, the
collective electrostatic force interpolated at macroparticle locations is
self-consistently coupled with the discretized Maxwell equations describing
the evolution of the (approximated) microscopic fields. The source terms of
the discretized Maxwell’s equations are obtained by depositing the electrical
charge and current of macroparticles onto the spatial grid such that the
Coulomb electrostatic fields generated by macroparticles are reduced at
distances smaller than the spatial grid spacing. If the minimum possible
distance between particles due to their binary electrostatic interaction (the
Landau length) is not resolved, then close-encounter binary Coulomb collisions
between macroparticles are underestimated in PIC codes. To understand the
statistical properties of PIC plasmas, we have applied the theory of
fluctuations in collisionless plasmas derived by [Rostoker, 1961, Klimontovich
and Silin, 1962, Dupree, 1963] starting from the statistical fluctuations
typically implemented by the random sampling used to initialize the velocities
of macroparticles. By neglecting numerical consistency errors from the
numerical schemes, we used the exact discrete Laplace-Fourier transform of the
PIC kinetic equations coupled with the exact discrete Laplace-Fourier
transform of self-consistent discretized Maxwell equations. We derived the
single-time autocorrelation of the electrostatic field and plasma densities
fluctuations, as well as the the single-time correlation between the
fluctuations of weighted macroparticle center phase-space density and
electrostatic field. This has allowed us to estimate the deviations from the
expected results of the simulation by applying the central limit theorem to
the statistical fluctuations in the limit of a large number of macroparticles
per cell as well as deriving the kinetic equations for each PIC simulation
plasma species. The former is useful to determine the spatial spectrum seeding
an eventual physical (or numerical) instability in a PIC simulation plasma or
to study a physical process of small amplitude that can be hidden in the
”numerical noise” (fluctuations), and the latter allowed us to emphasize that
”collisions” are indeed taken into account in a PIC plasma, where ”collisions”
mean here the the friction and diffusion of macroparticles in the fluctuating
component of the electrostatic field due to Landau damping and plasma waves
emission by macroparticles. The kinetic equations of the PIC plasma are indeed
described by a Vlasov-Lenard-Balescu-like equation with non-physical spatial
aliases in the ”collision” integral (that can be mitigated by using a
smoothing function with the cost of eventually losing a physical part of the
plasma spatial spectrum) and with time aliases (that can be mitigated by using
small time steps). Our theory recovers all results previously obtained
[Langdon, 1970b], confirming their approach of perturbing a discretized
Vlasovian plasma with a macroparticle test and then averaging the obtained
physical quantity over the initial macroparticle velocity distribution. Our
theory also extends their results to multiple-species multiple-weights PIC
plasmas. The underlying hypothesis and assumptions for the validity of the
theory have been systematically justified and/or estimated, closely connecting
this with the PIC plasma stability and the simulation time duration, with
estimates for these bounds provided as a function of the dimensionality of the
PIC plasma, the number of real electrons in the plasma Debye sphere and the
macroparticle weight factors. It is important to observe that the physical
statistical fluctuations are overestimated when weight factors $\delta N>1$
are used i.e. when a much-smaller-than-in-reality number of macroparticles is
chosen to represent the real plasma particles phase-space. We have used our
theoretical framework to explore the ”ensemble averaging technique” as a
possible technique to recover lower/physical statistical fluctuations. In this
method, a series of PIC simulations, representing identical statistical
realizations of the same macroscopic initial conditions, are averaged. The
equivalent statistical realizations, corresponding to different initial
microscopic conditions of the PIC plasma, are obtained from the random
sampling from the same distribution function of initial macroparticle
velocities. It is shown that the ensemble average of
$N_{\text{\Large\textcalligra{s}}}\,=\delta N$ electrostatic PIC simulations
is equivalent to the Vlasov-Lenard-Balescu/Maxwell set of equations
discretized in space and time. Even for macroparticle weights of unity, the
statistical fluctuations may hide theoretical features predicted in the Vlasov
limit of an infinite number of electrons in the Debye sphere. In this case,
the ensemble averaging technique can also be used to study in detail such
microphysics processes beyond the limits imposed by statistical fluctuations.
In order to illustrate the ”ensemble averaging technique”, we have applied it
to the study of linear Landau damping of one electrostatic wave perturbing a
PIC plasma in equilibrium by averaging the results obtained from
$N_{\text{\Large\textcalligra{s}}}\,=30$ simulations. The amplitude of the
statistical fluctuations is reduced by a factor
$\sqrt{N_{\text{\Large\textcalligra{s}}}\,}$, thus allowing the ensemble
averaging technique to identify short-lived very small amplitude perturbations
that would be undetected in a single simulation. The method shows good
agreement with the theory and the subtraction technique. This technique is
also particularly suited for modern computer architectures, since all
$N_{\text{\Large\textcalligra{s}}}\,$ simulations that will be averaged can
run in parallel, allowing for speed ups by a factor
$N_{\text{\Large\textcalligra{s}}}\,$ the global simulation (including the
individual simulations over which the ensemble averaging is going to be
performed), as compared with an equivalent simulation that uses
$N_{\text{\Large\textcalligra{s}}}\,\times N_{e,\text{mpc}}$ macroelectrons
per cell and a lower macroelectron weight $\delta
N_{e}/N_{\text{\Large\textcalligra{s}}}\,$. All
$N_{\text{\Large\textcalligra{s}}}\,$ simulations can be run in a sequential
way too. The ensemble averaging technique can thus also be used if the
equivalent simulation needs more RAM memory than what is available. The
flexibility of the technique, the theoretical understanding presented in this
paper, and the outstanding computational resources now available, open the way
to a broader and wider use of this technique in all areas of plasma physics
using PIC simulations.

## Appendix

### Appendix A
Derivation of the Laplace-Fourier transform of electrostatic field
fluctuations single-time autocorrelation

$\mathrm{For\,all\,functions\,g,}\displaystyle\left\\{\begin{array}[]{lll}&\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\displaystyle\left(2\pi\right)}^{3}}\displaystyle\sum_{p,q,r}g\left(\mathbf{k}_{p,q,r}\right)&=\displaystyle\int_{\mathbb{R}^{3}}\displaystyle\frac{d^{3}\mathbf{k}}{{\displaystyle\left(2\pi\right)}^{3}}g\left(\mathbf{k}\right)\cr&\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\sum_{m}g\left(\omega_{m}\right)&=\displaystyle\int_{\iota\nu-\infty}^{\iota\nu+\infty}\displaystyle\frac{d\omega}{2\pi}g\left(\omega\right)\end{array}\right..$
(142)

Therefore, by substituting (70) into (67), one can perform the integration
over $\mathbf{k}^{\prime}$ by using the $\mathbf{k}_{g}$-periodicity property
of $\left\\{\mathbf{k}\right\\}_{\mathbf{r}}$, $\varepsilon_{L}$ and
$\delta\widehat{f}_{a_{c}}$ and the changes of indices
$\left(p^{\prime},q^{\prime},r^{\prime}\right)$ into
$\left(p^{\prime}-p,q^{\prime}-q,r^{\prime}-r\right)$. One gets

$\mathbb{E}\displaystyle\left\\{\displaystyle\frac{{\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}^{2}}{8\pi}\right\\}=\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)$
(143)

where we have introduced the spectral energy density of the electrostatic
fluctuations

${\begin{array}[]{lcl}W_{\omega,\mathbf{k}}\displaystyle\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\frac{2\pi}{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)\varepsilon_{L}\displaystyle\left(\omega^{\prime},-\mathbf{k}\right)}\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)\overset{\frown}{\mathrm{S}}\left(-\mathbf{k}_{p,q,r}\right)\displaystyle\sum_{a}\delta
N_{a}\bar{n}_{a}{q_{a}}^{2}\cr&\times&\displaystyle\sum_{m}\displaystyle\sum_{m^{\prime}}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}^{\prime}\displaystyle\frac{F_{a_{0}}\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]}{\displaystyle\left(\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}\right)\displaystyle\left({\mathbf{k}}_{p,q,r}\cdot\mathbf{v}^{\prime}+{\omega^{\prime}}_{m^{\prime}}\right)}\exp{\displaystyle\left[-\iota\displaystyle\left(\omega+\omega^{\prime}\right)t_{n}\right]}.\end{array}}$
(144)

Since we have assumed that the PIC simulation plasma is stable which means
that all poles of the longitudinal PIC simulation plasma permittivity (63)
have a negative imaginary part, the second term in (144) $\propto
F_{a0}\left(\mathbf{v}\right)F_{a0}\left(\mathbf{v}^{\prime}\right)$ has
similarly poles $\omega^{\prime}$ in the lower half of the
$\omega^{\prime}$-planes, only. Thus, the undamped part of the electrostatic
field statistical fluctuations autocorrelation arises only from the
contribution of the first term on large time scale $\omega_{p}t_{n}\gg 1$. We
thus find after integration over $\mathbf{v}^{\prime}$ and by using again the
integral-sum property (142)

${\begin{array}[]{lcl}W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)&\underset{\omega_{p}t_{n}\gg
1}{=}&\displaystyle\int_{\iota\nu^{\prime}-\infty}^{\iota\nu^{\prime}+\infty}\displaystyle\frac{d\omega^{\prime}}{2\pi}\displaystyle\frac{2\pi}{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)\varepsilon_{L}\displaystyle\left(\omega^{\prime},-\mathbf{k}\right)}\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)\overset{\frown}{\mathrm{S}}\left(-\mathbf{k}_{p,q,r}\right)\displaystyle\sum_{a}\delta
N_{a}\bar{n}_{a}{q_{a}}^{2}\cr&\times&\displaystyle\sum_{m}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}\displaystyle\frac{F_{a_{0}}\displaystyle\left(\mathbf{v}\right)}{\displaystyle\left(\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}\right)\displaystyle\left({\mathbf{k}}_{p,q,r}\cdot\mathbf{v}+{\omega^{\prime}}\right)}\exp{\displaystyle\left[-\iota\displaystyle\left(\omega+\omega^{\prime}\right)t_{n}\right]}.\end{array}}$
(145)

Then, since

$\displaystyle\frac{1}{\left(\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}\right)\left(\mathbf{k}_{p,q,r}\cdot\mathbf{v}+\omega^{\prime}\right)}=\displaystyle\frac{1}{\omega_{m}+{\omega^{\prime}}}\displaystyle\left[\displaystyle\frac{1}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}}-\displaystyle\frac{1}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}+{\omega^{\prime}}}\right],$
(146)

the only one contributions that are not damped comes from the residue at poles
${\omega^{\prime}}=-\omega_{m}$. In this sense, the factor
$1/\left(\omega_{m}+\omega^{\prime}\right)$ is to be interpreted as $-\iota
2\pi\delta\left(\omega_{m}+\omega^{\prime}\right)$ leading to

${\begin{array}[]{lcl}W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)&\underset{\omega_{p}t_{n}\gg
1}{=}&-\iota\displaystyle\frac{2\pi}{{\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}\cr&\times&\displaystyle\sum_{m}\displaystyle\lim_{\nu^{\prime}\rightarrow
0^{+}}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}F_{a0}\displaystyle\left(\mathbf{v}\right)\displaystyle\left[\displaystyle\frac{1}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}+\iota\nu^{\prime}}-\displaystyle\frac{1}{\mathbf{k}_{p,q,r}\cdot\mathbf{v}-\omega_{m}-\iota\nu^{\prime}}\right]\cr&\underset{\omega_{p}t_{n}\gg
1}{=}&\displaystyle\frac{4\pi^{2}}{{\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}\cr&\times&\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}}$
(147)

Here, we have used the particle shape parity property
$\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)=\overset{\frown}{\mathrm{S}}\left(-\mathbf{k}\right)$
and the $\left(\omega_{g},\,\mathbf{k}_{g}\right)$-periodicity property of the
longitudinal PIC simulation plasma permittivity to simplify the equation. As a
conclusion, one may express

$\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)$
(148)

in the expression of electrostatic field statistical fluctuations
autocorrelation (67) where

$\widehat{\widehat{{\delta\mathbf{E}}^{2}}}\left(\omega,\mathbf{k}\right)=8\pi\,W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)$
(149)

reads (72).

### Appendix B
Derivation of the single-time autocorrelation of electrostatic field
fluctuations in PIC simulation plasmas at equilibrium

According to (156), one finds back the fluctuation-dissipation theorem derived
by [Callen and Welton, 1951]

$W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)=\displaystyle\frac{\widehat{\widehat{{\delta\mathbf{E}}^{2}}}\left(\omega,\mathbf{k}\right)}{8\pi}\underset{\omega_{p}t_{n}\gg
1}{=}-\delta
N\displaystyle\frac{k_{B}T}{\omega}\displaystyle\frac{\mathrm{Im}\displaystyle\left\\{\varepsilon_{L}\displaystyle\left(\omega,\,\mathbf{k}\right)\right\\}}{{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)}^{2}}$
(150)

only in the particular case of PIC simulations of a plasma at equilibrium for
which $\delta N_{e}=\delta N_{i}=\delta N$ and using a spectral solver with a
smoothing function filtering all spatial frequency aliases $p,q,r\neq 0,0,0$.
Indeed, only in this particular case, $\delta N$ can be extracted from the sum
over all species in (72) simplifying the electrostatic fluctuations energy
spectrum into

$W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right)\underset{\omega_{p}t_{n}\gg
1}{=}\displaystyle\frac{4\pi^{2}\delta
N}{{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)}^{2}}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{2}}{{\mathbf{k}}^{2}}\displaystyle\sum_{a}\bar{n}_{a}{q_{a}}^{2}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}F_{a0}\displaystyle\left(\mathbf{v}\right)\delta\displaystyle\left(\omega-\mathbf{k}\cdot\mathbf{v},\omega_{g}\right)$
(151)

while the imaginary part of PIC simulation electrostatic plasmas longitudinal
permittivity (156) reduces to

$\mathrm{Im}\displaystyle\left\\{\varepsilon_{L}\displaystyle\left(\omega,\,\mathbf{k}\right)\right\\}=-\displaystyle\frac{4\pi^{2}\omega}{k_{B}T}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{2}}{\mathbf{k}^{2}}\sum_{a}\bar{n}_{a}{q_{a}}^{2}\displaystyle\int_{\mathbb{R}^{3}}F_{a0}\displaystyle\left(\mathbf{v}\right)\delta\displaystyle\left(\omega-\mathbf{k}\cdot\mathbf{v},\omega_{g}\right)d^{3}\mathbf{v}$
(152)

thus providing the fluctuation-dissipation theorem (150) by recognizing (152)
in the numerator of (151). Still in this particular case, one can get a simple
estimate of the electrostatic field fluctuations energy spatial spectrum

$W_{\mathbf{k}}\left(\mathbf{k}\right)=\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}W_{\omega,\mathbf{k}}\displaystyle\left(\omega,\,\mathbf{k}\right).$
(153)

It reads

$W_{\mathbf{k}}\left(\mathbf{k}\right)\underset{\omega_{p}t_{n}\gg 1}{=}\delta
N\displaystyle\frac{k_{B}T}{2}\displaystyle\frac{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{2}}{{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)}^{2}+{\displaystyle\left(\mathbf{k}\lambda_{D}\right)}^{2}}$
(154)

by performing the $\omega$-integration using the residue theorem, assuming
$\Delta_{t}\leq\Delta_{x}/v_{T_{e}}$ and where we have introduced the Debye
screening length $\lambda_{D}=v_{T_{e}}/\omega_{p_{e}}$. Assuming then
$a_{x},\,a_{y}\,\mathrm{and}\,a_{z}\ll\lambda_{D}$ and
$\Delta_{x}=\Delta_{y}=\Delta_{z}\ll\lambda_{D}$ for simplicity, one may
consider $\overset{\frown}{\mathrm{S}}\left(\mathbf{k}\right)\sim 1$ and
estimate

$\mathbb{E}\displaystyle\left\\{\displaystyle\frac{{\delta\mathbf{E}\left(\mathbf{r}_{i,j,k},\,t_{n}\right)}^{2}}{8\pi}\right\\}=\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}W_{\mathbf{k}}\left(\mathbf{k}\right).$
(155)

While diverging in 2D and 3D collisionless plasmas for which $\delta N=1$,
$N_{e}\rightarrow\infty$ and
$\Delta_{x}=\Delta_{y}=\Delta_{z}=\Delta_{t}\rightarrow 0$, we find (73) for
PIC simulation plasmas at equilibrium by replacing
$d^{3}\mathbf{k}/{\left(2\pi\right)}^{3}$ with $dk_{x}/2\pi$ for 1D
simulations and with $dk_{x}dk_{y}/{\left(2\pi\right)}^{2}$ for 2D simulations
as well as the domain of integration $V_{\mathbf{k}_{g}}$ with the closed disk
$\\{(k_{x},\,k_{y})\,\mathrm{such}\,\mathrm{that}\,{k_{x}}^{2}+{k_{y}}^{2}\leq{(k_{g}/2)}^{2}\\}$
for 2D simulations and with the closed sphere
$\\{(k_{x},\,k_{y},\,k_{z})\,\mathrm{such}\,\mathrm{that}\,{k_{x}}^{2}+{k_{y}}^{2}+{k_{y}}^{2}\leq{(k_{g}/2)}^{2}\\}$
for 3D simulations instead of $V_{\mathbf{k}_{g}}=[-k_{g}/2,\,k_{g}/2]$, the
closed square $[-k_{g}/2,\,k_{g}/2]^{2}$ and the closed cube
$[-k_{g}/2,\,k_{g}/2]^{3}$ for 1D, 2D and 3D simulations respectively. An
empirical factor $1/2$ for 2D simulations and $2/3$ for 3D simulations has
been chosen to account for this latter approximation.

### Appendix C
Derivation of the ”collision operator”

According to (63) and the Plemelj formula, one may express the imaginary part
of PIC simulation electrostatic plasmas longitudinal permittivity according to

${\begin{array}[]{l}\mathrm{Im}\displaystyle\left\\{\varepsilon_{L}\displaystyle\left(\omega,\mathbf{k}\right)\right\\}=\displaystyle\frac{\pi}{{\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_{a}{\omega_{p_{a}}}^{2}\displaystyle\int_{\mathbb{R}^{3}}d^{3}\mathbf{v}\displaystyle\left\\{\mathbf{k}\right\\}_{\mathbf{r}}\cdot\displaystyle\frac{dF_{a0}}{d\mathbf{v}}\delta\displaystyle\left(\omega-\mathbf{k}_{p,q,r}\cdot\mathbf{v},\,\omega_{g}\right).\end{array}}$
(156)

Therefore, since the collision operator (81) is necessarily real, the first
term in (80) gives

${\begin{array}[]{lll}&&\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\sum_{p,q,r}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}\displaystyle\frac{8\pi^{2}\bar{n}_{a}q_{a}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}}{{\varepsilon_{L}\left(\omega,\,\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{2}}\delta
N_{a}\cr&&\hskip
120.00018pt\times\mathrm{Im}\displaystyle\left\\{\varepsilon_{L}\left(\omega,\,\mathbf{k}\right)\right\\}F_{a0}\left(\mathbf{v}_{a}\right)\delta\left(\omega-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a},\omega_{g}\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\sum_{p,q,r}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}\displaystyle\frac{32\pi^{4}\bar{n}_{a}q_{a}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}}{{\varepsilon_{L}\left(\omega,\,\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{4}}\delta
N_{a}\cr&&\hskip
120.00018pt\times\displaystyle\sum_{p^{\prime},q^{\prime},r^{\prime}}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\right)}^{2}\displaystyle\sum_{b}\displaystyle\frac{\bar{n}_{b}{q_{b}}^{2}}{m_{b}}\displaystyle\int_{\mathbb{R}^{3}}{\left\\{\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\right\\}_{\mathbf{r}}}\cdot\displaystyle\frac{dF_{b0}}{d\mathbf{v}_{b}}\cr&&\hskip
120.00018pt\times\;\delta\left(\omega-\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\cdot\mathbf{v}_{b},\,\omega_{g}\right)d^{3}\mathbf{v}_{b}F_{a0}\left(\mathbf{v}_{a}\right)\delta\left(\omega-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a},\,\omega_{g}\right).\end{array}}$
(157)

For the same reason, the contribution of the second term in (80) comes from
its imaginary part that can be simplified into

${\begin{array}[]{l}\hskip
10.00002pt-\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\mathrm{Im}\displaystyle\Big{\\{}\displaystyle\sum_{p,q,r}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}\displaystyle\frac{32\pi^{2}\bar{n}_{a}q_{a}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}}{{\varepsilon_{L}\left(\omega,\,\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{4}}\displaystyle\frac{\pi}{m_{a}}\displaystyle\sum_{p^{\prime},q^{\prime},r^{\prime}}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\right)}^{2}\cr\hskip
160.00024pt\times\displaystyle\sum_{m}\displaystyle\frac{{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}}{\omega_{m}-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a}}\cdot\displaystyle\frac{dF_{a0}}{d\mathbf{v}_{a}}\displaystyle\sum_{b}\delta
N_{b}\bar{n}_{b}{q_{b}}^{2}\cr\hskip
160.00024pt\times\displaystyle\int_{\mathbb{R}^{3}}F_{b0}\displaystyle\left(\mathbf{v}_{b}\right)\delta\displaystyle\left(\omega-\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\cdot\mathbf{v}_{b},\,\omega_{g}\right)d^{3}\mathbf{v}_{b}\Big{\\}}\cr=-\displaystyle\int_{V_{\mathbf{k}_{g}}}\displaystyle\frac{d^{3}\mathbf{k}}{{\left(2\pi\right)}^{3}}\displaystyle\int_{\iota\nu-\omega_{g}/2}^{\iota\nu+\omega_{g}/2}\displaystyle\frac{d\omega}{2\pi}\displaystyle\sum_{p,q,r}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}\displaystyle\frac{32\pi^{4}\bar{n}_{a}q_{a}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p,q,r}\right)}^{2}}{{\varepsilon_{L}\left(\omega,\,\mathbf{k}\right)}^{2}{\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}^{4}}\displaystyle\frac{1}{m_{a}}\displaystyle\sum_{p^{\prime},q^{\prime},r^{\prime}}{\overset{\frown}{\mathrm{S}}\left(\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\right)}^{2}\cr\hskip
140.00021pt\times\,\delta\left(\omega-\mathbf{k}_{p,q,r}\cdot\mathbf{v}_{a},\omega_{g}\right){\left\\{\mathbf{k}_{p,q,r}\right\\}_{\mathbf{r}}}\cdot\displaystyle\frac{dF_{a0}}{d\mathbf{v}_{a}}\displaystyle\sum_{b}\delta
N_{b}\bar{n}_{b}{q_{b}}^{2}\cr\hskip
140.00021pt\times\displaystyle\int_{\mathbb{R}^{3}}F_{b0}\displaystyle\left(\mathbf{v}_{b}\right)\delta\displaystyle\left(\omega-\mathbf{k}_{p^{\prime},q^{\prime},r^{\prime}}\cdot\mathbf{v}_{b},\,\omega_{g}\right)d^{3}\mathbf{v}_{b}\end{array}}$
(158)

according to the Plemelj formula. By performing the integration over $\omega$
using the integral-sum properties (142), one finally deduces (82).

### Appendix D - Landau damping BEPS input deck

pinput2
IDRUN = 1
INDX = 8, INDY = 1, NPX = 8388608, NPY = 2
NPXB = 0, NPYB = 0, INORDER = 2,
POPT = 2, DOPT = 2, DJOPT = 2
NTW = 1, NTP = 0, NTD = 10, NTA = 0, NTV = 0, NTS = 0
NUSTRT = 1, NTR = 0, PSOLVE = 1, RELATIVITY = 0
TEND = 500.0, DT = 0.200, CI = 0.1
QME = -1.000, VTX = 1.000, VTY = 1.000, VTZ = 1.000
VX0 = 0.000, VY0 = 0.000, VZ0 = 0.000
VDX = 0.000, VDY = 0.000, VDZ = 0.000
VTDX = 0.000, VTDY = 0.000, VTDZ = 0.000
AX = 0.866667, AY = 0.866667
NSRAND = 0, NDPROF = 0 AMPDX = 0.0, SCALEDX = 0.0, SHIFTDX = 0.0
AMPDY = 0.0, SCALEDY = 0.0, SHIFTDY = 0.0
MOVION = 0, NPXI = 8192, NPYI = 2048
NPXBI = 0, NPYBI = 0
QMI = 1.000, RMASS = 100.0
RTEMPXI = 20.000, RTEMPYI = 20.000, RTEMPZI = 20.000
VXI0 = 0.000, VYI0 = 0.000, VZI0 = 0.000
VDXI = 0.000, VDYI = 0.000, VDZI = 0.000
RTEMPDXI = 1.000, RTEMPDYI = 1.000, RTEMPDZI = 1.000
NSRANDI = 0, NDPROFI = 0
AMPDXI = 0.0, SCALEDXI = 0.0, SHIFTDXI = 0.0
AMPDYI = 0.0, SCALEDYI = 0.0, SHIFTDYI = 0.0
SORTIME = 100, SORTIMI = 250
NPLOT = 0, IDPAL = 1, NDSTYLE = 1, SNTASKS=1,
ntw=1,ntp=50

pinput2 jf
driver select= 1
amp=0.005,wavemode=13,wavew=1.18
rise = 0.,flat = 0.,fall=0.,yrise fall=1000.0
timerise=1.2,timeflat=1.2,timefall=1.2,
center1 = 200, center2 = 0, phase offset = 4
bow = 0., bow power = 0.
superGauss = 2

fvxmax=12., fvymax=12.
nphbx=150,nphby=150
nphxx=20,nphyx=20
nphxy=20,nphyy=2

## Acknowledgments

This work was partially sponsored by NNSA, NSF under Grant ACI-1339893 and the
European Research Council through the project ERC AdG InPairs no. 695088.

## References

  * [Abramowitz and Stegun, 1965] Abramowitz, M. and Stegun, I. (1965).  Handbook of Mathematical Functions.  Dover, New York.
  * [Balescu, 1960] Balescu, R. (1960).  Phys. of Fluids, 3(1):52–63.
  * [Birdsall and Fuss, 1969] Birdsall, C. K. and Fuss, D. (1969).  Journal of Computational Physics, 3(4):494.
  * [Birdsall and Langdon, 1991] Birdsall, C. K. and Langdon, A. B. (1991).  Plasma Physics via Computer Simulation.  The Adam Hilger Series on Plasma Physics. Adam Hilger Bristol, Philadelphia and New York.
  * [Bogoliubov, 1946] Bogoliubov, N. N. (1946).  JETP, 16(8):691–702.
  * [Born and Green, 1946] Born, M. and Green, H. S. (1946).  Proc. Roy. Soc. A, 188:10–18.
  * [Buneman, 1967] Buneman, O. (1967).  Journal of Computational Physics, 1(4):517 – 535.
  * [Burger et al., 1965] Burger, P., Dunn, D. A., and Halsted, A. S. (1965).  The Physics of Fluids, 8(12):2263–2272.
  * [Callen and Welton, 1951] Callen, H. B. and Welton, T. A. (1951).  Irreversibility and generalized noise.  Phys. Rev., 83:34–40.
  * [Dawson, 1964] Dawson, J. M. (1964).  Thermal relaxation in a one‐species, one‐dimensional plasma.  The Physics of Fluids, 7(3):419–425.
  * [Dawson, 1983] Dawson, J. M. (1983).  Rev. Mod. Phys., 55(2).
  * [de Boer and Uhlenbeck, 1962] de Boer, J. and Uhlenbeck, G. E. (1962).  N. N. Bogoliubov, in Studies in statistical mechanics.  eds. north-holland, amsterdam, 1962 edition.
  * [Decyk, 1982] Decyk, V. K. (1982).  Energy conservation theorem for electrostatic systems.  The Physics of Fluids, 25(7):1205–1206.
  * [Decyk, 1987] Decyk, V. K. (1987).  Simulation of microscopic processes in plasma.  Proc. 1987 International Conference on Plasma Physics, Kiev, USSR, April 1987, ed. A G Sitenko [World Scientific, Singapore, 1987] vol. II, p. 1075.
  * [Decyk, 1995] Decyk, V. K. (1995).  Skeleton pic codes for parallel computers.  Computer Physics Communications, 87(1):87–94.  Particle Simulation Methods.
  * [Decyk, 2007] Decyk, V. K. (2007).  Upic: A framework for massively parallel particle-in-cell codes.  Computer Physics Communications, 177(1):95–97.  Proceedings of the Conference on Computational Physics 2006.
  * [Decyk and Singh, 2014] Decyk, V. K. and Singh, T. V. (2014).  Particle-in-cell algorithms for emerging computer architectures.  Computer Physics Communications, 185(3):708–719.
  * [Dupree, 1963] Dupree, T. H. (1963).  The Physics of Fluids, 6(12):1714–1729.
  * [Esirkepov, 2001] Esirkepov, T. (2001).  Computer Physics Communications, 135(2):144 – 153.
  * [Grismayer et al., 2011] Grismayer, T., Fahlen, J. E., Decyk, V. K., and Mori, W. B. (2011).  Plasma Physics and Controlled Fusion, 53(7):074011.
  * [Guernsey, 1960] Guernsey, R. (1960).
  * [Hockney and Eastwood, 1981] Hockney, R. and Eastwood, J. (1981).  Computer Simulation Using Particles.  McGraw-Hill, New York.
  * [Hockney, 1965] Hockney, R. W. (1965).  J. ACM, 12(1):95–113.
  * [Hockney, 1966] Hockney, R. W. (1966).  The Physics of Fluids, 9(9):1826–1835.
  * [Hockney, 1971] Hockney, R. W. (1971).  Measurements of collision and heating times in a two-dimensional thermal computer plasma.  J. Comput. Phys. 8: No. 1, 19-44 (Aug 1971).
  * [Hubbard, 1961] Hubbard, J. (1961).  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 260(1300):114–126.
  * [Joglekar et al., 2014] Joglekar, A. S., Thomas, A. G. R., Fox, W., and Bhattacharjee, A. (2014).  Phys. Rev. Lett., 112:105004.
  * [Kirkwood, 1946] Kirkwood, J. G. (1946).  Journal of Chem. Phys., 14(3):180.
  * [Klimontovich, 1958] Klimontovich, Y. L. (1958).  JETP, 7(1):119.
  * [Klimontovich and Silin, 1962] Klimontovich, Y. L. and Silin, V. P. (1962).  JETP, 15(1).
  * [Landau, 1937] Landau, L. D. (1937).  JETP, 7:203.
  * [Landau, 1946] Landau, L. D. (1946).  Yad. Fiz., 10:25.
  * [Landau and Lifshitz, 1981] Landau, L. D. and Lifshitz, E. M. (1981).  Physical kinetics, volume 10 of A Course of Theoretical Physics.
  * [Langdon, 1970a] Langdon, A. (1970a).  Journal of Computational Physics, 6(2):247–267.
  * [Langdon, 1970b] Langdon, A. B. (1970b).  In Proc. Fourth Conf. Num. Sim. Plasmas, Naval Res. Lab., Washington D. C., pages 467–495.
  * [Langdon, 1979] Langdon, A. B. (1979).  Phys. fluids, 22:163.
  * [Langdon and Birdsall, 1970] Langdon, A. B. and Birdsall, C. K. (1970).  Theory of plasma simulation using finite‐size particles.  The Physics of Fluids, 13(8):2115–2122.
  * [Lenard, 1960] Lenard, A. (1960).  Annals of Physics, 10(3):390–400.
  * [Lindman, 1970] Lindman, E. (1970).  Journal of Computational Physics, 5(1):13 – 22.
  * [Liouville, 1838] Liouville, J. (1838).  Journ. de Math., 3(349).
  * [Montgomery and Nielson, 1970] Montgomery, D. and Nielson, C. W. (1970).  Thermal relaxation in one‐ and two‐dimensional plasma models.  The Physics of Fluids, 13(5):1405–1407.
  * [Nanbu, 1997] Nanbu, K. (1997).  Phys. Rev. E, 56:7314–7314.
  * [Okuda and Birdsall, 1970] Okuda, H. and Birdsall, C. K. (1970).  Collisions in a plasma of finite‐size particles.  The Physics of Fluids, 13(8):2123–2134.
  * [Rostoker, 1961] Rostoker, N. (1961).  Nuclear Fusion, 1(2):101.
  * [Rostoker, 1964] Rostoker, N. (1964).  Test particle method in kinetic theory of a plasma.  The Physics of Fluids, 7(4):491–498.
  * [Shanny et al., 1967] Shanny, R., Dawson, J. M., and Greene, J. M. (1967).  One‐dimensional model of a lorentz plasma.  The Physics of Fluids, 10(6):1281–1287.
  * [Thompson and Hubbard, 1960] Thompson, W. B. and Hubbard, J. (1960).  Rev. Mod. Phys., 32:714–718.
  * [Touati et al., 2014] Touati, M., Feugeas, J.-L., Nicolai, P., Gremillet, L., and Tikhonchuk, V. (2014).  New Jour. Phys., 16(073014).
  * [Tzoufras et al., 2011] Tzoufras, M., Bell, A., Norreys, P., and Tsung, F. (2011).  J. Comput. Phys., 230:6475–6494.
  * [Villasenor and Buneman, 1992] Villasenor, J. and Buneman, O. (1992).  Rigorous charge conservation for local electromagnetic field solvers.  Computer Physics Communications, 69(2):306–316.
  * [Virtamo and Tuomisto, 1979] Virtamo, J. and Tuomisto, H. (1979).  Verification of a simple collision operator for one‐dimensional plasma by simulation experiments.  The Physics of Fluids, 22(1):172–175.
  * [Vlasov, 1938] Vlasov, A. A. (1938).  JETP, 8(3):291.
  * [Yee, 1966] Yee, K. (1966).  IEEE Transactions on Antennas and Propagation, 14(3):302–307.
  * [Yu et al., 1965] Yu, S. P., Kooyers, G. P., and Buneman, O. (1965).  Journal of Applied Physics, 36(8):2550–2559.
  * [Yvon, 1935] Yvon, J. (1935).  Actual. Sci. et Indust., (203).


# Machine Learning Interatomic Potential for Simulations of Carbon at Extreme
Conditions

Jonathan T. Willman  University of South Florida, Tampa, Florida 33620, USA
Kien Nguyen-Cong  University of South Florida, Tampa, Florida 33620, USA
Ashley S. Williams  University of South Florida, Tampa, Florida 33620, USA
Anatoly B. Belonoshko  Royal Institute of Technology, 106691 Stockholm, Sweden
Stan G. Moore  Sandia National Laboratories, Albuquerque, New Mexico 87185,
USA  Aidan P. Thompson  Sandia National Laboratories, Albuquerque, New Mexico
87185, USA  Mitchell A. Wood  Sandia National Laboratories, Albuquerque, New
Mexico 87185, USA  Ivan I. Oleynik<EMAIL_ADDRESS>University of South
Florida, Tampa, Florida 33620, USA

###### Abstract

A Spectral Neighbor Analysis (SNAP) machine learning interatomic potential
(MLIP) has been developed for simulations of carbon at extreme pressures (up
to 5 TPa) and temperatures (up to 20,000 K). This was achieved using a large
database of experimentally relevant quantum molecular dynamics (QMD) data,
training the SNAP potential using a robust machine learning methodology, and
performing extensive validation against QMD and experimental data. The
resultant carbon MLIP demonstrates unprecedented accuracy and transferability
in predicting the carbon phase diagram, melting curves of crystalline phases,
and the shock Hugoniot, all within $3\,\mathrm{\%}$ of QMD. By achieving
quantum accuracy and efficient implementation on leadership class high
performance computing systems, SNAP advances frontiers of classical MD
simulations by enabling atomic-scale insights at experimental time and length
scales.

Carbon at extreme pressures and temperatures is a topic of great scientific
interest for several disciplines including planetary science (Ross, 1981;
Ancilotto  _et al._ , 1997; Benedetti  _et al._ , 1999; Kraus  _et al._ ,
2017; Madhusudhan  _et al._ , 2012) and inertial confinement fusion (ICF)
research (Ross  _et al._ , 2015; Millot  _et al._ , 2018; Hopkins  _et al._ ,
2019; Kritcher  _et al._ , 2021). Methane ice at megabar pressures and
temperatures of thousands of kelvins is predicted to convert to solid or
liquid carbon in the cores of giant planets (Ancilotto  _et al._ , 1997;
Benedetti  _et al._ , 1999; Kraus  _et al._ , 2017). A successful suppression
of hydrodynamic instabilities seeded by solid and liquid carbon phases
appearing upon strong compression of the outer diamond ablation shell of an
ICF capsule (Casey  _et al._ , 2021) was the key for achieving a record-
breaking, near-threshold fusion energy ignition at the National Ignition
Facility (Tollefson, 2021).

The exploration of carbon’s behavior at extreme conditions is challenging for
both theory and experiment. Shock and ramp compression experiments using
powerful lasers (Duffy and Smith, 2019), pulsed power (Sinars  _et al._ ,
2020) and bright X-ray sources (McMahon, 2020; Seddon  _et al._ , 2017)
uncovered complicated mechanisms of inelastic deformations (Pavlovskii, 1971;
Kondo and Ahrens, 1983; McWilliams  _et al._ , 2010; Lang  _et al._ , 2018;
Katagiri  _et al._ , 2020; Winey  _et al._ , 2020), anomalous strength of
diamond (McWilliams  _et al._ , 2010; Lang  _et al._ , 2018; Katagiri  _et
al._ , 2020; Winey  _et al._ , 2020; MacDonald  _et al._ , 2020), unusual
melting (Brygoo  _et al._ , 2007; Knudson  _et al._ , 2008; Eggert  _et al._ ,
2010; Gregor  _et al._ , 2017; Millot  _et al._ , 2018, 2020) and liquid
carbon properties (Knudson  _et al._ , 2008; Eggert  _et al._ , 2010), as well
as extreme metastability of diamond well beyond the pressure-temperature range
of its thermodynamic stability (Lazicki  _et al._ , 2021). Molecular dynamics
(MD) simulations can provide a fundamental understanding of these phenomena,
but to be of experimental relevance, it must accurately describe interatomic
interactions in a system consisting of a large number of atoms.

Previous simulations of carbon at extreme conditions were predominantly
performed using quantum molecular dynamics (QMD) based on density functional
theory (DFT) (Scandolo  _et al._ , 1996; Grumbach and Martin, 1996; Wang  _et
al._ , 2005; Correa  _et al._ , 2006, 2008; Knudson  _et al._ , 2008; Sun  _et
al._ , 2009; Martinez-Canales  _et al._ , 2012; Benedict  _et al._ , 2014).
Due to high computational cost, QMD simulations are limited to several hundred
atoms for up to tens of picoseconds, which is insufficient for uncovering non-
equilibrium processes at experimental time (ns) and length
($\mathrm{\mathrm{\mu}}m$) scales. In principle, these scales can be accessed
by classical MD simulations on massively parallel computers (Thompson  _et
al._ , 2022). However, empirical interatomic potentials developed for carbon
at near ambient conditions (Brenner  _et al._ , 2002; Stuart  _et al._ , 2000;
Los and Fasolino, 2003; Pastewka  _et al._ , 2008; Srinivasan  _et al._ ,
2015) singularly fail upon extension to high pressures and temperatures
(Oleynik  _et al._ , 2008; Perriot  _et al._ , 2013), thus compromising
predictive power of atomistic simulations.

The advent of machine learning interatomic potentials (MLIPs) (Behler and
Parrinello, 2007; Bartók  _et al._ , 2010) opens up exceptional opportunities
for achieving a classical description of chemical bonding with quantum
accuracy (Friederich  _et al._ , 2021). Although numerous MLIPs have been
recently introduced and successfully applied to modeling properties of
materials at ambient conditions (Behler and Parrinello, 2007; Bartók  _et al._
, 2010; Thompson  _et al._ , 2015; Wood and Thompson, 2018; Shapeev, 2015;
Drautz, 2019; Zhang  _et al._ , 2018), including carbon (Khaliullin  _et al._
, 2010; Rowe  _et al._ , 2020), their exceptional power in describing diverse
and complex atomic environments at megabar pressures and tens of thousand of
kelvins has yet to be demonstrated.

This letter reports a significant advance in development of a quantum-accurate
Spectral Neighbor Analysis Potential (SNAP) for simulations of carbon at
extreme pressure-temperature (P-T) conditions. This includes construction of
an experimentally relevant training database, implementation of a robust SNAP
machine learning training methodology, and extensive validation against QMD
and experimental data. The end result is the first MLIP that delivers
unprecedented accuracy in simulating carbon over a remarkably wide range of
pressures (from 0 to 50 Mbar) and temperatures (up to 20,000 K).

Figure 1: SNAP training database: a) pressure-temperature map of QMD and
static DFT simulations included in the database, each represented by a P-T
point on carbon phase diagram sampling FC8 (diamond), BC8, SC solid and liquid
phases (total number of structures - $636$); b) pressure-temperature-density-
energy/atom distribution.

In general, MLIPs fingerprint a unique local atomic environment around each
atom by a set of descriptors. SNAP’s descriptors are bispectrum components
$\\{\mathbf{B^{\mathit{i}}}\\}$ of the local neighbor density projected onto a
basis of hyper-spherical harmonics in four dimensions (Thompson  _et al._ ,
2015; Wood and Thompson, 2018). Other successful MLIPs – Neural Network
Potentials (NNP) (Behler and Parrinello, 2007; Zhang  _et al._ , 2018),
Gaussian Approximation Potential (GAP) (Bartók  _et al._ , 2010), the Moment
Tensor Potential (MTP) (Shapeev, 2015) – employ mathematically different, but
physically similar descriptors. All of them, including SNAP can be mapped into
a general Atomic Cluster Expansion (ACE) descriptor framework (Drautz, 2019).

Herein, the total potential energy of the system of $N$ atoms is written as a
sum of atomic energies $E^{i}$, which are quadratic functions of the
bispectrum coefficients $\mathbf{B^{\mathit{i}}}$ (Wood and Thompson, 2018)

$E_{tot}(\\{\mathbf{r}^{N}\\})=\sum_{i}E^{i};\qquad
E^{i}=\mathbf{\boldsymbol{\boldsymbol{\beta}}\cdotp
B^{\mathit{i}}+}\frac{1}{2}\mathbf{B^{\mathit{i}}\cdotp\boldsymbol{\alpha}\cdotp
B^{\mathit{i}}}.$ (1)

Machine-learning techniques are used to determine the symmetric matrix
$\boldsymbol{\boldsymbol{\alpha}}$ and the vector
$\boldsymbol{\boldsymbol{\beta}}$, the unknown parameters of SNAP, to
reproduce potential energy, atomic forces and stress tensor for each structure
in the DFT training database.

SNAP displays a good balance between computational cost and accuracy (Zuo  _et
al._ , 2020). Both are controlled by SNAP hyperparameters – the cutoff radius
$r_{cut}$ and the integer angular momentum $J$. The former specifies the
number of neighbor atoms participating in the fingerprinting of atomic
environment around each atom $i$ and the latter refers to the dimensionality
of the descriptor space, i.e. the number of bispectrum
descriptors$\\{\mathbf{B^{\mathit{i}}}\\}$: $(J+1)(J+3/2)(J+2)/2$ for each
atom $\mathit{i}$.

The SNAP development includes: (i) construction of a robust training database
of first-principles QMD data; (ii) machine-learning training; and (iii)
extensive validation against QMD and experiment. These three tightly connected
steps constitute a single development cycle. Several of such cycles are
executed to improve upon deficiencies observed in previous iterations, see the
development workflow in supplemental Fig. S2 (S_m, ). For example, during the
validation, two-phase SNAP MD simulations produce melting lines of several
high-pressure carbon phases in a substantial disagreement with QMD. The
problem has been traced back to SNAP inaccuracies in calculation of enthalpies
of solid and liquid phases along the melting line, which, according to
Clausius–Clapeyron relation, define its slope $dT/dP$. Therefore, separate
liquid and solid phases were added to the QMD training database to complement
original combined two-phase solid-liquid structures. This update resulted in
substantial improvements in SNAP accuracy upon execution of a new training
cycle.

Figure 2: Validation of SNAP against QMD: (a) carbon phase diagram, including
melting lines of diamond, BC8 and SC crystalline phases of carbon at pressures
up to $5,000$ GPa and hydrostatic shock Hugoniot; (b) Radial distribution
functions $g(r)$ for diamond, BC8 and liquid phase at diamond-bc8-liquid
triple point ($840$ GPa, $7510$ K). (MD trajectories are averaged over $20$ ps
time interval); c) density difference $\Delta\rho=\rho_{s}-\rho_{l}$ between
solid (s) and liquid (l) phases (top panel) and average coordination number of
carbon atoms (bottom panel) as a function of pressure along the melting lines
of diamond and BC8. $\Delta\rho=0$ corresponds to melting line maxima: diamond
– at $500\,\mathrm{GPa}$, BC8 – at $1,300\,\mathrm{GPa}$.

One of the highlights of our MLIP development is the dedicated construction of
a SNAP training database of experimentally relevant QMD of large systems (by
QMD standards) – up to 1,024 atoms. These are individual frames (1-3 frames
per simulation) from QMD production simulations of physical properties of
carbon (melting lines, hydrostatic and uniaxial isotherms, shock Hugoniot)
performed within a wide range of pressures (from $0$ to $5$ TPa), temperatures
(from $0$ to $40,000$ K) and densities (from $2.9$ to $13.6$
$\mathrm{g/cm^{3}}$). For example, MD frames for each $(P,T)$ point along the
melting lines of several crystalline carbon phases were taken from production
two-phase QMD simulation (Belonoshko, 1994). These complex cells contain both
liquid and solid parts separated by a realistic solid-liquid interface, which
add an additional complexity to the SNAP database. A series of complimentary
QMD simulations of liquid and solid phases were also included in the database.
The QMD data is supplemented by static DFT calculations of binding energy
curves, point and extended defects, and metastable carbon structures obtained
from dedicated crystal structure searching (Williams  _et al._ , 2022). Fig. 1
displays the range of pressures, temperatures and energies covered by QMD
simulations. The SNAP database, consisting of $636$ structures, samples
124,907 unique atomic environments. Each structure of $N$ atoms contributes
$1$ total energy, $6$ stress components and $3N$ atomic forces resulting in a
total training complexity (the number of regression equations to fit) of
$388,077$. Additional information on database composition and its generation
is provided in the Supplemental Material (S_m, ).

Once the database is constructed, the SNAP model for the total potential
energy, stress tensor components and atomic forces is trained using machine
learning techniques to determine SNAP parameters $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$ through minimization of an objective function \- a sum of
the normalized squared differences between DFT and SNAP energies, stresses and
atomic forces (Fig. S2) (S_m, ). The training is performed in a series of
iterations. For a given set of weights, the SNAP parameters
$\boldsymbol{\boldsymbol{\alpha}}$ and $\boldsymbol{\beta}$ are determined
through weighted linear regression as implemented in FitSNAP package (fit, ).
The weights are then optimized to minimize the objective function using a
genetic algorithm (GA) within DAKOTA software package (Adams, 2021). At every
step of GA minimization FitSNAP is called with the current set of weights to
determine new $\boldsymbol{\boldsymbol{\alpha}}$ and
$\mathbf{\boldsymbol{\beta}}$, which are then fed back to update the objective
function being minimized (Fig. S2) (S_m, ). The iterations stop when the GA
minimization converges to a final set of weights. In addition to group
weights, the SNAP hyperparameters, $r_{cut}$ and $J$ are optimized in an
offset fashion to find a right balance between SNAP accuracy and computational
efficiency (Fig. S2 (S_m, )). The final values for the SNAP hyperparameters
are $r_{cut}=2.7\,\mathrm{\text{\r{A}}}$ and $J=4$. The resultant quality of
SNAP training is discussed in Supplemental Material (S_m, ).

The critical part of SNAP potential development is a thorough validation of
SNAP MD results against QMD and experimental data. Only one to three QMD
frames per (P,T) state were included in the SNAP training database. Therefore,
simulating these states with SNAP using much larger simulations cells and
performing thermodynamic averaging of atomic trajectories containing tens of
thousands of frames to obtain stresses, densities, internal energies, as well
as radial distribution functions is considered as a rigorous validation test
against QMD. Further, these validation simulations are extended to sample a
variety of pressures and temperatures not included in the database and
compared to both QMD and experiment to demonstrate SNAP transferability.

The important validation test is concerned with carbon phase diagram,
including melting lines of several high-pressure phases at pressures up to
$5,000$ GPa, see Fig. (2a). A series of isobaric-isothermal NPT two-phase MD
simulations were run at a given pressure but varying temperature to determine
the P-T value of the phase coexistence (Belonoshko, 1994). SNAP melting lines
are in excellent agreement with those from QMD, the average temperature error
being $\sim 200\,\mathrm{K}$ or $3\%$ in a pressure interval from $0$ to
$5,000$ GPa. Fig. 2(b) displays SNAP and QMD radial distribution functions
$g(r)$ at diamond-BC8-liquid triple point: they are almost indistinguishable
from each other.

A remarkable property of carbon at extreme conditions is the negative slope
($dT/dP$) of diamond and BC8 melting lines at high pressures (Grumbach and
Martin, 1996; Wang  _et al._ , 2005; Correa  _et al._ , 2006, 2008; Benedict
_et al._ , 2014; Eggert  _et al._ , 2010). This is because liquid carbon
becomes denser than the corresponding solid phase upon increase of pressure
above $\sim 500$ GPa for diamond and $\sim 1,300$ GPa for BC8, see Fig. 2(c).
This can be traced back to a significant increase of carbon packing in the
liquid as carbon atom coordination changes from less than 4 to higher values,
see Fig. 2(d). SNAP accurately predicts this subtle change in density upon
increase of pressure as well as corresponding pressure-dependent evolution of
the average coordination of carbon atoms ( Fig. 2(d)), which is in an
excellent agreement with QMD. The latter is a result of “snap-on” agreement
between QMD and SNAP radial distribution functions for both diamond and BC8
over a large range of pressures (Fig. S4 (S_m, )).

Figure 3: Carbon shock Hugoniot calculated by QMD and SNAP and compared with
experimental data. (a) $U_{s}-U_{p}$ Hugoniot (b) pressure-density Hugoniots.
Points correspond to experimental data.

Another validation test of great experimental importance is the prediction of
the carbon shock Hugoniot, which passes through both solid and liquid parts of
the phase diagram (Fig. 2(a)). The Hugoniot points are calculated in a
hydrodynamic approximation by ignoring crystalline anisotropy in a series of
MD simulations at a given pressure $P$ but varying temperature $T$ to satisfy
the Hugoniot condition of conservation of mass, momentum, and energy:
$\frac{1}{2}(P+P_{0})(V-V_{0})=(E-E_{0})$, where $P$, $V,$ and $E$ are the
pressure, volume, and internal energy at a given point on the Hugoniot, and
$P_{0}$, $V_{0}$, and $E_{0}$ are those at ambient conditions of $0$ GPa and
$300$ K. The Hugoniot in $(P,T)$ space is shown in Fig. 2(a), in pressure-
density $(P-\rho)$ space – in Fig. 3(b) and particle velocity $U_{p}$ – shock
velocity $U_{s}$ space – in Fig. 3(a). The SNAP Hugoniots (red lines in Figs
2(a) and 3) are in a very good agreement with those from QMD (black lines in
Figs 2(a) and 3). Visible differences in temperature of the _P–T_ Hugoniot at
very high temperatures (Fig. 2(a)) are due to electronic entropy effects,
which are not captured by any classical interatomic potential.

Both SNAP and QMD Hugoniots are also in a very good agreement with multiple
experiments in a pressure range from $300$ to $1,500\,\mathrm{GPa}$ (Knudson
_et al._ , 2008; McWilliams  _et al._ , 2010; Gregor  _et al._ , 2017;
MacDonald  _et al._ , 2020; Millot  _et al._ , 2020; Winey  _et al._ , 2020;
Katagiri  _et al._ , 2020). Some difference between SNAP/QMD and experiment at
higher pressures is due to experimental uncertainty in density, which is not
measured directly, but rather determined using an impedance matching method
(Gregor  _et al._ , 2017; Millot  _et al._ , 2020). Differences at low
pressures between experiment (McWilliams  _et al._ , 2010; MacDonald  _et al._
, 2020; Winey  _et al._ , 2020) and SNAP/QMD (Fig. 3(a)) are due to
crystalline anisotropy and strength effects, which are not well described by a
hydrodynamic approximation (Winey  _et al._ , 2020). To make a proper
comparison with experiment in this split elastic-inelastic shock wave regime,
explicit large-scale SNAP MD simulations of piston-driven shock waves are
required, which will be the focus of future work.

Figure 4: (a) Comparison of diamond melting curves calculated by three
methods: QMD 2-phase (solid line), SNAP 2-phase (dashed line), and direct
melting of 1 million atom polycrystalline sample using SNAP (stars); (b) Time
progression of polycrystalline sample melting at $600$ GPa and $7750$ K. The
sample is initially composed of crystalline diamond regions (blue) separated
by grain boundaries (light blue). Liquid regions (red) emerge at the grain
boundaries, grow in size, and eventually consume the isolated diamond
crystallites.

To demonstrate SNAP’s ability to attack problems that are impossible to solve
with QMD, we simulated melting of polycrystalline diamond at pressures between
50 and 1,200 GPa using a 1 million atom sample (Fig. 4). For each pressure
point, a series of NPT simulations is performed to determine the temperature
at the onset of melting. In addition to validating SNAP, this simulation also
validates the two-phase melting curve calculation. The presence of grain
boundaries suppresses superheating characteristic of single crystals: the
melting starts at the most weakly bonded defective regions of the sample,
followed by the growth of liquid fraction at the expense of crystalline
grains, which gradually transform to shrinking round crystallites embedded in
the liquid carbon (Fig. 4).

This work represents a major step towards solving extremely challenging but
fundamentally important problem of predictive atomic-scale simulations of
carbon at extreme pressure-temperature conditions at experimental time and
length scales. The quantum-accurate SNAP is the first MLIP that describes the
properties of carbon at extreme pressures up to $5$ TPa and temperatures up to
$20,000$ K, including the phase diagram, melting curves of diamond, BC8 and
simple cubic phases of carbon and shock Hugoniots with unprecedented accuracy
within $3\,\mathrm{\%}$ of QMD results. SNAP’s linear scaling with number of
atoms, and its efficient implementation within the LAMMPS MD simulation
package (Thompson  _et al._ , 2022) allows billion atom simulations on
leadership class high performance computing systems (Nguyen-Cong  _et al._ ,
2021). By advancing frontiers of classical MD simulations, SNAP enables new
insights by uncovering fundamental atomic-scale mechanisms of materials
response which are difficult or even impossible to observe in experiment.

The work at USF is supported by DOE/NNSA (grant DE-NA- 0003910). The work at
KTH is supported by Swedish Scientific Council (VR) (grant 2017-03744). Sandia
National Laboratories is a multi-mission laboratory managed and operated by
National Technology and Engineering Solutions of Sandia, LLC, a wholly owned
subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s (DOE) National Nuclear Security Administration under Contract No. DE-
NA0003525. The computations were performed using leadership-class HPC systems:
OLCF Summit at Oak Ridge National Laboratory (ALCC and INCITE awards MAT198)
and TACC Frontera at University of Texas at Austin (LRAC award DMR21006).

## References

  * Ross (1981) M. Ross, Nature 292, 435 (1981).
  * Ancilotto  _et al._ (1997) F. Ancilotto, G. L. Chiarotti, S. Scandolo, and E. Tosatti, Science 275, 1288 (1997).
  * Benedetti  _et al._ (1999) L. R. Benedetti, J. H. Nguyen, W. A. Caldwell, H. Liu, M. Kruger, and R. Jeanloz, Science 286 (1999).
  * Kraus  _et al._ (2017) D. Kraus, J. Vorberger, A. Pak, N. J. Hartley, _et al._ , Nat. Astron. 1, 606 (2017).
  * Madhusudhan  _et al._ (2012) N. Madhusudhan, K. K. M. Lee, and O. Mousis, Astrophys. J. 759, L40 (2012).
  * Ross  _et al._ (2015) J. S. Ross, D. Ho, J. Milovich, T. Döppner, J. McNaney, _et al._ , Phys. Rev. E 91, 1 (2015).
  * Millot  _et al._ (2018) M. Millot, P. M. Celliers, P. A. Sterne, L. X. Benedict, A. A. Correa, S. Hamel, S. J. Ali, _et al._ , Phys. Rev. B 97, 144108 (2018).
  * Hopkins  _et al._ (2019) L. B. Hopkins, S. LePape, L. Divol, A. Pak, E. Dewald, _et al._ , Plasma Phys. Control. Fusion 61, 014023 (2019).
  * Kritcher  _et al._ (2021) A. L. Kritcher, A. B. Zylstra, D. A. Callahan, O. A. Hurricane, _et al._ , Phys. Plasmas 28, 072706 (2021).
  * Casey  _et al._ (2021) D. T. Casey, B. J. MacGowan, J. D. Sater, A. B. Zylstra, O. L. Landen, J. Milovich, O. A. Hurricane, A. L. Kritcher, _et al._ , Phys. Rev. Lett. 126, 25002 (2021).
  * Tollefson (2021) J. Tollefson, Nature 597, 163 (2021).
  * Duffy and Smith (2019) T. S. Duffy and R. F. Smith, Front. Earth Sci. 7, 1 (2019).
  * Sinars  _et al._ (2020) D. B. Sinars, M. A. Sweeney, C. S. Alexander, D. J. Ampleford, T. Ao, _et al._ , Phys. Plasmas 27, 070501 (2020).
  * McMahon (2020) M. I. McMahon, in  _Synchrotron Light Sources Free. Lasers_ , edited by E. J. Jaeschke, S. Khan, J. R. Schneider, and J. B. Hastings (Springer International Publishing, 2020) pp. 1857–1896.
  * Seddon  _et al._ (2017) E. A. Seddon, J. A. Clarke, D. J. Dunning, C. Masciovecchio, C. J. Milne, F. Parmigiani, D. Rugg, J. C. H. Spence, N. R. Thompson, K. Ueda, S. M. Vinko, J. S. Wark, and W. Wurth, Reports Prog. Phys. 80, 115901 (2017).
  * Pavlovskii (1971) M. N. Pavlovskii, Sov. Phys. - Solid State 13, 741 (1971).
  * Kondo and Ahrens (1983) K.-I. Kondo and T. J. Ahrens, Geophys. Res. Lett. 10, 281 (1983).
  * McWilliams  _et al._ (2010) R. S. McWilliams, J. H. Eggert, D. G. Hicks, D. K. Bradley, P. M. Celliers, D. K. Spaulding, T. R. Boehly, G. W. Collins, and R. Jeanloz, Phys. Rev. B 81, 27 (2010).
  * Lang  _et al._ (2018) J. M. Lang, J. M. Winey, and Y. M. Gupta, Phys. Rev. B 97, 104106 (2018).
  * Katagiri  _et al._ (2020) K. Katagiri, N. Ozaki, Y. Umeda, T. Irifune, N. Kamimura, K. Miyanishi, T. Sano, T. Sekine, and R. Kodama, Phys. Rev. Lett. 125, 185701 (2020).
  * Winey  _et al._ (2020) J. M. Winey, M. D. Knudson, and Y. M. Gupta, Phys. Rev. B 101, 184105 (2020).
  * MacDonald  _et al._ (2020) M. J. MacDonald, E. E. McBride, E. Galtier, M. Gauthier, E. Granados, D. Kraus, A. Krygier, A. L. Levitan, A. J. MacKinnon, I. Nam, W. Schumaker, P. Sun, T. B. van Driel, J. Vorberger, Z. Xing, R. P. Drake, S. H. Glenzer, and L. B. Fletcher, Appl. Phys. Lett. 116, 234104 (2020).
  * Brygoo  _et al._ (2007) S. Brygoo, E. Henry, P. Loubeyre, J. Eggert, M. Koenig, B. Loupias, A. Benuzzi-Mounaix, and M. Rabec Le Gloahec, Nat. Mater. 6, 274 (2007).
  * Knudson  _et al._ (2008) M. D. Knudson, M. P. Desjarlais, and D. H. Dolan, Science 322, 1822 (2008).
  * Eggert  _et al._ (2010) J. H. Eggert, D. G. Hicks, P. M. Celliers, D. K. Bradley, R. S. McWilliams, R. Jeanloz, J. E. Miller, T. R. Boehly, and G. W. Collins, Nat. Phys. 6, 40 (2010).
  * Gregor  _et al._ (2017) M. C. Gregor, D. E. Fratanduono, C. A. McCoy, D. N. Polsin, A. Sorce, J. R. Rygg, G. W. Collins, T. Braun, P. M. Celliers, J. H. Eggert, D. D. Meyerhofer, and T. R. Boehly, Phys. Rev. B 95, 144114 (2017).
  * Millot  _et al._ (2020) M. Millot, P. A. Sterne, J. H. Eggert, S. Hamel, M. C. Marshall, and P. M. Celliers, Phys. Plasmas 27, 102711 (2020).
  * Lazicki  _et al._ (2021) A. Lazicki, D. McGonegle, J. R. Rygg, D. G. Braun, D. C. Swift, _et al._ , Nature 589, 532 (2021).
  * Scandolo  _et al._ (1996) S. Scandolo, G. Chiarotti, and E. Tosatti, Phys. Rev. B 53, 5051 (1996).
  * Grumbach and Martin (1996) M. P. Grumbach and R. M. Martin, Phys. Rev. B 54, 15730 (1996).
  * Wang  _et al._ (2005) X. Wang, S. Scandolo, and R. Car, Phys. Rev. Lett. 95, 1 (2005).
  * Correa  _et al._ (2006) A. A. Correa, S. A. Bonev, and G. Galli, Proc. Natl. Acad. Sci. U. S. A. 103, 1204 (2006).
  * Correa  _et al._ (2008) A. Correa, L. Benedict, D. Young, E. Schwegler, and S. A. Bonev, Phys. Rev. B 78, 25 (2008).
  * Sun  _et al._ (2009) J. Sun, D. D. Klug, and R. Marton̂ák, J. Chem. Phys. 130, 194512 (2009).
  * Martinez-Canales  _et al._ (2012) M. Martinez-Canales, C. J. Pickard, and R. J. Needs, Phys. Rev. Lett. 108, 45704 (2012).
  * Benedict  _et al._ (2014) L. X. Benedict, K. P. Driver, S. Hamel, B. Militzer, T. Qi, A. A. Correa, A. Saul, and E. Schwegler, Phys. Rev. B 89, 224109 (2014).
  * Thompson  _et al._ (2022) A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comput. Phys. Commun. 271, 108171 (2022).
  * Brenner  _et al._ (2002) D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J. Phys. Condens. Matter 14, 783 (2002).
  * Stuart  _et al._ (2000) S. J. Stuart, A. B. Tutein, and J. A. Harrison, J. Chem. Phys 112, 6472 (2000).
  * Los and Fasolino (2003) J. H. Los and A. Fasolino, Phys. Rev. B 68, 024107 (2003).
  * Pastewka  _et al._ (2008) L. Pastewka, P. Pou, R. Pérez, P. Gumbsch, and M. Moseler, Phys. Rev. B 78, 161402 (2008).
  * Srinivasan  _et al._ (2015) S. G. Srinivasan, A. C. T. van Duin, and P. Ganesh, J. Phys. Chem. A 119, 571 (2015).
  * Oleynik  _et al._ (2008) I. I. Oleynik, A. C. Landerville, S. V. Zybin, M. L. Elert, and C. T. White, Phys. Rev. B 78, 180101 (2008).
  * Perriot  _et al._ (2013) R. Perriot, X. Gu, Y. Lin, V. V. Zhakhovsky, and I. I. Oleynik, Phys. Rev. B 88, 064101 (2013).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 1 (2007).
  * Bartók  _et al._ (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
  * Friederich  _et al._ (2021) P. Friederich, F. Häse, J. Proppe, and A. Aspuru-Guzik, Nat. Mater. 20, 750 (2021).
  * Thompson  _et al._ (2015) A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles, and G. J. Tucker, J. Comput. Phys. 285, 316 (2015).
  * Wood and Thompson (2018) M. A. Wood and A. P. Thompson, J. Chem. Phys. 148, 241721 (2018).
  * Shapeev (2015) A. V. Shapeev, Multiscale Model. Simul. 14, 1 (2015).
  * Drautz (2019) R. Drautz, Phys. Rev. B 99, 1 (2019).
  * Zhang  _et al._ (2018) L. Zhang, J. Han, H. Wang, R. Car, and E. Weinan, Phys. Rev. Lett. 120, 143001 (2018).
  * Khaliullin  _et al._ (2010) R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler, and M. Parrinello, Phys. Rev. B 81, 18 (2010).
  * Rowe  _et al._ (2020) P. Rowe, V. L. Deringer, P. Gasparotto, G. Csányi, and A. Michaelides, J. Chem. Phys. 153, 034702 (2020).
  * Zuo  _et al._ (2020) Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler, G. Csányi, A. V. Shapeev, A. P. Thompson, M. A. Wood, and S. P. Ong, J. Phys. Chem. A 124, 731 (2020).
  * (56) See Supplemental Material at http://link.aps.org/ supplemental/ for details on (1) SNAP training database; (2) SNAP machine learning development workflow; (3) Accuracy assessment of SNAP training; (4) SNAP validation: Radial Distribution Functions of solid and liquid phases along diamond and BC8 melting lines.
  * Belonoshko (1994) A. B. Belonoshko, Geochimica et Cosmochimica Acta 58, 4039 (1994).
  * Williams  _et al._ (2022) A. S. Williams, K. Nguyen-Cong, J. T. Willman, and I. I. Oleynik, Phys. Chem. Lett. , submitted for publication (2022).
  * (59) FitSNAP: Software for generating SNAP machine learning interatomic potentials: http://github.com/FitSNAP.
  * Adams (2021) B. Adams,  _Dakota, A Multilevel Parallel Object-Oriented Framework for Design Optimization, Parameter Estimation, Uncertainty Quantification, and Sensitivity Analysis: Version 6.14 User’s Manual_ (2021).
  * Nguyen-Cong  _et al._ (2021) K. Nguyen-Cong, J. T. Willman, S. G. Moore, A. B. Belonoshko, R. Gayatri, E. Weinberg, M. A. Wood, A. P. Thompson, and I. I. Oleynik, Proc. Int. Conf. High Perform. Comput. Networking, Storage Anal. , 1 (2021).