Pierre-Simon De Laplace

D F Escande

D Bénisti

Y Elskens
email: yves.elskens@univ-amu.fr

‡ D Zarzoso

F Doveil

Basic microscopic plasma physics from N -body mechanics A tribute to

Keywords: N -body dynamics, Debye shielding, Landau damping, wave-particle interaction, spontaneous emission, Coulomb scattering PACS 2010 45.50.Jf Few-and many-body systems 45.50.Tn Collisions 52.35.Fp Electrostatic waves and oscillations (e.g., ion-acoustic waves) 52.20.Fs Electron collisions

Computing is not understanding. This is exemplified by the multiple and discordant interpretations of Landau damping still present after seventy years. For long deemed impossible, the mechanical N -body description of this damping, not only enables its rigorous and simple calculation, but makes unequivocal and intuitive its interpretation as the synchronization of almost resonant passing particles. This synchronization justifies mechanically why a single formula applies to both Landau growth and damping. As to the electrostatic potential, the phase mixing of many beam modes produces Landau damping, but it is unexpectedly essential for Landau growth too. Moreover, collisions play an essential role in collisionless plasmas. In particular, Debye shielding results from a cooperative dynamical self-organization process, where "collisional" deflections due to a given electron diminish the apparent number of charges about it. The finite value of exponentiation rates due to collisions is crucial for the equivalent of the van Kampen phase mixing to occur in the N -body system. The N -body approach incorporates spontaneous emission naturally, whose compound effect with Landau damping drives a thermalization of Langmuir waves. O'Neil's damping with trapping typical of initially large enough Langmuir waves results from a phase transition. As to Coulomb scattering, there is a smooth connection between impact parameters where the two-body Rutherford picture is correct, and those where a collective description is mandatory. The N -body approach reveals two important features of the Vlasovian limit: it is singular and it corresponds to a renormalized description of the actual N -body dynamics.

I. INTRODUCTION

This review deals with the microscopic physics of plasmas, mainly collisionless ones. Its main purpose is to improve the foundations of this physics by laying out, in a pedagogic and elementary manner, a systematic exposition of its approaches using N -body classical mechanics. These approaches enable in particular: (i) a short, yet explicit and rigorous, derivation of Landau damping and growth, (ii) an associated intuitive, yet rigorous, interpretation of these phenomena, (iii) another derivation showing Landau damping to result from phase mixing, as originally proved by van Kampen in a Vlasovian setting, (iv) a third derivation recovering the usual Vlasovian dielectric function and unifying Landau damping and Debye shielding (or screening), (v) unveiling how the microscopic mechanism of Debye shielding is intimately connected to collisions, (vi) a unification of the derivations of spontaneous emission and of Landau damping, (vii) proving the depletion of nonlinearity when there is a plateau in the tail distribution function, (viii) proving that damping with trapping results from a phase transition, (ix) a calculation of Coulomb scattering describing for the first time correctly all impact parameters with no ad hoc cut-off.

Item (ii) is important, since the lack of mechanical interpretation prevented the plasma community from accepting the reality of Landau damping from its publication in 1946 [START_REF] Landau | On the vibrations of the electronic plasma[END_REF] till 1964 when Malmberg and Wharton's celebrated experiment proved its existence (see e.g. [115] for a historical sketch). Unfortunately, more than seven decades after Landau's publication, textbooks still propose discordant physical interpretations of the effect, invoking trapping, surfing, or synchronization of almost resonant passing particles. As to item (v), the new approach elucidates a longstanding mystery: how can a given particle be shielded by all other particles, while contributing to their individual shieldings? Finally, we stress that the mechanical approach of this review does not aim at challenging the Vlasovian one, as far as efficiency or convenience is concerned to perform kinetic calculations in plasma physics.

For macroscopic classical systems, the N -body description by classical mechanics was deemed impossible. This led to the development of thermodynamics, of fluid mechanics, and of kinetic equations to describe various macroscopic systems made up of particles like electrons, gas atoms or molecules, stars, or microorganisms. When plasma physicists had to address the microscopic description of their state(s) of matter, they did not consider using N -body classical mechanics, but directly derived kinetic analogues of the Boltzmann and Liouville equations, in particular Vlasov with his celebrated equation. This trend has been dominant till nowadays. However, the theoretical approach presented in this review shows that, actually, N -body classical mechanics is a practical tool both for solving important linear and nonlinear problems in microscopic plasma physics: Laplace's dream 1 was not a mere utopia, since the calculation of classical orbits starting from prescribed initial conditions can genuinely describe and explain non-trivial aspects of the macroscopic dynamics of a many-body system. This is the reason why this review is "A tribute to Pierre-Simon de Laplace". A second reason is the important role of the Laplace transform in the present N -body approach.

While deemed impossible to use, N -body classical mechanics is generally considered as the ultimate reference for the description of the microscopic physics of classical plasmas. However, if attempted, the corresponding ultimate description may a priori run into difficulties. Indeed, as shown in Chibbaro, Rondoni, and Vulpiani's book "Reductionism, Emergence and Levels of Reality" [START_REF] Chibbaro | Reductionism, emergence and levels of reality[END_REF] with various examples in statistical mechanics, in chaotic systems, and in chemistry, the relations between different levels of description of physical phenomena are not simple. As stated by Michael Berry in the foreword to this book, "It is far from straightforward to derive the formula relating the object and image of a simple lens by starting from the field operators of quantum optics supposedly the deepest of our current pictures of light. [...] The resolution of these difficulties starts from the observation that the theories of physics are mathematical, and relations between them involve limits as some parameter vanishes: wave optics 'reduces to' geometrical optics when the wavelength is negligibly small, quantum physics 'reduces' to classical physics when Planck's constant can be neglected, etc. Therefore understanding relations between levels must involve the study of limits, that is, mathematical asymptotics. And the central reason why 'reduces to' is so problematic is the fact that the limits involved are usually singular." Plasma physics would be exceptional if it would escape this problem. It does not. 1 The sentence defining what was called afterwards Laplace's demon is well known : "Une intelligence qui, pour un instant donné, connaîtrait toutes les forces dont la nature est animée et la situation respective des êtres qui la composent, si d'ailleurs elle était assez vaste pour soumettre ces données à l'analyse, embrasserait dans la même formule les mouvements des plus grands corps de l'univers et ceux du plus léger atome : rien ne serait incertain pour elle, et l'avenir, comme le passé, serait présent à ses yeux." [Essai philosophique sur les probabilités (1814) [START_REF] Laplace | Essai philosophique sur les probabilités[END_REF]] English translation : "An intellect, which at a certain moment would know all forces that set nature in motion, and all positions of all items of which nature is composed, if this intellect were also vast enough to submit these data to analysis, it would embrace in a single formula the movements of the greatest bodies of the universe and those of the tiniest atom; for such an intellect nothing would be uncertain and the future just like the past would be present before its eyes." However, the genuine Laplace's dream is reasonable, since a few sentences later he states: "Tous ses efforts dans la recherche de la vérité tendent à rapprocher [l'esprit humain] sans cesse de l'intelligence que nous venons de concevoir, mais dont il restera toujours infiniment éloigné." English translation : "All its efforts in the quest of truth tend at moving [the human spirit] closer to the intelligence we have just conceived, but from which it will always stay infinitely distant".

Indeed, this issue is already present in the derivation of the kinetic equations of plasmas. Deriving the Vlasov equation by the BBGKY hierarchy assumes the two-point correlation function to be small. Because of the singularity of the Coulomb force at vanishing distances, this is not true when two particles come close to each other, and the singularity must be cured by imposing a short-range cut-off, but this works only for a uniform plasma, without waves in particular (see section 5.2 of [START_REF] Nicholson | Introduction to plasma theory[END_REF]). The usual way to obtain a vanishing correlation function is the "pulverization procedure" : one cuts each particle into M equal pieces, and one lets M and the number of particles in the Debye sphere (defined in section II A) go to infinity (see section 4.3 of [START_REF] Nicholson | Introduction to plasma theory[END_REF]). However, this entails the loss of the actual correlations provided by the genuine two-point correlation function of the plasma of interest. Furthermore, the issue of the short-range cut-off is important to evaluate Coulomb scattering. The general practice is to take this cut-off equal to the classical distance of closest approach for the lower range of temperatures, and the quantum uncertainty on the electron's position for temperatures high enough for the latter to be larger than the former.

The singularity of the Coulomb force at vanishing distances is also a problem (see [118] and references therein for a discussion) in the mean-field derivation of the Vlasov equation [START_REF] Boers | On mean field limits for dynamical systems[END_REF][START_REF] Braun | The Vlasov dynamics and its fluctuations in the 1/N limit of interacting classical particles[END_REF][START_REF] Dobrushin | Vlasov equations[END_REF][START_REF] Elskens | Vlasov equation and N -body dynamics : How central is particle dynamics to our understanding of plasmas ?[END_REF][START_REF] Jabin | Mean field limit and propagation of chaos for Vlasov systems with bounded forces[END_REF]101,[START_REF] Neunzert | Die Approximation der Lösung von Integro-Differentialgleichungen durch endliche Punktmengen[END_REF]118]. This derivation requires a smoothing of this singularity at a length scale about N -α , for some α > 0 [START_REF] Hauray | Particle approximation of Vlasov equations with singular forces : Propagation of chaos[END_REF][START_REF] Kiessling | The microscopic foundations of Vlasov theory for jellium-like Newtonian Nbody systems[END_REF]. However, even so, the equation is proved to be accurate only over a time of the order of the inverse of the largest Lyapunov exponent of the N -body system. This derivation of Vlasov equation deals with a distribution describing a single realization of the plasma, while the BBGKY derivation deals with a function describing an ensemble of plasmas.

These two cases indicate that the issue of singular limits is likely to come about in approaches using N -body classical mechanics too. It does, as we will see.

This paper was initially meant as a mere review of results already published on the description of microscopic plasma physics by N -body mechanics. However, this endeavour rapidly ran into a difficulty: a part of these results were derived for three-dimensional plasmas described as sets of particles coupled by the Coulomb force, while another part came from a one-dimensional waveparticle description. This induced an artificial separation in the presentation of concepts and opposed logical continuity. This was an incentive to reformulate some of the one-dimensional results in three dimensions. When doing so for the proof of the average synchronization of particles at work in Landau damping, it proved useful to extend somewhat the calculation with respect to the one-dimensional case, to make it more intuitive. Unexpectedly, this eventually led to a very simple and rigorous description of Langmuir waves analogous to that published by Kaufman in a Vlasovian setting in 1972. Because of its simplicity, this new description is the first presented in this review 2 . Paradoxically, going to the more realistic three dimensions brought a simplified theory! The new description is followed by a reformulation in three dimensions of a previous van Kampen-Dawson-like description of Langmuir waves. This review is organized as follows. Section II introduces the basic scales, the equations of motion of the present N -body approach, and a class of granular distributions close to being spatially uniform. Then, the first half of this review provides three different derivations of Landau damping.

First, section III deals with Langmuir waves by using a technique introduced by Kaufman in a Vlasovian setting, which provides a very short, though rigorous, derivation of Landau damping together with the interpretation of this effect as an average synchronization of particles. Then, section IV introduces a fundamental equation for the electrostatic potential which is reminiscent of the one obtained by the calculation à la Landau starting with the Vlasov equation. This equation is the basis of the next two derivations of Landau damping. First, section V deals with Langmuir waves in a way reminiscent of both van Kampen's and Dawson's works: phase mixing is at work in Landau damping. Second, section VI takes a singular limit of the fundamental equation for the electrostatic potential; it recovers simultaneously Debye shielding and the classical expression for Langmuir waves obtained by dealing à la Landau with Vlasov equation. Section VII shows the Vlasovian limit to be a singular one providing a renormalized description of the plasma.

The second half of this paper reviews a series of previously published results. In section VIII a calculation using Picard technique shows that the acceleration of a particle due to another one is mediated by all other particles. This unveils the microscopic mechanism of Debye shielding to be intimately connected to collisions, and explains how a given particle can be shielded by all other particles, while contributing to their individual shieldings. Section IX is devoted to waveparticle interaction. First one derives from the N -body dynamics a Hamiltonian describing the self-consistent motion of M waves with N particles in the tail of the velocity distribution. On this basis, a simple statistical calculation derives both spontaneous emission and Landau damping, and the corresponding quasilinear friction and diffusion coefficients. Then, the saturation of the weak warm beam instability and the dynamics with a single wave are discussed. Finally section X deals with a problem completely out of reach of a Vlasovian description: it provides a calculation of Coulomb scattering describing correctly all impact parameters b, with a convergent expression reducing to Rutherford scattering for small b. Nonlinear issues are dealt with in sections IX and X.

Finally, section XI discusses in a self-contained way the new physical picture of basic microscopic plasma physics provided by the N -body approach and the new insight into the Vlasovian limit resulting from this approach. The reader might benefit from reading this last section first.

II. ONE COMPONENT PLASMA MODEL

This section defines the basic scales and equations of the One Component Plasma model used up to and including section IX, and introduces a class of granular distributions close to being spatially uniform. Since each new charged species brings new elements of complexity to the system, for the sake of simplicity we focus on the dynamics due to electrons, while neglecting that of ions because of their higher inertia. We also neglect any magnetic effect. More specifically, we deal with the One Component Plasma (OCP) model [START_REF] Abe | Giant cluster expansion theory and its application to high temperature plasma[END_REF][START_REF] Baus | Statistical mechanics of simple Coulomb systems[END_REF]116], which considers the plasma as infinite with spatial periodicity L in three orthogonal directions with coordinates (x, y, z), and made up of N electrons in each elementary cube with volume L 3 . Ions are present only as a uniform neutralizing background, enabling periodic boundary conditions.

A. Basic scales

The number density is n = N/L 3 . It enables the definition of a fundamental length scale of the system,

λ id = n -1/3 , (1)
which is the average inter-particle distance. At this scale, the plasma looks granular. The number density also enables the definition of a fundamental frequency of the system, the plasma (angular)

frequency

ω p = e 2 n ε 0 m 1/2 , ( 2
)
where m is the electron mass, e is the proton charge, and ε 0 is the vacuum permittivity. In solid state physics, since ω p T , where T is the plasma temperature (T stands for k B T , where k B is the Boltzmann constant), it is very hard to excite oscillations at the plasma frequency; these oscillations are quantized and are called plasmons. For the plasmas of interest here, the density and temperature are such that the opposite ordering holds. Therefore, there are many plasmons already at the thermal level, and they may be described through their coherent states, i.e. classical electrostatic waves.

If the plasma has a temperature T measured in units of energy, two new fundamental length scales may be defined. First, the classical distance of closest approach

λ ca = e 2 4π ε 0 T (3)
is the minimum distance of two electrons in a Rutherford collision, since e 2 /(4π ε 0 λ ca ) = mv 2 T , where v T = [T /m] 1/2 is their thermal velocity 3 .

The second length scale involves both the density and the temperature, and is the Debye length

λ D = v T ω p = ε 0 T ne 2 1/2 . ( 4
)
This length will turn out to be the typical distance where the Coulomb potential of a point charge is shielded out (see section VI). For the plasma to be quasineutral, we impose L λ D .

We define the plasma parameter

Λ = nλ 3 D = λ D λ id 3 = λ id 4πλ ca 3/2 , (5)
which is about one fourth of the number of particles in the Debye sphere, a sphere of radius λ D .

A large value of Λ implies λ ca λ id λ D . This ordering is represented in figure 1 By definition, λ ca is the typical distance where the electrostatic potential of a pair of electrons balances their kinetic energy. Since for Λ 1, λ ca λ id , particles typically interact weakly with one another, which corresponds to weakly coupled plasmas.

B. Equations of motion

Because of the spatial periodicity of our One Component Plasma model, we write the potential created by the N particles as a Fourier series

ϕ(r) = 1 L 3 m, kmb smooth ≤1 φ(m) exp(ik m • r), (6)
where k m = 2π L m, and where we restrict the Fourier expansion to k m 's such that k m b smooth ≤ 1, where

b smooth λ D , (7)
3 Many textbooks multiply this expression by √ 2.

φ(m) = - e ε 0 k 2 m N j=1 exp[-ik m • r j (t)], (8)
where r j (t) is the position at time t of particle j acting as a source, and k m = k m . Since the selffield due to the smoothed Coulomb potential vanishes, it is not necessary to exclude self-interaction in the equations of motion.

Newton's second law of motion for particle j reads rj = e m ∇ϕ(r j ).

The N -body dynamics is defined by equations ( 6) and ( 9) written for the N particles.

C. Spatially uniform granular distribution of particles

We introduce a class of spatially uniform granular distributions, which will be perturbed in the next sections. In the case of a cold plasma, such a distribution can be obtained by setting particles with a vanishing velocity over a cubic array. For a multi-velocity distribution, we take a set of monokinetic beams where each beam is a simple cubic array of particles whose elementary cube has its edges along the three orthogonal directions with coordinates (x, y, z). Such a set is called multi-beam-multi-array. Figure 2 displays such a distribution of particles for a one-dimensional plasma.

We first focus on a given beam. It is convenient to consider the index j of its particles as a three-dimensional vector p whose each integer components run from 1 to n edge . Since its particles have the same velocity v j , in equation ( 9) combined with equation ( 8), the part of the sum due to the particles of this beam bears on exp[ik m • r j0 ] only, where r j0 = L n edge p is the initial position of particle j (L/n edge is the edge length of the elementary cube). Due to the periodicity of the r j0 's, the corresponding sum vanishes unless the three components of m are on the simple cubic lattice (n edge Z) 3 with mesh length n edge . Therefore, the sum vanishes for k m = 0, if we assume

πb smooth n edge /L > 1. ( 10
)
This condition can be satisfied provided the inter-particle distance is such that it fulfills the condition

λ id b smooth . (11)
This and condition [START_REF] Balescu | Statistical mechanics of charged particles[END_REF] imply λ id λ D , which implies a large value of Λ as a consequence of equation [START_REF] Baker | Ion-wave echoes[END_REF]. As the contribution of each beam vanishes, the total Coulomb force for this set of beams vanishes identically. Therefore, the system of beams is force-free, and their distribution is invariant in time. We notice that, for this set of beams, both collisions and Debye shielding do not work. As is shown in section VIII, the fact that these two mechanisms fail simultaneously is natural, because the former produces the latter. For kinetic/collisionless plasmas, the position distribution corresponding to this set of beams is very atypical, and Debye shielding can generically set in.

In contrast, cold plasmas take on a crystalline structure where particles vibrate about positions corresponding to the nodes of perfect lattices, and do not experience shielding.

In the following, it will be useful to consider continuous limits of multi-beam-multi-arrays with bounded velocities related to very large numbers of particles in the Debye sphere (Λ 1). We consider the beam velocity distributions to be a granular approximation of a spatially uniform smooth velocity distribution f 0 (v). One can deal in this limit with a large number of beams whose velocities are on a three-dimensional grid becoming tighter and tighter in this limit. These tight grids may be taken as cubic arrays with an orientation arbitrary with respect to the spatial frame. When dealing with wavevector k m , the calculations are simpler when the velocity cubic array has an edge parallel to k m and one of its points is at the origin of velocities. Then, we consider the beam velocity distribution to be a granular approximation of a smooth velocity distribution f 0 (v) whose integral perpendicular to k m is g(v) (we omit the index m to simplify notations). We

assume f 0 (v)d 3 v = g(v)dv = 1. When Λ increases, f 0 (v) and g(v)
are split over a growing number of beams whose velocities lie on grids with a mesh size going to 0.

III. LANGMUIR WAVES À LA KAUFMAN

This section deals with Langmuir waves by using a technique introduced by Kaufman in a

Vlasovian setting [START_REF] Kaufman | Reformulation of quasi-linear theory[END_REF]. The calculation provides a very short and rigorous derivation of Landau damping together with the intuitive interpretation of this effect5 .

A. Landau damping

