# Benchmark between N-body, PIC, and semi-Lagrangian simulations of Landau-damped Langmuir wave

Maxime Lesur, Jérôme Moritz, Etienne Gravier, Thomas Drouot

► To cite this version:

Maxime Lesur, Jérôme Moritz, Etienne Gravier, Thomas Drouot. Benchmark between N-body, PIC, and semi-Lagrangian simulations of Landau-damped Langmuir wave. Physics of Plasmas, 2025, 32 (1), pp.013901. 10.1063/5.0240349. hal-04835129

HAL Id: hal-04835129

https://hal.univ-lorraine.fr/hal-04835129v1

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Benchmark between N-body, PIC, and semi-Lagrangian
simulations of Landau-damped Langmuir wave

M. Lesur¹,², J. Moritz¹, E. Gravier¹, T. Drouot¹

¹Institut Jean Lamour, Université de Lorraine,

54506 Vandoeuvre-lès-Nancy Cedex, France and

²Institut Universitaire de France (IUF), 75005 Paris, France*

(Dated: December 10, 2024)
Abstract

This work presents a benchmark study comparing three distinct numerical methods – Particle-In-Cell (PIC), semi-Lagrangian, and N-body simulations – for analyzing the damping of Langmuir waves in a one-dimensional Vlasov-Poisson plasma system. Each approach has unique advantages in terms of accuracy, resolution, and computational cost. The comparison aims to discriminate between numerical artifacts and physical phenomena, identifying the contribution of finite particle numbers and boundary conditions in both linear and nonlinear regimes. The study demonstrates strong agreement between the PIC and semi-Lagrangian methods in both regimes. N-body simulations, while requiring a specific method to overcome statistical noise, agree in the limit of many bodies (> 500). Crucial subtleties regarding initial and boundary conditions are discussed throughout.

PACS numbers: 07.05.Tp,41.20.Cv,47.27.E-,47.27.ek,52.25.Dg,52.35.Fp,52.35.Mw,52.65.Ff

# I. INTRODUCTION

The 1D Vlasov-Poisson model is a fundamental tool in plasma physics, particularly for studying the behavior of collisionless plasmas. It is extensively used to analyze the propagation and interaction of plasma waves, such as Langmuir waves and ion-acoustic waves, and to understand phenomena like Landau damping [24]. This model is also crucial in exploring beam-plasma interactions, where it helps to investigate instabilities (in both linear and nonlinear terms [45]), and nonlinear phase space structures like electron holes [4, 11, 25, 33, 52].

In fusion research, this model is applied to understand the interaction between particle beams and plasmas, which is essential for both inertial confinement and magnetic confinement fusion systems. It also plays a significant role in the design and analysis of plasma-based devices, such as particle accelerators and plasma thrusters. Its applications extend to space and astrophysical plasmas, aiding in the study of kinetic instabilities and the dynamics of charged particles in environments like the solar wind and Earth's magnetosphere. Additionally, the model is the foundation for many numerical simulations, and a tool for verification of both numerical simulation codes and new theoretical frameworks [5, 7, 9, 16, 44, 46, 48].

<footnote>* maxime.lesur@polytechnique.org</footnote>
Most high-temperature plasmas of interest do not strictly adhere to the Vlasov (or collisionless Boltzmann) equation, even notwithstanding effects of finite particle numbers, and higher-order correlations [17]. They are influenced by collisions, sources and sinks, energetic, chemical, atomic, and nuclear processes, and possibly external forces. However, in many practical cases, it is possible to focus on scales (spatial and temporal) where the Vlasov equation [41, 59] serves as a reasonable approximation or starting point.

In this paper, we focus on the damping of Langmuir waves, which can be caused by Coulomb collisions, and/or by collisionless Landau damping. Landau damping is due to the resonant interaction between the wave and the particles. This phenomenon has been extensively studied both analytically and numerically (recent examples: [28, 38, 57, 63], and review: [49]), as it plays a crucial role in many plasma processes.

In this study, we investigate the Landau-damped Langmuir wave using three different numerical methods: a Particle-In-Cell (PIC) method, a semi-Lagrangian approach, and a N-body approach. Each of these methods offers distinct advantages and challenges in simulating the Vlasov-Poisson system. The objective is to benchmark these numerical approaches against each other to determine their accuracy, and suitability for different plasma conditions and investigations. The main point of this paper is, by analysing fine details of a Landau-damped wave using these 3 fundamentally different numerical methods, to discriminate without ambiguity, among the observed phenomena, which are numerical artifacts and which are physical. Then, among the physical phenomena, we can discriminate which are due to the small number of particles, which depend on boundary conditions, which depend on fluctuation amplitude (we investigate both linear and nonlinear regimes), and which are universal.

This study is most of all an effort to clarify numerical artifacts (e.g., noise in PIC and statistical effects in N-body) versus physical phenomena (e.g., nonlinear Landau damping and BGK modes), and to provide insights into cross-code validation, with crucial subtleties regarding initial and boundary conditions.

While our results are presented for a 1D electrostatic plasma, the implications extend to more complex systems. For instance, understanding the noise characteristics and resolution limits of PIC codes in this simplified setting provides a foundation for their application to
2D/3D electromagnetic problems.

The remainder of this paper is structured as follows: Section II describes the 1D Vlasov-Poisson model and physical plasma setup. Section III details the numerical methods employed, including the PIC, semi-Lagrangian, and N-body approaches. Section IV discusses the boundary conditions, including initial conditions, and their implementation in each method. Section V provides a refresher on Landau damping, including the high-amplitude regime. Section VI presents the results of the benchmarks, comparing the different methods in both linear and nonlinear regimes. Finally, Section VII provides conclusions and insights drawn from the comparative analysis.

## II. MODEL

### A. 1D Vlasov-Poisson model

We consider a uniform, unmagnetized 1D plasma, which is relevant when the physics under consideration is essentially one-dimensional. This is the case in contexts of e.g. solar wind plasmas, laser/plasma interactions, and even magneto-hydrodynamics modes in strongly magnetized fusion plasmas (where the fast nonlinear wave/particle interactions can be modeled by an electrostatic 1D model [3, 13, 27, 29, 60]).

We further assume that ions form a fixed, neutralizing background. This precludes the investigation of ion-acoustic-type waves, but is usually a good approximation for electrostatic plasma (Langmuir) waves.

We adopt the 1D Vlasov-Poisson model, which yields the self-consistent evolution of the distribution function of electrons $f(x, v, t)$ and the electric field $E(x, t)$. The Vlasov equation is cast as

$$ \frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} - \frac{eE}{m} \frac{\partial f}{\partial v} = 0 \qquad (1) $$

where $e$ is the elementary charge, $m$ is the electron mass, and $E$ is the electric field. Here collisions are neglected.

Self-consistency is ensured by the Maxwell-Gauss-Poisson (hereafter, Poisson for concision) equation,

$$ \frac{\partial E}{\partial x} = \frac{e}{\epsilon_0} \left[ n_0 - \int f dv \right] \qquad (2) $$
where $n_0$ is a constant density which accounts for the neutralizing background of ions (and is equal to the ion density if ions are singly-charged).

This is a well-known classical model, which has been investigated by many authors, both analytically (e.g. [24, 43, 58]) and numerically (e.g. [1, 6, 10, 21–23, 28, 42]). It is extensively used to study the propagation and evolution of Langmuir waves (electron plasma waves) and ion-acoustic waves. It describes resonant phenomena like Landau damping and beam-plasma interactions, as well as nonlinear phenomena like wave steepening, soliton formation, and phase-space structures (e.g. electron holes).

The subtleties lie in the treatment of boundaries, and in the numerical method.

## B. Plasma setup

