Direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction

# Direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction

D F Escande, Yves Elskens and F Doveil Aix-Marseille Université, CNRS, PIIM, UMR 7345,
case 321, campus Saint-Jérôme, FR-13013 Marseille, France
###### Abstract

The derivation of Debye shielding and Landau damping from the -body description of plasmas is performed directly by using Newton’s second law for the -body system. This is done in a few steps with elementary calculations using standard tools of calculus, and no probabilistic setting. Unexpectedly, Debye shielding is encountered together with Landau damping. This approach is shown to be justified in the one-dimensional case when the number of particles in a Debye sphere becomes large. The theory is extended to accommodate a correct description of trapping and chaos due to Langmuir waves. Shielding and collisional transport are found to be two related aspects of the repulsive deflections of electrons, in such a way that each particle is shielded by all other ones while keeping in uninterrupted motion.

PACS numbers :
52.20.-j Elementary processes in plasmas
52.35.Fp Plasma : electrostatic waves and oscillations
45.50.-j Dynamics and kinematics of a particle and a system of particles
05.20.Dd Kinetic theory

Keywords : basic plasma physics, Debye shielding, Landau damping, wave-particle interaction, spontaneous emission, amplitude equation, N-body dynamics

## I Introduction

For macroscopic classical systems, the -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 the -body description by classical mechanics, but directly derived kinetic analogues of the Boltzmann equation, in particular the Vlasov equation. This trend has been dominant till nowadays.

However, for plasmas where transport due to short range interactions is weak, -body classical mechanics yields useful results. As will be recalled in section VI, it already enabled a description of wave-particle interaction making it more intuitive, incorporating modern chaotic dynamics, and unifying particle and wave evolutions, as well as collective and finite- physics AEE (); EZE (); EEbook (); Houches (). The present paper makes an even more thorough use of -body mechanics by working directly with Newton’s second law for this system. It shows, in particular, that basic phenomena like Debye shielding and Landau damping can be easily derived while avoiding kinetic and statistical calculations altogether. Furthermore, the new derivation brings together Debye shielding and Landau damping, a totally unexpected fact in view of presentations in classical textbooks.

## Ii Main results and paper outline

Here are the main results of this paper and its organization :

1. Section III analyses a perturbation from ballistic orbits of an infinite plasma made up of the periodic replication of electrons coupled by Coulomb forces in a cubic volume with a neutralizing ionic background (One Component Plasma (OCP) model Salp (); Abe (); BH (); Kie14 ()). By using the Fourier and Laplace transforms in a way similar to that in the Vlasovian derivation of Landau damping, an equation (Eq. (14)) is derived for a version of the electrostatic potential linearized from ballistic orbits. This equation has the form , where is a linear operator, acting on the infinite dimensional array whose components are all the Doppler shifted Fourier-Laplace components of the potential. Both and the source term are sums over the particles.

2. In section IV.1, the discrete sums in are substituted with integrals over a continuous uniform distribution function . Then becomes diagonal, which brings Eq. (20), , where is the classical plasma dielectric function, familiar in a Vlasovian context. Then, the new approximate potential turns out to be the sum of two parts. The first one is the sum of the shielded Coulomb potentials of the individual particles (Eq. (22)). Such potentials were first computed by a kinetic approach in section II.A of Ref. Gasio () and later on in Bal (); Rost (). So, Debye shielding is computed for a single mechanical realization of the plasma.

3. Section IV.2 describes the second part of . It is the sum of the contributions to the excitation of Langmuir waves by the individual particles (Eq. (23)). Substituting again the discrete sums over particles with integrals over a continuous distribution function yields the total amplitude of this wave in a compact way (Eq. (24)), which works also for distributions that are non-smooth in (for instance a two-stream one).

4. In order to recover the usual expressions of textbooks, section IV.3 introduces before performing the inverse Laplace transform providing . This yields Eqs (20) and (25) which are the expressions including initial conditions in Landau contour calculations of Langmuir wave growth or damping, obtained from linearizing Vlasov equation and using Fourier-Laplace transform, as described in many textbooks (see for instance Refs HW (); Nicholson (); BoydSan ()).

5. Section IV.5 shows that the previous analysis can be performed for a smoothed Coulomb potential whose Fourier expansion is bounded to wavenumbers smaller than , where . Then the previous linearization is justified. An argument developed in B shows that the presence of many particles in the Debye sphere justifies the passage to a continuous distribution function in the one-dimensional case.

6. In section V, Picard iteration technique (one of the standard methods to prove the existence and uniqueness of solutions to first-order equations with given initial conditions) is applied to the equation of motion of a particle due to the Coulomb forces of all other ones. It stresses that a part of the effect on particle of another particle is mediated by all other particles (Eq. (31)) and reduces the direct part. Indeed, particle modifies the motion of all other particles, implying that the action of the latter ones on particle is affected by particle .