We now consider the case where all particles have their initial positions slightly perturbed with respect to those of a multi-beam-multi-array; this choice will provide useful cancellations in the following. Let r j0 be the initial position of the unperturbed beam particle with index j, and v j be its velocity. Let δr j be the smooth amplitude of modulation at t = 0 of the positions of the particles belonging to the same beam as particle j, and ∆r j (t) = r j (t)-r j0 -v j t be the mismatch of particle j with respect to its ballistic position. We now assume that the beam particles have their initial positions sinusoidally modulated with the wavevector k m according to ∆r j (0) = δr j sin(k m • r j0 ), while keeping the same velocity 6 . Setting these positions in the expression of φ(l) provided by equation ( 8), shows that if the δr j 's are of order ε, because of condition [START_REF] Bekefi | Collective emission processes in plasmas[END_REF], φ(l) is of order ε 2 or higher for l = ±m. In contrast, φ(m) is 7 of order ε. This limits the number of terms to be taken into account when computing φ(m, t) to this order. However, the excitation of other wavevectors at order ε 2 is important physically, since it corresponds to spontaneous emission enabling the excitation of a wave with an initially vanishing amplitude. This emission will be described in a statistical setting in section IX B.

We consider the linearized motion of the OCP about the multi-beam-multi-array distribution, and we look for a solution of the type φ(m, t) = A(t) exp(-iωt), where ω is real and A(t) is a complex amplitude of order ε varying on the scale τ A 1/|ω| (its dependence on m is not shown to simplify the notations). Retaining only the contribution of the wave of interest, the equation of

motion of particles is rj = αk m A(t) exp[i(k m • r j (t) -ωt)] + c.c., (12)
where "c.c." means complex conjugate, and

α = ie/(mL 3 ). ( 13
)
Integrating formally twice in time the equation of motion for particle j up to first order in ε, and changing the order of integrations yields

∆r j1 (t) = δr j sin(k m • r j0 ) + αk m t 0 τ A(t -τ ) exp[i(Ω j (t -τ ) + k m • r j0 )]dτ + c.c. , (14)
6 An arbitrary phase should be present in the sine, but it is set to 0, since it is not important for the derivation.

The case with a concomitant small velocity modulation can be dealt with in a similar manner, but with longer expressions. Then a ∆vj t sin(km • rj0 + ψj) contribution must be added to ∆rj(t), where ψj and ∆vj are respectively the phase and the amplitude of modulation of the velocities of the particles belonging to the same beam as particle j. 7 The same property would hold as well for other types of interactions than Coulombian ones.

where Ω j = k m • v j -ω, and both contributions of φ(m, t) and of its complex conjugate φ(-m, t) are taken into account; we stress that the resonant component of the particle position (and velocity) is along the considered wavevector. Setting this into the definition (8) of φ(m), and approximating exp[-ik m • ∆r j (t)] by 1 -ik m • ∆r j (t) yields to order ε

A(t) = e 2ε 0 k 2 m N j=1 k m • δr j exp(-iΩ j t) -2αk 2 m ∂ ∂Ω j t 0 A(t -τ ) exp(-iΩ j τ )dτ , (15)
where we took advantage of cancellations resulting from condition [START_REF] Bekefi | Collective emission processes in plasmas[END_REF] and from the fact that beam particles sit on arrays. In the spirit of the continuous limits of multi-beam-multi-arrays introduced in section II C, we now substitute the discrete sum over particles with the integral over a smooth distribution function f (v) whose integral perpendicular to the k m axis is a function g(v) of the type introduced in section II C (this substitution is justified below equation ( 22)). This yields

A(t) = f (v) N e 2ε 0 k 2 m k m • δr(v) exp[-i(k m • v -ω)t] d 3 v -i ω 2 p g(v) ∂ ∂Ω t 0 A(t -τ ) exp(-iΩτ )dτ dv, (16)
where δr(v) is the equivalent of δr j in the continuous limit,

Ω = k m v -ω with k m = k m , g(v) = f (v)d 2 v ⊥
, with v and v ⊥ the components of v respectively parallel and perpendicular to k m . Then, a Taylor expansion of A(t -τ ) to first order in τ yields

A(t) = N e 2ε 0 k 2 m h(v) exp(-iΩt) + ω 2 p A(t)g(v) ∂ ∂Ω exp(-iΩt) -1 Ω dv -i ω 2 p Ȧ(t) g(v) ∂ 2 ∂Ω 2 exp(-iΩt) -1 Ω dv, (17)
where

h(v) = f (v) k m •δr(v) d 2 v ⊥ .
This truncated Taylor expansion makes sense since, as shown in Appendix A, A(t) is an entire function, and the successive contributions in the expansion decay like powers of (τ A k m v T ) -1 . After integrating by parts, this becomes

A(t) = N e 2ε 0 k 2 m h(v) exp(-iΩt) + ω 2 p k m A(t)g (v) 1 -cos(Ωt) + i sin(Ωt) Ω dv -i ω 2 p k 2 m Ȧ(t) g (v) cos(Ωt) -i sin(Ωt) -1 Ω dv, (18)
where g (v) and g (v) are respectively the first and the second derivative of g(v).

For t = 0, this equation defines A(0). When t grows, the first exponential and the cosine terms in equation ( 18) produce a phase mixing 8 leading asymptotically to the expression 9

A(t) = ω 2 p k m A(t)P g (v) 1 Ω dv + i πω 2 p k 2 m g ( ω k m )A(t) + i ω 2 p k 2 m Ȧ(t)P g (v) 1 Ω dv, (19)
where P stands for the Cauchy principal value, and only the dominant real and imaginary terms are kept; the other ones are at most of order (k m v T τ A ) -1 , where τ A 1/|ω| is the time scale of variation of A(t). The typical time scale for phase mixing to become strong is

τ mix = (k m v T ) -1 .
Equation ( 19) is a homogeneous, linear ordinary differential equation with constant coefficients.

Its general solution thus reads A(t) = A 0 exp(γt), where γ is the solution to the algebraic equation

r (m, ω) -i πω 2 p k 2 m g ( ω k m ) = γ i ω 2 p k 2 m P g (v) 1 Ω dv ( 20
)
where we define

r (m, ω) = 1 - ω 2 p k m P g (v) 1 Ω dv. (21)
The unique solution of the problem is obtained by requiring r (m, ω) = 0, which provides the usual definition of ω mr , the real part of the frequency in the Landau calculation of Langmuir waves.

Assuming g is regular enough (e.g. with a uniformly continuous second derivative), integrating by parts under the principal value shows that the right-hand side in equation ( 20) is -iγ ∂ ω r , and equation [START_REF] Boers | On mean field limits for dynamical systems[END_REF] reduces to γ = γ L (k m ) where

γ L (k m ) = πω 2 p k 2 m ∂ r ∂ω g ( ω mr k m ) ( 22
)
8 Phase mixing is a classical concept in the theory of kinetic plasma waves, and especially in the van Kampen-Case approach to Landau damping [START_REF] Case | Plasma oscillations[END_REF][START_REF] Van Kampen | On the theory of stationary waves in plasmas[END_REF][START_REF] Van Kampen | The dispersion equation for plasma waves[END_REF]. Intuitively, it corresponds to the idea that the integral of a rapidly oscillating function is close to zero. Mathematically, it is grounded on the fact that dνF (ν) exp(-iνt) is the Fourier transform of F (ν), which decays for large t's in various cases. This occurs in particular, at least on average, if F (ν) is an L 2 function; also if F (ν) has an integrable derivative of order at least one. Having one of such properties is natural for h(v), g (v), and g (v), especially if g(v) is analytic, as assumed in Landau's derivation of Landau damping [START_REF] Lifshitz | Landau and Lifshitz Course of Theoretical Physics 10 : Physical Kinetics[END_REF]. 9 Phase mixing works here in the following way: because of the integration over v, the first exponential and the cosine terms in equation [START_REF] Besse | On the validity of quasilinear theory[END_REF] produce the Fourier transforms of h(v), g (v) and g (v). For t large, the vanishing values of the tails of these transforms are involved, and thus neglected. For t large, in the integrals involving sin(Ωt)/Ω, this factor has non vanishing values over a vanishing domain in v. This enables extracting g ( ω km ) out of the integral. The similar contribution involving g (v) is neglected, since it is of higher order in (τAkmvT) -1 .

We stress that we do not use Plemelj formula, in contrast with what was done in reference [START_REF] Kaufman | Reformulation of quasi-linear theory[END_REF].

is the Landau damping/growth rate to lowest order in γ L (k m )(k m v T ) -1 , whose sign is 10 the one of g ( ωmr km ). The above calculation is performed to all orders in Appendix A, and yields exactly Landau's complete formula if g(v) is analytic. The latter property is not necessary for the present calculation.

In practice, for the Landau effect to be observable, ωmr km must stay close to v T . Then the relevant mixing time is close to a plasma period. If the considered multi-beam-multi-array corresponds in v to a series of adjacent equidistant monokinetic beams with a mismatch δ in velocity, the passage to a smooth distribution function in equation ( 16) makes sense provided k m δt 1. Therefore, to be observable, Landau damping requires

11 k m δ |γ L (k m )| -1 1. Also γ L (k m ) k m v T for phase
mixing to occur faster than the damping time. This condition is well satisfied in a Maxwellian plasma, when using the Bohm-Gross dispersion relation 12 . We notice that, in the present derivation of Landau damping, phase mixing corresponds to the transients related to the terms complementing those of Landau poles when inverting the Laplace transform in the usual Vlasovian derivation.

The relevance of transients, depending on the way the self-consistent electrostatic field has been generated, is further discussed in Appendix B.

Finally, since A(t) is entire, the above Taylor expansion is converging for all τ 's. Therefore we may compute the real part of the dispersion relation and the damping rate to arbitrary order in

γ L (k m )(k m v T ) -1
. This only requires g(v) to be infinitely differentiable, not analytic. The next order terms in equation [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical confirmation[END_REF] give perturbatively 13 a modified r (m, ω)

r2 (m, ω) = r (m, ω) + πγ L (k m )ω 2 p k 3 m g ( ω mr k m ), (23)
and a more precise expression of Landau damping

γ L2 (k m ) = γ L (k m )[1 - γ 2 L (k m ) 2k 2 m g ( ωmr km ) g ( ωmr km ) ]. (24)
If g(v) is a Maxwellian, this becomes

γ L2 (k m ) = γ L (k m )[1 + γ 2 L (k m ) 2(k m v T ) 2 (3 - ω 2 mr (k m v T ) 2 )] ( 25
)
10 At this point, we notice that the derivation of Landau damping in the mechanical N -body setting is more accessible to students than Landau's derivation: the number of pages is divided by three (see for instance sections 6.3 to 6.5

of [START_REF] Nicholson | Introduction to plasma theory[END_REF]), the mathematics is elementary, and there is no need to introduce Vlasov equation. 11 Because of the Floquet exponents introduced in Appendix A and computed in Appendix C, the actual condition [START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF] is somewhat stronger. 12 The Bohm-Gross dispersion relation is ω 2 = ω 2 p + 3k 2 v 2 T . 13 Using again phase mixing, but going to next order in (kmvTτA) -1 in equation [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical confirmation[END_REF]

γ L2 (k m ) = γ L (k m )[1 - (γ L (k m ) ω p ) 2 2(k m v T ) 4 ] (26)
by using the Bohm-Gross dispersion relation.

B. Average synchronization of particles with waves

After introducing Landau damping, we now consider the corresponding dynamics of particles.

To this end, we compute the sum of the rj 's for the particles considered at the beginning of the previous subsection whose position is slightly perturbed with respect to monokinetic beams. This sum vanishes, since the system is isolated (multiplying this sum by m corresponds to the time derivative of the total momentum of the system, which is conserved). Using equation ( 14) to express r j (t) in equation [START_REF] Bénisti | Origin of diffusion in Hamiltonian dynamics[END_REF] to second order in A, and changing the order of integrations, we obtain

0 =-αA(t)k m N j=1 k m • δr j exp(iΩ j t) + c.c. + N j=1 2i|α| 2 k 2 m A(t)k m t 0 τ A * (t -τ ) exp(iΩ j τ ) dτ + c.c. , (27)
where A * (t) is the complex conjugate of A(t), α is defined in equation [START_REF] Bénisti | Finite range of large perturbations in Hamiltonian dynamics[END_REF], and where we took again advantage of cancellations resulting from the fact that beam particles sit on arrays and from condition [START_REF] Bekefi | Collective emission processes in plasmas[END_REF]. For t τ A , the integral term and its complex conjugate combine to produce a term scaling like |A(t)| 2 t 0 τ sin(Ω j τ )dτ . For |Ω j |t < π, it has the sign of Ω j , which implies an acceleration with the opposite sign. This corresponds to an average synchronization of the particles with the wave. Such a synchronization was observed experimentally for the particles of a monokinetic beam in a travelling wave tube [START_REF] Doveil | Experimental observation of nonlinear synchronization due to a single wave[END_REF]. Equation [START_REF] Cary | Enhancement of the velocity diffusion in longitudinal plasma turbulence[END_REF] may also be written

0 =-αA(t)k m N j=1 k m • δr j exp(iΩ j t) + c.c. + N j=1 2|α| 2 k 2 m A(t)k m ∂ ∂Ω j t 0 A * (t -τ ) exp(iΩ j τ ) dτ + c.c. . (28)
Using again Taylor expansion in τ , and following the procedure leading to equation [START_REF] Berndtson | Analysis of velocity diffusion of electrons with Vlasov-Poisson simulations[END_REF], equation [START_REF] Case | Plasma oscillations[END_REF] becomes

0 =-αA(t)k m N j=1 k m • δr j exp(iΩ j t) + c.c. + N j=1 4|α| 2 k 2 m k m |A(t)| 2 ∂ ∂Ω j sin(Ω j t) Ω j + A(t) Ȧ * (t) ∂ 2 ∂Ω 2 j exp(iΩ j t) -1 2Ω j + c.c. . (29)
As in the previous subsection, we use the continuous limit of multi-beam-multi-arrays introduced in section II C to introduce a smooth distribution function, and we integrate again by parts, which yields

0 =-N αA(t)k m h(v) exp[iΩ j t] dv + c.c. -4N |α| 2 k 2 m k m |A(t)| 2 k m g (v) sin(Ωt) Ω dv - A(t) Ȧ * (t) k 2 m g (v) cos(Ωt) -1 + i sin(Ωt) 2Ω dv + c.c. . (30)
For t = 0, the three contributions vanish identically. Moreover, since cos(Ωt)-1 Ω vanishes when Ω goes to 0, this factor cannot provide a contribution of nearly-resonant particles (those with Ω 0).

On the contrary, since sin(Ωt) Ω does not vanish when Ω goes to 0, this factor provides a contribution of nearly-resonant particles. However, this factor eventually provides a vanishing contribution of the third integral, because it involves the factor A(t) Ȧ * (t) -A * (t) Ȧ(t), which vanishes since γ was found to be real in the previous subsection. Therefore, when t grows, equation [START_REF] Chibbaro | Reductionism, emergence and levels of reality[END_REF] becomes

asymptotically 0 = -4πN |α| 2 k m g ( ω k m )|A(t)| 2 + 2N |α| 2 k m k 2 m ω 2 p ∂ r ∂ω d|A(t)| 2 dt , (31)
where the first term corresponds to the contribution of particles nearly-resonant with the wave, while the second one is the contribution of non-resonant particles. Multiplying by the electron mass m yields

dP res dt + dP wave dt = 0, (32)
where

dP res dt = - 4πε 0 ω 2 p L 3 g ( ω k m )|A(t)| 2 k m ( 33
)
may be interpreted as the derivative of the momentum of particles nearly-resonant with the wave, and

P wave = 2ε 0 k 2 m L 3 ∂ r (m, ω) ∂ω |A(t)| 2 k m ( 34
)
is the total wave momentum in the volume L 3 (we recover the conservation of the total momentum indicated at the beginning of this subsection). Equation [START_REF] Dobrushin | Vlasov equations[END_REF] can also be written as

d|A(t)| 2 dt = 2γ L (k m )|A(t)| 2 , (35)
giving again Landau growth or damping to lowest order.

We now focus on the term of equation ( 29) involving sin(Ωt) Ω , calling it rsec1 . As indicated after equation [START_REF] Dawson | Plasma oscillations of a large number of electron beams[END_REF], it corresponds to the contribution of particles close to being resonant with the wave.

Since the derivative in Ω of sin(Ωt) The previous discussion is illustrated in figure 3, where the behaviour of B as a function of t/Ω is plotted for a given t. The left panel shows that B indeed increases with |Ω| and therefore also does the synchronization. For |Ω| large enough, the behaviour B ∼ t/Ω is also illustrated in the magnified view, which highlights the oscillations with decreasing amplitude between positive and negative values of B.

Ω is B = [Ωt cos(Ωt) -sin(Ωt)]Ω -2 ,
The average synchronization for small times can be intuitively understood by computing the velocity modulation of the particles having the same parallel velocity to the wave at t = 0. Indeed, those whose absolute differential velocity to the wave is increasing at t = 0 (blue line in figure 4) have an average velocity which is higher than that of particles with an initially decreasing absolute differential velocity (red line in figure 4). Since the former have a smaller velocity modulation than the latter, the desynchronization of the former is smaller than the synchronization of the latter, which yields an average synchronization.

As just shown, synchronization is mostly experienced by the particles whose absolute differential velocity to the wave scales as |γ L |/k m . When |γ L |/k m is small enough, these particles are quasiresonant. Consequently, there is a net loss of particles momentum when g ( ωmr km ) > 0, and a net gain of momentum in the opposite case, in agreement with equation [START_REF] Doveil | Two regimes of self-consistent heating of charged particles[END_REF]. Equation [START_REF] Dobrushin | Vlasov equations[END_REF] shows that this momentum is exchanged with the wave momentum.

If the initial distribution of particles is meant to approximate a smooth distribution g(v), this means that the length L of the system is taken large enough to fulfill the condition "there are many particles with parallel velocities in the range (ω

-|γ L |)/k < v < (ω + |γ L |)/k". However, whatever
L, for ω k large enough, this condition is unfulfilled. Therefore, one should not use r (m, ω), but its discrete counterpart computed with the N -body distribution function: Landau damping will not occur, but only the beam modes to be described in section V B.

Several textbooks, when trying to provide an intuitive explanation of Landau damping computed by a Vlasovian approach, make mechanical calculations similar to the above ones, in particular

Nicholson's one in section 6.7 of [START_REF] Nicholson | Introduction to plasma theory[END_REF]. Here, the calculation is performed with the N -body mechanical model used to compute Landau damping itself. This calculation involves a smooth velocity distribution g(v) at the end of the calculation only, while the Vlasovian approach deals with the dynamics of g(v) from the outset.

The fact that both Landau damping and instability result from the same synchronization mechanism of particles with waves, is the fundamental reason why there is a single formula for the rates of Landau growth and damping. While this comes naturally in Kaufman's derivation and in section III A, growth and damping look very different in the van Kampen-Case [START_REF] Case | Plasma oscillations[END_REF][START_REF] Van Kampen | On the theory of stationary waves in plasmas[END_REF][START_REF] Van Kampen | The dispersion equation for plasma waves[END_REF] approach and in the germane one of section V, since growth involves an eigenmode and damping only a phasemixing instead; very different also in Landau's derivation where damping requires an analytic continuation, but not growth.

IV. FUNDAMENTAL EQUATION FOR THE ELECTROSTATIC POTENTIAL

This section derives a fundamental equation for the electrostatic potential [START_REF] Escande | Corrigendum: Direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction (2015 plasma phys. control. fusion 57 025017)[END_REF][START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF][START_REF] Escande | N -body description of Debye shielding and Landau damping[END_REF], which is reminiscent of the one obtained by the calculation à la Landau starting with the Vlasov equation.

This equation is introduced in a specific section, because it is the basis for the description of Langmuir waves of the next two sections. The derivation parallels completely the one in textbooks using the Vlasov-Poisson system: the motion is linearized with respect to the ballistic motion of particles of the system with a uniform density, as well as the electrostatic potential, and the Laplace transform in time is applied to the linearized equations.

