{'fraction_non_alphanumeric': 0.03807078672290482, 'fraction_numerical': 0.010526315789473684, 'mean_word_length': 4.50383528461849, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
N-body description of Debye shielding and Landau damping

This paper brings further insight into the recently published N-body description of Debye shielding and Landau damping [Escande D F, Elskens Y and Doveil F 2014 Plasma Phys. Control. Fusion 57 025017]. Its fundamental equation for the electrostatic potential is derived in a simpler and more rigorous way. Various physical consequences of the new approach are discussed, and this approach is compared with the seminal one by Pines and Bohm [Pines D and Bohm D 1952 Phys. Rev. 85 338--353].


I. INTRODUCTION
In order to discuss microscopic plasma physics, most textbooks start with invoking the N -body description of plasmas, but they deem it too difficult to tackle, and they rapidly switch to kinetic equations. The kinetic approach has proved extremely powerful to compute a wealth of plasma effects, but it often hides some aspects of their mechanical nature ; for instance, after the Vlasovian derivation of Landau damping, most textbooks use simple mechanical models to provide physical interpretations of this damping.
However, the N -body description of plasmas has proved for more than six decades to be very powerful at describing analytically effects like Debye shielding, Landau damping, spontaneous emission of Langmuir waves, and at revealing the corresponding mechanical behavior of both particles and waves. This started in 1952 with a seminal paper by Pines and Bohm [27]. From the nineties on, there was a revival of the N -body approach to describe analytically wave-particle interaction in plasmas [2,12].
Last year, reference [15] provided a direct and simultaneous derivation of Debye shielding and Landau damping by an N -body description of the plasma. This derivation works directly with Newton's second law and with a version of the electrostatic potential linearized from ballistic orbits of the particles ; it does not use probabilistic arguments or partial differential equations. However, in retrospect, the formulation of this paper now appears to be unnecessarily intricate, because it was still too close to the heuristic derivation of the theory ; it also required an erratum [16].
The heuristic derivation substituted the discrete summations over the particles, which are intrinsic to the N -body approach, with integrals over a continuous velocity distribution, and was presented at the beginning of reference [15]. In retrospect, this substitution turns out to be an efficient shortcut to the simultaneous derivation of Debye shielding and Landau damping, but this non rigorous step can be avoided by considering the initial N -body distribution as close to a set of monokinetic arrays of particles. This other approach was presented in a second part of reference [15]. In strong resonance with Dawson's seminal multi-beam fluid theory [7], this second approach shows that Landau damping is due to phase mixingà la van Kampen [15], while the first approach leads to the Vlasovian expression introduced by Landau that hides this fact.
Sections II and III (and the appendix) provide a shorter and more rigorous version of the approach of reference [15] using a smooth distribution function. The calculation reduces to a mere exercise using Fourier series and Laplace transform. Section IV discusses the various physical consequences of the present N -body approach. Section V is devoted to the comparison of the present N -body approach with the 1952 one by Pines and Bohm [27]. This latter work is often overlooked by plasma physics textbooks and was overlooked also by the present authors in reference [15]. The main technical difference of the present theory with reference [27] is the use of the Laplace transform. It turns out to have the same heuristic power as in the Vlasovian approach, and enables a compact and straightforward simultaneous derivation of Debye shielding and Landau damping without any a priori physical intuition.

II. SIMPLER DERIVATION OF THE FUNDAMENTAL EQUATION FOR THE ELECTROSTATIC POTENTIAL
As in reference [15], we deal with the One Component Plasma (OCP) model [1,4,29], 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. We now provide a shorter and more rigorous version of the derivation of reference [15] using a smooth distribution function. The actual calculation involves sums over the N particles and over all Fourier components of the potential. However, the principle of this calculation can be explained more simply by performing it formally for the unphysical case of a single electron acted upon by a single Fourier harmonic of its own field. This is done in the next subsection, while the complete derivation is in appendix.
A. Single particle acted upon by a single Fourier harmonic of its own field While the Fourier decomposition of the periodic delta distribution corresponding to the charge density of an electron entails an infinity of harmonics, in this subsection we keep just a single one with wavevector k which provides through Poisson equation the single harmonic of the electrostatic potential where k = k , ǫ 0 is the vacuum permittivity, and r(t) is the position at time t of the particle acting as a source with charge −e. Let r 0 and v be respectively the initial position and velocity of the particle, and let ∆r(t) = r(t) − r 0 − vt. We now compute a perturbative solution to the dynamics. To this end, setting r(t) = r 0 + vt + ∆r(t) in equation (1), we replaceφ with its expansion to first order in ∆r(t) We now use the time Laplace transform which maps a function g(t) to g(ω) = ∞ 0 g(t) exp(iωt)dt (with ω complex). Since the arguments of functions are spelled explicitly, from now on we omit diacritics for the Laplace (or Fourier) transformed quantities. The time Laplace transform of equation (2) is where the Doppler shift −k · v comes from the linear dependence on t of the exponent of equation (2). To compute ∆r(ω − k · v), we use Newton's equation, which writes r = e m e ∇ϕ lin (r) (4) for the particle with mass m e . The single harmonic of the electrostatic potential defined by equation (1) implies where r now is the position of the particle as a massive object subjected to a force. In order to derive our perturbative solution, we consider k ∆r to be small in equation (2), so that the linearized particle dynamics defined by equation (4) reduces to Since ∆r(0) = 0 and ∆ṙ(0) = 0, the time Laplace transform of equation (6) is where the Doppler shift k · v comes from the linear dependence on t in the exponent of equation (6). Computing ∆r(ω − k · v) in equation (3) from the expression of ∆r(ω) given by equation (7) yields where ϕ (bal) is the time Laplace transform of the ballistic potential obtained by setting ∆r(t) = 0 in equation (2). We keep on purpose the terms (k − k) · r 0 and k · v − k · v due to the two Doppler shifts, in order to exhibit the structure of the equation when the full set of wavevectors is taken into account.

B. Fundamental equation for the electrostatic potential
As derived in the appendix, when the full set of wavevectors and particles is taken into account, equation (8) where ω p = [e 2 n/(m e ǫ 0 )] 1/2 is the plasma frequency (with n = N/L 3 ), where ϕ lin (m, ω) is the time Laplace transform of the Fourier component of the electrostatic potential with wavevector k m = 2π L m, and m = (m x , m y , m z ) is a vector with three integer components (we comment below on the restriction k n b smooth ≤ 1 making k n run over a bounded domain). Note that g(m, ω) * = g(−m, −ω * ) for real m and complex ω if g(r, t) is real-valued. We also defined with the ballistic potential related to particle j, i.e. the unshielded Coulomb potential of an electron with position r j (t) = r j0 + v j t.
In this complete derivation, for any finite value of the ∆r 's due to the various particles, the linearization in the ∆r's is justified for finite values of the wavenumber moduli k m = k m only. Therefore, as in section 4.5 of reference [15], we restrict the Fourier expansion of the Coulomb potential to k m 's such that the Debye length, with k B the Boltzmann constant, and T the temperature ; let us denote by M smooth the set of such wavevectors. Then, we may consider finite ∆r 's, with ∆r ≪ b smooth . The above truncation of the Fourier expansion smoothes the Coulomb potential. We notice that a similar smoothing of the 1/r singularity of this potential (or some penalization of very close approaches [22]) is performed in the mean-field derivation of the Vlasov equation [6,8,11,23,24,30]. In any case, the Vlasov equation cannot describe scales smaller than the inter-particle distance. Since the self-field due to the smoothed Coulomb potential vanishes, it is not necessary to exclude self-interaction in the complete derivation [16].
Equation (10) is the fundamental equation of this paper. This equation is of the type Eϕ lin = source term, where E is a linear operator, acting on the infinite dimensional array whose components are all the Doppler shifted ϕ lin (m, ω)'s.

III. SMOOTHED LINEAR OPERATOR OF THE POTENTIAL EQUATION
We now approximate the granular (or empirical, Kolmogorov) distribution F (r, v) = N j=1 δ(r − r j0 ) δ(v − v j ) with a position and velocity distribution function f (r, v) that is continuous in r, such that distribution f − F yields a negligible contribution when applied to space dependent functions which vary slowly on the scale b smooth . This means that we use a version of F coarse-grained on this scale. On replacing the discrete sums over particles with integrals over the smooth distribution function f (r, v) in its left hand side, equation (10) becomes where ϕ cg is the smoothed version of ϕ lin , andf is the spatial Fourier transform of f . We consider almost uniform initial spatial distributions of particles. So, f is a spatially uniform distribution function of the order of ϕ cg . In linearizing equation (13) for which vanishes for k = 0. Therefore, operator E becomes diagonal with respect to both m and ω (a complex quantity), and equation (13) becomes with The validity of the above coarse-graining was discussed in section 4.5 of reference [15]. Here we only recall the main steps of the discussion. First, the presence of many particles in the Debye sphere justifies neglecting the nondiagonal terms in equation (10). Second, Dawson's work [7] justifies the Vlasovian expression (15) for ǫ(m, ω) for a one-dimensional plasma. The generalization of Dawson's calculation to the three-dimensional case, though tedious, does not a priori involve any new conceptual difficulty. These results back up the above coarse-graining.
However, we now show these results come with conditions on the parameters of our problem. Invoking Dawson's calculation means that, when interested in a wave with phase velocity v w (with v w = v w ), there are many particles with parallel velocities in the 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 wavenumber. Indeed, this range is the one where the phase mixingà la van Kampen is occurring, and where the synchronization of particles with the wave brings the change of particle momentum inducing Landau damping or growth of this wave [9,10,12,13,17]. If the initial distribution of particles is meant to approximate a smooth distribution f 0 (v), this means that the length L of the system is taken large enough to fulfill the condition "many particles with parallel velocities in the range [v w − |γ L |/k, v w + |γ L |/k]". However, whatever L, for v large enough, one should not use ǫ(m, ω), but its discrete counterpart computed with the N -body distribution function.

IV. DEBYE SHIELDING, LANGMUIR WAVES AND LANDAU DAMPING
To keep this paper self-contained, we now discuss how Debye shielding and Landau damping were derived from equations (14)- (15) in reference [15]. These equations are obtained by introducing the continuous velocity distribution after particle dynamics is taken into account, and not before, as occurs when kinetic equations are used. However, they look very similar to those including initial conditions in Landau contour calculations of Langmuir wave growth or damping, usually obtained by linearizing Vlasov equation and using Fourier-Laplace transform, as described in many textbooks. Indeed, the smoothed self-consistent potential ϕ cg is determined by the response function ǫ(m, ω), viz. the classical plasma dielectric function, whose Vlasovian expression involving ∂f 0 /∂v obtains by a mere integration by parts if f 0 is differentiable. For a cold plasma, ǫ(m, ω) = 1 − ω 2 p /ω 2 , as expected (however this case is ruled out since the present theory is correct only for a large number of particles in the Debye sphere, as recorded in the previous section).
In equation (14), ϕ lin (m, ω) is a discrete sum over the N particles. As computed in [3,19,28], in the inverse Laplace transform of ϕ cg (m, ω), the pole ω = k m ·ṙ j (0) of ϕ (bal) lin (m, ω) brings a contribution to ϕ cg (m, t) which is the shielded Coulomb potential of particle j. This potential is often computed in textbooks by adding a test particle to a Vlasovian plasma (see for instance chapter 9 of [25]). Since Debye shielding is effective over distances larger than the Debye length, this discussion is relevant under the assumption that our box size is large in comparison, viz. L ≫ λ D .
As in the Vlasovian approach, the part of ϕ cg (m, t) due to dominant poles of 1/ǫ(m, ω) corresponds to one Langmuir wave propagating in the direction of k m and one propagating in the opposite direction, with the corresponding Landau growth or damping. Therefore, the full potential in the plasma turns out to be the sum of two contributions only : one due to Langmuir waves and one due to all the shielded Coulomb potentials of individual particles.
Equation (14) becomes the Vlasovian one on replacing in equation (12) the discrete sums over particles with integrals over the smooth distribution function f (r, v). This yields the Vlasovian formula for the ballistic potential. The N -body description reveals that the contribution of ϕ (bal) cg (m, ω) to ϕ cg (m, t) in the Vlasovian approach is the continuous limit of the sum of all the Debye-shielded potential of the particles. In this approach, this contribution is only a transient [21]. Though its physical meaning is unclear, it does not correspond to a robust supplementary contribution to the full electrostatic potential in the plasma which would add to the potential of Langmuir waves and to the shielded Coulomb potentials of particles.
To afford the description of nonlinear effects in wave-particle dynamics, section 6.1 of reference [15] generalizes the linear analysis analogous to that in section II of the present paper by applying it to the particles which are not resonant with Langmuir waves, while keeping the exact contribution to the electrostatic potential of the remaining N ′ tail particles that may experience trapping or chaotic dynamics. This eventually leads to equations describing the self-consistent dynamics of M Langmuir waves and of N ′ tail particles (these equations could be also obtained with the formalism of the present paper). For a one-dimensional plasma, the linear theory of Langmuir waves can be formulated on this basis, with the bonus that the information on particle dynamics is provided in parallel with the waves' [10,12,13,17]. This approach also provides spontaneous emission [12,15,17], and some nonlinear effects like damping with trapping [18], and the decoupling of the waves from the strongly chaotic motion of the particles when there is a plateau in the particle velocity distribution function (see sec. 2.2 of [5]).

V. COMPARISON WITH PINES AND BOHM 1952
In 1952, Pines and Bohm showed that an N -body approach can bring a lot of insight into basic plasma effects [27]. We now summarize the part of their results overlapping with our theory. Pines and Bohm already considered a One Component Plasma in a periodic box and wrote the equations of motion of its N electrons acted upon by their electrostatic potential. Working with the physical time t, they got an equation (number (8) in their paper) for the time evolution of a Fourier component ρ k of the charge density. This equation involves a contribution of all Fourier components whose structure is similar to equation (10), but without linearizing particle dynamics. In order to neglect the contribution of all other Fourier components, they introduced the now celebrated "random phase approximation" by showing that the particle positions may be considered as nearly random for all times if there are many particles in a Debye sphere of a plasma in thermal equilibrium.
As a first step, they considered the collective part of ρ k , and showed that, for small wavenumbers k , the charge density oscillates at the plasma frequency, and that plasma oscillations exist for k ≪ λ −1 D . Inductively, they derived a harmonic equation for non-vanishing k 's with an oscillation frequency ω assumed implicitly real, and satisfying a discrete version of the Vlasovian dispersion relation (their equation (23)). They recovered the Bohm-Gross dispersion relation for small k 's.
As a second step, in their section III they split ρ k into a collective and a granular part. The evaluation of the latter for thermal particles yielded correctly the Yukawa version of Debye shielding.
In contrast, as those presented in most textbooks, our approach to Debye shielding and Landau damping is a linear theory. Laplace transform can then be applied readily. As already demonstrated by Landau's theory, the Laplace transform is a powerful heuristic tool. In the N -body approach, it shows the full potential in the plasma to split naturally in two contributions related to two families of poles : one due to Langmuir waves subjected to Landau growth and damping, and one due to all the shielded Coulomb potentials of individual particles, whatever be their velocity. By considering the initial N -body distribution as close to a set of monokinetic arrays of particles, it was shown in Appendix B of reference [15] that linearization provides naturally the cancellation of non-diagonal terms in equation (10) ; it is not necessary to use the random phase approximation, or other statistical arguments.
Particle deflection and Gauss' theorem were used by Pines and Bohm to interpret their Yukawa version of Debye shielding. Section 5 of reference [15] brings more mechanical insight to this interpretation. There, instead of applying the Laplace transform to linearized particle dynamics, one applies Picard iteration technique (a standard method to prove the existence and uniqueness of solutions to first-order equations with given initial conditions) to the equation of motion of a particle, say P , due to the Coulomb forces of all other ones without linearizing particle dynamics. This calculation shows that particle P acts on another particle P ′ , directly by its Coulombian force, and indirectly by deflecting the orbits of all other particles whose force on particle P ′ is modified by this deflection. The direct Coulombian force is reduced by the indirect action, since the Coulombian deflections due to particle P diminish the number of other particles in a sphere surrounding it, which diminishes the apparent charge of particle P according to Gauss' theorem. This shielding effect becomes dominant at a distance on the order of λ D . 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 . This is also the time for the equilibrium pair-correlation function to settle at its time-asymptotic values. Therefore, shielding results from the accumulation of almost independent repulsive deflections with the same qualitative impact on the effective electric field of particle P . If point-like ions were present, the attractive deflection of charges with opposite signs would have the same effect. The Coulombian deflections induce a self-organization of the plasma : it becomes a dielectric. So, these deflections do not produce only the disorder leading to collisional transport, but also the order corresponding to shielding and to the behaviour as a dielectric. Paradoxically, "collisionless" Landau damping turns out to occur because of what is usually called "collisions". Since the global deflection of particles includes the contributions of many other ones, the density of the electrons does not change, at variance with what occurs for the shielding at work next to a Langmuir probe. Gauss' theorem was used in section 2.2.1 of [26] to provide an intuitive explanation of the shielding of such a probe.

VI. CONCLUSION
This paper provided further insight into the unified N -body description of Debye shielding and Landau damping. First, by deriving the fundamental equation for the electrostatic potential of reference [15] in a simpler and more rigorous way. Second, by stressing that the term from initial conditions in Landau's calculation of Langmuir waves does not bring new physics, but is a mere damped continuous version of the Debye shielded potentials of all particles whose physical meaning is unclear in the Vlasov picture. Third, by discussing more precisely the condition for approximating the N -body distribution by a smooth one in the dispersion relation. Fourth, by comparing the present approach with the 1952 one by Pines and Bohm [27] which already brought a lot of insight into Debye shielding and Langmuir waves. In particular, we have shown how the Coulombian deflections of particles bring order to the plasma by producing Debye shielding and by making the plasma behave as a dielectric. Furthermore, the Laplace transform confirms to be a powerful heuristic tool providing a compact and straightforward simultaneous derivation of Debye shielding and Landau damping. Linearization was shown to provide naturally the cancellation of non-diagonal terms in the fundamental equation for the electrostatic potential. This avoids using the random phase approximation or statistical arguments.
The N -body approach provides a derivation of Debye shielding for a single mechanical realization of the plasma. This was used in an N -body calculation of collisional transport of the One Component Plasma (OCP) incorporating all impact parameters with no ad hoc cut-off [14]. This calculation explains why a two-body calculation yields a correct estimate of collisional transport, 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.
When generalizing this theory to the collisional transport of a plasma with granular electrons and ions, the singularity of head-on collisions of particles with opposite charges can be avoided by excluding the case of a vanishing impact parameter which has probability 0. Then, the deflection of particles with same sign is just the opposite of the deflection of particles with opposite signs, so that the classical Rutherford cross-section is the same for charges of like sign and opposite signs. However, one might wonder whether a quantum mechanical treatment would not modify the scattering cross-section. Yet, it does not. This result, traditionally derived with the Born approximation, was obtained more recently with a much simpler calculation using Fermi golden rule that makes intuitive why the classical and quantum results agree [20].
The theory of collisional transport of reference [14] and the approach of the present paper show that N -body dynamics, which has always been the ultimate reference in plasma textbooks, is now also a practical tool. Furthermore, Laplace's dream was not a mere utopia, since classical mechanics can genuinely describe non trivial aspects of the macroscopic dynamics of a many-body system.
One of us (DFE) thanks the Program Committee of the 42st EPS Conference on Plasma Physics and G.L. Delzanno to have led him to look for a more pedagogical derivation of the results of reference [15] in order to prepare an invited talk and a course ; this was the origin of section II of this paper. He also thanks J. Daligault for recording him the pioneering 1952 contribution of Pines and Bohm [27] to the N -body approach to basic plasma physics, and J. Daligault and S. Khrapak for a critical reading of our manuscript. We thank the two anonymous referees for leading us to formulate our theory in a more rigorous way, and to write a clearer and more self-contained paper.
We now consider ϕ(r), the potential created by the N particles at any point where there is no particle. Its discrete Fourier transform is readily obtained from the Poisson equation, and is given byφ(0) = 0, and for m = 0 bỹ where r j (t) is the position at time t of particle j acting as a source,φ(m) = ϕ(r) exp(−ik m ·r) d 3 r, with k m = 2π L m, and k m = k m .
Let r j0 and v j be respectively the initial position and velocity of particle j, and let ∆r j (t) = r j (t) − r j0 − v j t. We now compute a perturbative solution to the full N -body dynamics and the resulting Fourier components of the potential for m running in the finite domain M smooth defined in subsection II B, such that k m b smooth ≤ 1.
To this end, setting r j = r j0 + v j t + ∆r j (t) in equation (A1), we replaceφ with its expansion to first order in the ∆r j (t)'sφ whose time Laplace transform is where the Doppler shift −k m · v j comes from the linear dependence on t of the exponent of equation (A2).
To compute ∆r j (ω − k m · v j ), we use Newton's equation for the particles r j = e m e ∇ϕ lin (r j ), where m e is the electron mass. When 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 (A1), which yields ϕ lin (r j ) = 1 L 3 m, kmb smooth ≤1φ lin (m) exp(ik m · r j ), where r j now is the position of particle j as a massive object subjected to a force. Using equation (A2), the linearized particles dynamics defined by equation (A4) is then given by ∆r j = ie L 3 m e n, knb smooth ≤1 k nφlin (n, t) exp[ik n · (r j0 + v j t)], Since ∆r j = 0 and ∆ṙ j = 0 at t = 0, the time Laplace transform of equation (A6) is − ω 2 ∆r j (ω) = ie L 3 m e n, knb smooth ≤1 k n exp(ik n · r j0 ) ϕ lin (n, ω + k n · v j ), where the Doppler shift k n · v j comes from the linear dependence on t in the exponent of equation (A6). Computing ∆r(ω − k · v) in equation (A3) from the expression of ∆r(ω) given by equation (A7) yields the fundamental equation (10) for the potential. Note that the assumption k m b smooth ≤ 1, with b smooth ≪ λ D , defining the set M smooth , excludes scales which are irrelevant to Debye shielding and Landau damping, since these phenomena involve scales larger than,