7. This calculation yields the following interpretation of shielding. At , consider a set of uniformly distributed particles, and especially particle . At a later time , the latter has deflected all particles which made a closest approach to it with a typical impact parameter , where is the thermal velocity. This part of their global deflection due to particle reduces the number of particles inside the sphere of radius about it. Therefore, according to Gauss’ theorem, the effective charge of particle as seen out of is reduced : the charge of particle is shielded due to these deflections. This shielding effect increases with , and thus with the distance to particle . It becomes complete at a distance on the order of . 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 APiel ()).

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. , where 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 (if point-like ions were present, the attractive deflection of charges with opposite signs would have the same effect). So, shielding and collisional transport are two aspects of the same two-body repulsive process.

8. In section VI, in the spirit of Refs OWM (); OLMSS (); AEE (); EEbook (), to accommodate a correct description of trapping or chaos due to Langmuir waves, the set of particles is split into bulk and tail, where the bulk is the set of particles which cannot resonate with Langmuir waves. Repeating for the bulk particles the analysis leading to Eq. (14), the same equation is recovered with an additional source term due to the tail particles (Eq. (36)).

9. Using the fact that the number of tail particles is small with respect to the bulk one, and a technique introduced in Refs OWM (); OLMSS (), an amplitude equation is derived for any Fourier component of the potential where tail particles provide a source term (Eq. (41)).

10. This equation, together with the equation of motion of the tail particles, enables one to show that, in the linear regime, the amplitude of a Langmuir wave is ruled by Landau growth or damping, and by spontaneous emission (Eq. (43)), a generalization to 3 dimensions of the one-dimensional result of Refs EZE (); EEbook ().

## Iii Fundamental linear equation for the potential

This paper deals with the One Component Plasma (OCP) model Salp (); Abe (); BH (), which considers the plasma as infinite with spatial periodicity in three orthogonal directions with coordinates , and made up of electrons in each elementary cube with volume . Ions are present only as a uniform neutralizing background, enabling periodic boundary conditions. This choice is made to simplify the analysis which focuses on , the potential created by the particles at any point where there is no particle. The discrete Fourier transform of , readily obtained from the Poisson equation, is given by , and for by

 ~φ(m)=−eϵ0k2m∑j∈Sexp(−ikm⋅rj), (1)

where is the electron charge, is the vacuum permittivity, is the position of particle , , , with a vector with three integer components running from to , , and . Reciprocally,

 φ(r)=1L3∑m~φ(m)exp(ikm⋅r). (2)

The dynamics of particle follows Newton’s equation

 ¨rl=eme∇φl(rl), (3)

with the electron mass, and the electrostatic potential acting on particle , i.e. the one created by all other particles and by the background charge. Its Fourier transform is given by Eq. (1) with the restriction in order to exclude self-interaction. Let

 r(0)l=rl0+vlt (4)

be a ballistic approximation to the motion of particle , and let . We leave some freedom to this approximation : and may be respectively the initial position and velocity of particle , or they may be slightly shifted from these values, for instance by low amplitude Langmuir waves. We now define the ballistic approximation to which is computed from Eq. (1) on setting for all ’s in the latter and excluding the -th term. We define the two mismatches to ballistic values

 Δξl(t) = ξl(t)−ξl(0)−˙ξl(0)t, (5) Δ~φl(m,t) = ~φl(m,t)−ϕ(bal)l(m,t). (6)

Until the end of section IV, we consider cases where all the ’s are small. So we approximate by its expansion to first order in the ’s (Approximation 1, discussed in A)

 Δ~ϕl(m,t)=∑j∈S;j≠lδ~ϕj(m,t), (7)

with

 δ~ϕj(m,t)=ieϵ0k2mexp[−ikm⋅r(0)j(t)] km⋅Δξj, (8)

the contribution of particle . This provides a first order approximation in the ’s to , namely

 ~ϕl(m,t)=ϕ(bal)l(m,t)+Δ~ϕl(m,t). (9)

We further consider to be small, and the ’s to be of the order of (Approximation 2). At lowest order, the particles dynamics defined by Eq. (3) is then given by

 (10)

We now use the time Laplace transform which maps a function to (with complex). Since the arguments of functions are spelled explicitly, from now on we omit diacritics for the Fourier-Laplace transformed quantities. Note that for real and complex if is real-valued. The Laplace transform of Eq. (10) is

 ω2ξl(ω)=−ieL3me∑nknexp(ikn⋅rl0) ϕl(n,ω+ωn,l)+iωξl(0)−˙ξl(0), (11)

or

 ω2Δξl(ω)=−ieL3me∑nknexp(ikn⋅rl0) ϕl(n,ω+ωn,l), (12)