We consider a given multi-beam-multi-array defined in section II C. Let r j0 be the initial position of the unperturbed beam particle with index j, and v j be its velocity, and let ∆r j (t) = r j (t)r j0 -v j t be the mismatch of the actual position of particle j with respect to the unperturbed beam particle with the same index. We define the ballistic approximation φ(bal) (m, t) to φ(m, t) which is computed from equation ( 8) on setting r j (t) = r j0 + ∆r j (0) + [v j + ∆ṙ j (0)]t for all j's in the latter

φ(bal) (m, t) = - e ε 0 k 2 m N j=1 exp -ik m • [r j0 + ∆r j (0) + (v j + ∆ṙ j (0) ) t] . (36)
If all ∆r j (0)'s and ∆ṙ j (0)'s are small, so is φ(bal) (m, t) since φ(m, t) vanishes identically for the particles of a multi-beam-multi-array. We define the two mismatches to ballistic values

δr j (t) = ∆r j (t) -∆r j (0) -∆ṙ j (0)t, (37)
δ φ(m, t) = φ(m, t) -φ(bal) (m, t). (38)
We now compute a linearized solution of the full N -body dynamics about the chosen multibeam-multi-array 15 , and the resulting Fourier components of the potential. We assume ∆r j (t)

b smooth for all j's. To this end, we replace δ φ with its expansion to first order in the ∆r j (t)'s

δ φ(m, t) = i N j=1 e ε 0 k 2 m exp[-ik m • (r j0 + v j t)] k m • δr j (t), (39)
whose Laplace transform in time 16 is 17

δϕ(m, ω) = i N j=1 e ε 0 k 2 m exp[-ik m • r j0 ] k m • δr j (ω -k m • v j )], (40)
where the Doppler shift -k m •v j comes from the linear dependence on t of the exponent of equation [START_REF] Elskens | Finite-N dynamics admit no travelling-waves solutions for the hamiltonian XY model and single-wave collisionless plasma model[END_REF].

To compute δr j (ω -k m • v j ), we use Newton's equation ( 9) for the particles. Since m runs over the finite domain such that k m b smooth ≤ 1, the self-field due to ϕ vanishes, and it is not necessary to exclude self-interactions. Therefore, one may use the harmonics of the electrostatic potential due to all particles defined by equation ( 8), which yields

ϕ(r j , t) = 1 L 3 m, kmb smooth ≤1 φ(m, t) exp(ik m • r j ), (41)
where φ(m, t) = φ(bal) (m, t) + δ φ(m, t). Using equation ( 39), the linearized particles dynamics defined by equation ( 9) is then given by

δr j = ie L 3 m n, knb smooth ≤1 k n φ(n, t) exp[ik n • (r j0 + v j t)]. (42)
The Laplace transform in time of this equation is

-ω 2 δr j (ω) = ie L 3 m n, knb smooth ≤1 k n exp(ik n • r j0 ) ϕ(n, ω + k n • v j ), ( 43
)
15 This derivation is close to that in [START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF], but takes advantage of the simplification of [START_REF] Escande | Corrigendum: Direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction (2015 plasma phys. control. fusion 57 025017)[END_REF]. The derivation of [START_REF] Escande | N -body description of Debye shielding and Landau damping[END_REF] was convenient only when taking right away the singular limit of section VI. 16 The Laplace transform in time maps a function g(t) to g(ω) = ∞ 0 g(t) exp(iωt)dt (with ω complex). 17 Since the arguments of functions are spelled explicitly, from now on we omit diacritics for the Laplace (or Fourier) transformed quantities.

where the Doppler shift k n • v j comes from the linear dependence on t in the exponent of equation [START_REF] Elskens | Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes[END_REF], and where we take into account that δr j (0) = δ ṙj (0) = 0. Computing δr j (ω -k • v) in equation ( 40) from the expression of δr(ω) given by equation ( 43) yields

k 2 m ϕ(m, ω) - ω 2 p N n, knb smooth ≤1 k m • k n × N j=1 ϕ(n, ω + k n • v j -k m • v j ) (ω -k m • v j ) 2 exp[i(k n -k m ) • r j0 ] = k 2 m ϕ (bal) (m, ω), (44)
with

ϕ (bal) (m, ω) = N j=1 ϕ (bal) j (m, ω), (45)
where

ϕ (bal) j (m, ω) = - ie ε 0 k 2 m exp[-ik m • r j (0)] ω -k m • ṙj (0) , (46)
is the ballistic potential related to particle j. Note that assumption k m b smooth ≤ 1 (equation ( 6)), with b smooth λ D (equation ( 7)) excludes scales which are irrelevant to Debye shielding and Landau damping, since these phenomena involve scales larger than, or of the order of λ D .

Because of equation [START_REF] Bekefi | Collective emission processes in plasmas[END_REF], if n = m, for each beam the corresponding values of exp

[i(k n -k m )•r j0 ]
are uniformly distributed on the unit circle, and their global contribution to the non-diagonal term in equation ( 44) vanishes. Therefore, equation [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF] becomes

d (m, ω) ϕ(m, ω) = ϕ (bal) (m, ω), (47)
where

d (m, ω) = 1 - ω 2 p N N j=1 1 (ω -k m • v j ) 2 . ( 48
)
This is the fundamental equation for the potential, which is going to be used for the description The truncated Coulomb potential cannot correctly describe the encounters between particles with impact parameters smaller than b smooth , which makes our description of the dynamics relevant for times shorter than the collision time τ coll = 3(2π) 3/2 Λ/(ω p ln Λ), as happens for the Vlasovian description. Both descriptions are relevant for Langmuir waves, since conditions Λ 1 and

ω p τ coll 1 are equivalent.
As a result of equations [START_REF] Elskens | Diffusion limit for many particles in a periodic stochastic acceleration field[END_REF][START_REF] Elskens | Vlasov equation and N -body dynamics : How central is particle dynamics to our understanding of plasmas ?[END_REF][START_REF] Escande | Description of Landau damping and weak Langmuir turbulence through microscopic dynamics[END_REF][START_REF] Escande | Large scale structures in kinetic plasma turbulence[END_REF], the part of ϕ(m, ω) generated by particle j is

δϕ j (m, ω) = ϕ (bal) j (m, ω)/ d (m, ω). (49)
By inverse Fourier-Laplace transform, after some transient discussed later, the potential due to particle j is the sum of two parts : one due to the zeros of d (m, ω) and one to the pole ω = k m • ṙj (0). These two contributions are computed in the next two sections: "Landau damping à la van Kampen-Dawson" and "Langmuir waves and Debye shielding à la Landau".

V. LANDAU DAMPING À LA VAN KAMPEN-DAWSON

This section deals anew with Langmuir waves by keeping the discrete summation in the expression of d (m, ω). Such a dielectric function was considered by Dawson in 1960 for a one-dimensional plasma made up of many fluid monokinetic beams [START_REF] Dawson | Plasma oscillations of a large number of electron beams[END_REF] whose velocities are successive multiples of a small velocity δ. He showed that d (m, ω) brings two beam modes per beam. Their eigenfrequencies are pairs of complex conjugate values for ω, whose imaginary parts tend to vanish when δ decreases : this makes these modes analogous to the van Kampen modes.

Here, Landau damping is recovered by phase mixing of these modes by following a procedure similar to that used in section 3.8 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]. The spacing of the beam velocities is kept finite, but the limit of a vanishing spacing is used at the end of the calculation to approximate finite sums by integrals. Most of this section deals with Φ j m (r, t), the part of the potential due to particle j, provided by the zeros of d (m, ω) in equation (49). A final subsection shows that the part of the potential coming from the pole ω = k m • ṙj (0) vanishes.

A. Granular dielectric function

The calculation of the zeros of d (m, ω) is simpler when the cubic array has an edge parallel to k m and one of its points is at the origin of velocities. In this section, we make this choice, and show that, in this limit, the contribution of these zeros to Φ j m (r, t) corresponds to two Langmuir waves, and that the chosen discretization of the beams is absent in the final expression of Φ j m (r, t). One of these waves propagates in the direction of k m and the other one in the direction of k -m = -k m .

In the following, we consider only the first type of wave, and we associate the other one to k -m .

Let δ be the mesh size of the grid, 1 be the unit vector parallel to k m , and 2 and 3, be the other two unit vectors. We index a beam by three integers (σ, τ, υ) corresponding to its position on the grid, so that its velocity is w σ,τ,υ = (σ 1 + τ 2 + υ 3)δ. Let N σ,τ,υ be its number of particles. Then equation [START_REF] Escande | Large scale structures in kinetic plasma turbulence[END_REF] becomes

d (m, ω) = 1 - ω 2 p N σ,τ,υ N σ,τ,υ (ω -σk m δ) 2 = 1 - ω 2 p N σ M σ (ω -σk m δ) 2 , ( 50
)
where

M σ = τ,υ N σ,τ,υ .
We consider again the beam velocity distribution to be a granular approximation to a spatially uniform smooth velocity distribution whose integral perpendicular to the 1 axis is g(v). Then M σ is taken equal to the integer part of N g(σδ)δ. When the number of particles in the Debye sphere increases, g(v) is split over an increasing number of beams (∼ δ -1 ) whose velocities lie on grids with mesh size δ going to 0. The dielectric function of the granular plasma is now expressed with the one-dimensional distribution function g(v) and reads

d1 (m, ω) = 1 -ω 2 p σ g(σδ)δ (ω -σk m δ) 2 . ( 51
)
Subsection V C shows that at most two zeros of d1 (m, ω) have a finite imaginary part when δ goes to zero. The other zeros have an imaginary part vanishing in the limit where δ goes to zero. Dawson introduced in 1960 a very clever technique to compute them, when he considered a one-dimensional plasma made up of many fluid beams [START_REF] Dawson | Plasma oscillations of a large number of electron beams[END_REF]. He decomposed the sum in equation ( 51) into a regular part converging to

1 (m, ω) = 1 -ω 2 p g(v) (ω -k m v) 2 dv = 1 + ω 2 p k m g (v) (ω -k m v) dv, (52)
and a singular part, which is summed exactly using classical summation formulas for trigonometric functions (see Appendix C). The conjugate zeros of d1 (m, ω) have real parts between the frequencies of the beams σk m δ, and their imaginary part scales like δ | ln(δ/v T ) | (see equations ( 123)-( 124)

of Appendix C). We write these zeros as ν σ,µ = α σ + µiβ σ , where µ = ±1.

B. Phase mixing

We now follow a procedure similar to that used in section 3.8.3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF] (see also [START_REF] Elskens | Irreversible behaviours in Vlasov equations and many-body hamiltonian dynamics: Landau damping, chaos and granularity[END_REF]). The zero of d1 (m, ω) with index (σ, µ) brings to Φ j m (r, t) a contribution

Φ j m σ µ (r, t) = - e ε 0 k 2 m L 3 σ µ exp[i(k m • (r -r j (0)) -ν σ,µ t)] ν σ,µ -k m • ṙj (0) + c.c. , (53)
where

σ,µ = ∂ d1 ∂ω (m, ν σ,µ ) = 2ω 2 p σ g(σk m δ)δ (ν σ,µ -σk m δ) 3 . ( 54
)
The sum over j of the Φ j m σ µ (r, t)'s yields

Φ m,σ,µ (r, t) = - exp ik m • r ε 0 k 2 m L 3 eN f (m, v) ν σ,µ -k m • v d 3 v exp(-iν σ,µ t) σ,µ + c.c. , (55)
where f (r, v) = N -1 j δ(r -r j (0)) δ(v -ṙj (0)), and f (m, v) is its spatial Fourier transform (δ(•) stands for the Dirac distribution).

Using again a decomposition into a regular part and a singular one, for small δ, the expression

of σ,µ is shown in Appendix C to converge toward σ,µ -iµ 2π k m δ 1 + P ∞ -∞ ω 2 p g (v) k m (α σ -k m v) dv -i πµω 2 p k 2 m g ( α σ k m ) . (56)
The sum of the first two terms in the bracket of this equation is the real part of 1 (m, ν σ,µ ) defined by equation [START_REF] Bos | Reduction of mean-square advection in turbulent passive scalar mixing[END_REF]. Since 1 (m, ω m ) = 0, this implies that the ν σ,µ 's bring a contribution to Φ m,σ,µ (r, t), which is resonant in the vicinity of ω m . Then, in the summation over (σ, µ), we may take

σ,µ - 2π k m δ m (γ L (k m ) + iµ(α σ -ω mr )), (57)
where m = ∂ 1 ∂ω (m, ω m ), and where the definition [START_REF] Boyd | The physics of plasmas[END_REF] of γ L (k m ) was used 18 . This yields

Φ m (r, t) = exp ik m • r ε 0 k 2 m L 3 m eN f (m, v) ω mr -k m • v d 3 v S + c.c. , (58)
where

S = k m δ 2π Σ σ,µ exp(-iν σ,µ t) γ L (k m ) + iµ(α σ -ω mr ) = k m δ π Σ σ exp(-iα σ t) γ L (k m ) cosh(β σ t) γ 2 L (k m ) + (α σ -ω mr ) 2 sgn(γ L (k m )) exp[-(|γ L (k m )| + iω mr )t], (59)
where the last expression makes cosh(β σ t) = 1, since β σ goes to zero in the continuous limit, and uses a Fourier transform identity in this limit. Equations ( 58) and ( 59) yields a contribution to the wave

Φ m (r, t) = sgn(γ L (k m )) exp[i(k m • r -ω m t) -|γ L (k m )|t] ε 0 k 2 m L 3 m × eN f (m, v) ω m -k m • v d 3 v + c.c. . (60)
We notice that the initially chosen discretization of the beams is absent in this expression.

C. Vlasovian zero of the granular dielectric function

When δ goes to 0, d1 (m, ω) as defined by equation (51) converges toward its continuous limit 1 (m, ω) defined by equation ( 52), provided that Im ω = 0 or g (v) vanishes in a finite domain about ω/k m . In ( 52) one recognizes the Vlasovian expression of the dielectric function.

We first look for the zeros ω mr + i ω mi of d1 (m, ω) having a non vanishing ω mi when δ goes to zero. Since |ω mi | does not vanish, the convergence of d1 (m, ω) toward 1 (m, ω) implies that the considered zero of the former converges to a zero of the latter. If the zero of the former has a positive ω mi when δ goes to zero, its limit value corresponds to a Vlasovian zero with the same property, which requires g (ω mr /k m ) > 0. Indeed, the corresponding Vlasovian growth rate is Landau's one and is given by equation [START_REF] Boyd | The physics of plasmas[END_REF]. Since d1 (m, ω) has zeros coming in conjugate pairs, there is also a zero with the damping rate -γ L (k m ). This zero will be important to recover the correct Landau instability when g (ω mr /k m ) > 0. This is surprising for people used to the Vlasovian setting, where only an unstable root is present.

If g (ω mr /k m ) < 0, d1 (m, ω) cannot have the corresponding Vlasovian root when δ goes to zero, since its existence would imply that of the complex conjugate zero, which itself would reduce this case to the previous one, implying a contradiction. In the usual Vlasovian approach [START_REF] Lifshitz | Landau and Lifshitz Course of Theoretical Physics 10 : Physical Kinetics[END_REF], finding this zero results from continuing analytically 19 g (v) outside the real axis, which is nonsense for a derivation involving the sum of Dirac distributions corresponding to the many-beam distribution.

Therefore, there is no zero of d1 (m, ω) having a non-vanishing imaginary part when δ goes to zero. We now show that the missing Vlasovian Landau damping is compensated by the damping resulting from the phase mixing of the many beam modes defined by the zeros of d1 (m, ω). 19 Actually, the analysis of Cauchy integrals, a more advanced topic in mathematics (section VIII.12 of [START_REF] Godement | Mathematical Analysis III[END_REF]), allows wilder g(v)'s. There, the analyticity of the integrals with respect to ω, not v, in the upper or lower complex halfplanes mirrors the use of Laplace transform. The relevant g(v)'s are such that their absolute values are integrable, as well as that of their Fourier transforms.

D. Landau damping

As just explained, if γ L (k m ) < 0, all zeros of d1 (m, ω) have an imaginary part vanishing in the limit where δ goes to zero. Therefore, the only contribution to the potential of the wave is provided by equation [START_REF] Escande | Application of Picard iteration technique to selfconsistent wave-particle interaction in plasmas[END_REF], which is the Vlasovian result for the Langmuir wave with wavevector k m propagating in the direction of the wavevector, except for the granular character of the initial distribution in space and velocity. Φ m (r, t) is exponentially damped, because of the phase mixing of the beam modes, which makes these modes analogous to the van Kampen modes [START_REF] Case | Plasma oscillations[END_REF][START_REF] Van Kampen | On the theory of stationary waves in plasmas[END_REF][START_REF] Van Kampen | The dispersion equation for plasma waves[END_REF] (see also section 6.14 of [START_REF] Nicholson | Introduction to plasma theory[END_REF]). However, while van Kampen modes and those in section 3.8.3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF] are eigenmodes, the Laplace transformed quantities of the present derivation are not. Finally, we stress that the phase mixing of this derivation is different from that in section III. Indeed the latter corresponds to the transients related to the terms complementing those of Landau poles when inverting the Laplace transform in the present derivation.

E. Landau instability

For γ L (k m ) > 0, one must add to the expression of Φ m (r, t) provided by equation ( 60) the contributions of the two roots with the finite imaginary parts 20 ±γ L (k m ). This yields the sum of three exponentials e -iωmrt e γ L (km)t + e -γ L (km)t -e -γ L (km)t = e -iωmrt+γ L (km)t ,

where the first two come from the two roots with finite imaginary parts, and the last one comes from equation [START_REF] Escande | Application of Picard iteration technique to selfconsistent wave-particle interaction in plasmas[END_REF]. Then, the above phase mixing term cancels the contribution of the damped root and one finally obtains again the Vlasovian expression (60) 21 . It is therefore of paramount importance that the single root of 1 (m, ω) yields two conjugate roots of d1 (m, ω). In section 3.8.3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF], a similar calculation was done by using normal modes, and not the Laplace transform. An equation similar to equation ( 61) was obtained where the third exponential was e -γ L (km)|t| instead;

this brings the cancellation of the second exponential for positive times, and of the first one for negative ones, because of the time-reversible character of the dynamics.

To illustrate the previous discussion, we have numerically calculated the zeros of the granular dielectric function [START_REF] Escande | Complexity and simplicity of plasmas[END_REF]. This was done for the following realistic plasma parameters: n e = 10 19 m -3 , 20 In the limit where δ vanishes, the absence of these two roots does not modify equation [START_REF] Escande | Application of Picard iteration technique to selfconsistent wave-particle interaction in plasmas[END_REF], since they bring an infinitesimal contribution to the result. 21 Without the phase mixing term, a wave with an initial amplitude 1 would have only half the amplitude e γ L (km)t after a long enough time. This remark by Dr A. Samain was at the origin of the development of the van Kampen-like calculation in section 3.8.3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF].

T = 1 keV and L = 30λ D . We performed the calculation for the mode number m = 1. To model the distribution of particles in velocity space, we considered a thermal population modelled by a Maxwellian distribution and an energetic tail modelled by a shifted Maxwellian

M σ ∼ f th e -v 2 σ 2 + f t e -(vσ -v 0 ) 2 2 (62)
where, for educational reasons in order to give a clear example of the Landau instability, we have chosen the values f t = 0.15, f th = 0.85 and v 0 = 4. Velocities are normalised to the thermal velocity v th and frequencies in equation ( 50) are normalised to (2π/L)v th . The velocity space is discretised to form the grid v min ≤ v σ = δσ ≤ v max , where v min = -5 and v max = 7, both normalised to the thermal velocity, and σ is used to parametrise the distribution of particles as in the previous discussion. The zeros of the granular dielectric function are plotted in figures 5b and 5d for two values of δ, namely δ = 0.5 and δ = 0.1, respectively 22 . Figures 5a and5c represent the distribution of particles in velocity space, where the inversion of the slope is clearly visible around v ≈ 3 < v 0 = 4. It can also be observed that the imaginary part of the zeros decreases when δ decreases, and only two zeros emerge from the whole set of zeros, leading to the Landau instability.