Hereafter, the spatial extent (simulation box length) is set to $L = 100 \lambda_D$, where $\lambda_D$ is the Debye length. The equilibrium electron distribution function is a Maxwellian, centered at $v = 0$, with homogeneous density $n_0$ and temperature $T_e$ (or equivalently, thermal velocity $v_{Te} = (T_e/m_e)^{1/2}$). Both density and temperature are normalised out of the system $^1$, so that the results presented in terms of normalised quantities (e.g. time in terms of $\omega_{pet}$, where $\omega_{pe}$ is the electron plasma frequency) are valid for a large range of reasonable values of density and temperature. In this sense, to characterize the equilibrium, the only remaining input parameter is $L/\lambda_D$.

On top of this equilibrium, we introduce at $t = 0$ a monochromatic wave of amplitude $A$ and wavenumber $k$. The way this is implemented slightly depends on the method, as described in subsection IV A.

## III. NUMERICAL METHODS

In this paper, we focus on 3 types of numerical approaches, two of which take a Lagrangian point-of-view: Particle-In-Cell (PIC) and N-body types. The third type of numerical approach, called as semi-Lagrangian, takes a hybrid point-of-view: quantities are sampled on

<footnote n="1">¹ Noting normalised quantities with a tilde, density and temperature (or equivalently, thermal velocity) can be normalised out of the system by setting $\tilde{t} = \omega_{pet}t$, $\tilde{x} = x/\lambda_D$, $\tilde{v} = v/v_{Te}$, $\tilde{f} = v_{Tef}/n_0$, $\tilde{E} = \lambda_D eE/T_e$.</footnote>
a Eulerian grid, but are advected by tracing back Lagrangian trajectories. Note that we do not include any fully-Eulerian approach in our investigation. For a comparison between PIC and Eulerian approaches, we refer to [50].

A. Particle-in-Cell

The transformation of the Vlasov equation into the PIC method involves discretizing the plasma distribution function in phase-space and approximating it using a finite number of "superparticles" or "macroparticles". Each macroparticle represents a large number of real particles. The collective motion of macroparticles, and their interaction with the fields are used to simulate the behavior of the entire plasma. The connection between the Vlasov equation and the PIC method can be formalized by introducing weight functions and translating the continuous distribution function into a discrete form. The macroparticles move through a computational grid that discretizes the spatial domain. The fields are solved on this grid, and the forces acting on the particles are computed by interpolating the fields at the macroparticles' positions. The particles are then advanced in time according to these forces, and the resulting charge and current densities are used to update the fields through Maxwell's equations.

To transform the Vlasov equation into a form suitable for the PIC method, we first approximate the distribution function as a sum of weighted Dirac delta functions, each representing a macroparticle. In our simple 1D model,

$$f(x, v, t) \approx \sum_{p} w_p \delta[x - x_p(t)] \delta[v - v_p(t)] \quad (3)$$

where $p$ indexes the macroparticles, $w_p$ is the weight of the $p$-th macroparticle, $x_p(t)$ and $v_p(t)$ are the position and velocity of the $p$-th macroparticle at time $t$, and $\delta$ is the Dirac delta function.

To reduce noise, and consistently with the statistical meaning of the Vlasov equation, instead of using Dirac delta functions, a weight function $S(x - x_p)$ is introduced to spread the contribution of each macroparticle over nearby grid points. The distribution function then becomes:

$$f(x, v, t) \approx \sum_{p} w_p S[x - x_p(t)] \delta[v - v_p(t)] \quad (4)$$
The weight function $S(x - x_p)$ is typically a shape function, such as a linear or quadratic function, which determines how the charge density is distributed onto the grid from the macroparticles.

The charge density $\rho(x)$ on the grid is calculated by summing the contributions from all macroparticles, weighted by the shape function:

$$ \rho(x) = \sum_{p} q_p S[x - x_p(t)] \qquad (5) $$

where $q_p = w_p q$ is the charge of the p-th macroparticle. The Poisson equation is solved on the grid to update the electric field $E(x)$ based on the charge densities. The electric field is interpolated from the grid $x_i$ to the positions of the macroparticles using the shape function,

$$ E_p = \sum_i E(x_i) S(x - x_p) \qquad (6) $$

The equations of motion for each macroparticle $p$ are solved to update their positions and velocities, based on the interpolated electric field:

$$ \frac{dx_p}{dt} = v_p \qquad (7) $$

$$ \frac{dv_p}{dt} = \frac{q_p}{m_p} E_p \qquad (8) $$

This step is called as particle pushing.

The process of charge density assignment, field solving, force interpolation, and particle pushing is repeated iteratively to simulate the time evolution of the plasma.

In summary, by using weight functions to discretize the distribution function and assign charge and current densities to the grid, the Vlasov equation is effectively converted into a set of equations which model plasma dynamics by tracking a manageable number of macroparticles and their interaction with fields on a computational grid.

The specific PIC code used in this article is relatively versatile and can be adapted to different physical models. For instance, it has been used to simulate the plasma-wall transition in strongly magnetized plasmas [36, 37] or to self-consistently calculate the temperature of a wall facing a high heat load, taking into account thermionic emission and its impact on the sheath structure [35].

In the present study, two types of boundary conditions for electrons were employed, i.e. periodic or mirror. The macroparticles were initially uniformly distributed on the grid at
the start of the calculation, with a possible modulation if necessary (in the case described in subsection IV A of initial potential wave). The electron velocities were randomly drawn from a Maxwellian distribution, with either zero drift velocity or a drift depending on position (in the case described in subsection IV A of initial velocity drift). Finally, the Poisson equation was solved iteratively with $\phi(0) = \phi(L) = 0$.

## B. Semi-Lagrangian

The semi-Lagrangian method is designed to solve the Vlasov equation by following the characteristics of the phase-space flow backward in time. The updated value of $f$ on any fixed grid point is calculated by tracing back the trajectory (characteristic) from the current grid point at time $t = t_n$ to an earlier time $t = t_{n-1} = t_n - \Delta t$. This backward tracing determines the origin of the characteristics that arrive at the grid point. In other words, although $f$ is calculated on a Eulerian grid, we use the Lagrangian fact that the distribution function is conserved along particle trajectories,

$$f[x_p(t_n), v_p(t_n), t_n] = f[x_p(t_{n-1}), v_p(t_{n-1}), t_{n-1}] \quad (9)$$

where $(x_p(t), v_p(t))$ is the Lagrangian trajectory of a virtual particle $p$ which coincides at time $t = t_n$ with a point $(x_i, v_j)$ of the Eulerian grid.

Once the origin $(x_p(t_{n-1}), v_p(t_{n-1}))$ of the characteristic is found, the distribution function $f$ at this origin point is interpolated from the grid values of $f$ at the previous time step. This interpolated value is then used to update the grid point at the current time step. The interpolation can be performed using various methods, such as linear, cubic, or higher-order schemes, depending on the desired accuracy. This approach is said to combine the strengths of Eulerian (relatively noise-free compared to PIC [50]) and Lagrangian approaches.

In previous works [28], we developed a 1D semi-Lagrangian Vlasov code, based on the Cubic-Interpolated-Propagation (CIP) scheme [40] and the splitting method [6], which enabled accurate simulations of the Vlasov-Poisson system.

All quantities such as $f$ are discretized on uniform Eulerian grids with $N_x$ and $N_v$ grid points in the $x$ and $v$ directions, respectively, covering the computational domain $\{(x, v) | 0 \le x < L, -v_{\max} \le v \le v_{\max}\}$. Periodic boundary conditions are applied in the $x$-
direction, while fixed boundary conditions are used in the v-direction. To achieve second-order accuracy in time, we employ the Strang splitting method [54] as described in [56]. Each time-step, of duration $\Delta t$, involves the following sequence of operations:

1. Perform advection according to $\partial_t f + v \partial_x f = 0$ for a duration of $\Delta t/2$.

2. Update the electric field $E$ by solving Poisson equation in Fourier space.

3. Perform advection according to $\partial_t f + (qE(x)/m) \partial_v f = 0$ for a duration of $\Delta t$.

4. Repeat step 1.

Steps 1, 3 and 4 involve solving a 1D hyperbolic advection equation of the form:

$$ \partial_t F + u \partial_s F = 0, \qquad (10) $$

where $u$ is a constant in the $s$ direction, $s$ represents a generalized advection coordinate, and $F$ is a function of $s$ and $t$.

In the present paper, we use the COBBLES code, where this advection equation is treated in a conservative form, using the CIP-CSL algorithm [61], extended to the $(x,v)$ phase-space [39]. The key idea is that, in order to keep a phase-space density flux balance between neighboring cells, the algorithm advects not only the distribution function $f$, but integrated quantities $\int f dx$, $\int f dv$, and $\int f dx dv$ as well, where integrals are between two successive grid points. This allows long-time accurate simulations [26], with minimal loss of information [30].

Hereafter, we describe simulations performed with $N_x \times N_v = 1024 \times 1024$ grid points, and a time-step-width such that CFL < 0.5, where CFL is the Courant-Friedrichs-Lewy number $CFL \equiv v_{max} \Delta t N_x/(2 L)$.

### C. N-body (infinite charged sheets)

A key assumption in the Vlasov equation is that the particles interact through mean fields, ignoring direct particle-particle collisions (collisionless regime). The N-body method, on the other hand, is a direct numerical approach that simulates the plasma dynamics by tracking the motion of individual particles. In this method, the equations of motion for each particle are solved explicitly, accounting for the forces due to interactions with all other
particles in the system. The equations of motion are the same as for the PIC method, but they are applied to each real particle rather than to macroparticles, and the main difference is in the way the force is calculated.

The N-body method is more detailed but computationally intensive. In the most basic version of this method, each of the *N* particles interacts with every other particle, making the computational cost scale as $O(N^2)$ in general. Techniques [2] like tree codes [19, 47, 51] or fast multipole methods [15, 18] are often used to reduce computational complexity.

In high-density regimes where collisions become important, extensions of the Vlasov equation, like the Vlasov-Boltzmann equation, might be used to incorporate some of the effects that are naturally captured by the N-body method. However, recent studies [12, 55, 62], based on weak turbulence calculations derived from the Klimontovich equation (incorporating all wavenumbers), have raised open questions about the collision operator. These studies propose a modification of the collision operator for the Vlasov distribution function. To clarify this issue, N-body simulations provide crucial information. For example, in a previous article [14], we estimated collision rates by conducting computer simulations using a 1D N-body plasma model, enabling comparisons between theoretical predictions and measured collision rates.

The N-body code used in this paper is the same as in the latter reference [14]. The approach is inspired by Dawson's sheet model [8], which modeled a 1D plasma using an ensemble of infinite charged sheets interacting via their self-consistent electric field. All sheets are negatively charged (representing electrons) and evolve in a neutralizing background (homogeneous positive charge density). Similarly, we model a 1D plasma as a collection of infinite plane sheets, but consider sheets of both electrons and ions, confined between two reflecting walls. The motion of these sheets is governed by Newton's laws, where the electric field is constant between the sheets, leading to uniform acceleration unless a sheet crosses a neighbor. This allows Newton's equations to be integrated analytically between two sheet-crossing times. The code computes the collision time for each sheet with its neighbor and evolves the system up to the smallest crossing time. The main computational challenge is the identification of this minimum time, which is efficiently handled using a heap-based algorithm. The scaling of the algorithm with the number of sheets $N_s$ is $N_s \log N_s$, which allows
for the simulation of systems with relatively large $N_s$ (up to the order of $10^4$) with relatively low computational cost. This is an important improvement from Dawson's approach, which was limited in terms of scalability due to $O(N_s^2)$ computational cost. Another important improvement is the way we incorporate ensemble averaging to mitigate noise, as described thereafter.

Both PIC and semi-Lagrangian approaches are based on the collisionless Vlasov equation, which, in some sense, corresponds to a limit of infinite number of particles (for a given, fixed density). Consistently, our PIC simulations typically include a large number of macro-particles, of the order of several millions. On the other hand, the N-body simulations performed here describe only a relatively small number of particles: $N_s$ in the range $10^2-10^4$. For this limited number of particles, Coulomb collisions can play important roles. This is discussed further in Appendix B.

# IV. BOUNDARY CONDITIONS

In this section, let us discuss the boundary conditions in $x$, $v$, and $t$ (the initial conditions), for the three approaches described in Sec. III.

## A. Initial conditions

To investigate Landau damping, an initial perturbation corresponding to a monochromatic wave of amplitude $A$ and wavenumber $k$ is imposed at $t = 0$. The wavenumber $k = m k_0$, where $m$ is an integer, is always a multiple of the smallest non-zero wavenumber $k_0 = 2\pi/L$.

In the semi-Lagrangian code COBBLES, the initial distribution function is set as

$$f(x, v, t = 0) = f_M(v) \left[ 1 - A (\lambda_D k)^2 \cos(kx) \right] \quad (11)$$

where $f_M(v)$ is a Maxwellian distribution centered at $v = 0$. Substituting into Eq. (2) yields the initial electric potential,

$$\phi(x, t = 0) = \frac{e A \lambda_D^2 n_0}{\epsilon_0} \cos(kx) \quad (12)$$

Note that the normalized initial electric potential is

$$\frac{e \phi}{T_e} = A \cos(kx) \quad (13)$$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Early-times effect of an initial sine velocity drift, assuming free-streaming. Top: contours of constant distribution function. Bottom: electron density.</figcaption></figure>

so that in this sense, A is a good dimensionless measure of fluctuation amplitude.

In the PIC code, there are two ways to set the initial perturbation. The first method is to initialize macro-particles according to the distribution Eq. (11) (and according to a uniform distribution in space). Given that through Poisson equation it is equivalent to Eq. (12), we refer to this method as "initial potential wave". The second method is to impose a drift velocity,

$$v_{\text{drift}} = A v_{Te} \sin(kx) \qquad (14)$$

which corresponds to an initial distribution of macro-particles

$$f(x, v, t = 0) = f_M[v - v_{\text{drift}}(x)] \qquad (15)$$

This initial velocity "kick" induces an electron density wave and a corresponding electric wave. We refer to this method as "initial velocity drift".

If A is small enough, we can estimate early-times distribution function based on the free-streaming ($E = 0$) solution,

$$f(x, v, t) \approx f(x - vt, v, t = 0) = f_M[v - v_{\text{drift}}(x - vt)] \qquad (16)$$
<figure><img src="image_3.png" /><figcaption>FIG. 2. Early-times effect of an initial sine velocity drift: comparison with a PIC simulation with initial velocity drift ($A = 0.5, m = 5$), and with an ensemble of N-body simulations ($N_s = 8192, N_r = 8$) with the same parameters. Top: spatial distribution of density for 4 different times. Bottom: time-evolution of $m = 5$ component and zoom on the early times.</figcaption></figure>

Based on this free-streaming approximation, Fig. 1 illustrates how an initial velocity drift with $m = 3$ and $A = 0.1$ induces after a short time a density wave. Fig. 2 confirms that the early-time dynamics in both PIC and N-body simulations is consistent with the free-streaming picture. Disagreement for $\omega_{pe} t > 1$ is due to significant electric field already built
up by that time. Although we do not show the figures for concision, harmonics ($m=10, 15, 20$, etc.) are also excited, also in agreement with the free-streaming solution for $\omega_{pe} t \ll 1$, but they remain small. The $m = 10$ harmonic remains 9 times smaller than the $m = 5$ component, and the following harmonics are negligible.

In this paper, we mention for each PIC simulation which of the two methods are used: initial potential wave, or initial velocity drift.

In the N-body code, it is always the latter method, initial velocity drift, that is applied: the drift velocity Eq. (14) is applied to the initial random distribution of sheets.

## B. Spatial boundary conditions

Two types of spatial boundary conditions are investigated.

* Periodic boundary conditions: any $x + L$ is the same point as any $x$.

* Mirror boundary conditions: particles are reflected at two walls $x = 0$ and $x = L$.