where comes from the time dependence of , as defined by Eq. (4), in the exponent of Eq. (10). Equation (11) and the Laplace transform of Eqs (7) and (8) then yield

 k2mΔϕl(m,ω)=ω2pN∑nkm⋅kn∑j∈S;j≠lϕj(n,ω+ωn,j−ωm,j)(ω−ωm,j)2exp[i(kn−km)⋅rj0], (13)

where comes from the time dependence of in the exponent of Eq. (8), and is the plasma frequency.

Summing Eq. (13) over , using Eq.  (9), and dividing by yields

 (14) = k2mϕ(bal)(m,ω),

where and are respectively and complemented with the missing -th term. More explicitly,

 ϕ(bal)(m,ω)=∑j∈Sδϕ(bal)j(m,ω), (15)

where

 (16)

The division by and not by is a direct consequence of the exclusion of self-interaction terms in the ’s. Equation (14) is the fundamental equation of this paper. This equation is of the type source term, where is a linear operator, acting on the infinite dimensional array whose components are all the Doppler shifted ’s.

## Iv Debye shielding, Langmuir waves and Landau damping

We rewrite Eq. (14) by separating its diagonal terms and its non-diagonal terms as

 ϵd(m,ω)ϕ(m,ω) (17) = ϕ(bal)(m,ω)+ω2pN∑n≠mkm⋅knk2m ∑j∈Sϕ(n,ω+ωn,j−ωm,j)(ω−ωm,j)2exp[i(kn−km)⋅rj0],

where

 ϵd(m,ω)=1−ω2pN∑j∈S1(ω−ωm,j)2. (18)

### iv.1 Shielded Coulomb potential

We now approximate the granular distribution with a position and velocity distribution function that is continuous in , such that distribution yields a negligible contribution when applied to space dependent functions which vary slowly on the scale of the inter-particle distance. We further assume to be spatially uniform, and we write it (normalized with ).

We replace the discrete sums over particles with integrals over (Approximation 3). This has two consequences : operator becomes diagonal with respect to both and (which is complex), and becomes

 ϵ(m,ω)=1−ω2pL3N∫f0(v)(ω−km⋅v)2 d3v. (19)

Therefore, Eq. (17) becomes

 (20)

where is the smoothed version of (approximating ) resulting from Approximations 1 to 3.

This shows that the smoothed self-consistent potential is determined by the response function , viz. the classical plasma dielectric function. One readily checks this for a cold plasma : then , where is now computed with the plasma density . The classical expression involving obtains by a mere integration by parts if is differentiable.

As a result of Eq. (15), the part of generated by particle is . By inverse Fourier-Laplace transform, after some transient discussed later, the potential due to particle is the sum of two parts : one due to the pole , and one to the poles of .

The first part is the shielded Coulomb potential Gasio (); Bal (); Rost ()

 δΦj(r,t)=δΦ(r−rj(0)−˙rj(0)t,˙rj(0)), (21)

where

 δΦ(r,v)=−eL3ϵ0∑m≠0exp(ikm⋅r)k2mϵ(m,km⋅v+iε) (22)

with the usual prescription resulting from inverting the Laplace transform, as the integral in Eq. (19) is undefined for the real-valued . Therefore, after this transient, except for possible Langmuir waves, the dominant contribution to the full potential in the plasma turns out to be the sum of the shielded Coulomb potentials of individual particles located at their ballistic positions computed with their initial positions and velocities. This property enables a calculation of the collisional transport for all impact parameters EED14 ().

Let be the Debye length, where is the Boltzmann constant and the temperature. The wavenumbers resolving scale are such that . Shielding involves scales on the order of . The transient is given by the zeros of . For shielding scales, these zeros correspond to a strong damping over time scales on the order of the plasma period. Therefore, the transient is damped after such a period, as estimated in statement II.(7).

If , the corresponding wavenumbers are such that . Therefore, there is no shielding for , since where is the thermal velocity.

### iv.2 Langmuir waves and Landau damping

We now consider the part of the potential due to particle provided by the two dominant (Landau) poles of . These two poles correspond to one wave propagating in the direction of and one propagating in the opposite direction. For simplicity, we retain only the first one (which we call ), and associate the other one to . By retaining this pole in the inverse Laplace transform of Eq. (20) and using Eqs (2) and (15)-(16), one finds that particle brings a contribution to the wave with frequency

 (23)

where “c. c.” means complex conjugate, and .

We now approximate the granular distribution with a smooth position and velocity distribution function in the sum over of the ’s. This sum yields the total amplitude of the wave with frequency

 Φm(r,t)=−e exp[i(km⋅r−ωmt)]ϵ0k2mL3ϵ′m∫f(m,v)ωm−km⋅v d3v+c.c., (24)