The wave echo experiment [START_REF] Baker | Ion-wave echoes[END_REF] proved beam modes to be the actual support of Landau damped Langmuir waves excited by grids in a magnetized plasma column. In such a plasma, they would exist in a Landau unstable case too, but the above result shows their phase mixing contribution would be cancelled by a damped eigenmode.

For a wave with phase velocity v w , as in section III B, the above passage to the smooth velocity distribution g(v) implicitly assumes that there are many particles with parallel velocities in the

typical range [v w -|γ L |/k, v w + |γ L |/k],
where γ L is the Landau growth or damping rate of the wave, and k the modulus of its wavevector. Indeed, this range is the one where the phase mixing à la van Kampen is occurring, as shown in section V B, and where the synchronization of particles with the wave brings the change of particle momentum inducing Landau damping or growth of this wave, as shown in section III B. If the number of particles in the range

[v w -|γ L |/k, v w + |γ L |/k] is
not large enough, the system behaves as a multi-beam in this range.

F. Collisions at work

We now comment on the nature of the unstable modes of the above many-beam system. When a cold beam with a high enough density is added to a Maxwellian plasma, it triggers the cold beam- 22 These two figures represent actually the poles of 1/| d (m, ω)|, which are nothing else but the zeros of d (m, ω).

Identifying the poles of a function of complex variable is easily done by identifying the closed contours in the complex plane. plasma instability whose growth rate is defined by the beam and plasma densities (see for instance

[107]). When the beam becomes warmer and warmer, the beam-plasma instability becomes kinetic and its growth rate is given by Landau's formula [107]. At some moment, the velocity spreading of the beam becomes so large that the slope of the global velocity distribution of the beam-plasma system becomes negative everywhere: the instability is quenched. Therefore, in the case of a kinetically stable distribution, the unstable modes of the above many-beam system cannot be interpreted as beam-plasma ones. Indeed, their growth rate (124) is mainly defined by global features of the plasma like the slope of g(v), and depends weakly on the beam density through a logarithm of g(nδ)/δ.

Actually, the instability is due to what is traditionally called collisions, and will be described in section X. Indeed, section II C showed the distribution of the pure many-beam system is invariant in time. Therefore, it is insensitive to collisions. However, a typical perturbation to this distribution evolves because of collisions, and exponentially diverges from the many-beam distribution with the many-beam growth rates. What is traditionally called collisions involves scales up to the Debye length. This explains why the present derivation keeps a "collisional" contribution, while The whole derivation of section V B heavily relies on the fact that the β σ 's are larger than k m δ, i.e. that collisions are strong enough. Therefore, not only collisions make Landau damping irreversible [START_REF] Callen | Coulomb collision effects on linear Landau damping[END_REF], but they are necessary for the effect to exist. Furthermore, in the second equality of equation ( 59), for a finite δ, the exp(-iα σ t)'s might produce a recurrence of S after a time (k m δ) -1 .

However, the cosh(β σ t)'s grow on the smaller time-scale [k m δ | ln(δ/v T ) |] -1 , which prevents the recurrence 24 . The last equality of equation ( 59) makes sense over a time |γ

L (k m )| -1 only if condition k m δ | ln(δ/v T )γ L (k m ) | -1 1 ( 63
)
is fulfilled.

G. Vanishing contribution of the ballistic poles

We now consider the part of the potential due to particle j coming from the pole ω = k m • ṙj (0).

It is

Φ j d (r, t) = Φ d (r -r j (0) -ṙj (0)t, ṙj (0)), (64)
where 23 It was recently noted that the collisional damping rate of Langmuir waves is much smaller than the one provided by the collision frequency [121,[START_REF] Yoon | Weak turbulence theory for collisional plasmas[END_REF]. 24 This property is the signature of the Floquet exponents mentioned in Appendix A.

Φ d (r, u) = - e L 3 ε 0 m =0 exp(ik m • r) k 2 m d (m, k m • u + iε) (65)
with the usual iε prescription resulting from inverting the Laplace transform, as d (m, ω) diverges for some real ω's. We again compute the limit of the contribution of component k m in this expression by using the convenient grid of the previous subsection. By taking the limit ε → 0 for a fixed δ, the first term in the bracket of equation ( 118) diverges when δ goes to zero. This implies that d (m, k m • v + iε) becomes infinite, which cancels the contribution of the pole ω = k m • ṙj (0) in equation (49). Therefore, apart from transients, δϕ j (r, t) is a superposition of waves only.

VI. LANGMUIR WAVES AND DEBYE SHIELDING À LA LANDAU

The Vlasovian limit recovers the Landau damping/growth predicted by the N -body dynamics, by substituting d (m, ω) in equation ( 47) with (m, ω) defined by

(m, ω) = 1 -ω 2 p f 0 (v) (ω -k m • v) 2 d 3 v = 1 + ω 2 p k 2 m k m • ∇ v f 0 (v) (ω -k m • v) d 3 v, (66)
where f 0 (v) is the smooth approximation of the velocity distributions of the multi-beam-multiarrays defined in section II C. This section considers the solutions of equation ( 47) when this substitution is made, which corresponds to a singular limit of the problem. It becomes

(m, ω) ϕ(m, ω) = ϕ (bal) (m, ω). (67)
This singular limit was taken in [START_REF] Escande | Corrigendum: Direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction (2015 plasma phys. control. fusion 57 025017)[END_REF][START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF][START_REF] Escande | N -body description of Debye shielding and Landau damping[END_REF]. Then, Φ j m (r, t), the part of the potential due to particle j obtained by inverse Fourier-Laplace transform, was shown to comprise of two parts: one provided by the zeros of (m, ω) corresponding to the usual Vlasovian expression of Langmuir waves, and one provided by the pole ω = k m • ṙj (0) corresponding to the shielded Coulomb potential of the particle [START_REF] Balescu | Statistical mechanics of charged particles[END_REF][START_REF] Gasiorowicz | Dynamics of ionized media[END_REF]114] (see for instance section 9.2 of [105]) 25 . For particle j, this potential is

Φ j m (r, t) = Φ(r -r j (0) -ṙj (0)t, ṙj (0)), (68)
where

Φ(r, v) = - e L 3 ε 0 m =0 exp(ik m • r) k 2 m (m, k m • v + iε) (69)
with the usual iε prescription resulting from inverting the Laplace transform, as the integral in equation ( 52) is undefined for the real-valued ω = k m • v. The Vlasovian derivation of shielding requires a test particle to be added to the Vlasovian plasma (see for instance section 9.2 of [START_REF] Nicholson | Introduction to plasma theory[END_REF]). Similarly, the N -body approach requires d (m, ω)

to be substituted with its Vlasovian analog. In the standard Vlasovian calculation, ϕ (bal) (m, ω) is replaced by its continuous limit for Im ω finite

Φ (bal) (m, ω) = - ie ε 0 k 2 m f (m, v) ω -k m • v d 3 v, (70)
which is the smoothed version of the actual shielded potential in the plasma. Unfortunately, this quantity has no obvious interpretation. Therefore, generally textbooks do not dwell upon it, and just use it as the term of initial condition for the calculation of Langmuir waves. The N -body description reveals that it is the continuous limit of a granular source term bringing not only the excitation of Langmuir waves, but also the Debye-shielded potential of the particles.

The Debye shielded potential is one of the simplest examples of a renormalized potential [START_REF] Mccomb | Renormalization methods: a guide for beginners[END_REF].

This potential is a mean-field potential produced by the Coulomb deflections of the particles, as will be explained in section VIII. Vlasov equation deals also with a mean-field potential, as obvious from its mean-field derivation recalled in the introduction. This is also the case for the BBGKY derivation, because of its statistical aspect. Therefore, (m, ω) may be interpreted as a renormalized version of d (m, ω). When using (m, ω) = 0 as the dispersion relation, one Langmuir wave is actually the renormalized version of a set of beam modes of the N -body system. As will be shown in section VIII, similarly the Debye shielded potential of a particle j is its Coulomb potential dressed with all the modifications of the potential of the other particles due to their deflections by j. Analytically, it is the renormalized potential provided by equations ( 68) and ( 69) by using the renormalized dielectric function (m, ω). Therefore, it is natural to recover both the Vlasovian Langmuir waves and the Debye shielded potentials of the particles, when using (m, ω) = 0 as dispersion relation.

Equation ( 67) can be introduced intuitively by the following simple argument. An almost uniform distribution of particles, written as a perturbation of a multi-beam-multi-array distribution related to a given grid, may be also written as a perturbation for many other nearby multi-beammulti-array distributions; actually a continuum set. While d (m, ω) depends on the choice of the grid, the time evolution of the un-linearized dynamics is independent of it, and all the linearized dynamics corresponding to these various grids stay close to it for some time. This almost invariance of the linearized N -body dynamics under the above continuous set of distributions suggests approximating this dynamics by a coarse-grained version of it, where d (m, ω) is substituted with (m, ω), the Vlasovian expression of the dielectric function, which is present in equation [START_REF] Firpo | Long-time discrete particle effects versus kinetic theory in the self-consistent single-wave model[END_REF].

Appendix D describes a coarse-graining leading to this equation. As shown in Appendix E, the shielded Coulomb potential can also be recovered by a singular limit of the many-beam description.

VII. THE VLASOVIAN LIMIT: A SINGULAR ONE

The singularity of the Vlasovian limit is mathematically visible in the dielectric function. Indeed, finding the zeros of d (m, ω) corresponds to finding those of a polynomial of degree 2n b with real coefficients (n b is the number of beams). While when increasing n b , the number of zeros of d (m, ω) keeps increasing, this number is fixed for (m, ω). Even more singular is the transition from a finite increasing number of poles for d (m, ω) to a cut for (m, ω) 26 .

The results of sections III and V show that the Vlasovian limit is also singular from a physical point of view. Indeed, the mechanical description of the plasma as a granular system shows that, as to the electrostatic potential, Landau damping is due to a phase mixing, and that Landau growth involves a phase mixing too. As already mentioned, the wave echo experiment [START_REF] Baker | Ion-wave echoes[END_REF] proved the existence of the beam modes in a genuine (granular) plasma. In contrast, in the usual Vlasovian setting, the derivation of Landau damping requires an analytic continuation, which is, though powerful, far from physically intuitive. While section III B shows that Landau growth and damping are due to the same average synchronization of particles nearly-resonant with the wave, the fate of individual particles is usually eluded in a Vlasovian setting 27 . As recalled in the introduction of this paper, this lack of physical intuition made it hard for the plasma community to accept the reality of this damping, and is the origin of still discordant physical interpretations of the effect, some invoking trapping or surfing of the particles, which is wrong mechanically.

Adding a test particle to the N -body system, induces a mere modification of the calculations of section IV: one more term is added to the ballistic potential [START_REF] Doveil | Experimental observation of nonlinear synchronization due to a single wave[END_REF], but the remaining of section IV stays the same. Therefore, equation (47) stays the same except for the contribution of the test particle in the ballistic term. The calculation of section V G also applies to this contribution, which

does not provide the shielded potential of this test particle. This is in contrast with the Vlasovian case, another symptom of the singularity of the Vlasovian limit.

The reason, why the derivation of Landau damping in a Vlasovian setting does not correspond to an eigenmode, is deeply rooted in mechanics. Indeed, imagine Landau damping would correspond to a zero of d (m, ω): since the above polynomial has real coefficients, its complex solutions come in conjugate pairs, and the damped eigenmode would come with a growing companion, which would dominate for a typical initial condition of the system. This contradicts observation. Therefore, Landau damping cannot correspond to a zero of d (m, ω), and section V B showed that, for the genuine granular system, Landau damping occurs by phase mixing of beam modes.

As a matter of fact, this reductio ad absurdum can be done without any calculation. Indeed, since the N -body system is a Hamiltonian one, its Lyapunov exponents come in pairs: consistently with the preservation of phase space areas, a damped eigenmode comes with a growing companion.

Therefore, while Landau damping corresponds to a zero of (m, ω), it cannot for d (m, ω). This issue does not exist for the unstable case.

Furthermore, while introducing the van Kampen modes was a creative contribution in the Vlasovian case, they are spontaneously present in the N -body approach, because the multi-beammulti-array distribution is a natural way to deal with a "uniform" plasma. The van Kampen modes are beam modes recovered with (m, ω), which is a renormalized dielectric function, as was shown in section VI. As was explained in section V F, the instability of the beam modes is the signature of Coulomb collisions. If a collision operator is added in Vlasov equation, the continuous spectrum of the Van Kampen-Case modes is eliminated and replaced by a discrete spectrum, even in the limit of zero collision, and these discrete eigenmodes form a complete set of solutions [START_REF] Ng | Kinetic eigenmodes and discrete spectrum of plasma oscillations in a weakly collisional plasma[END_REF]104]. This is similar to the N -body case. However, the Landau-damped solutions are recovered as true eigenmodes, which they are neither in the collisionless Vlasovian theory, nor in the N -body approach. Finally, we notice that in the N -body approach, Landau damping occurs without requiring g(v) to be analytically continuable like in the usual Vlasovian approach, which is better for the experimental observability of the phenomenon.

VIII. MEDIATED INTERACTIONS IMPLY DEBYE SHIELDING

In this section 28 , Picard's iteration technique 29 is applied to the full equation of motion of a particle j due to the Coulomb forces of all other ones for any Λ, i.e. by using equation ( 6) without the restriction k m b smooth ≤ 1. It stresses that, for Λ 1, a part of the effect on particle j of another particle j is mediated by all other particles (equation ( 75)) and reduces the direct part.

Indeed, particle j modifies the motion of all other particles, implying that the action of the latter ones on particle j is affected by particle j . This will be shown to provide an intuitive picture of Debye shielding.

We consider the dynamics in real time of the particle with index l defined by equation ( 9) with the full OCP Coulomb potential of equation [START_REF] Baldwin | Magnetized plasma kinetic theory. II. derivation of Rosenbluth potentials for a uniform magnetized plasma[END_REF]. To this end we use Picard's iteration technique.

From equation ( 9), r

l , the n-th iterate for r l , is computed from

r(n) l = e m ∇ϕ (n-1) (r (n-1) l ), (71)
where ϕ (n-1) is obtained by the inverse Fourier transform of equation ( 8) with the r j 's substituted with the r

m 's and ξ(0) m 's vanish identically for all m's. It is convenient to write equation [START_REF] Godement | Mathematical Analysis III[END_REF] as

ξ(n) l = j∈S;j =l ξ(n) lj , (72)
with S denoting the set of integers from 1 to N labeling particles, and

ξ(n) lj = a C (r (n-1) l -r (n-1) j ) ( 73
)
and

a C (r) = ie 2 ε 0 mL 3 m =0 k -2 m k m exp(ik m • r). ( 74
)
28 For completeness, most of this section follows closely section 5 of [START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF]. 29 Picard's iteration technique is one of the standard methods to prove the existence and uniqueness of solutions to first-order equations with given initial conditions. It uses the fact that the exact solution of equation dX dt = f (X) is the fixed point of the iterative process starting from n = 0, and providing Xn+1 from Xn by equation

dX n+1 dt = f (
Xn) with any choice of X0(t). This iteration technique is very convenient, in particular to alleviate the algebra of perturbation calculations.

Let ξ

(n) lj = t 0 t 0 ξ(n) lj (t ) dt dt = t 0 (t -t ) ξ(n) lj (t ) dt . For n ≥ 2, one finds ξ(n) l = j∈S;j =l [ ξ(1) lj + M (n-1) lj + 2∇a C (r (0) l -r (0) j ) • ξ (n-1) lj ] + O(a 3 ), ( 75
)
where a is the order of magnitude of the total Coulombian acceleration, and

M (n-1) lj = ∇a C (r (0) l -r (0) j ) • [ξ (n-1) l -ξ (n-1) j -2ξ (n-1) lj ] (76)
= ∇a C (r

(0) l -r (0) j ) • i∈S;i =l,j (ξ (n-1) li + ξ (n-1) ij ), (77)
where the second expression takes into account that a C (r) is anti-symmetrical in r. The latter expression displays ξ

(n-1) ij
which is the deflection of particle i by particle j. It shows how the bare Coulomb acceleration of particle l due to particle j is modified by the following process : particle j modifies the motion of all other particles, so that the action of the latter ones on particle l is modified by particle j. Therefore M

(n-1) lj is the acceleration of particle l due to particle j mediated by all other particles. The last term in the bracket in equation ( 75) accounts for the fact that both particles j and l are shifted with respect to their ballistic positions. Both M

(n-1) lj and this last term are anti-symmetrical with respect to the labels j and l, since ∇a C (r) is an even function of r.

The full Coulombian dynamics of the plasma includes the relaxation to a thermal state. If the corresponding temperature is low enough, the particles have a vanishing mean velocity. In order to describe such an equilibrium, it is advisable to take v m = 0 for all m's for a faster convergence of the iterative process toward the equilibrium solution.

Since the shielded potential of section VI was found by first order perturbation theory, it is felt in the acceleration of particles computed to second order. This acceleration is provided by equation ( 75) for n = 2. Therefore, its term in brackets is the shielded acceleration of particle l due to particle j. As a result, though the summation runs over all particles, its effective part is only due to particles j typically inside the Debye sphere (with radius λ D ) about particle l. Starting from the third iterate of the Picard scheme, the effective part of the summation in equation [START_REF] Hauray | Particle approximation of Vlasov equations with singular forces : Propagation of chaos[END_REF] ranges inside this Debye sphere, since the ξ (n-1) lj 's are then computed with a shielded acceleration.

In equation ( 77), the compound effect of the ξ (n-1) ij 's, the deflections of particle i by particle j, is to diminish the negative charge inside a sphere centered on particle j.

This calculation yields the following interpretation of shielding. At t = 0, consider a set of uniformly distributed particles, and especially particle j. At a later time t, the latter has deflected all particles which made a closest approach to it with a typical impact parameter b v T t, where v T is the thermal velocity. This part of their global deflection due to particle j reduces the number of particles inside the sphere S j (t) of radius v T t about it. Therefore, according to Gauss' theorem, the effective charge of particle j as seen out of S j (t) is reduced : the charge of particle j is shielded due to these deflections. This shielding effect increases with t, and thus with the distance to particle j. It becomes complete at a distance on the order of λ D . Since the global deflection of particles includes the contributions of many other ones, the density of the electrons does not change, at variance with the shielding at work next to a probe (see e.g. section 2.2.1 of [START_REF] Piel | Plasma physics: an introduction to laboratory, space, and fusion plasmas[END_REF]). This interpretation explains how a given particle can be shielded by all other particles, while contributing to their individual shieldings. The Debye shielding of a test particle can be computed by using explicitly the just described Coulombian deflections [START_REF] Meyer-Vernet | Aspects of Debye shielding[END_REF].

When starting from random particle positions, the typical time-scale for shielding to set in is the time for a thermal particle to cross a Debye sphere, i.e. ω -1 p , where ω p is the plasma frequency. Furthermore, shielding, though very fast a process, is a cooperative dynamical one, not a collective (viz. coherent) one : it results from the accumulation of almost independent repulsive deflections with the same qualitative impact on the effective electric field of particle j (if point-like ions were present, the attractive deflection of charges with opposite signs would have the same effect). So, shielding and Coulomb scattering are two aspects of the same two-body repulsive process. For Λ 1, this transport is negligible on fairly long time scales. However, collisions are of paramount importance to provide shielding over the plasma period ω -1 p . We now understand that Debye shielding cannot work for a multi-beam-multi-array, because it does not experience collisions, as anticipated in section II C. Finally, in contrast to what occurs if electron j is substituted with a Langmuir probe, this does not change the density of charges of the plasma, because the deflections due to j are compensated by those of the other electrons. This approach clarifies the mechanical background of the calculation of shielding using the equilibrium pair correlation function, which shows shielding to result from the correlation of two particles occurring through the action of all the other ones (see e.g. section 12.3 of [START_REF] Boyd | The physics of plasmas[END_REF]). As discussed in section VI, the derivation of Debye shielding using Vlasov equation plus a test particle takes advantage of the renormalization involved in the Vlasovian dielectric function. However, rigorously speaking, the derivation of Debye shielding using fluid equations is not justified in a collisionless plasma, since these equations are not justified, especially in a non magnetized plasma (see chapter 3 of [START_REF] Hazeltine | The framework of plasma physics[END_REF]).