The concept of our N-body code imposes mirror conditions, while the implementation of our semi-Lagrangian code COBBLES imposes periodic boundary conditions. Our PIC code, on the other hand, can be switched between both, which is useful to compare the behavior of a same plasma under these two boundary conditions.

We conducted 2 PIC simulations with the same initial perturbation ($m = 5$, initial velocity drift $A = 0.5$) but different boundary conditions. Fig. 3 shows the location in phase-space of PIC macro-particles near the resonant velocity $v = \omega/k$, side-by-side between periodic boundary conditions and mirror boundary conditions. At $t = 20 \omega_{pe}^{-1}$, there are some noticeable differences between the periodic case and the mirror case, however they may be attributed to an effect of chaos, given the random initial positioning of macro-particles. There is no evidence from this data that the periodic or mirror conditions have any significant impact on overall (rather than individual) phase-space dynamics in this case where there is a kind of "anti-symmetry" between $v > 0$ and $v < 0$ half-planes (i.e. the only difference is the direction of propagation).

Fig. 4 shows the time-evolution of the main velocity wave for both simulations. The discrepancy remains below 1.5%. This confirms that the periodic or mirror conditions do
<figure><img src="image_4.png" /><figcaption>FIG. 3. Snapshot of the location in phase-space of macro-particles near the resonant velocity $v = \omega/k$, in a PIC simulation with typical parameters ($m = 5$, initial velocity drift $A = 0.5$). Left: periodic boundary conditions. Right: mirror boundary conditions. Top: $t = 0$, center: $t = 10 \omega_{pe}^{-1}$, bottom: $t = 20 \omega_{pe}^{-1}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation in two PIC simulations with typical parameters ($m = 5$, initial velocity drift $A = 0.5$) and either periodic or mirror boundary conditions.</figcaption></figure>

not have any significant impact on overall dynamics, at least for the timescales we are investigating in this paper.

The fact that periodic and mirror solutions are essentially the same may appear surprising at first (at least they did surprise us). However, for our simple setup where the equilibrium distribution is symmetric in $v$, the linear dispersion is symmetric<fnref n="2" /> in $k$. In this case, both forward- and backward-propagating waves are solutions with the same damping rate, resulting in an overall damped standing wave. Even the nonlinear solution appears to be a sum of standing waves. Then, it is intuitive to realize that a standing solution would not depend on periodic or mirror solutions.

Let us think about the problem from the other point-of-view – that of propagating waves. First, consider a case with only one wave, propagating forward. For example, if a peak of the forward-propagating wave arrives to a right boundary, coming from the left, it will lead in the case of periodic boundaries to a peak propagating forward from the left boundary, and in the case of mirror boundaries to a peak propagating backward from the right boundary. In that case we expect the type of boundary to matter. However, since there is a similar wave, propagating backward, then the following happens simultaneously: a peak of the

<footnote n="2">² With the notations of Eq 18, we have $D(k,p) = D(-k,p)$, as can be obtained with a change of variables $v \to -v$.</footnote>
backward-propagating wave arrives to a left boundary, coming from the right, which leads in the case of periodic boundaries to a peak propagating backward from the right boundary, and in the case of mirror boundaries to a peak propagating forward from the left boundary. Generalizing, the mirror version of the backward-propagating wave plays the role of the periodic version of the forward-propagating wave, and vice-versa. Overall, this explains again why the solution should not depend on periodic or mirror boundaries.

### C. Velocity boundary conditions

In the velocity direction, there is in practice no boundary for the PIC and N-body approaches. Macro-particles and charged sheets are free to take any velocity value (there is a limitation due to floating point representation, but in double precision that limit is never reached in practice). For the semi-Lagrangian approach, the distribution function is computed within a bounded range $\{-v_{max} \le v \le v_{max}\}$. However, we always choose $v_{max}$ large enough (many thermal velocities) that the values of the distribution function near the velocity boundaries are negligible. In this sense, velocity boundaries are not a limitation for the semi-Lagrangian approach either. Therefore, we do not expect, nor did we observe, any difference between the 3 approaches stemming from boundary conditions in the velocity direction.

## V. LANDAU DAMPING AND NONLINEAR REGIME

Before describing novel results, let us review known facts about Landau damped Langmuir waves. For weak-enough perturbation, linear theory shows that the Fourier component $m$ of the electric potential, $\hat{\phi}_m$, behaves as

$$ \hat{\phi}_m = C \exp(pt) \qquad (17) $$

for large enough $t$, where $C$ is a complex constant, and $p = \gamma - i\omega$ is the complex root of a dispersion relation $D(k, p) = 0$ with the largest real value (least negative $\gamma$). This is true for $t \gg |\gamma|^{-1}$. For $t \ll |\gamma|^{-1}$, it is not necessarily the case, as several poles in the complex $p$-plane can significantly contribute to $\hat{\phi}_m$.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Time-evolution of the absolute value of the real part of Fourier component $m = 5$ of potential fluctuation, for various initial amplitudes of potential wave $m = 5$, in semi-Lagrangian (COBBLES) simulations. To facilitate comparison, each simulation is renormalised so that they all start from the value 1. Dashed curve: $|\cos(\omega t)|\exp(\gamma t)$, renormalised to fit the simulation result for $A \to 0$ and $t \sim |\gamma|^{-1}$. Inset: zoom on the second “rebound”.</figcaption></figure>

The linear dispersion relation is

$$D(k, p) = 1 - \frac{\iota \omega_{pe}}{k n_0} \int_P \frac{\partial f_0(v)}{\partial v} \frac{dv}{p + \iota k v} = 0 \quad (18)$$

where $P$ is a Landau-modified integral path, and $f_0(v)$ is the analytic continuation in the complex $v$-plane of the initial velocity distribution. This dispersion relation can be solved analytically for limiting cases, and numerically in general. Here, to avoid introducing any unnecessary approximation, we solve it numerically by locating the zeros of $D(p)$ in the complex $p$-plane. With $L = 100 \lambda_D$, for mode $m = 5$, on which we focus hereafter, we obtain $\omega = 1.17667 \omega_{pe}$ and $\gamma = -0.0176382 \omega_{pe}^3$.

We performed a series of semi-Lagrangian (COBBLES) simulations for various amplitudes $A$ of an initial monochromatic wave $m = 5$. Fig. 5 shows the time-evolution of the

<footnote n="3">³ Note that $\omega$ is close to the well-known long-wavelength approximation $\omega/\omega_{pe} \approx 1 + 3(\lambda_D k)^2/2 = 1.148$, but this level of agreement is fortuitous since $\omega$ is further from the finer approximation $\omega/\omega_{pe} \approx [(1 + (1 + 12\lambda_D^2 k^2)^{1/2})/2]^{1/2} = 1.1131$ from which the former approximation stems. As for $\gamma$, it is rather far from the well-known approximation $\gamma \approx \pi\omega^3/(2n_0k^2) \partial_v f_0|_{v=\omega/k} = -0.0348\omega_{pe}$.</footnote>
Fourier component $m = 5$ of potential fluctuation. As expected, as $A \to 0$, the numerical solution approaches linear theory for $t \sim |\gamma|^{-1}$. For $A = 0.03$, the difference is barely noticeable on this plot. In contrast, for $A$ of the order of $10^{-1}$ and above, nonlinear effects are evident. We can list 4 effects. In the short term, visible in particular in the first few bounces of the potential amplitude (for $\omega_{pe}t < 10$), we observe, as $A$ increases: 1. a decrease of the oscillation period; and 2. an enhancement of damping compared to linear Landau damping. In a longer timeframe (for $10 < \omega_{pe}t < 20$), we observe, as $A$ increases: 3. an increase of the apparent oscillation period (the oscillation has no longer a single, fixed frequency); and 4. a recovery of the amplitude. This behavior [31] is linked to the evolution of a BGK structure in phase-space$^4$. The evolution for further time is shown afterwards in Fig. 8.