which is the Vlasovian result for this wave.

Since Eqs (19) and (24) do not involve derivatives of the distribution functions, they also enable computing Langmuir waves induced by an initial perturbation in the case where they are non-differentiable (for instance a two-stream one).

The other poles of correspond to a strong damping, and bring again a transient which vanishes over time scales on the order of the plasma period. We notice that converges to 1 when goes to infinity. Therefore, the inverse Laplace transform of has a singular part, the Dirac distribution . As a result, the convolution involved in the inverse Laplace transform of Eq. (20) brings the ballistic potential which dominates transiently, as expected.

### iv.3 Vlasovian formula for the ballistic potential

We now recover the usual textbook expressions by introducing before performing the inverse Laplace transform providing . On neglecting the ’s and the ’s to lowest order in Eq. (16), this yields a whose Laplace transform is

 (25)

the smoothed version of the actual shielded potential in the plasma. This shows that this second continuous approximation makes Eq. (20) to become the expression 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. In general, textbooks do not dwell upon . The -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.

### iv.4 Comparison with kinetic approaches

It is interesting to compare the above derivation with that presented by classical textbooks when they start with the -body description to derive both Debye shielding and the combination of Eqs (20) and (25). Debye shielding is exhibited in the equilibrium pair correlation function computed after deriving the first two equations of the BBGKY hierarchy (see e.g. ch. 12 of BoydSan ()). The combination of Eqs (20) and (25) is obtained independently by linearizing Vlasov equation about a uniform velocity distribution function, and by using the Fourier-Laplace transform. A prerequisite is the derivation of Vlasov equation by either of two approaches : a mean-field derivation NeunzertWick (); Neunzert84 (); Dobru (); BraunHepp (); Spohn (); Kie14 (); ElsVla () that does not involve statistical arguments, or the BBGKY hierarchy that involves statistical arguments and starts with the Liouville equation (see e.g. Nicholson ()). Rigorously speaking, these two derivations do not give the same definition to the velocity distribution function.

The present derivation Laplace-transforms in time the linearized dynamics of a single realization of the -body system. This yields Eq. (14) which keeps the full graininess of the system. A first continuous approximation involving a velocity distribution function yields Eqs (21)-(22), and a second one yields Eq. (25) combined with Eq. (20). This provides a short connection between these equations and the underlying -body problem. In this derivation, the continuous velocity distribution is introduced after particle dynamics has been taken into account, and not before, as occurs when kinetic equations are used. This avoids addressing the issues of the exact definition of the continuous distribution for a given realization of the plasma, and of the uncertainty as to the way the continuous dynamics departs from the actual -body one NeunzertWick (); Neunzert84 (); Dobru (); BraunHepp (); Spohn (); Kie14 (); ElsVla (). Introducing the continuous velocity distribution after particle dynamics has been taken into account may bring important new physics : it exhibited that the true (granular) plasma goes beyond the Vlasovian saturation for the growth of a single wave due to a warm beam Firpo ().

### iv.5 More rigorous derivation and range of validity

Up to this point, linearization was applied without questioning its range of validity. However, for any finite value of , the linearization of Eq. (8) is justified for finite values of only. Fortunately, we notice that the above introduction of Debye shielding and of Langmuir waves involves scales larger than, or of the order of . As a result, we do not need a correct description of the dynamics due to the large wavenumbers, i.e. due to the singularity of the Coulomb potential. Therefore, we may restrict the expansion of Eq. (2) to ’s such that , where . Then the summation in Eq. (14) is bounded similarly, and we may consider finite ’s, much smaller than though. We notice that a similar smoothing of the singularity of the Coulomb potential (or some penalization of very close approaches Kie14 ()) is performed in the mean-field derivation of the Vlasov equation NeunzertWick (); Neunzert84 (); Dobru (); BraunHepp (); Spohn (); ElsVla ().

Using the smoothed version of Coulomb potential recovers the shielded Coulomb potential of Eqs (21)-(22) down to a distance from the particle of interest. Since , a part of the dependence of the genuine Coulomb potential is recovered, and can be matched with the central dependence for closer distances.

B considers a velocity distribution corresponding to a large number of monokinetic beams. It shows that, if the inter-particle distance verifies condition

 d≪bsmooth, (26)

in Eq. (17) the non-diagonal terms vanish, which is the main simplification brought by the summation over the continuous distribution function at the beginning of section IV.1. Since , then too, which is equivalent to . This shows that the presence of many particles in the Debye sphere justifies neglecting the non-diagonal terms in Eq. (17).