IX. WAVE-PARTICLE INTERACTION: LINEAR AND NONLINEAR EFFECTS

This section introduces a wave-particle description of Langmuir waves. This is done by splitting the set of particles into a bulk, which cannot resonate with Langmuir waves, and a tail. Then, amplitude equations are derived for the Fourier components of the potential where tail particles provide a source term. These equations, together with the equations of motion of the tail particles, provide a Hamiltonian description of the dynamics of tail particles coupled with Langmuir waves. This description is well fitted to statistical descriptions unifying Landau damping and spontaneous emission, or showing that the transition from Landau damping to O'Neil's damping with trapping is a second order phase transition. This description is also amenable for studying the nonlinear saturation of various regimes of the beam-plasma instability, and to prove the nonexistence of Bernstein-Greene-Kruskal (BGK) modes corresponding to traveling-wave solutions. Finally, this description enables the use of modern tools of nonlinear dynamics and chaos available for finite dimensional systems. This section first derives the Hamiltonian ruling the self-consistent dynamics of tail particles with Langmuir waves. On this basis, a statistical analysis shows Langmuir waves evolve under the compound action of Landau damping and of spontaneous emission. Then, these two effects are shown to have their counterpart in the diffusion and friction coefficients of the quasilinear Fokker-Planck equation ruling the dynamics of particles in a broad spectrum of waves. This sets the ground to face the controversial issue of the saturation of the weak warm beam instability.

The remaining of the section deals with the single wave Hamiltonian both in the unstable and the damped case.

A. Self-consistent Hamiltonian

Till now, we described Langmuir waves by a fully linear theory. Following the analysis of section 6 of [START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF], we now generalize the analysis of section IV to afford the description of both linear and nonlinear effects in wave-particle dynamics. Indeed, resonant particles may experience trapping or chaotic dynamics, which imply k m • ∆r j 's of the order of 2π or larger for wave k m 's.

For such particles, it is not appropriate to make the linearizations leading to equations ( 39) and [START_REF] Elskens | Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes[END_REF]. However, these linearizations may still be justified for non-resonant particles over times where trapping and chaos show up for resonant ones. In order to keep the capability to describe the latter effects, we now split the set of N particles into bulk and tail, in the spirit of references [START_REF] Antoni | Explicit reduction of N-body dynamics to self-consistent particle-wave interaction[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF]108,109]. The bulk is defined as the set of particles which are not resonant with Langmuir waves. We then perform the analysis of section IV for these N bulk particles, while keeping the exact contribution of the remaining N tail particles to the electrostatic potential. To this end, we number the tail particles from 1 to N tail , the bulk ones from N tail + 1 to N = N bulk + N tail , and we call these respective sets of integers S tail and S bulk . Then, equation ( 44) is substituted with

d (m, ω) ϕ(m, ω) = ϕ (bal) (m, ω) + U (m, ω), (78)
where d (m, ω) is defined by equation ( 48), and ϕ

(bal)
bulk (m, ω) is defined by equation [START_REF] Elskens | Diffusion limit for many particles in a periodic stochastic acceleration field[END_REF], where the sums over j runs over S bulk only, and where

U (m, t) = - e ε 0 k 2 m j∈S tail exp(-ik m • r j (t)), (79)
for N bulk 1.

Since Langmuir waves are not resonant with the bulk, the singularities of d (m, ω) do not show up, and in the limit of large numbers of particles in the Debye sphere, this quantity may be approximated by (m, ω), or more precisely by bulk (m, ω) defined by

bulk (m, ω) = 1 -ω 2 p f 0 (v) (ω -k m • v) 2 d 3 v = 1 + ω 2 p k 2 m k m • ∇ v f 0 (v) (ω -k m • v) d 3 v, (80)
where f 0 (v) is the bulk distribution only, and ω p is computed with the bulk particles only. Then,

equation (78) becomes bulk (m, ω) ϕ(m, ω) = ϕ (bal) bulk (m, ω) + U (m, ω). (81)
For the scales much larger than λ D , the electric potential for the bulk is a superposition of Langmuir waves. The presence of tail particles slightly modifies these waves. Therefore, as shown in section 6.1 of [START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF], one can derive an amplitude equation for the potential ϕ(m, t) of the wave with wavevector k m in a way similar to references [108, 109]

dϕ(m, t) dt + iω m ϕ(m, t) = ie ε 0 k 2 m ∂ bulk ∂ω (m, ω m ) j∈S tail exp(-ik m • r j ), (82)
where ω m is the eigenfrequency solving bulk (m, ω m ) = 0 corresponding to the wave propagating in the direction of k m ; this frequency is real, since it is not resonant with the support of the bulk distribution function (indeed, we assumed f 0 (v) = 0 for all v's such that k m • v = ω m ). We notice that bulk (m, ω m ) = 0 is a Bohm-Gross type dispersion relation associated with plasma oscillations of the bulk.

The self-consistent dynamics of M Langmuir waves and of the tail particles is ruled by equation ( 82) written for each wave and by the equation of motion of these particles due to the M waves,

rj = ie L 3 m n∈M k n ϕ(n, t) exp(ik n • r j ), ( 83
)
where M is the set of the indices of the M waves, and the tail-tail interactions were neglected owing to the low density of the tail particles. These two sets of equations generalize to three dimensions the self-consistent dynamics defined in references [START_REF] Antoni | Explicit reduction of N-body dynamics to self-consistent particle-wave interaction[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]100] for the one-dimensional N -body dynamics.

We now write ϕ(n, t) as a n √ 2I n exp(-iθ n ), where I n and θ n are real, and where

a n = L 3 2ε 0 k 2 n ∂ bulk ∂ω (n, ω n ) 1/2 . ( 84
)
Then equations ( 82)-( 83) can be cast as the canonical equations related to the self-consistent

Hamiltonian

H sc = N tail j=1 p 2 j 2m + n∈M ω n I n -ε N tail j=1 n∈M k -1 n κ n 2I n cos(k n • r j -θ n ) (85)
where κ n = [∂ bulk (n, ω n )/∂ω] -1/2 , and ε = ω p [2m/N ] 1/2 is the coupling parameter ruling the intensity of the wave-particle interaction. The conjugate variables for H sc are (p j , r j ) for the particles and (I n , θ n ) for the waves. H sc is the sum of the kinetic energy of particles, of the energy of waves (harmonic oscillators described in action-angle variables), and of a coupling term. It may be useful to write H sc by using the "Cartesian" coordinates of the harmonic oscillators instead of their intensity-phase (or action-angle) components. To this end we write ϕ(n, t) as a n (X n + iY n ),

where X n and Y n are real. This yields

H sc = N tail j=1 p 2 j 2m + n∈M ω n X 2 n + Y 2 n 2 -ε N tail j=1 n∈M k -1 n κ n [X n cos(k n • r j ) -Y n sin(k n • r j )]. ( 86
)
For N tail /N fixed, the coupling parameter ε scales like N -1/2 tail when N goes to infinity. This keeps constant the typical size of the coupling term for random phases and/or particle positions in this limit. It can be checked that, on top of the total energy E sc = H sc , the total momentum

P sc = N tail j=1 p j + n∈M I n k n , (87)
is conserved, which is natural since the system is isolated. The total momentum is made of a wave contribution and of a particle one, while the total energy has a coupling contribution on top of these two. This invalidates the interpretations of Landau damping using energy conservation, while neglecting this coupling. We further notice that waves are now explicit degrees of freedom, while the electric potential is slaved to the particles in the One Component Plasma.

In references [START_REF] Antoni | Explicit reduction of N-body dynamics to self-consistent particle-wave interaction[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF] the one-dimensional analog of H sc was obtained by a direct mechanical reduction of degrees of freedom starting with the one-dimensional N -body problem 30 . This derivation had been introduced in [START_REF] Escande | Description of Landau damping and weak Langmuir turbulence through microscopic dynamics[END_REF][START_REF] Escande | Large scale structures in kinetic plasma turbulence[END_REF]. In chapter 2 of reference [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF], this derivation was performed rigorously by making error estimates involving three small parameters. These parameters control the slow evolution of the wave amplitudes, the slow evolution of the average velocity of a given bulk particle in presence of the waves, and ensure that the bulk particles are in the linear regime

B. Spontaneous emission and Landau damping

For the sake of brevity, we do not develop here the full generalization of the one-dimensional analysis in Refs [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF] ; it is lengthy, but straightforward. However, since this analysis unifies spontaneous emission with Landau growth and damping, we provide the generalization to three dimensions of the result of [START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF] and of section 4.1.4 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF] ruling the evolution of the amplitude of a Langmuir wave. It is provided by perturbation calculations where the right hand sides of equations ( 82)-( 83) are considered small (of first order). This is natural for equation [START_REF] Van Kampen | The dispersion equation for plasma waves[END_REF] since N tail N bulk , and for equation [START_REF] Kaufman | Reformulation of quasi-linear theory[END_REF] if the Langmuir waves have a low amplitude. Let J(m, t) = ϕ(m, t)ϕ(-m, t) , where the average (mathematical expectation) is over the random initial positions of the tail particles (their distribution being spatially uniform and pairwise independent). Then a calculation to second order in ϕ yields dJ(m, t)

dt = 2γ L (k m )J(m, t) + S m spont , (88)
where γ L (k m ) is the Landau growth or damping rate defined by equation ( 22), and S m spont is given by32

S m spont = 2πN e 2 ε 2 0 [ ∂ bulk ∂ω (m, ω m )] 2 k 5 m g ω mr k m . ( 89
)
S m spont corresponds to the spontaneous emission of waves by particles and induces an exponential relaxation of the waves to the thermal level in the case of Landau damping (the three-dimensional analogue of what was found in [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF]). This spontaneous emission was already present for a given realization of the N -body system in the beginning of section III A, by terms enabling the excitation of a wave with an initially vanishing amplitude. It was implicitly present in the sum of the shielded potentials of equation [START_REF] Firpo | Equilibrium statistical mechanics for single waves and wave spectra in Langmuir wave-particle interaction[END_REF]. Indeed, the space-time average of the square of the corresponding electric field provides the estimate of spontaneous emission [START_REF] Bekefi | Collective emission processes in plasmas[END_REF].

C. Diffusion and friction coefficients

This subsection and the next ones deal with the one-dimensional case and review, in particular, the main results of [START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF] and of chapters 4, 8, and 9 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]. The corresponding Hamiltonians are given by equations ( 85) and ( 86), where p j , r j , n, k n are scalars. In order to go to the limit of a continuous wave spectrum, we define an interpolating function J(k) such that

J(k n ) = I n L 2π . ( 90
)
Then the continuous spectrum limit reads

j I j0 • → •I j L 2π dk = •J(k)dk. (91)
In the case of a broad spectrum of Langmuir waves, a second order perturbation calculation for the particles (see [START_REF] Escande | Intuitive and rigorous derivation of spontaneous emission and Landau damping of Langmuir waves through classical mechanics[END_REF] and sections 4.1.3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]), similar to the above one for the waves, yields the diffusion and friction coefficients of a Fokker-Planck equation ruling the particle dynamics. This equation can be written in the compact way

∂f ∂t = ∂ ∂p D QL (p) ∂f ∂p - ∂ ∂p (F c f ), (92)
with

D QL (p) = πk(p)ε 2 κ 2 (k(p)) ω p J(k(p)), (93)
F c (p) = - k 2 (p)Lε 2 κ 2 (k(p)) 4ω p , (94)
where κ(k) and k(p), are respectively κ n and k n with n chosen such that respectively |k -k n | and |p -ω n /k n | are minimum. The Fokker-Planck equation displays two terms on the right hand side.

The first one involves only the wave amplitude through the diffusion coefficient, and the second one is a friction term due to the previously introduced spontaneous, or Cherenkov, emission. The former is the so-called quasilinear diffusive term [START_REF] Drummond | Nonlinear stability of plasma oscillations[END_REF]125]. The latter accounts for the dynamical friction due to the spontaneous emission responsible for S j in equation [START_REF] Laplace | Essai philosophique sur les probabilités[END_REF].

Equation ( 92) is coupled to a wave evolution equation, which is the analogous of equation ( 88)

in one dimension dJ(k, t) dt = 2γ L (k)J(k, t) + S spont (k), (95)
with

S spont (k) = N e 2 κ 2 (k) ε 0 L 2 k 3 g ω r (k) k . ( 96
)
The second order perturbation calculations leading to the quasilinear results make sense only over time scales ∆t such that τ c ∆t τ spread , where τ c and τ spread are respectively the correlation time of the wave potential and the time over which the particle positions spread over k -1 , where k is a typical wave number, because of the diffusion of velocities. These two times are defined as

τ c = (k∆u) -1 (97)
with ∆u the spread in particle and phase velocities, and

τ spread = 4 k 2 D QL0 -1/3 , (98)
with D QL0 a typical value of D QL (p).

When particles diffuse, the local momentum conservation underlying Landau damping and growth no longer corresponds to the synchronization with a single wave. Therefore, this effect corresponds to a coherent mechanism for isolated waves and to a (quasilinear) diffusive one for a broad spectrum.

The coupled equations ( 88) and ( 92) may be used to describe the nonlinear evolution of the kinetic beam-plasma instability, as long as the above perturbation theory remains correct. They are called the quasilinear equations in the literature when the Cherenkov terms are omitted, which occurs naturally in a Vlasovian setting.

The latter equations imply that the initial bump (corresponding to the beam) in the distribution function flattens into a plateau (so that ∂f /∂v = 0 over the corresponding velocity interval). The presence of the Cherenkov term implies a relaxation toward a distribution function with a negative slope on a longer time scale. This is the prelude to the relaxation toward a new thermal equilibrium for the whole plasma. 33 The system ( 92), [START_REF] Liang | Weak turbulence theory of Langmuir waves : a reconsideration of the validity of quasilinear theory[END_REF], without the terms due to spontaneous emission, was first introduced in 1962 in the framework of the Vlasov-Poisson description of the saturation of the weak warm beam instability, by Vedenov,Velikhov,and Sagdeev [125], and by Drummond and Pines [START_REF] Drummond | Nonlinear stability of plasma oscillations[END_REF]. These two works dealt with the plasma in a quasilinear way, as they neglected mode coupling for the wave growth and considered the particle evolution to be close to a ballistic one 34 . This enabled the introduction of the quasilinear diffusion coefficient long before chaos became fashionable in physics. Furthermore, this introduction was first done for the self-consistent problem.

D. Saturation of the weak warm beam instability

The theory of the saturation of the weak warm beam instability is at the origin of the N -body approach of this review paper. This is why we briefly recall the corresponding historical background in the next paragraph. The following one summarizes results about the chaotic transport of a particle in a prescribed spectrum of waves, which are needed to tackle the full self-consistent problem. The next one shows there is a depletion of nonlinearity when the distribution is a plateau, and the last one summarizes the results of numerical simulations of the instability.

a. Historical background Quasilinear theory shows that the weak warm beam instability saturates by the formation of a plateau in the distribution function [START_REF] Drummond | Nonlinear stability of plasma oscillations[END_REF]125]. This agrees with experimental observation [START_REF] Roberson | Experimental test of the quasilinear theory of the gentle bump instability[END_REF].

During the above defined spreading time τ spread , the velocity diffuses by the amount ∆v spread = 4(D QL0 /k) 1/3 , whose size is defined by the turbulent wave spectrum, which contains an energy bounded by the weak beam energy. Therefore, ∆v spread is typically a quantity much smaller than the width of the plateau. This implies that the saturation time is much longer than τ spread , which turns out to be also the (Lyapunov) time of separation of nearby orbits in the chaos induced by the waves (see for instance section 6.8.2 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]). Therefore, the formation of the plateau can be proved without quasilinear theory: it comes from the chaos induced by the unstable Langmuir waves among the resonant particles, whatever be the precise description of the corresponding chaotic transport. 33 This further relaxation cannot be described by the self-consistent Hamiltonian, since the latter corresponds to a given bulk. 34 Rigorously speaking, the quasilinear equations had already been mentioned in 1961 by Romanov and Filippov [112]. This reference makes the Ansatz of a Fokker-Planck equation for particle evolution and computes the corresponding diffusion coefficient; it describes the evolution of the Langmuir wave amplitude as the result of spontaneous and stimulated emission of quanta and estimates the corresponding coefficients.

When the motion of particles is chaotic, the perturbative approach used in the usual derivations, in particular in section IX C, of the quasilinear equations cannot be justified. Indeed, waves scatter the particle positions far away from their ballistic value. At time τ spread , when the corresponding spreading of positions becomes on the order of the wavelength, the perturbative approach fails.

Since leveling out the beam distribution function needs a time much longer than τ spread , one may doubt the validity of quasilinear equations to describe the saturation of the instability.

The validity of quasilinear theory was first questioned in 1979 by Adam, Laval, and Pesme [START_REF] Adam | Reconsideration of quasilinear theory[END_REF] when accounting for nonlinear wave coupling, but the importance of this coupling was denied by Galeev et al. in 1980 [69]. In 1981, Krivoruchko et al. performed a beam-plasma experiment in a plasma column, which confirmed the quasilinear predictions to be correct, but displaying the non-quasilinear onset of field correlations (formation of coherent packets) and particle-motion correlation (formation of tails of accelerated particles and acceleration of the energy exchange between the waves and the particles) [START_REF] Krivoruchko | Experimental investigations of correlation phenomena in the relaxation of velocity-spread beam in a plasma[END_REF]. However, in 1983, Laval and Pesme proved the existence of a renormalization of the growth rate and of the diffusion coefficient during the growth of the instability [START_REF] Laval | Breakdown of quasilinear theory for incoherent 1-D Langmuir waves[END_REF], and the inconsistency of quasilinear theory due to mode coupling [START_REF] Laval | Inconsistency of quasilinear theory[END_REF]. In 1984, they proposed a new Ansatz to substitute the quasilinear one, and predicted that whenever γ Landau τ spread 1 both the wave growth rate and the velocity diffusion coefficient should be renormalized by a factor 2.2 [START_REF] Laval | Self-consistency effects in quasilinear theory: a model for turbulent trapping[END_REF] 35 .

This motivated Tsunoda, Doveil, and Malmberg to perform a new experiment with a traveling wave tube in order to have a much lower noise than in a magnetized plasma column [START_REF] Tsunoda | Experimental test of the quasilinear theory of the interaction between a weak warm electron beam and a spectrum of waves[END_REF][123][124].

It brought a surprising result: quasilinear predictions looked right, while quasilinear assumptions were completely violated. Indeed no renormalization was measured, but mode-mode coupling was not negligible at all, while it is neglected in the quasilinear approach. These results were further documented in an extension of the experiment [START_REF] Hartmann | Measurements of the weak warm beam instability[END_REF]. This sets the issue: is there a rigorous way to justify quasilinear estimates in the chaotic regime of the beam-plasma instability?

Tackling this issue in a Vlasovian setting sounded formidable. However, the theory of chaos for finite number of degrees of freedom Hamiltonian systems had been developing in the plasma physics community for more than a decade [START_REF] Escande | Contributions of plasma physics to chaos and nonlinear dynamics[END_REF][START_REF] Escande | From thermonuclear fusion to hamiltonian chaos[END_REF], and this was an incentive to address the weak warm beam-plasma instability by generalizing [START_REF] Escande | Large scale structures in kinetic plasma turbulence[END_REF]119] a model originally introduced for the numerical simulation of the cold beam-plasma instability [108, 109]. There the beam was described as a set of particles, while the wave was present as a harmonic oscillator. If one considers a wave-particle interaction occurring in a finite range of velocities [v min , v max ], then it is sufficient to include in 35 More information about this controversy can be found in [START_REF] Doxas | Numerical observation of turbulence enhanced growth rates[END_REF][START_REF] Laval | Controversies about quasi-linear theory[END_REF].