Given these nonlinear effects as fluctuation amplitude approaches 10%, we argue that comparison with linear theory is not always enough to confirm the validity of a numerical model to simulate Landau damping. It is important to benchmark both linear and nonlinear regimes [20, 32, 34]. This is the main objective of Sec. VI.

# VI. BENCHMARK

Cross-code comparison is an essential part of the development of new algorithms, especially if they involve new numerical methods (recent example: [53]), or simply to measure the strengths and weaknesses of various approaches. Here, there is an additional objective: taking advantage of the possibilities of the PIC code to switch between periodic and mirror boundary conditions, and between initial potential wave and initial velocity drift, we aim to clarify what information can be compared, and to what accuracy, between setups with subtly different initial and boundary conditions.

## A. semi-Lagrangian vs PIC

We conducted PIC simulations with the same conditions as in the COBBLES simulations presented in Fig. 5, that is with an initial potential wave $m = 5$ ($\lambda_D k = 5 \times 2\pi/100$), with periodic boundary conditions, and with varying initial amplitudes of fluctuation.

<footnote n="4">⁴ One piece of evidence is that the time of half-turn of the phase-space structure, $T_h = \pi/\omega_b$ where $\omega_b = k(e\phi/m)^{1/2}$ is the bounce frequency, does correspond to the observed evolution of the envelope as shown in Fig. 8: $T_h/\omega_{pe} \approx 45$ for $e\phi/T_e \approx 0.05$ corresponding to the case $A = 0.1$, and $T_b/\omega_{pe} \approx 82$ for $e\phi/T_e \approx 0.015$ corresponding to the case $A = 0.03$.</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Benchmark of PIC vs semi-Lagrangian: time-evolution of the absolute value of the real part of Fourier component $m = 5$ of potential fluctuation. The vertical axis is renormalized such that fluctuation amplitude starts from unity. Note that unlike Fig. 5, the vertical scale is linear. Very large initial amplitude of potential wave, $A = 0.5$. Inset: zoom on the $120 < \omega_{pe}t < 150$ section, with additional curve (green, bold) corresponding to a PIC simulation with 10 times more macroparticles (200 million).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Convergence test on the spatial step width $\Delta x$ in PIC simulations. Two typical peaks around $\omega_{pe}t \approx 60$ – instead of the whole time-evolution – are shown. The legend indicates the value of $\Delta x/\Lambda_D$. Here the Fourier component of the density wave is shown, but this is equivalent to that of the potential wave.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8. Benchmark of PIC vs semi-Lagrangian: time-evolution of the absolute value of the real part of Fourier component $m = 5$ of potential fluctuation. The vertical axis is renormalized such that fluctuation amplitude starts from unity. Note that unlike Fig. 5, the vertical scale is linear. Top: large initial amplitude of potential wave, $A = 0.1$. Bottom: smaller initial amplitude, $A = 0.03$ (dashed: renormalised $A = 0.1$ for comparison).</figcaption></figure>
Fig. 6 shows the time-evolution of the $m = 5$ component of potential fluctuation for both PIC and COBBLES simulations, for large initial amplitude $A = 0.5$. For $\omega_{pe}t < 100$, there is accurate quantitative agreement, including basic linear behavior and nonlinear effects (in terms of frequencies, phases, and amplitudes). The agreement starts to deteriorate after about 20 oscillations. Increasing the number of PIC macroparticles from 20 to 200 millions postpones the time where discrepancies grow, but only by about 3 oscillation periods.

Here, the PIC simulations use 1000 spatial grid points on the box of length $L = 100\Lambda_D$. In other words, the spatial step width is $\Delta x = 0.1\Lambda_D$. To check that this is small enough, we perform a test of convergence. Fig. 7 shows two typical peaks around $\omega_{pe}t \approx 60$ of density fluctuations. The discrepancy between the case $\Delta x = 0.1\Lambda_D$ and the case $\Delta x = 0.05\Lambda_D$ remains below 1%, which confirms that the simulation is converged in terms of spatial step width.

Fig. 8 shows the time-evolution of the $m = 5$ component of potential fluctuation for both kinds of simulations again, but for smaller initial amplitudes: $A = 0.1$ and $A = 0.03$. (These are less stringent tests that for $A = 0.5$, because nonlinear particle trapping plays less significant roles). We stop the PIC simulations around $\omega_{pe}t < 60$, which is enough to observe the nonlinear behavior. This figure confirms that the quantitative agreement for $\omega_{pe}t < 60$ remains of very good quality for lower amplitudes.

## B. N-body vs PIC

A straightforward approach to test the capability of an N-body code to reproduce collisionless Landau damping is to increase the number of particles as much as possible. This aims to approach the collisionless regime, where collective effects dominate over collisional ones due to the small number of particle-particle interactions, in order to compare simulation results against linear theory. Calculation cost limits how high the number of particles, or sheets in our case, can be increased. However, Ref. [14] confirms that for $N_s = 8192$, which is tractable, collisional effects are negligible compared to Landau damping of reasonably sized ($\lambda_D k \ll 1$) waves.

The problem with this approach is that a wave, whose amplitude must be small enough for linear theory to be valid, is easily drowned in the noise due to the limited number of particles. Fig. 9 shows the time-evolution of re-normalized amplitude of the $m = 5$ wave for
<figure><img src="image_10.png" /><figcaption>FIG. 9. N-body ($N_s = 8192$) simulations compared to linear theory: time-evolution of re-normalized amplitude of the $m = 5$ wave for varying initial amplitudes of velocity drift. Ensemble averaging is not applied; each curve is an individual simulation. The dashed curve corresponds to collisionless linear theory.</figcaption></figure>

$N_s = 8192$ and varying initial amplitudes, against collisionless linear theory. None of the 4 simulations matches linear theory, as can be seen for example by focusing on the end of the second oscillation, near $\omega_{pe} t = 11$. The 3 simulations for $A \le 0.2$ suffer too much from noise, and the 3 simulations for $A \ge 0.1$ are too far from the linear regime. Therefore there is no range of $A$ where a match with linear theory can be confirmed.

A method to extract low-amplitude waves by averaging ensembles of noisy N-body simulations is described in appendix A. Signal-to-noise ratio can be increased by an order-of-magnitude. However, it turns out that, when $A$ is small enough to neglect nonlinear effects, this is still not enough to confidently make accurate comparisons between simulation results and linear theory.

Therefore, to test the capability of our N-body code to accurately simulate Landau damping, we switch from a verification against linear theory, to a benchmark against PIC simulations with mirror boundary conditions, which can be done for higher values of $A$.

Though collisionless PIC and collisional N-body simulations are fundamentally different, we expect them to match for a large enough number $N_s$ of sheets and early enough times
<figure><img src="image_11.png" /><figcaption>FIG. 10. Benchmark of N-body ($N_s = 8192$, $N_r = 8$) vs PIC: time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation for large initial amplitude of velocity drift, $A = 0.5$.</figcaption></figure>

– before collisional effects had time to affect the evolution. As is described in Appendix B, collisional effects are much slower than collisionless damping for $N_s > 500$. Here, we compare a PIC simulation and an ensemble of $N_r = 8$ N-body simulation with a relatively large number of sheets $N_s = 8192$, both with mirror boundary conditions and initial velocity drift of amplitude $A = 0.5$. Fig. 10 shows the time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation. PIC and N-body solutions agree very well for $t < 18 \omega_{pe}^{-1}$. We assume that the discrepancy which gradually appears after $t = 18 \omega_{pe}^{-1}$ can be attributed to Coulomb collisions.