Then, this equation describes the response of a large number of monokinetic beams to the initial perturbation defined by . This problem was considered by Dawson in 1960 for a one-dimensional plasma Dawson60 (). He showed that brings two beam modes per beam. Their eigenfrequencies are pairs of complex conjugate values for , whose imaginary parts vanish when the spacing of the beam velocity decreases : this makes these modes analogous to the van Kampen modes. Landau damping is recovered by phase mixing of these modes. Dawson’s analysis can be directly used by restricting the plasma of interest to a one-dimensional OCP. Therefore, in this case the use of the continuous distribution function at the beginning of section IV.1 is justified when the number of particles in a Debye sphere becomes large. The generalization of Dawson’s calculation to the three-dimensional case a priori does not display any new conceptual difficulty, but is more involved and will not be given here.

The truncated Coulomb potential cannot correctly describe the encounters between particles with impact parameters smaller than , which makes the description of the dynamics relevant for times shorter than the collision time , as happens for the Vlasovian description. Both descriptions are relevant for Langmuir waves, since conditions and are equivalent.

## V Mediated interactions imply Debye shielding

In the above derivation of Debye shielding, using the Laplace transform of the particle positions does not provide an intuitive picture of this effect. We now show that such a picture is obtained directly from the mechanical description of microscopic dynamics with the full OCP Coulomb potential of Eq. (1) in real time. To compute the dynamics, we use Picard iteration technique. From Eq. (3), , the -th iterate for , is computed from

 ¨r(n)l=eme∇φ(n−1)l(r(n−1)l), (27)

where is computed by the inverse Fourier transform of Eq. (1) with the ’s substituted with the ’s. The iteration starts with the ballistic approximation of the dynamics defined by Eq. (4), and the actual orbit of Eq. (3) corresponds to . Let be the mismatch of the position of particle with respect to the ballistic one at the -th iterate, viz. the effect of Coulomb interactions to that order of iterations ; we assume all ’s and ’s to vanish identically. It is convenient to write Eq. (27) as

 ¨ξ(n)l=∑j∈S;j≠l¨ξ(n)lj, (28)

with

 ¨ξ(n)lj=aC(r(n−1)l−r(n−1)j) (29)

and

 aC(r)=ie2ϵ0meL3∑m≠0k−2mkmexp(ikm⋅r). (30)

Let . For , one finds

 ¨ξ(n)l=∑j∈S;j≠l[¨ξ(1)lj+M(n−1)lj+2∇aC(r(0)l−r(0)j)⋅ξ(n−1)lj]+O(a3), (31)

where is the order of magnitude of the total Coulombian acceleration, and

 M(n−1)lj = ∇aC(r(0)l−r(0)j)⋅[ξ(n−1)l−ξ(n−1)j−2ξ(n−1)lj] (32) = ∇aC(r(0)l−r(0)j)⋅∑i∈S;i≠l,j(ξ(n−1)li+ξ(n−1)ij), (33)

where the second expression takes into account that is anti-symmetrical in . The latter expression displays which is the deflection of particle by particle . It shows the bare Coulomb acceleration of particle due to particle is modified by the following process : particle modifies the motion of all other particles, so that the action of the latter ones on particle is modified by particle . Therefore is the acceleration of particle due to particle mediated by all other particles. The last term in the bracket in Eq. (31) accounts for the fact that both particles and are shifted with respect to their ballistic positions. Both and this last term are anti-symmetrical with respect to the labels and as is an even function of .

Since the shielded potential of the previous paragraph was found by first order perturbation theory, it is felt in the acceleration of particles computed to second order. This acceleration is provided by Eq. (31) for . Therefore its term in brackets is the shielded acceleration of particle due to particle . As a result, though the summation runs over all particles, its effective part is only due to particles typically inside the Debye sphere (with radius ) about particle . Starting from the third iterate of the Picard scheme, the effective part of the summation in Eq. (31) ranges inside this Debye sphere, since the ’s are then computed with a shielded acceleration. 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 BoydSan ()). In Eq. (33), the compound effect of the ’s, the deflection of particle by particle , is to diminish the negative charge inside a sphere centered on particle . This yields the interpretation of shielding given in statement II.(7) using Gauss’ theorem.

## Vi Wave-particle dynamics

Section IV.2 enables the calculation of Langmuir waves excited by a given initial perturbation after a continuous approximation of the true granular distribution. To describe Langmuir waves with discrete particles, we still consider that the ’s are uniformly distributed, and we allow for non-zero ’s and ’s for the ’s of section III. Therefore, in the formulas of section IV, the ’s and ’s are slightly shifted with respect to the ’s and ’s, due to Langmuir waves.