the Hamiltonian the waves with phase velocities in this interval, which defines their number M when the length L of the system is defined. This leads to the self-consistent Hamiltonian [START_REF] Kraichnan | Depression of nonlinearity in decaying isotropic turbulence[END_REF]. In retrospect, in a system where the transport due to short range interactions ("collisions") is weak, it is natural to think about plasma dynamics by working directly with classical mechanics, and by taking into account that the collective field dominates over the graininess field.

Before embarking in the study of the chaotic dynamics of the self-consistent Hamiltonian corresponding to the weak warm beam instability, three preliminary investigations sounded useful:

deriving the self-consistent Hamiltonian from the underlying N -body dynamics, recovering Landau damping/growth from this Hamiltonian, and studying the chaotic transport of a particle in a prescribed spectrum of waves, which was propaedeutic to the self-consistent case. The first investigation led to the derivations of the one-dimensional version of the self-consistent Hamiltonian recalled at the end of section IX A. The second one led to the one-dimensional version of the van Kampen-Dawson approach presented in section V. The third one brought a series of results concerning chaotic transport in a prescribed spectrum of waves, which are now briefly recalled 36 .

b. Chaotic transport While for uncorrelated phases, it is natural to expect the diffusion coefficient to converge to its quasilinear estimate from below when the resonance overlap of the waves increases, for intermediate values of the overlap, unexpectedly the diffusion coefficient turns out to exceed its quasilinear value by a factor about 2.5 [START_REF] Cary | Non quasilinear diffusion far from the chaotic threshold[END_REF]. This further suggested the possibility of a renormalization for the self-consistent case. Nevertheless, the diffusive picture for the chaos due to waves needed to be substantiated. It was shown to be right, provided adequate averages are performed on the dynamics; however, this picture is wrong if one averages only over the initial positions of particles with the same initial velocity [START_REF] Bénisti | Origin of diffusion in Hamiltonian dynamics[END_REF] (see also section 6.2 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF], and [START_REF] Escande | When can the Fokker-Planck equation describe anomalous or chaotic transport?[END_REF][START_REF] Escande | When can the Fokker-Planck equation describe anomalous or chaotic transport? Intuitive aspects[END_REF]).

These results were completed by mathematical results: individual diffusion and particle decorrelation were proved for the dynamics of a particle in a set of waves with the same wavenumber and integer frequencies if their electric field is gaussian [START_REF] Elskens | Diffusion limit for many particles in a periodic stochastic acceleration field[END_REF], or if their phases have enough randomness [START_REF] Elskens | Gaussian convergence for stochastic acceleration of N particles in the dense spectrum limit[END_REF]. The randomness of amplitudes also affects significantly the overal diffusion [START_REF] Elskens | Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes[END_REF].

The intuitive reason for the validity of the diffusive picture is given in [START_REF] Bénisti | Origin of diffusion in Hamiltonian dynamics[END_REF]: it is due to the locality in velocity of the wave-particle interaction, which makes the particle to be acted upon by a series of uncorrelated dynamics when experiencing large scale chaos. This locality of the wave-particle interaction was rigorously proved by Bénisti in [START_REF] Bénisti | Finite range of large perturbations in Hamiltonian dynamics[END_REF]. On taking into account that the 36 Experimentally, studying the chaotic transport of particles in a prescribed spectrum of waves was propaedeutic to the self-consistent case too. This led to the experimental observation of resonance overlap [START_REF] Doveil | Experimental observation of resonance overlap responsible for Hamiltonian chaos[END_REF], and of the transition from stochastic diffusion in a large set of waves to slow chaos associated to a pulsating separatrix [START_REF] Doveil | Two regimes of self-consistent heating of charged particles[END_REF].

Nonlinear resonances excited by injected waves were both observed as a "devil's staircase" [START_REF] Macor | Electrons climbing a "devils staircase" in wave-particle interaction[END_REF] and cancelled to build a barrier to transport [START_REF] Chandre | Channeling chaos by building barriers[END_REF].

effect of two phases on the dynamics is felt only after a long time when there is strong resonance overlap, it can be approximately proved that the diffusion coefficient is larger than quasilinear, but converges to this value when the resonance overlap goes to infinity [START_REF] Bénisti | Origin of diffusion in Hamiltonian dynamics[END_REF][START_REF] Escande | Proof of quasilinear equations in the chaotic regime of the weak warm beam instability[END_REF] (see also [START_REF] Escande | Quasilinear diffusion for the chaotic motion of a particle in a set of longitudinal waves[END_REF][START_REF] Escande | Proof of quasilinear equations in the strongly nonlinear regime of the weak warm beam instability[END_REF], and section 6.8.2 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]).

c. Depletion of nonlinearity when the distribution is a plateau When the plateau forms in velocity during the saturation of the weak warm beam instability, density becomes also almost uniform spatially in this range of velocities. Indeed, chaos tends at equidistributing particles all over the chaotic domain in (r, v)-space. Actually, KAM tori, bounding the chaotic domain defined by a prescribed spectrum of waves, experience a sloshing motion due to the waves. This brings a small spatial modulation to the particle density which provides a source term for the Langmuir waves. However, if the plateau is broad, the source term in equation ( 82) almost vanishes, since the particles are equidistributed spatially, and the waves keep a fixed amplitude : the self-consistency of equations ( 82)-( 83) is quenched and the wave spectrum is frozen, even when particle dynamics is strongly chaotic in the plateau domain 37 (see sec. 2.2 of [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical confirmation[END_REF]). Clumps of particles may experience a strong temporary trapping, but the distribution function stays flat. As a result, the plateau dynamics is almost the same as in a prescribed field of Langmuir waves 38 . This is an instance where nonlinear effects increase the symmetry of the system, and lead to a depletion of nonlinearity 39 .

Then, it is possible to use the tools of 1.5 degree-of-freedom Hamiltonian chaos mentioned in the above paragraph "Chaotic transport" to compute the diffusion of particle velocities.

In a Vlasovian description, the bump-on-tail instability saturates with the previous plateau substituted with a very jagged distribution in both space and velocity, resulting from the chaotic stretching and bending of the initial beam-plasma distribution. Indeed, since the initial distribution f 0 (v) is conserved along particle motion, there is no finite range in velocity where f 0 (v) has the amplitude of the plateau. This plateau can be obtained only by coarse-graining (local averaging) of the Vlasovian distribution of the saturated state.

37 Equations ( 82)-( 83) are used here to avoid writing their one-dimensional analog, which is the relevant one for this discussion. 38 For a plateau with a finite width, the small remaining source brings a further evolution of the wave-particle system toward a state where the wave spectrum collapses toward small wavelengths together with the escape of initially resonant particles towards low bulk plasma thermal speeds [START_REF] Firpo | Equilibrium statistical mechanics for single waves and wave spectra in Langmuir wave-particle interaction[END_REF]. This corresponds to a further step toward a new thermal equilibrium of the N -body system corresponding to the initial beam-plasma system. The description of the subsequent steps toward thermal equilibration require to use a full N -body model. 39 This phenomenon, also called depression of nonlinearity, was introduced in fluid mechanics [START_REF] Kraichnan | Depression of nonlinearity in decaying isotropic turbulence[END_REF]. In Navier-Stokes turbulence, the mean-square value of the nonlinear term of the equation was found significantly depressed, i.e. smaller than the same quantity in the Gaussian field with the same energy distribution. This was identified to result from the emergence of long-lived vortices where the enstrophy cascade is inhibited. It also exists in systems with quadratic nonlinearities [START_REF] Bos | Reduction of mean-square advection in turbulent passive scalar mixing[END_REF][START_REF] Kraichnan | Depression of nonlinearity in decaying isotropic turbulence[END_REF].

In the above paragraph "Chaotic transport", we recalled that, in the chaotic regime, D/D QL may cover a large range of values [START_REF] Cary | Non quasilinear diffusion far from the chaotic threshold[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]. In particular D D QL is obtained for random phases of the waves and strong resonance overlap [START_REF] Cary | Non quasilinear diffusion far from the chaotic threshold[END_REF][START_REF] Elskens | Quasilinear limit for particle motion in a prescribed spectrum of random waves[END_REF][START_REF] Elskens | Nonquasilinear evolution of particle velocity in incoherent waves with random amplitudes[END_REF][START_REF] Elskens | Gaussian convergence for stochastic acceleration of N particles in the dense spectrum limit[END_REF][START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Elskens | Diffusion limit for many particles in a periodic stochastic acceleration field[END_REF]. The plateau regime corresponds to γ L = 0 and therefore to γ L τ spread = 0. Since D/D QL may cover a large range of values in this regime, γ L τ spread 1 does not imply per se any renormalization or non-renormalization of D/D QL (nor of γ/γ L by wave-particle momentum conservation). This contradicts previous works using γ L τ spread 1 to try and prove the validity of quasilinear theory [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF][START_REF] Escande | Proof of quasilinear equations in the chaotic regime of the weak warm beam instability[END_REF][START_REF] Liang | Revisiting the validity of quasilinear theory[END_REF][START_REF] Liang | Weak turbulence theory of Langmuir waves : a reconsideration of the validity of quasilinear theory[END_REF] along with the "turbulent trapping" Ansatz aiming at the contrary [START_REF] Laval | Self-consistency effects in quasilinear theory: a model for turbulent trapping[END_REF]. The value of D/D QL in the plateau regime of the bump-on-tail instability depends on the kind of wave spectrum the beam-plasma system reaches during the saturation of the instability, and not only on condition γ L τ spread 1, as assumed by these works.

d. Numerical simulations

The difficulty of dealing analytically with the strongly nonlinear regime of the Vlasov-Poisson system led from 1989 to the development of the finite number of degrees of freedom approach using the self-consistent Hamiltonian (85). This enabled numerical simulations of the self-consistent dynamics to be performed [START_REF] Cary | Enhancement of the velocity diffusion in longitudinal plasma turbulence[END_REF][START_REF] Doxas | Numerical observation of turbulence enhanced growth rates[END_REF].

Let τ w be the typical time for the growth of the wave amplitudes. In the regime τ w τ spread the system (94)-( 96) was found to be correct, but for one realization the wave spectrum appeared to be jagged with respect to the average one, each wave having a temporal behaviour strongly marked by nonlinear wave coupling. This confirmed the behaviour found by Theilhaber, Laval and Pesme

[120], and Berndtson [START_REF] Berndtson | Analysis of velocity diffusion of electrons with Vlasov-Poisson simulations[END_REF] with Vlasovian codes. The simulation of Doxas and Cary [START_REF] Doxas | Numerical observation of turbulence enhanced growth rates[END_REF] indicated a possible renormalization, but the wave spectrum was not dense enough for this to be completely convincing. Furthermore, the renormalization factor was much smaller than the one proposed by [START_REF] Laval | Self-consistency effects in quasilinear theory: a model for turbulent trapping[END_REF].

In 2011, taking advantage of the increased power of computers, more precise numerical simulations were performed using a semi-Lagrangian code for the Vlasov-wave model [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical confirmation[END_REF]. This model is the mean-field limit of the granular dynamics defined by the self-consistent Hamiltonian: waves are still present as M harmonic oscillators, but particles are described by a continuous distribution function (which is discretized in the numerical scheme, though). The simulations were benchmarked in various ways. In particular, the conservation laws were checked, as well as the above depletion of nonlinearity when the distribution is a plateau. They were repeated for a large number of random realizations of the initial wave phases for a fixed initial spectrum of amplitudes.

As shown by previous simulations, the final wave spectrum was found to be quite jaggy, and not smooth as that predicted by QL theory [START_REF] Drummond | Nonlinear stability of plasma oscillations[END_REF]125].

For each of the realizations, one computed the spreading of the velocities of test particles when acted upon by the final set of waves. The first four even moments of this spreading were compared with those of the solution to the quasilinear Fokker-Planck equation for velocity diffusion, using the velocity-dependent diffusion constant D QL computed with the final wave spectrum. The agreement was found to be excellent: the plateau verified the predictions of QL theory. However, as found in previous numerical simulations and experimentally, mode-mode coupling was found to be very strong during the saturation, which invalidated the QL assumptions. Similar results were obtained in [START_REF] Volokitin | Velocity diffusion in plasma waves excited by electron beams[END_REF] with a code using the dynamics of the self-consistent Hamiltonian.

At this point, the validity of QL predictions while QL assumptions are wrong sounded still like a mystery. However, the simulations brought an unexpected clue to elucidate it: the variation ∆φ(t) = φ(t) -φ(0) of the phase φ(t) of a given wave, was found to be almost non fluctuating with the random realizations of the initial φ(0)'s of the waves [START_REF] Besse | On the validity of quasilinear theory[END_REF] 40 . Since ∆φ(t) does not depend on φ(0), the randomness of the final wave phases is the same as that of initial phases. As a result, the self-consistent dynamics was shown to display an important ingredient for the validity of a quasilinear diffusion coefficient for the dynamics in a prescribed spectrum. Some analytic support was brought to this finding by using a third order Picard iterate of the dynamics [START_REF] Escande | Application of Picard iteration technique to selfconsistent wave-particle interaction in plasmas[END_REF].

E. Dynamics with a single wave

The single wave case is relevant to two kinds of situation. First, in the unstable beam-plasma system, one wave may have a larger growth rate than the others, so that it soon dominates over these and the single wave approximation is reasonable : this is the case for a cold beam. Second, the physical device of interest, because of its finite length, may have a single wave being resonant with a beam, even a warm one.

If a single wave interacts with all particles, the locality in velocity of the wave-particle interaction is particularly important. In the nonlinear regime of the beam-wave system, it appears that only particles with velocities close to the wave phase velocity, up to the resonance width, interact significantly with the wave, and that the relevant velocities thus depend on the instantaneous wave intensity (see chapter 8 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]).

Numerical simulations considered the growth and saturation of a single wave whose phase velocity lies in the range with positive slope of the velocity distribution function of a warm beam 40 This is reminiscent of the fact that the initial correlations were not disturbed in the course of the relaxation in the beam-plasma experiment [START_REF] Krivoruchko | Experimental investigations of correlation phenomena in the relaxation of velocity-spread beam in a plasma[END_REF]. [START_REF] Doveil | Trapping oscillations, discrete particle effects and kinetic theory of collisionless plasma[END_REF][START_REF] Firpo | Long-time discrete particle effects versus kinetic theory in the self-consistent single-wave model[END_REF]. They showed that the dynamics with a finite number N of particles leads to a long-time evolution of a single wave qualitatively distinct from the numerical integration of the kinetic Vlasov model. Indeed, while the latter integration leads the amplitude to saturate with an amplitude A, the dynamics with a finite number N of particles brings a slow second growth of the amplitude after an apparent saturation at the value A.

This occurs despite the fact that the limit N → ∞ formally reduces the many-particle evolution equation to the Vlasov kinetic equation (see section G.1.2 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]), but there is no contradiction.

Kinetic limit theorems show that the limit N → ∞ commutes with the dynamics over any finite time interval 0 ≤ t ≤ t * , but the discrepancy between both evolutions may well diverge exponentially as t * → ∞, especially because the wave-particle system has dynamical instabilities [START_REF] Firpo | Kinetic limit of N -body description of wave-particle self-consistent interaction[END_REF]. Then, in the long term, given a finite number of particles, there is a time beyond which the smooth solution to the kinetic model (approximating the finite N system initially) may evolve significantly differently from the physical finite-N system : the limits t → ∞ and N → ∞ need not commute [START_REF] Escande | Description of Landau damping and weak Langmuir turbulence through microscopic dynamics[END_REF].

In short, in its evolution, the plasma eventually reminds the physicist of its microscopic granular nature, so that the Vlasov equation cannot accurately describe it over long times [START_REF] Escande | Description of Landau damping and weak Langmuir turbulence through microscopic dynamics[END_REF], another symptom of the singular Vlasovian limit. In the above single wave simulations, the granular nature of the plasma induces fluctuations of the width of the separatrix of the single wave in the apparent saturated state. This enables the wave to exchange almost trapped particles with average velocities smaller and higher than the phase velocity of the wave. However, when the wave reaches amplitude A, the beam distribution is not globally flattened, and there are still more faster particles than slower particles. Therefore, the exchange of almost trapped particles provides momentum to the wave, i.e. further growth.

Stating this result in a different way, the BGK equilibrium corresponding to the saturated Vlasovian mode does not exist because of the fluctuations induced by the granular nature of the plasma. In reality, this problem is for all Vlasovian travelling-wave BGK equilibria, as shown in [START_REF] Elskens | Finite-N dynamics admit no travelling-waves solutions for the hamiltonian XY model and single-wave collisionless plasma model[END_REF] : they cannot exist for the dynamics ruled by the single wave finite-N tail self-consistent Hamiltonian, because they would come with singularities in complex time, which are inconsistent with the travelling-wave assumption.

F. Damping with trapping results from a phase transition

As was shown by O'Neil in 1965, when the initial amplitude of a Landau damped Langmuir wave is increased, there is a threshold above which the wave amplitude enters an oscillatory regime after a small initial damping [START_REF] Neil | Collisionless damping of nonlinear plasma oscillations[END_REF]. This oscillatory regime results from the trapping of particles inside the wave. The description of this phenomenon by the self-consistent dynamics enables proving that this transition is a second order phase transition [START_REF] Firpo | Phase transition in the collisionless damping regime for wave-particle interaction[END_REF]. This is shown in a series of steps:

(i) In the many-particle limit N → ∞, the self-consistent system is described by an invariant distribution over (r, v)-space, and explicit forms can be found for several physical observables in a Gibbsian setting. (ii) For a single wave, the partition function can be computed analytically, and one distinguishes two regimes separated by a phase transition. The two thermodynamic phases are linked to the different dynamical evolutions of the system. In particular, the two regimes of wave damping (Landau damping and damping with trapping) are recovered in connection with the two phases 41 .

A striking property of the one-dimensional One Component Plasma is that it undergoes no phase transition. A phase transition in the wave-particle model was thus unexpected. However, the two systems are very different thermodynamically. Indeed, the thermodynamic treatment of the binary interaction assumes that all particles reach a "global" equilibrium, whereas the thermodynamic treatment of the wave-particle system applies to a reduced system, in which only the faster, more efficient interactions are taken into account.

G. Dynamics with two waves

The dynamics of H sc defined by equation ( 86) are extremely rich, and do not at all reduce to those of the single or of the many waves cases described above. In particular, a work with the one-dimensional analog of this Hamiltonian involving two waves coupled to many particles displays the emergence of long-lived quasi-stationary states (QSS) [START_REF] Carlevaro | Beam-plasma instability and fast particles: the Lynden-Bell approach[END_REF]. Motivated by the problem of α-particle thermalization in a burning thermonuclear plasma, this work focuses on the case where ω n 1 for the two waves, which endows them with vanishing phase velocities, and yields a dramatic importance to the coupling term. Spatially, one of the two waves is the second harmonic of the other one. In numerical simulations, at t = 0 particles are spread uniformly spatially and in velocity in the interval [-p 0 , p 0 ]. A threshold value of p 0 is found, such that below this threshold the system responds to the beam injection by the emergence of long-lived QSS: the clustering of the beam particles into resonant rotating clumps. This threshold is also present for the case of a single wave. The values of the two coupling constants select a leading harmonic, and the passage 41 The envelope equation of an electron plasma wave has a sudden variation when going from the linear to the trapping regimes, in a way similar to a first order phase transition [START_REF] Bénisti | Self-consistent theory for the linear and nonlinear propagation of a sinusoidal electron plasma wave. application to stimulated raman scattering in a non-uniform and non-stationary plasma[END_REF]. Furthermore, there are other aspects of non collisional damping for a wave having trapped electrons [START_REF] Bénisti | The various manifestations of collisionless dissipation in wave propagation[END_REF].

from a leading one to the other one corresponds to a first-order phase transition. This transition can be described analytically by applying Lynden-Bell's technique to the Vlasovian limit of the problem [START_REF] Carlevaro | Beam-plasma instability and fast particles: the Lynden-Bell approach[END_REF].