Fig. 11 shows how the solution depends on the number of sheets, albeit for a smaller amplitude of initial perturbation, and for shorter time. The solutions for $N_s = 512$, 1024 and 2048 are very similar. This indicates that $N_s = 8192$ is enough to be well into the collisionless limit.
<figure><img src="image_12.png" /><figcaption>FIG. 11. Time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation for moderate initial amplitude of velocity drift, $A = 0.1$, and for varying number of sheets $N_s$ (and $N_r \ge 80$.</figcaption></figure>

VII. CONCLUSION

The 1D Vlasov-Poisson model is a powerful tool in plasma physics, providing essential insights into wave dynamics, particle interactions, and instabilities in a wide range of physical systems, from laboratory plasmas to astrophysical environments. In this paper, we conducted a benchmark study comparing three numerical methods – Particle-In-Cell (PIC), semi-Lagrangian, and N-body simulations – for the analysis of Landau-damped Langmuir waves in a 1D Vlasov-Poisson system. Although semi-Lagrangian and N-body simulations have different boundary conditions, and slightly different initial conditions, they can be bridged via PIC simulations which handle both types of boundary conditions and both types of initial conditions.

The results showed that, for $\omega_{pe}t < 100$, our PIC code and the semi-Lagrangian COBLES code produce near-identical solutions for both linear and nonlinear Landau damping in collisionless plasmas. Specifically, for initial wave amplitudes $A \le 0.03$, the damping rates observed in PIC and semi-Lagrangian simulations were consistent with linear Landau damping theory. Nonlinear effects, such as oscillation frequency shifts and amplitude recovery, were observed for higher amplitudes $A \ge 0.1$, aligning with theoretical predictions of phase-space structure evolution.

The N-body method, which models charged particles as infinite sheets and explicitly
tracks particle-particle interactions, was shown to accurately reproduce Landau damping in the collisionless regime when the number of sheets $N_s$, was sufficiently large (typically $N_s > 500$). The N-body simulations exhibited much higher noise levels compared to PIC and semi-Lagrangian methods, but the method presented in Appendix A is successful in extracting the relevant information, namely the time-evolution of a plasma-velocity wave.

The quantitative agreements obtained in the two benchmarks (semi-Lagrangian vs PIC; and N-body vs PIC) serve to give confidence in the 3 codes, which have been developed separately by 3 different persons or groups of persons.

Both PIC and Vlasov codes are quite efficient at capturing the physics of Landau damping, offering reasonable computation times and good time resolution compared to the N-body approach. However, the N-body approach has the advantage of accurately describing Coulomb collisions (from first principles), as is described in Appendix B.

The exploration of boundary conditions (periodic vs. mirror) and initial perturbation types (potential wave vs. velocity drift) underscores the importance of setup-specific effects on simulation outcomes. These findings are not only relevant to 1D electrostatic systems but also provide a foundation for benchmarking and improving numerical methods in more complex scenarios, such as 2D/3D electromagnetic plasmas.

## ACKNOWLEDGMENTS

This work was funded by the Institut Universitaire de France (chaire Fondamentale, Maxime LESUR), and by the Agence Nationale de la Recherche for the Project GRANUL (ANR-19-CE30-0005). This work was granted access to the HPC resources of EXPLOR (Project No. 2017M4XXX0251). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
<figure><img src="image_13.png" /><figcaption>FIG. 12. Time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation in N-body simulations with $N_s = 8192$, for varying initial amplitudes of velocity drift. The equivalent COBBLES simulations are shown for comparison, though boundary and initial conditions differ. Note that the vertical axis is renormalized to $A$ so that all solutions are expected to start from the value 1 at $t = 0$.</figcaption></figure>

## Appendix A: Methodology for measuring waves in N-body simulations

The output of one N-body simulation, at a given time $t$, is the list of location in phase-space of all $N_s$ sheets. The first step in our method is to group this data in $N_h$ boxes, regularly spaced, to obtain the sheet density $n_e$ and the average sheet velocity $v_e$ as functions of $x$. We choose a number of boxes that is not a multiple of the mode number $m$ under consideration, to avoid an artefact due to aliasing. For $m = 5$, we arbitrarily focused on multiples of 11 and 17 ($N_h = 11, 17, 22, 33, 34, 44, 51, 55, ...$). $N_h$ must be large enough to identify without ambiguity a given $m$-components in the Fourier spectrum, but small enough that a box
<figure><img src="image_14.png" /><figcaption>FIG. 13. Histogram ($N_h = 44$) of sheet velocity $v_e(x)$ at $t = 0$, with $N_s = 1024$ sheets and initial amplitude of velocity drift $A = 0.1$, without ($N_r = 1$) and with ensemble average ($N_r = 16, 32, 80$).</figcaption></figure>

includes a lot of sheets, thus mitigating noise due to the random initialisation. Scanning different values of $N_h$, we found $N_h = 88$ to be a good compromise.

Both $n_e$ and $v_e$ are then transformed to Fourier space using a FFT algorithm. Alternatively, we use a direct sine transform. It turns out that $v_e$ shows clearer trends compared to $n_e$. Fig. 12 shows the time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation in N-body simulations with $N_s = 8192$, for varying initial amplitudes. For this large number of sheets, and for high fluctuation amplitude $A = 0.5$, the data appears to not suffer much from noise. The problem with noise becomes evident as $A$ decreases. For $A = 0.03$, the noise level is such that even at $t = 0$ the measured amplitude of the wave is twice larger than the imposed one. Although this is not illustrated here, noise causes even more evident issues as the number of sheets is reduced.

Therefore, in general the histogram obtained from one simulation is too noisy for a meaningful analysis. To solve this issue, we group in one histogram a number $N_r$ of realisations of the same simulation (where only the seed of the random number generator used to initialize sheets was varied). When the number of sheets $N_s$ is small, a large number $N_r$ of simulations is required to reach a low level of noise. Typically, we use $N_r = 10$ for $N_s = 8192$, and up to $N_r = 640$ for $N_s = 128$. Fig. 13 shows an example of the velocity histogram at $t = 0$, without and with ensemble averaging. The $m = 5$ wave becomes evident as $N_r$ approaches
<figure><img src="image_15.png" /><figcaption>FIG. 14. Initial relative error in N-body simulations (with initial velocity drift $A = 0.1, m = 5$) against the number of realisations $N_r$ for various numbers of sheets $N_s$. Top: relative error in $v_e(x)$, accumulated in $x$. Bottom: relative error in the norm of $m = 5$ Fourier component. The legend is shared.</figcaption></figure>

~ 10. When the number of sheets $N_s$ is small, a large number $N_r$ of simulations is required to reach a low level of noise. Typically, we use $N_r = 10$ for $N_s = 8192$, and up to $N_r = 640$ for $N_s = 128$.

We can quantify the signal-to-noise ratio by focusing on initial time $t = 0$, where we expect the histogram of $v_e$ to follow $v_{\text{drift}} = Av_{Te} \sin(kx)$ (and $|\hat{v}_{e,5}| = A$) in the limit of infinite sheets. Fig. 14 shows the relative error in $v_e$, defined as the average over $x$ of
<figure><img src="image_16.png" /><figcaption>FIG. 15. Time evolution of mode $m = 5$ of velocity fluctuations in N-body simulations (with initial velocity drift $A = 0.1, m = 5$), for various numbers of realisations $N_r$, and for various numbers of sheets $N_s$.</figcaption></figure>

$|v_e - v_{\text{drift}}|$, normalized to $A$. As expected, the error decreases following a square root law with increasing number of realisations. Deviations from the square root law, and crossing of some curves, are due to mere randomness. For $N_s = 128$, averaging over 256 realisations reduces the relative error from 447% to 32%. The latter error is still quite large, but it is the result of accumulating absolute values. In terms of the wave itself, which is what we aim to measure, the relative error is much smaller. Fig. 14 includes the relative error in $\hat{v}_{e,5}$, defined as $|A_5 - A|/A$, where $A_5 = |\hat{v}_{e,5}|$. Here, the typical error after a reasonable number of realisations is less than 10%. The few points with very small relative error result from mere coincidences.

Finally, Fig. 15 shows the time evolution of mode $m = 5$ of velocity fluctuations in N-body simulations (with initial velocity drift $A = 0.1, m = 5$), for various numbers of realisations
<figure><img src="image_17.png" /><figcaption>FIG. 16. Time-evolution of the absolute value of the imaginary part of Fourier component $m = 5$ of velocity fluctuation in ensemble-averaged N-body simulations with initial amplitude of velocity drift $A = 0.1$, and varying number of sheets $N_s$ (noted $N$ in the legend).</figcaption></figure>

$N_r$. Each sub-figure corresponds to a different number of sheets $N_s$ from 128 to 1024. The smoothness of the curve for high-enough number of realisations can be interpreted as further confirmation of sufficient statistical quality of the data.

## Appendix B: Collisionless and collisional Landau damping

Here we describe research that was enabled by comparing results between N-body and COBBLES simulations, which was first reported in Ref. [14], and we extend on it by providing comparison with PIC simulations. The benchmark presented in Sec. VI, which was actually completed before we submitted the latter reference, served to build confidence in the simulation codes and in the methodology. Hereafter we describe the methodology in more details.

The goal was to measure collisional damping in N-body simulations by subtracting collisionless Landau damping from measured total damping. The initial perturbation is a
<figure><img src="image_18.png" /><figcaption>FIG. 17. Damping rate measured in N-body simulations with varying number of sheets and two choices of mode number *m*. Two horizontal lines show the values of damping measured in corresponding COBBLES simulations.</figcaption></figure>

velocity drift with either $m = 5$ or $m = 6$. By trial and error, we found that the choice $A = 0.1$ gave a good compromise between the difficulty of extracting the wave from the noise for small $A$, and the complications due to nonlinear effects for large $A$ (In the sense that for this value $A = 0.1$, the measured value of damping agrees between PIC, COBBLES, and N-body simulations with high-enough $N_s$).

Fig. 16 shows the time-evolution of the imaginary part of Fourier component $m = 5$ of velocity fluctuation in ensemble-averaged N-body simulations (method described in appendix A) with initial amplitude $A = 0.1$, and varying number of sheets $N_s$. Qualitatively, the enhancement of damping as $N_s$ decreases is clearly visible in this figure. Quantitatively, to measure the total damping rate more accurately, we further increase the number of realisations.

Our method was finally applied to a collection of more than 2000 N-body simulations. Fig. 17 shows the damping rate measured in N-body simulations with varying number of sheets and two choices of mode number $m$. The figure includes the damping rate measured in PIC simulations (initial velocity drift, $A = 0.1$, mirror boundary conditions), and in
COBBLES simulations (initial potential wave, periodic boundary conditions, and *A* adjusted such that velocity fluctuations are similar).

A simplified version was presented in Fig. 7 of Ref. [14]. Here, additionally, we include the comparison with PIC simulations (which was done, after publication of the latter reference, for the present paper), and for both N-body and PIC results we include vertical bars to indicate the uncertainty, estimated from an ensemble of plausible fits of the local maxima of fluctuation amplitude to $\exp(\gamma t)$. The additional comparison with the PIC code is important because it removes ambiguity due to different initial conditions and boundary conditions. We can confirm that for $N_s$ of the order of $10^3$ and above, the effect of Coulomb collisions appears to be negligible compared to Landau damping. The consequent separation of Landau and collisional timescales is what makes possible the benchmark described in subsection VIB. For smaller values of $N_s$, the difference between measured damping and the corresponding horizontal line can be attributed to Coulomb collisions. In Ref. [14], this technique was used to obtain the rate of collisional damping as a function of the plasma parameter *g*.

[1] TP Armstrong, RC Harding, G Knorr, and D Montgomery. Solution of vlasov's equation by transform methods. Technical report, Univ. of Kansas, Lawrence, 1970.

[2] Jasjeet Singh Bagla. Cosmological n-body simulation: Techniques, scope and status. *Current Sc.*, 88(7):1088-1100, 2005.

[3] H. L. Berk and B. N. Breizman. Saturation of a single mode driven by an energetic injected beam. iii. Alfvén wave problem. *Phys. Fluids B*, 2(9):2246-2252, 1990.

[4] R. H. Berman, D. J. Tetreault, and T. H. Dupree. Simulation of phase space hole growth and the development of intermittent plasma turbulence. *Phys. Fluids*, 28:155-176, January 1985.

[5] Edwin Chacon-Golcher, Sever A Hirstoaga, and Mathieu Lutz. Optimization of particle-in-cell simulations for vlasov-poisson system with strong magnetic field. *ESAIM: Proc. and Surveys*, 53:177-190, 2016.

[6] C. Z. Cheng and G. Knorr. The integration of the vlasov equation in configuration space. *J. Comput. Phys.*, 22(3):330-351, November 1976.

[7] G.-H. Cottet and P.-A. Raviart. Particle methods for the one-dimensional vlasov-poisson equations. *SIAM J. Num. Analysis*, 21(1):52-76, 1984.
[8] John Dawson. One-dimensional plasma model. *Phys. Fluids*, 5(4):445-459, 1962.

[9] Blanca Ayuso de Dios and Soheil Hajian. High order and energy preserving discontinuous galerkin methods for the vlasov-poisson system, 2012.

[10] Jacques Denavit and WL Kruer. Comparison of numerical solutions of the vlasov equation with particle simulations of collisionless plasmas. *Phys. Fluids*, 14(8):1782, 1971.

[11] T. H. Dupree. Growth of phase-space density holes. *Phys. Fluids*, 26(9):2460-2481, 1983.

[12] Dominique F Escande, Didier Bénisti, Yves Elskens, David Zarzoso, and Fabrice Doveil. Basic microscopic plasma physics from n-body mechanics: A tribute to pierre-simon de laplace. *Rev. Mod. Plasma Phys.*, 2:1-68, 2018.

[13] A. Fasoli, B. N. Breizman, D. Borba, R. F. Heeter, M. S. Pekker, and S. E. Sharapov. Nonlinear splitting of fast particle driven waves in a plasma: observation and theory. *Phys. Rev. Lett.*, 81:5564-5567, Decembre 1998.

[14] Etienne Gravier, Thomas Drouot, Maxime Lesur, Alejandro Guillevic, Guillaume Lo-Cascio, Jérôme Moritz, Dominique Escande, and Giovanni Manfredi. Collision rates estimated from exact n-body simulations of a one-dimensional plasma. *Phys. Plasmas*, 30(1):012102, 2023.

[15] Leslie Greengard and Vladimir Rokhlin. A fast algorithm for particle simulations. *J. Comput. Phys.*, 73(2):325-348, 1987.

[16] Richard B Horne and Mervyn P Freeman. A new code for electrostatic simulation by numerical integration of the vlasov and ampère equations using maccottack’s method. *J. Comput. Phys.*, 171(1):182-200, 2001.

[17] Setsuo Ichimaru. Kinetic equations for turbulent plasmas. i. quasilinear theory. *Phys. Rev.*, 174(1):289, 1968.

[18] Alexander T Ihler. An overview of fast multipole methods. *MIT Area*, pages 1-20, 2004.

[19] Muhammad Isnaenda Ikhsan. Performance comparison between direct method and tree-code used in n-body simulation using python. *J. Multidisc. Ac.*, 4(2):122-126, 2020.

[20] MB Isichenko. Nonlinear landau damping in collisionless plasma and inviscid fluid. *Phys. Rev. Lett.*, 78(12):2369, 1997.

[21] Glenn Joyce, Georg Knorr, and Homer K Meier. Numerical integration methods of the vlasov equation. *J. Comput. Phys.*, 8(1):53-63, 1971.

[22] Alexander J Klimas and William M Farrell. A splitting algorithm for vlasov simulation with filamentation filtration. *J. Comput. Phys.*, 110(1):150-163, 1994.
[23] G Knorr. Zur lösung der nicht-linearen vlasov-gleichung. Z. Naturforsch, 18a(12):1304–1315, 1963.

[24] L. Landau. On the vibrations of the electronic plasma. J. Phys., 10:25–34, 1946.

[25] M Lesur, P H Diamond, and Y Kosuga. Nonlinear current-driven ion-acoustic instability driven by phase-space structures. Plasma Phys. Control. Fusion, 56(7):075005, 2014.

[26] M. Lesur, Y. Idomura, and X. Garbet. Fully nonlinear features of the energetic beam-driven instability. Phys. Plasmas, 16(9):092305–+, September 2009.

[27] M. Lesur, Y. Idomura, K. Shinohara, X. Garbet, and the JT-60 Team. Spectroscopic determination of kinetic parameters for frequency sweeping Alfvén eigenmodes. Phys. Plasmas, 17(12):122311, 2010.

[28] M. Lesur, Y. Idomura, and S. Tokuda. Kinetic simulations of electrostatic plasma waves using cubic-interpolated-propagation scheme. JAEA-Research, 2006-089:1–29, 2007.

[29] M. Lesur, K. Itoh, T. Ido, M. Osakabe, K. Ogawa, A. Shimizu, M. Sasaki, K. Ida, S. Inagaki, S.-I. Itoh, and the LHD Experiment Group. Nonlinear excitation of subcritical instabilities in a toroidal plasma. Phys. Rev. Lett., 116:015003, Jan 2016.

[30] Maxime Lesur. Method- and scheme-independent entropy production in turbulent kinetic simulations. Comput. Phys. Commun., 200:182–189, 2016.

[31] Maxime Lesur. Nonlinear features of instabilities, turbulence and transport in hot plasmas. HDR thesis, Université de Lorraine, 2020.

[32] Zhiwu Lin and Chongchun Zeng. Small bgk waves and nonlinear landau damping. Commun. Math. Phys., 306(2):291–331, 2011.

[33] Debraj Mandal and Devendra Sharma. Nonlinearly interacting trapped particle solitons in collisionless plasmas. Phys. Plasmas, 23(2):022108, 2016.

[34] Giovanni Manfredi. Long-time behavior of nonlinear landau damping. Phys. Rev. Lett., 79(15):2815, 1997.

[35] J Moritz, S Heuraux, N Lemoine, M Lesur, E Gravier, F Brochard, L Marot, and P Hiret. Thermionic emission of a tungsten surface in high heat flux plasma: Pic simulations. Phys. Plasmas, 30(8):083514, 2023.

[36] J Moritz, M Lesur, E Faudot, S Devaux, S Heuraux, and J Ledig. The plasma-wall transition with collisions and an oblique magnetic field: Reversal of potential drops at grazing incidences. Phys. Plasmas, 26(1):013507, 2019.
[37] Jérôme Moritz, Stéphane Heuraux, Etienne Gravier, Maxime Lesur, Frédéric Brochard, L De Poucques, E Faudot, and Nicolas Lemoine. Sheath size and child-langmuir law in one dimensional bounded plasma system in the presence of an oblique magnetic field: Pic results. *Phys. Plasmas*, 28(8):083501, 2021.

[38] Clément Mouhot and Cédric Villani. On Landau damping. *Acta Mathematica*, 207(1):29–201, 2011.

[39] T. Nakamura, R. Tanaka, T. Yabe, and K. Takizawa. Exactly conservative semi-lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique. *J. Comput. Phys.*, 174(1):171–207, 2001.

[40] T. Nakamura and T. Yabe. Cubic interpolated propagation scheme for solving the hyperdimensional vlasov-poisson equation in phase space. *Comput. Phys. Commun.*, 120:122, 1999.

[41] Helmut Neunzert. An introduction to the nonlinear boltzmann-vlasov equation. In *Kinetic Theories and the Boltzmann Equation: Lectures given at the 1st 1981 Session of the Centro Internazionale Matematico Estivo (CIME) Held at Montecatini, Italy, June 10–18, 1981*, pages 60–110. Springer, Berlin, Heidelberg, 2006.

[42] Jürgen Nührenberg. A difference scheme for vlasov’s equation. *Zeitschrift für angewandte Mathematik und Physik ZAMP*, 22(6):1057–1076, 1971.

[43] T. O’Neil. Collisionless damping of nonlinear plasma oscillations. *Phys. Fluids*, 8:2255–2262, 1965.

[44] C. Pagliantini, G. Delzanno, and S. Markidis. Physics-based adaptivity of a spectral method for the vlasov-poisson equations based on the asymmetrically-weighted hermite expansion in velocity space. *J. Comput. Phys.*, 488:112252, 2022.

[45] Gerhard Rein. Non-linear stability for the vlasov–poisson system—the energy-casimir method. *Math. Methods App. Sc.*, 17(14):1129–1140, 1994.

[46] Thomas Respaud and E. Sonnendrücker. Analysis of a new class of forward semi-lagrangian schemes for the 1d vlasov poisson equations. *Numerische Mathematik*, 118:329–366, 2009.

[47] Vladimir Rokhlin. Rapid solution of integral equations of classical potential theory. *J. Comput. Phys.*, 60(2):187–207, 1985.

[48] J. Rossmanith and David C. Seal. A positivity-preserving high-order semi-lagrangian discontinuous galerkin scheme for the vlasov-poisson equations. *J. Comput. Phys.*, 230:6203–6232, 2010.
[49] DD Ryutov. Landau damping: half a century with the great discovery. *Plasma Phys. Control. Fusion*, 41(3A):A1, 1999.

[50] V. Saini, S. Pandey, P. Trivedi, and R. Ganesh. Coherent phase space structures in a 1d electrostatic plasma using particle-in-cell and vlasov simulations: A comparative study. *Phys. Plasmas*, 25:092107, 2018.

[51] John K Salmon and Michael S Warren. Fast parallel tree codes for gravitational and fluid dynamical n-body problems. *Int. J. Supercomp. App. High Perf. Comp.*, 8(2):129-142, 1994.

[52] H Schamel. Theory of electron holes. *Phys. Scr.*, 20(3-4):336, 1979.

[53] Mohamad Shalaby, Avery E Broderick, Philip Chang, Christoph Pfrommer, Astrid Lamberts, and Ewald Puchwein. Sharp: a spatially higher-order, relativistic particle-in-cell code. *Astr. J.*, 841(1):52, 2017.

[54] G. Strang. On the construction and comparison of difference schemes. *SIAM J. Numerical Analysis*, 5:506-517, 1968.

[55] SF Tigik, LF Ziebell, and PH Yoon. Collisional damping rates for plasma waves. *Phys. Plasmas*, 23(6):064504, 2016.

[56] R. G. L. Vann, R. O. Dendy, G. Rowlands, T. D. Arber, and N. d'Ambrumenil. Fully nonlinear phenomenology of the berk-breizman augmentation of the vlasov-maxwell system. *Phys. Plasmas*, 10(3):623-630, 2003.

[57] Cédric Villani. Particle systems and nonlinear landau damping. *Phys. Plasmas*, 21(3):030901, 2014.

[58] A Vlasov. On the kinetic theory of an assembly of particles with collective interaction. *Russ. Phys. J.*, 9:25-40, 1945.

[59] A A Vlasov. On vibration properties of electron gas. *J. Exp. Theor. Phys.*, 8(3):291, 1938. (in Russian).

[60] H. V. Wong and H. L. Berk. Growth and saturation of toroidal Alfvén eigenmode modes destabilized by ion cyclotron range of frequency produced tails. *Phys. Plasmas*, 5(7):2781-2796, 1998.

[61] F. Xiao, T. Yabe, X. Peng, and H. Kobayashi. Conservative and oscillation-less atmospheric transport schemes based on rational functions. *J. Geophys. Research (Atmospheres)*, 107:4609, November 2002.

[62] Peter H Yoon, Luiz Fernando Ziebell, EP Kontar, and Reinhard Schlickeiser. Weak turbulence
theory for collisional plasmas. *Phys. Rev. E*, 93(3):033203, 2016.

[63] Tie Zhou, Yan Guo, and Chi-Wang Shu. Numerical study on landau damping. *Phys. D: Nonlinear Phenomena*, 157(4):322-333, 2001.