Up to this point, we described Langmuir waves by a fully linear theory. We now generalize the analysis of section III to afford the description of nonlinear effects in wave-particle dynamics. Indeed, resonant particles may experience trapping or chaotic dynamics, which imply ’s of the order of or larger for wave ’s. For such particles, it is not appropriate to make the linearizations leading to Eqs (8) and (10). 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 particles into bulk and tail, in the spirit of Refs OWM (); OLMSS (); AEE (); EZE (); EEbook (). The bulk is defined as the set of particles which are not resonant with Langmuir waves. We then perform the analysis of section III for these particles, while keeping the exact contribution of the remaining particles to the electrostatic potential. To this end, we number the tail particles from 1 to , the bulk ones from to , and we call these respective sets of integers and . For , we now substitute Eq. (9) with

 ϕl(m,t)=ϕ(bal)bulk(m,t)+∑j∈Sbulk;j≠l δϕj(m,t)+Nbulk−1NbulkU(m,t), (34)

where

 U(m,t)=−eNbulkϵ0k2m(Nbulk−1)∑j∈Stailexp(−ikm⋅rj), (35)

and from now on, all quantities with subscript “” are those of previous sections computed with the bulk particles only. In the r.h.s. of Eq. (34), the third term vanishes if . We now perform the calculation of section III on substituting the previous summations with index running from 1 to by ones where the index runs over , while keeping the exclusion of where indicated. The previous division by preceding Eq. (14) is now a division by . This yields

 (36) = k2mϕ(bal)bulk(m,ω)+k2mU(m,ω).

On performing for the bulk particles the same approximations as in sections IV.1 and IV.2, Eq. (20) becomes

 ϵbulk(m,ω)Φ(m,ω)=ϕ(bal)bulk(m,ω)+U(m,ω). (37)

### vi.1 Amplitude equations

For the scales much larger than , the electric potential for the bulk is a superposition of Langmuir waves. The presence of tail particles slightly modifies these waves. We now derive an amplitude equation for the potential of the wave with wavenumber in a way similar to Refs OWM (); OLMSS (). On the l.h.s. of Eq. (37), the dielectric function defines a Bohm-Gross type dispersion relation associated with plasma oscillations of the bulk, and on the r.h.s. the tail particles provide with contributions oscillating in the Langmuir wave spectrum. Therefore, the frequency of interest in Eq. (37) for wavevector is close to the eigenfrequency solving corresponding to the wave propagating in the direction of  ; this frequency is real, since it is not resonant with the support of the bulk distribution function ( for all ’s such that ).

Let be the solution of Eq. (37) computed for , so that

 ϵbulk(m,ω)(Φ(m,ω)−Φbulk(m,ω))=U(m,ω). (38)

We are now looking for a which is close to the potential of the Langmuir wave , corresponding to , where is a constant equal to the sum for of

 Aj=−eϵ0k2mϵ′m exp[−ikm⋅rj(0)]ωm−km⋅˙rj(0). (39)

in Eq. (23) is . Let with . Therefore , which together with Eq. (38) yields

 Aϵbulk(m,ωm+ω′)[g(m,ω′)−iω′]=U(m,ωm+ω′), (40)

where . If , is a slowly evolving complex amplitude, and the dominant part of is concentrated near zero. This justifies Taylor-expanding about in Eq. (40), which yields to lowest order. Setting this into Eq. (40) and performing the inverse Laplace transform finally yields an amplitude equation for

 ∂Φ(m,t)∂t+iωmΦ(m,t)=ieNbulkϵ0k2m(Nbulk−1)∂ϵbulk∂ω(m,ωm)∑j∈Stailexp(−ikm⋅rj) (41)

similar to those in Refs OWM (); OLMSS (). The self-consistent dynamics of Langmuir waves and of the tail particles is ruled by Eq. (41) written for each wave and by the equation of motion of these particles due to the waves,

 ¨rj=ieL3me∑knΦ(n,t)exp(ikn⋅rj), (42)

where the summation runs over the indices of the 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 Refs MK (); AEE (); EEbook (). The study of this dynamics enables recovering Vlasovian linear theory with a mechanical understanding (see E013 (); SenFest () for a synthetic presentation). In particular, the reason why Landau damping cannot be a damped eigenmode is shown to be rooted deeply in Hamiltonian mechanics : a damped eigenmode must exist along with an unstable one, which is going to dominate with probability 1. Landau damping is recovered as an analogue of van Kampen phase-mixing effect. This phase-mixing in turn plays an essential role in the calculation of Landau instability in order to cancel the damped eigenmode (section 3.8.3 of Ref. EEbook ()). The self-consistent dynamics comes with an important bonus : it provides the information on particle dynamics in parallel with the wave’s. In particular, it reveals that both Landau damping and instability result from the same synchronization mechanism of particles with waves. This explains why there is a single formula for the rates of growth and damping, though growth involves an eigenmode and damping a phase-mixing instead EZE (); EEbook (); EFields (); Houches (). This synchronization mechanism was indeed evidenced experimentally DovEsMa (). In contrast, as yet the Vlasovian derivations of Landau damping do not provide the description of the corresponding individual evolution of particles, which forces textbooks to come up with complementary mechanical models. We point out that in Refs AEE (); EEbook () the equivalent of Eqs (41)-(42) was obtained without using any continuous approximation, but by a direct mechanical reduction of degrees of freedom starting with the -body problem.