X. COULOMB SCATTERING

This section reviews the theory developed in [START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF], which shows that an N -body model where electrons interact through their Debye shielded Coulomb potential, enables the calculation of Coulomb scattering describing correctly all impact parameters b, with a convergent expression reducing to Rutherford scattering for small b. This derivation explains why a two-body calculation yields a correct estimate of Coulomb scattering, while most of this transport is due to the simultaneous action of many particles with impact parameters between the inter-particle distance and the Debye length.

The word "collision" is borrowed from the physics of gases where it qualifies the close encounter of two particles. However, the interaction of particles in a plasma corresponds seldom to twobody collisions, even when taking into account Debye shielding : in the plasmas considered here, where the interparticle distance λ id is much smaller than the Debye length λ D , a particle j feels the simultaneous unshielded short-range action of many particles. Except for those particles very close to j, this action produces a slow and simultaneous deflection of j. Rigorously speaking, one should speak about "short range induced interactions", "unshielded Coulomb interactions", or so, and not "collisions". Such a designation is a result of the development of plasma physics after that of gases, which made natural for the former to borrow concepts and tools from the latter. In particular, the unshielded interactions of particles in kinetic plasmas were considered as collisions.

These interactions, occurring at scales λ D or smaller, set a bound to the Vlasovian description, and require another specific one.

About sixty years ago, two groups at UC Berkeley's Radiation Laboratory simultaneously studied transport due to collisions in non-magnetized plasmas, and they quoted each other's results

in their respective papers : one in 1956 by Gasiorowicz, Neuman and Riddell [START_REF] Gasiorowicz | Dynamics of ionized media[END_REF] and, in 1957, one by Rosenbluth, MacDonald and Judd [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF]. The second group of authors used the Rutherford picture of two-body collisions, while the first group of authors dealt with the mean-field part of the interaction by using perturbation theory in electric field amplitude 42 . Later on, within the 42 Their derivation suggested to the first author of the present review that a direct N -body approach might be possible.

same approximations as [START_REF] Gasiorowicz | Dynamics of ionized media[END_REF], a more elegant derivation of the Coulomb scattering coefficients was provided, in a "post-Vlasovian" approach, by taking the limit "infinite number of particles in the Debye sphere" of the Balescu-Lenard equation (see section 8.4 of [START_REF] Balescu | Statistical dynamics: matter out of equilibrium[END_REF] and sections 7.3 and 7.4 of [START_REF] Hazeltine | The framework of plasma physics[END_REF]), though the rigorous foundation of this equation is still a challenge [START_REF] Lancellotti | From Vlasov fluctuations to the BGL kinetic equation[END_REF]118].

Each of the above works on Coulomb scattering has a difficulty in describing the interactions at distances of the order of the typical interparticle distance λ id . Indeed, the mean-field approach cannot describe the graininess of these scales, and the Rutherford picture cannot describe the simultaneous collisions with several particles. Consequently, the mean-field approach is suited to describing scales larger than λ id , and should be used with a corresponding ultraviolet cutoff, while the Rutherford picture holds for scales smaller than λ id , and should be used with a corresponding infrared cutoff. Fortunately, in both approaches the transport coefficients depend only logarithmically on these cutoffs. Furthermore, forgetting about the latter ones, and considering in both cases the scales typically between λ ca and λ D , the two results are found to agree [START_REF] Gasiorowicz | Dynamics of ionized media[END_REF][START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF]. This provided confidence in these complementary extrapolations, which were for long the basis of the description of Coulomb scattering in plasmas, as presented in many plasma physics textbooks.

However, till 2015 a calculation of the contribution of scales about λ id to Coulomb scattering had been missing, and no theory provided a calculation of this transport covering all scales between λ ca and λ D . This gap was filled by reference [START_REF] Escande | Uniform derivation of Coulomb collisional transport thanks to Debye shielding[END_REF], which computes the trace T D of the velocity diffusion tensor of a given particle by a convergent expression including the particle deflections for all impact parameters. The main ideas of the new derivation are (i) the substitution of the Coulomb potential of a particle with its Debye-shielded potential, i.e. the substitution of the bare potential with its "dressed" one defined by equation ( 68), (ii) the computation of Coulombian deflections by first order perturbation theory in the total electric field, except for those due to close encounters, (iii) the contribution to T D of the former ones is matched with that of the latter ones computed by [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF]. The detailed matching procedure includes the scale of the inter-particle distance, and is reminiscent of that of [START_REF] Hubbard | The friction and diffusion coefficients of the Fokker-Planck equation in a plasma[END_REF], however without invoking the cancellation of three infinite integrals. It leads to the same expression as [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF], except for the Coulomb logarithm which is modified by a velocity dependent quantity of the order of 1. More precisely, the computation of the deflection of particle l is performed in four steps.

The first step uses first order perturbation theory in the electrostatic potential, which shows the total deflection to be the sum of the individual deflections due to all other particles. Indeed, to this order δ ṙl (t) = j∈S;j =l δ ṙlj (t), [START_REF] Meyer-Vernet | Aspects of Debye shielding[END_REF] where S denotes the set of integers from 1 to N labeling particles, and δ ṙlj (t) is the contribution of particle j to the change of velocity of particle l. This yields

δ ṙl (t) 2 = j∈S;j =l δ ṙlj (t) 2 , ( 100
)
where the average is performed over the initial random positions of particles, which kills the terms involving the contributions of two different particles. Therefore, though being due to the simultaneous scattering of particle l with the many particles inside its Debye sphere, δ ṙl (t) 2 turns out to be the sum of individual two-body deflections for impact parameters such that first order perturbation theory is correct.

For an impact parameter b much smaller than λ D , the main contribution of the acceleration due to particle j to the deflection of particle l comes from times when this acceleration takes on its bare Coulombian value. Therefore, δ ṙlj (t) is a first order approximation of the effect on particle l of a Rutherford collision with particle j. The perturbative calculation is seen to be correct for λ D b λ ca . This explains why the contribution to δ ṙl (t) 2 of this range of b's can be computed as if it would result from successive two-body collisions, as was done by [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF] and in many textbooks.

The second step of the computation of the deflection of particle l proves that for a close encounter with particle j, the deflection of particle l is exactly the one it would undergo if the other N -2 particles were absent. The contribution of such collisions to δ ṙl (t) 2 was calculated by [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF]. Now, since the deflection of particle l due to particle j as computed by the above perturbation theory is an approximation to the Rutherford deflection for the same impact parameter, one may conversely approximate the perturbative deflection with the full Rutherford one, and obtain an obvious matching of the theories for λ id b ∼ λ ca and for λ D b ∼ λ id : one may thus use the estimate of [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF] in the whole domain b λ D .

The third step shows that the deflection for an impact parameter of the order of λ D is given by the Rutherford expression multiplied by some function of the impact parameter reflecting shielding.

This enables a good matching of the deflections for large impact parameters with those for smaller ones.

The final and fourth step obtains an analytic expression for deflection whatever the impact parameter, by taking advantage of the fact that the individual deflections due to impact parameters b exceeding λ D decay rapidly with b. One finds that the Coulomb logarithm ln(λ D /λ ca ) of the second equation ( 14) of [START_REF] Rosenbluth | Fokker-Planck equation for an inverse-square force[END_REF] becomes ln(λ D /λ ca )+C where C is of order unity. The same approach can provide the calculation of the other elements of the diffusion tensor and of dynamical friction, the latter requiring second order perturbation theory. The tensors corresponding to electron-ion collisions and to ion-ion collisions can be computed likewise.

As was shown in section VIII, Debye shielding results from the Coulombian deflections of particles usually called "collisions". In turn, the resulting Debye shielded effective potential yields a description of pair interaction which provides a direct calculation of particle deflections, viz. of Coulomb scattering. Shielding and collisions are thus intimately linked, and the present ability of a thorough calculation of Coulomb scattering rests on this link.

A startling aspect of collisions in plasmas is that, although each particle interacts simultaneously with many other ones on the Debye length scale (suggesting the need for a collective description), the transport effect of these interactions is well approximated by a sum of independent binary estimates, because the deflections are so weak that they can be treated perturbatively.

XI. DISCUSSION

A. New physical picture of microscopic plasma physics

The N -body approach brings a new physical picture of microscopic plasma physics. First, it

shows that collisions play an essential role in collisionless plasmas. In particular, Debye shielding is a direct consequence of collisions (section VIII). Indeed, the Coulombian deflections of electrons by a given electron P decrease the number of electrons about P, which decreases its apparent charge, according to Gauss' theorem. In contrast to what occurs if electron P is substituted with a Langmuir probe, this does not change the density of charges of the plasma, because the deflections due to P are compensated by those of the other electrons. Debye shielding results from a cooperative dynamical self-organization process, produced by the accumulation of almost independent Coulomb deflections over a time-scale ω -1 p . It is now clear how a given particle can be shielded by all other particles, while contributing to their individual shieldings. Unexpectedly on the basis of the implicit picture of most basic textbooks, shielding and Coulomb scattering are two aspects of the same two-body repulsive/attractive process. The N -body approach to Coulomb scattering shows there is a smooth connection between impact parameters where the two-body Rutherford picture is correct, and those where a collective description is mandatory (section X).

The essential role of collisions was also exemplified in section II C for distributions of particles corresponding to a set of monokinetic beams where each beam is a simple cubic array of particles: multi-beam-multi-arrays. Indeed, such distributions correspond to invariant states of the dynamics where neither collisions, nor waves are present. When perturbing the positions of particles of a multi-beam-multi-array with a dense enough set of velocities, they exponentially diverge from their initial positions, because of collisions (section V F). The finite value of the corresponding exponentiation rates is crucial for the equivalent of the van Kampen phase mixing to occur in the N -body system (section V B). Finally, Landau damping occurs after a phase mixing time (section III A), which is longer than ω -1 p , the time necessary for collisions to establish Debye shielding. Therefore, this damping occurs in a plasma already organized by collisions.

The N -body approach shows unequivocally that Landau damping results from the simultaneous average synchronization of almost resonant passing particles with the wave (section III B). It is because this synchronization is the same for Landau growth and damping, that a single formula applies to both phenomena (section III B). The phase mixing of many beam modes produces Landau damping, which cannot correspond to a damped eigenmode because of Hamiltonian mechanics (section VII). This phase mixing is also active for Landau growth (section V E). If there is a too low number of particles in a range of velocities proportional to the damping rate, about the phase velocity of a wave, the system behaves as a multi-beam in this range, and no Landau damping occurs (section V E). When particles diffuse, Landau damping and growth no longer correspond to the synchronization with a single wave, but result from this diffusion (see section IX C). If there is a plateau in the distribution function and the waves with phase velocities in the range of this plateau make particle dynamics chaotic in this range, these waves have a stationary amplitude because of a depletion of nonlinearity due to the chaotic dynamics of particles: there is no component of the electron density resonating with these waves (section IX D). This implies that the case with a plateau is very different from one where there is a vanishing number of particles in the same velocity interval, while Landau damping vanishes in both case in a Vlasovian setting.

The N -body approach incorporates spontaneous emission naturally (sections III A and IX B).

When this emission is taken into account, Landau damping is nothing but a relaxation mechanism to the thermal level of Langmuir waves (section IX B). This emission prevents a plateau on the tail of the velocity distribution to be stationary, and triggers the walk toward the eventual thermal distribution (section IX C).

The travelling-wave BGK equilibria of the Vlasovian case do not exist in the N -body one, because of the fluctuations induced by the granular nature of the plasma (section IX E). In particular, a beam-plasma system with a single unstable wave, after reaching the equivalent of the Vlasovian saturation, evolves on a longer time-scale toward a saturated state where the wave has a larger amplitude (section IX E). O'Neil's damping with trapping typical of initially large enough Langmuir waves results from a phase transition (section IX F).

B. New insight into the Vlasovian limit

The N -body approach reveals two important features of the Vlasovian limit: it is singular and it corresponds to a renormalized description of the actual N -body dynamics. The singularity of the Vlasovian limit shows up: (i) in the dielectric function: that of a multibeam-multi-array does not converge toward the Vlasovian expression when the density of the beam velocities increases, since both zeros and poles do not match (section VII), (ii) in the importance of phase mixing for Landau growth in the N -body approach (section V E); van Kampen's theory displays a view of the dynamics of a plasma closer to the genuine one than Landau's, (iii) in the impossibility to give a physical interpretation to the term of initial conditions (70) in Landau's calculation of Langmuir waves, which is nothing but the continuous limit of the sum of the ballistic potentials of the N electrons, (iv) in the fact that adding a test particle to the N -body system does not provide the shielded potential of this particle (section V G), in contrast with the Vlasovian case, (v) in the usual requirement of an analytically continuable velocity distribution function too (see section V C; also its footnote, which shows that a softer requirement is possible).

One of the simplest examples of a renormalized potential is the Debye shielded potential [START_REF] Mccomb | Renormalization methods: a guide for beginners[END_REF].

It is a mean-field potential produced by the Coulomb deflections of the particles (section VIII).

Both its mean-field and BBGKY derivations show that Vlasov equation deals with a mean-field potential (section VI). Therefore, the Vlasovian dielectric function is a renormalized version of that of a multi-beam-multi-array, and a Vlasovian Langmuir wave is the renormalized version of a set of beam modes of the N -body system (section VI). The renormalized dielectric function enables the calculation of the shielded potentials of the N particles of the granular plasma considered here, or of a test particle added to a Vlasovian plasma (section VI).

XII. CONCLUSION

Laplace's dream was not a mere utopia, since the calculation of classical orbits starting from prescribed initial conditions can genuinely describe and explain many phenomena of microscopic plasma physics : Landau damping and growth, Debye shielding, Coulomb scattering, etc... Many of the calculations can be done for a single realization of the plasma by using standard tools of elementary mechanics, calculus and no probabilistic setting. This provides a stronger ground to face the complexity of plasmas, a difficult issue [49,[START_REF] Escande | Complexity and simplicity of plasmas[END_REF]. An alternative title of this review paper might be "Mechanical foundations of microscopic plasma physics".

Vlasov equation and calculations à la Landau have proved to be very efficient tools of theoretical plasma physics. However, computing and understanding are very different, as shown by the almost two decades the plasma physics community took to accept the reality of Landau damping.

To introduce basic microscopic plasma phenomena, the N -body approach is short, since it avoids the introduction of a kinetic equation, and provides a clear-cut interpretation of the derived phenomena. So, Landau damping can be taught or given as an exercise to students knowing Newton's second law of motion, Fourier series, but neither Vlasov equation, nor Laplace transform; in particular, to students studying Newtonian mechanics. The additional knowledge of Laplace transform makes Debye shielding accessible too. Since it is rigorous and starts from first principles, the Nbody approach could make plasma physics more attractive to colleagues in other disciplines and to prospective students who are fundamentally minded -even more so, because of the issues of singular limits, of cooperative self-organization, of renormalization, and of the depletion of nonlinearity. Furthermore, basic microscopic plasma physics has both a description corresponding to a singular limit, the Vlasovian one, and a non-singular one, the N -body approach: exceptionally, reductionism works for this physics.

Plasma physicists might enjoy the intuitive mechanics, the unifications and simplifications of the N -body approach, and the new insight into the Vlasovian limit it provides. Retrospectively, Landau's derivation of his damping was the first, but is the less physically intuitive; indeed, the Vlasovian limit is a singular and renormalized one. The N -body approach is also the occasion to consider the many facets of what is usually called "the distribution function". Most theoretical calculations do not specify which aspect of this function is considered. In particular, is it thought in a statistical setting describing an ensemble of plasmas, or as an idealization of a single granular distribution [START_REF] Elskens | Vlasov equation and N -body dynamics : How central is particle dynamics to our understanding of plasmas ?[END_REF]? In the N -body approach, this distribution is introduced at the end of the mechanical calculations. Finally, plasma physicists might enjoy the full description of the contribution of all impact parameters to Coulomb scattering (section X). In particular, most of this contribution corresponds to a sum of independent binary estimates, while each particle interacts simultaneously with many other ones on the Debye length scale.