### vi.2 Spontaneous emission

For the sake of brevity, we do not develop here the full generalization of the analysis in Refs AEE (); EEbook () ; it is lengthy, but straightforward. However, since this analysis unifies spontaneous emission with Landau growth and damping, we recall the result ruling the evolution of the amplitude of a Langmuir wave provided by perturbation calculations where the right hand sides of Eqs (41)-(42) are considered small (of first order). This is natural for Eq. (41) since , and for Eq. (42) if the Langmuir waves have a low amplitude. Let , where the average is over the random initial positions of the tail particles (their distribution being spatially uniform). Then a calculation to second order in yields

 dJ(m,t)dt=2γmLJ(m,t)+Smspont, (43)

where is the Landau growth or damping rate given by

 γmL=αmdfreddv(ωmkm;m) (44)

with

 αm=πe2meϵ0k2m∂ϵ∂ω(m,ωm), (45)

and is the reduced smoothed distribution function where is the unit vector along and is the component of the velocity perpendicular to  ; is given by

 Smspont=2α2mπe2kmnfred(ωmkm;m), (46)

where is the plasma density. 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 analogue of what was found in EZE (); EEbook ()). The second order calculation for the particles yields the diffusion and friction coefficients of the Fokker-Planck equation ruling the tail dynamics. This equation corresponds to the classical quasilinear result, plus a dynamical friction term mirroring the spontaneous emission of waves by particles, as found in the one-dimensional case in Refs EZE (); EEbook ().

An important aspect of the self-consistent dynamics defined by Eqs (41)-(42) is that it enables to use the modern tools of nonlinear dynamics and chaos available for finite dimensional systems. Let us consider two examples. First, the van Kampen phase-mixing effect leading to Landau damping is nowadays a classical result of Vlasovian theory. However, one may wonder whether nonlinear effects do not destroy these linear modes and the corresponding phase mixing. Proving the innocuity of nonlinear effects is the equivalent of deriving a Kolmogorov-Arnold-Moser (KAM) theorem for a continuous system (the Vlasov-Poisson one), a tour de force which partly earned C. Villani the 2010 Fields medal MV (); Villani (). The same result for the above finite-dimensional self-consistent dynamics requires the standard KAM theorem only : it is much simpler to keep the genuine granularity of the plasma.

Second, consider a tail distribution function which is a plateau in both velocity and space (this occurs for instance at the saturation of the bump-on-tail instability in a particle description of the plasma). Then the source term in Eq. (41) vanishes, as well as mode coupling, and the waves keep a fixed amplitude : the self-consistency of Eqs (41)-(42) is quenched, even when particle dynamics is strongly chaotic in the plateau domain (see sec. 2.2 of BEEB ()). Then, it is possible to use the tools of 1.5 degree-of-freedom Hamiltonian chaos to compute the diffusion of particle velocities. In particular, if chaos is strong enough, one may use a quasilinear diffusion coefficient (see section 2.2 of EE2 ()). 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 (the initial distribution is conserved along particle motion) ; a plateau in velocity exists for the spaced-averaged distribution function only, and a plateau in space exists for the velocity-averaged distribution only.

Third, 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 O'N (). The description of this phenomenon by the above self-consistent dynamics enables showing that this corresponds to a second order phase transition FirEl2 ().

## Vii Conclusion

This paper provides a direct path from microscopic mechanics to Debye shielding, Landau damping, and wave-particle interaction. This is performed by using Newton’s second law for the -body description, and standard tools of calculus. In particular, sections III to IV.3 provide the explicit, yet very compact derivation of formulas requiring about twenty pages (e.g. ch. 4, secs 6.3–5 and sec. 9.2 in Ref. Nicholson () and sec. 12.3 in Ref. BoydSan ()) in classical textbooks proceeding also explicitly from the -body description. This occurs thanks to a considerable simplification of the mathematical framework, in particular because no probabilistic argument and no partial differential equation are used.

This paper also solves the following issue : How can each particle be shielded by all other ones, while all the plasma particles are in uninterrupted motion ? This dynamical shielding turns out to be a mere consequence of the almost independent deflections of particles due to the Coulomb force. Therefore, shielding and collisional transport are two related aspects of the repulsive deflections of electrons.

The new approach also shows that each particle brings two simultaneous contributions to the plasma potential : a short-range one, its shielded Coulomb potential, and a long-range one, the excitation of Langmuir waves. However, these contributions need a time on the order of the plasma period to set in. Therefore, the plasma behaves like a dielectric only after this time during which Coulomb deflections bring their cooperative contribution to this self-organization. Paradoxically, Landau damping turns out to be correct because of what is usually called “collisions”.

The -body dynamics has always been the ultimate reference in plasma textbooks. Here it becomes a practical tool, which contributes to conform more strictly to the rules of inference from first principles : classical mechanics describes and unifies non-trivial aspects of the macroscopic dynamics of a many-body system, and brings more insight in plasma self-organization.

One might contemplate applying the above mechanical approach to plasmas with more species, or with a magnetic field, or where particles experience trapping and chaotic dynamics. The first generalization sounds rather trivial, and the third one is under way, at least in one dimension (see a pedestrian introduction in Houches () and more specific results in BEEB (); BEEBEPS ()).

J. Callen, Ph. Choquard, L. Couëdel, M.-C. Firpo, W. Horton, P.K. Kaw, J.T. Mendonça, F. Pegoraro, Y. Peysson, A. Piel, H. Schamel, D. Zarzoso, J.-Z. Zhu and an anonymous referee are thanked for very useful comments and references. The Program Committee of the 41st EPS Conference on Plasma Physics is thanked for providing one of us (DFE) with an extensive discussion of some aspects of the present paper.

## Appendix A Physical meaning of the linearization of the Coulomb potential (Approximation 1)

Approximation 1 in section III corresponds to substituting the Coulomb potential in Eq. (3) with as defined by Eq. (9). Using the definition of the ballistic approximation and Eqs (7)-(8),

 ϕl(m,t) = −eϵ0k2m∑j∈S;j≠lexp[−ikm⋅r(0)j(t)](e−ikm⋅(ξj(0)+˙ξj(0)t)−ikm⋅Δξj). (47)

Taking into account that the ’s are small, we may write

 ϕl(m,t)≃∑j∈S;j≠lδψj(m,t), (48)

with

 δψj(m,t)=−eϵ0k2mexp(−ikm⋅r(0)j(t))(1−ikm⋅ξj(t)). (49)

The inverse Fourier transform of verifies

 limL→∞δψj(r,t)=−e4πϵ0∥r−r(0)j∥−eξj⋅(r−r(0)j)4πϵ0∥r−r(0)j∥3. (50)

The -th contribution to the approximate electric field turns out to be due to a particle located at instead of , and is made up of a Coulombian part and of a dipolar part with dipole moment . The cross-over between these two parts occurs for on the order of , i.e. when the distance to the ballistic particle is about the distance between the latter and the true particle . For larger values of , the dipolar component is subdominant. For smaller ones, it is dominant, but with a direction which is a priori uncorrelated with the Coulombian one ( is almost independent from ). Since the ’s are assumed small, for particles having an almost uniform spatial distribution, this dipolar part should be exceptionally dominant, as it corresponds to a very close encounter between particles. As a result, the approximate electric field stays dominantly of Coulombian nature, but with a small mismatch of the charge positions with respect to the actual ones.

## Appendix B Vanishing non-diagonal contribution

We consider a class of granular distributions close to being spatially uniform. 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 . In the following, we are interested in a sequence of such sets of beams converging for large toward a spatially uniform Maxwellian distribution. We now discretize this distribution into beams having the same number of particles. The edge length for each beam is with an integer. Therefore the number of particles in each beam is .

In order to discretize the beam velocities, we first define in the space of velocities a sequence of shells with increasing finite radii centered on the origin and containing each particles ; the spacing between these radii is not constant, and increases like as . This means the far tail of the Maxwellian is not described, but the range of described velocities increases with . We discretize the directions of the velocities into adjacent cones with a solid angle . These cones are chosen such that they fit into circular cones with a solid angle scaling like (a particular instance of the bases of such cones can be obtained by tesselating the sphere as follows : choose a polar axis on the sphere, cut the latter with equidistant meridians, and introduce a family of parallels defining equal annular areas on the sphere). Therefore, the intersections of these cones with the shells define velocity domains with each particles, whose extension in all directions vanishes when and go to infinity. Each domain is related to one beam whose velocity is the average value of the velocity in the domain. Then . The continuous distribution function is recovered in the limit where , , and all three go to infinity. In order to describe particle distributions close to being spatially uniform, we consider that the ’s and ’s are slightly shifted with respect to the ’s and ’s.

We now focus on a given beam. It is convenient to consider the index of its particles as a vector whose integer components run each from 1 to . The particles have the same value of . Therefore, and take on a single value each in Eq. (14), and the summation over the particles of the beam of interest bears on only, where . Due to the periodicity of the ’s, the corresponding sum vanishes unless the three components of are on the simple cubic lattice with mesh length . For the range of ’s considered in section IV.5, if