It would be useless to rewrite the whole of plasma physics in an N -body setting. However, some extensions might turn useful. In particular, to clarify the issue of Debye shielding in a magnetized plasma. Indeed, in the direction perpendicular to the magnetic field, the repulsive electric field due to a particle provides an E × B drift of the "colliding" ones. What is the corresponding shielding in the direction perpendicular to B? Is it the same as in a non-magnetized plasma where Coulomb deflections are present 43 ? Using equations ( 8) and (101), the linearized particle dynamics defined by equation ( 9

This defines a system of N linear differential equations whose coefficients are entire functions of t.

Therefore, the ∆r l (t)'s are entire functions. Through equation ( 101), this property is transferred to φ(m, t) and to the amplitude A(t) = φ(m, t) exp(iωt) for ω real : A(t) is an entire function.

If the k m • v j 's are multiples of a given number, equation ( 103) has coefficients with some period T , and belongs to the Floquet class of differential equations. Then its solutions are of the type

U (t) = V (t)e αt , (104)
where V (t) is a vector of period T , and α a complex number. The corresponding Floquet exponents are the β σ 's introduced at the end of section V C. A similar equation was met in the self-consistent wave-particle approach introduced in section IX A (section 3 of [START_REF] Elskens | Microscopic dynamics of plasmas and chaos[END_REF]). Its full solution turned out to be a superposition of wave-like and ballistic solutions. This is remarkable, since such equations are generally not explicitly solvable with elementary functions, even for the simplest one, the Mathieu equation.

b. Solution to all orders By expanding A(t -τ ) in Taylor series, the second term of equation

The time evolution of A(t), and the frequency ω, are then derived from

1 = ω 2 p A +∞ n=0 i n n!k n+1 d n A(t) dt n P g (n+1) (v) Ω dv + i π k g (n+1) (v φ ) . (106)
Looking for a solution A(t) = A(0)e νt with ν real, equation ( 106) reads

1 = ω 2 p +∞ n=0 i n n!k n+1 ν n P g (n+1) (v) Ω dv + i π k g (n+1) (v φ ) , (107)
Since g (n+1) (v) ∼ v -n T g (v), equations ( 108) and (109) may be seen as expansions in the small parameter ε ≡ ν/(kλ D ω p ). Then, these equations may be solved order by order. The first two orders were used in section III A.

It is also noteworthy that, when g(v) is analytic, equation (107) reads

1 = ω 2 p k P g (v + iν/k) Ω dv + i π k g (v φ + iν/k) , (110)
which is exactly the formula obtained by Landau.

provided that its complex amplitude did not change much during a time interval of the order of τ mix = (kv T ) -1 . More precisely, the typical wave growth rate, γ, had to be such that |γ|τ mix < 0.1.

We now further discuss the importance of the product |γ|τ mix in terms of the transients. For the sake of simplicity, the discussion is restricted to the situation when the electrons only feel the effect of the drive during a finite time interval, namely when 0 < t < τ d (see Ref. [START_REF] Bénisti | Nonlinear plasma response to a slowly varying electrostatic wave, and application to stimulated raman scattering[END_REF] for a calculation that does not make use of this hypothesis). We also assume that the total force, including the effect of the self-consistent field and the drive, derives from an effective potential which reads ϕ = A(t) exp[i(k m .r -ωt)] + c.c. (which has been proved to be correct in Ref. [START_REF] Bénisti | Nonlinear plasma response to a slowly varying electrostatic wave, and application to stimulated raman scattering[END_REF] when the plasma wave was laser driven). Because we include the effect of the drive, we can now integrate the electrons motion from the time when they are at equilibrium, r j = r j0 when t = 0. Then, Eq. ( 14)

is
which is the same as Eq. ( 18) except that the term proportional to h(v) no longer appears, since the calculation has been performed with δr j = 0. Hence, unlike in Eq. ( 18), no transient is expected from this term. Now, when t > τ d it is clear that A(t) is nothing but the amplitude of the selfconsistent potential, since the electrons no longer feel the effect of the drive. Moreover, because it is slowly driven, only when τ d τ mix may the self-consistent field reach a significant amplitude, and may effectively be Landau damped. Hence, when t > τ d , Eq. ( 112) leads to Eq. ( 19). This means that, if the wave is slowly driven, it is Landau damped just after the drive has been turned off, and there is no transient.

Let us now investigate the situation when the wave may no longer be considered as slowlyvarying when t < τ d . Then, Eq. ( 112) is no longer valid, and we have to pay a specific attention to the electron motion during the driving phase (t ≤ τ d ). To do so, we still assume that the electrons are at equilibrium when t = 0. Moreover, we use Eq. [START_REF] Roberson | Experimental test of the quasilinear theory of the gentle bump instability[END_REF], which is exact, to calculate the shift in their positions at t = τ d , induced by the drive and the self-consistent electric fields.

As noted in Section III A, the phase ψ 0 has no importance in the derivation, and we henceforth drop it. Using the result of Eq. ( 113), we now calculate the shift in position when t > τ d , which will let us conclude about the evolution of the wave amplitude. Using again Eq. ( 111), we find ∆r j 1 (t) = δr j sin(k m • r j0 ) + αk m

where we introduced

A(t) ≡ A(t) exp[i(Ω j t + k m • r j0 )]. (115)
In Section III A, the limit τ d → 0 was considered, which yields Eq. ( 14). In this Appendix, we specify how small τ d has to be for the results of Section III A to be valid. First, we want the third term in the right-hand side of Eq. (114) to be negligible, which is only true if t τ d and if A(t)

does not abruptly vanish within a time interval smaller than τ d . The first condition implies that the results of Section III A are only valid when t τ d . The second condition is only true provided that γ L τ d 1.

Neglecting the third term of Eq. ( 114) and following the same steps as in Section III A, one finds that Eq. ( 18) is changed into

A(t
Then, Eq. ( 18) is recovered only in the limit Ωτ d → 0. Hence, this equation is only valid provided that the wave could reach a significant amplitude, due to the external drive, during a time much smaller than τ mix and the plasma period. When the latter condition is fulfilled, Eq. ( 18) implies that, once the drive is turned off, one has to wait for a time of the order of τ mix before Landau damping is effective.

Therefore, we recover here the usual difference in a system's response to a sudden excitation compared to an adiabatic one. When a system is subjected to a sudden change, it usually rapidly oscillates before entering a stationary, or slowly-varying, regime. These transient oscillations do not exist under an adiabatic-like external force, and the system keeps on varying in a smooth way.

Note that we only derived the wave evolution in two opposite limits, either when the driving time, τ d , was much smaller than τ mix and the plasma period, or when it was much smaller than τ mix . The situation when τ d is either of the order of the plasma period or of τ mix would deserve further investigation.

Appendix C. INFINITE NUMBER OF BEAMS

This appendix provides calculations similar to those in [START_REF] Dawson | Plasma oscillations of a large number of electron beams[END_REF], but reformulated in a form suitable for the derivation of section V, and with the correction of several errors. We first focus on d1 (m, ω) defined by equation [START_REF] Escande | Contributions of plasma physics to chaos and nonlinear dynamics[END_REF]. In order to separate the regular and the singular parts of this quantity, we add and subtract to the right hand side of equation ( 51) the quantity

ω 2 p [ π 2 g(ω/k) k 2 δ sin 2 (πω/kδ) - 2πg (ω/k) k 2 cot(πω/kδ)], (117)
where k = k m . This yields

d1 (m, ω) =1 -ω 2 p π 2 g(ω/k) k 2 δ sin 2 (πω/kδ) - 2πg (ω/k) k 2 cot(πω/kδ) +Σ ∞ σ=-∞ ( [g(σδ) -g(ω/k)]δ (ω -σkδ) 2 + 2g (ω/k)ωδ k(ω 2 -(σkδ) 2 ) ) , (118)
where use has been made of the relations (see [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF]) 44 π 2 sin 2 (πx) = Σ ∞ σ=-∞ 1 (x -σ) 2 (119)

π cot(πx) = Σ ∞ σ=-∞ x (x 2 -σ 2 ) . ( 120
)
We notice that the function of σδ inside the summation in equation (118) has no real poles.

Therefore, the sum passes smoothly to an integral as δ goes to zero 45 . In this limit, equation (118

) becomes d1 (m, ω) =1 -ω 2 p π 2 g(ω/k) k 2 δ sin 2 (πω/kδ) - 2πg (ω/k) k 2 cot(πω/kδ) + ∞ -∞ ( g (v) k(ω -kv) + 2g (ω/k)ω k(ω 2 -k 2 v 2 ) )dv . ( 121
)
We now compute the zeros of d1 (m, ω) and write ω = α + iβ. We first consider those with β vanishing when δ goes to zero. If β vanished like or faster than δ, the first term in the bracket of equation (121) would diverge, while the second and third one would remain finite, which is impossible. Therefore, β vanishes slower than δ, which forces the cotangent to converge toward -µ i, where µ = ±1 is the sign of β. Then, in order to stay finite, the first term requires β to 44 Equation (120) corrects a typo in equation ( 33) of [START_REF] Dawson | Plasma oscillations of a large number of electron beams[END_REF]. 45 We assume g (v) continuous, |g(v)| to be integrable, and g(v) ≥ 0.

k 2 δ + iµπg (nδ) k 2 ] = 1 + P ∞ -∞ ω 2 p g (v) k 2 (nδ -v) dv, (122)
for nδ in the support of g, with α 1 = α -nkδ. Equation (121) provides two contributions to the term in g (nδ): one from the term in cotangent, and one from the pole of the term in g (v) in the integral, while the two poles of the term in g (ω/k) bring contributions cancelling each other.

Solving for α 1 and β yields

tan 2πα 1 kδ = - πω 2 p k 2 g (nδ)/[1 + P ∞ -∞ dv ω 2 p g (v) k 2 (nδ -v) ] , (123)
β = µ kδ 2π ln k 2 δ 4π 2 ω 2 p g(nδ) ×   1+ P ∞ -∞ dv ω 2 p g (v) k 2 (nδ -v) 2 + πω 2 p g (nδ) k 2 2   . (124)
Equation (123) yields 2πα 1 /kδ modulo π, and the right solution is obtained by requiring cos(2πµα 1 /kδ)

to have the opposite sign to the denominator of this equation.

Like the natural frequencies of the beams, the roots are spaced kδ apart in α. Therefore, the above zeros have real parts between these natural frequencies. There are two roots for each beam, since µ can be either positive or negative, and no root away form the support of g. Thus we obtain two modes for each beam, as required.

We now compute ∂ ∂ω (m, ω) defined in equation [START_REF] Escande | Proof of quasilinear equations in the chaotic regime of the weak warm beam instability[END_REF] for the case of a vanishing imaginary part of ω when δ goes to zero. Here again, we handle the singularity by adding and subtracting to the right hand side of equation ( 54 These terms can also be written in the form of sums by using again equations ( 119) and (120), and

(see [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF])

π 3 cos(πx) sin 3 (πx) = Σ ∞ σ=-∞ 1 (x -σ) 3 . ( 126
)
Using the latter expression, for δ small we find that σ,µ is given by equation ( 56).

  using a logarithmic scale, for a plasma with a density n = 10 19 m -3 and a temperature T = 1 keV. It rules out strongly coupled and/or degenerate plasmas.

FIG. 1 :

 1 FIG.1: View in logarithmic scale of the different scale lengths discussed in this section.

FIG. 2 :

 2 FIG.2: One-dimensional multi-beam-multi-array.

  and

FIG. 3 :FIG. 4 :

 34 FIG. 3: (Left) B as a function of t/Ω for t = 10 3 , showing the decreasing amplitude when |Ω| → 0. (Right) Magnified view of the rectangle indicated on the left panel, showing the linear dependence on t/Ω, as well as the decreasing oscillations when |Ω| grows.

  of the Langmuir waves à la van Kampen-Dawson and à la Landau. Since b smooth λ D , b smooth is an intermediate scale between λ id and λ D , which exists provided Λ is large enough (see section II A). Then, equation (47) describes the response of the plasma to the initial perturbation defined by ϕ (bal) (m, ω), and d (m, ω) is the dielectric function of the granular plasma.

FIG. 5 :

 5 FIG. 5: (Left) Distribution function. The position of the phase velocity is indicated by a vertical arrow. (Right) Contour plots of the modulus of the granular dielectric function. The calculations are done for δ = 0.5 (top panels) and δ = 0.1 (bottom panels). The poles identified as Van Kampen poles get closer to the real axis as δ decreases. Only two poles, identified as Landau poles and responsible for the Landau instability, exhibit a finite imaginary part when δ → 0.

  using a smoothed version of the Coulomb potential. Because of the cutoff at b smooth , the collision frequency 23 involves a Coulomb logarithm ln(λ D /b smooth ) instead of the usual one ln(λ D /λ ca ).

  Using the smoothed version of Coulomb potential recovers the shielded Coulomb potential of equations (68)-(69) down to a distance ∼ b smooth from the particle of interest. Since b smooth λ D , a part of the 1/r dependence of the genuine Coulomb potential is recovered, and can be matched with the central 1/r dependence for closer distances.

m

  The iteration starts with the ballistic approximation of the dynamics r (0) m = r m0 + v m t, where r m0 and v m are respectively the initial position and velocity of particle m. The actual orbit of equation (9) corresponds to n → ∞. Let ξ be the mismatch of the position of particle m with respect to the ballistic one at the n-th iterate, viz. the effect of Coulomb interactions to that order of iterations ; because of the initial conditions, ξ

  of oscillation. Such a Hamiltonian was first introduced by Mynick and Kaufman [100], and derived in a consistently Hamiltonian way from the Vlasov-Poisson system of equations by Tennyson et al. [119] 31 .

  k m φ(m, t) exp[ik m • (r j0 + v j t)] + c.c. .(102)Because of ordering[START_REF] Bekefi | Collective emission processes in plasmas[END_REF], for each beam the corresponding values of exp[ik m • r l0 ] are uniformly distributed on the unit circle, and their global contribution to the 1 factor in the last bracket of equation (101) vanishes. Therefore equation (102) yields the compact expression∆r j =m • [r j0 -r l0 + (v j -v l ) t ] k m • ∆r l (t).

( 16 )

 16 becomes (with k = k m and v φ = ω/k)) n+1 d n A(t) dt n (-1) n+1 g (n+1) (v) exp(-iΩt) -n+1 d n A(t) dt n g (n+1) (v) exp(-iΩt) -) (v φ ) .

(- 1 )

 1 p (2p)!k 2p+1 ν 2p P g (2p+1) (v) Ω dv -πν (2p + 1)k 2 g (2p+2) (v φ ) , p ν 2p (2p + 1)!k 2p νP g (2p+2) (v) Ω dv + π(2p + 1)g (2p+1) (v φ ) .

  changed into ∆r j1 (t) = αk m t 0 τ A(t -τ ) exp[i(Ω j (t -τ ) + k m • r j0 )]dτ + c.c. (111)If the wave is slowly driven, A(t) is a slowly-varying function, and one may stop the Taylor expansion of A(t -τ ) at first order, which yields ) g (v) cos(Ωt) -i sin(Ωt) -1 Ω dv,

  This yields∆r j1 (t) = αk m τ d 0 τ A(τ d -τ ) exp[i(Ω j (τ d -τ ) + k m • r j0 )]dτ + c.c. ≡ δr j sin(k m • r j0 + ψ 0 ).

  (t -τ ) -A(τ d -τ )] dτ + c.c. ,

  )g (v) e -iΩτ d -cos(Ωt) + i sin(Ωt) ) g (v) cos(Ωt) -i sin(Ωt) -e -iΩτ d Ω dv.

  scale like δ | ln(δ/v T ) |. With this in mind, and looking for solutions in the vicinity of nkδ, equation (121) requires ω 2 p [-4π 2 g(nδ) exp[2π(iµα 1 -|β|)/kδ]

π 3

 3 cos(πω/kδ)g(ω/k) k 3 δ 2 sin 3 (πω/kδ) -π 2 g (ω/k) k 2 δ sin 2 (πω/kδ) + πg (ω/k) 2k 3 cot(πω/kδ)].(125)

  , adds

	i	πω 2 p 2k 4 m	πω 2 p k 3 m km ) Ä(t) in the right hand side of this equation. The next orders can be obtained from equations (108)-g ( ω km ) Ȧ(t) -g ( ω
	(109).

  the component of rsec1 along k m has the sign of B. For |Ω|t 1, B -Ωt 3 /3 , which implies a stronger average synchronization when |Ω| grows at fixed t. This also shows this synchronization to vanish with |Ω|, which rules out any role of trapped particles in Landau damping or growth 14 . One easily sees that B keeps the sign of -Ω for Ωt of the order of a few units. Therefore, since Landau damping occurs over a time scale γ -1 L (k m ), synchronization occurs for |Ω| up to the order of |γ L |. For |Ω| |γ L |, B scales like t/|Ω| over the time scale γ -1 L (k m ). Consequently, the average synchronization is maximum for |Ω| of the order of |γ L |. More precisely, B performs decreasing oscillations between positive and negative values when |Ω| grows at fixed t.

This rigorous derivation is accessible to students knowing Newton's second law of motion and the Fourier transform, but neither analytic functions, nor the Laplace transform. It also provides a correction to the lowest order expression of Landau damping. In a second step, the calculation is extended to the corresponding particle dynamics, showing that it produces an average synchronization of almost resonant passing particles with the wave.

This smoothing is similar to the one performed in the mean-field derivation of the Vlasov equation, which was recalled in the introduction.

In the similar spirit of considering a wave with a slowly varying amplitude, but by using smooth distribution functions, references[START_REF] Bénisti | Nonlinear plasma response to a slowly varying electrostatic wave, and application to stimulated raman scattering[END_REF][START_REF] Bénisti | Nonlinear landau damping rate of a driven plasma wave[END_REF] derive not only Landau damping, but also its equivalent in the nonlinear regime.

Actually, a possible role of trapping is a priori excluded since the bounce period is unbounded in the linear regime of Langmuir waves.

We stress that the derivation of Landau damping in this section is completely independent of that in section III though.

Therefore, taking this singular limit is a way to bring to students a rapid derivation of both Landau damping and Debye shielding.

However, this has no practical consequence, since the time scale over which the system is observed, or measured, introduces a natural granularity in frequency space below which the actual continuum (or singular) limit and the discretized system should not be distinguishable.

In contrast, the global exchange of energy and momentum between waves and particles is easily available in a Vlasovian setting.

This was the starting point of the one-dimensional N -body approach.

We notice that the derivations of the self-consistent dynamics starting with a Vlasovian description[108, 109, 119] perform kind of a zigzag with respect to the N -body description, since they go back to a finite number of degrees of freedom after going through the continuous Vlasovian description.

This expression corrects equation (46) of[START_REF] Escande | Direct path from microscopic mechanics to Debye shielding, Landau damping and wave-particle interaction[END_REF].

ACKNOWLEDGMENTS D. F. E. is grateful to the members of Equipe Turbulence Plasma in Marseilles, since the theory reviewed in this paper is the result of three decades of collaboration with them. He thanks Professor M. Kikuchi for suggesting him to write this review. He also thanks Professor A. Sen for many useful suggestions, and Professor P. Huneman for pointing out to him the book "Reductionism, emergence and levels of reality" by Chibbaro et al. He thanks Drs F. Bonneau, M.-C. Firpo, and F. Sattin for helpful comments on the manuscript. Also D. F. G. Minenna who brought the precious views of a newcomer in the field. One of the authors (D. Z.) has been supported by the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the "Investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR).

Appendix A. RIGOROUS DERIVATION À LA KAUFMAN a. A(t) is entire The dynamics considered in section III corresponds to the linearized motion of the N electrons with respect to a given multi-beam-multi-array, when a single wave with wavevector k m is excited at t = 0. Let r j0 be the initial position of the unperturbed beam particle with index j, and v j be its velocity, and let ∆r j (t) = r j (t) -r j0 -v j t be the mismatch of the actual position of particle j with respect to the unperturbed beam particle with the same index.

Setting r j = r j0 + v j t + ∆r j (t) in equation ( 8), we replace φ with its expansion to first order in the ∆r

(101) 43 The conspicuous modification of relaxation processes by magnetic fields [START_REF] Baldwin | Magnetized plasma kinetic theory. II. derivation of Rosenbluth potentials for a uniform magnetized plasma[END_REF][START_REF] Hassan | Magnetized plasma kinetic theory. I. Derivation of the kinetic equation for a uniform magnetized plasma[END_REF][START_REF] Hassan | Magnetized plasma kinetic theory. III. Fokker-Planck coefficients for a uniform magnetized plasma[END_REF][START_REF] Ichimaru | Relaxation processes in plasmas with magnetic field. Temperature relaxations[END_REF][START_REF] Silin | On relaxation of electron and ion temperatures of fully ionized plasma in a strong magnetic field[END_REF], suggests that shielding might be modified too. They also suggest the importance of the relative ordering of spatial scales, such as Larmor radius versus the Debye length.

Appendix B. RELEVANCE OF TRANSIENTS BEFORE LANDAU DAMPING

In this Appendix, we discuss the relevance of transients that would occur before the selfconsistent electrostatic field may experience Landau damping. To do so, we have to specify how the self-consistent electrostatic field has actually been generated, an issue that is usually eluded.

This implies that we do account for the external drive (e.g. a laser, a polarized grid, electrodes. . . ) used to induce the self-consistent field in the plasma, when calculating the electron motion. The calculation is, therefore, slightly more general than that leading to Eqs ( 14)- [START_REF] Besse | On the validity of quasilinear theory[END_REF]. Such a generalization has already been performed in Ref. [START_REF] Bénisti | Nonlinear plasma response to a slowly varying electrostatic wave, and application to stimulated raman scattering[END_REF] when the self-consistent field was slowly driven.

An important point of Ref. [START_REF] Bénisti | Nonlinear plasma response to a slowly varying electrostatic wave, and application to stimulated raman scattering[END_REF] was to prove that a wave may be considered as slowly varying,

Appendix D. APPROACHING THE SINGULAR LIMIT BY COARSE-GRAINING

In the introduction, we recalled the pulverization procedure for deriving Vlasov equation from the BBGKY hierarchy. The singular limit can be obtained by a coarse-graining procedure, which is germane to the pulverization procedure: each particle is substituted by a continuum of particles with velocities close to its velocity, with a mismatch in velocity ∆v distributed with the continuous distribution P (∆v), instead of a discrete distribution in the case of the pulverization 46 . This procedure may be viewed as the counterpart in velocity of the quantum regularization of small spatial scales for collisions recalled in the introduction. Indeed, the coarse-graining in velocity may appear as a way to account for the quantum uncertainty on the particle velocities.

The calculation leading to equation ( 47) can be performed again, but the summation over particles now involves an integral over the nearby velocities of the coarse-grained system. Then

and

where n b is the number of beams. If the width of P is large with respect to the edge of an elementary cube of the velocity grid, f 0 is a smooth function, the grid may be taken as very tight, and for practical purposes P (u) may be taken as a Dirac distribution in equation [START_REF] Yoon | Weak turbulence theory for collisional plasmas[END_REF], which then becomes equation [START_REF] Elskens | Vlasov equation and N -body dynamics : How central is particle dynamics to our understanding of plasmas ?[END_REF].

Appendix E. SHIELDED COULOMB POTENTIAL BY A SINGULAR LIMIT OF THE MANY-BEAM DESCRIPTION

If in equation (65) we take first the limit δ → 0 and then the limit ε → 0 + , according to equation

This is nothing but the contribution of (m, k m • u + iε) in the same limit (see for instance equation (9.12) of [START_REF] Nicholson | Introduction to plasma theory[END_REF]). Therefore, for δ small enough, equation [START_REF] Escande | N -body description of Debye shielding and Landau damping[END_REF] becomes equation [START_REF] Firpo | Equilibrium statistical mechanics for single waves and wave spectra in Langmuir wave-particle interaction[END_REF], and equation 46 More precisely, the pulverization leads to infinitesimal beamlets whose summation corresponds to a coarse-graining of the previous multi-beam-multi-arrays.