Basic microscopic plasma physics from N-body mechanics

Abstract

Computing is not understanding. This is exemplified by the multiple and discordant interpretations of Landau damping still present after 70 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.

Appendices

Appendix 1: Rigorous derivation à la Kaufman

A(t) is entire

The dynamics considered in Sect. 3 corresponds to the linearized motion of the N electrons with respect to a given multi-beam-multi-array, when a single wave with wavevector $\mathbf {k}_{\mathbf {m}}$ is excited at $t=0$. Let $\mathbf {r}_{j0}$ be the initial position of the unperturbed beam particle with index j, and $\mathbf {v}_{j}$ be its velocity, and let $\varDelta {\mathbf {r}}_j (t)= \mathbf {r}_j (t) - \mathbf {r}_{j0} - \mathbf {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 $\mathbf {r}_j = \mathbf {r}_{j0} + \mathbf {v}_j t + \varDelta {\mathbf {r}}_j (t)$ in Eq. (8), we replace $\tilde{\varphi }$ with its expansion to first order in the $\varDelta {\mathbf {r}}_j (t)$’s

$$\begin{aligned} \tilde{\varphi } (\mathbf {m},t) = - \sum _{l = 1}^N \frac{ e}{\varepsilon _0 k_{\mathbf {m}}^2} \exp [- {\mathrm {i}}\mathbf {k}_{\mathbf {m}} \cdot (\mathbf {r}_{l0} + \mathbf {v}_{l} t)] \ [1 - {\mathrm {i}}\mathbf {k}_{\mathbf {m}} \cdot \varDelta {\mathbf {r}}_l(t)]. \end{aligned}$$

(101)

Using Eqs. (8 and 101), the linearized particle dynamics defined by Eq. (9) is then given by

$$\begin{aligned} \varDelta \ddot{\mathbf {r}}_j = \frac{{\mathrm {i}}e}{L^3 m} \mathbf {k}_{\mathbf {m}} \ \tilde{\varphi }(\mathbf {m},t) \exp [{\mathrm {i}}\mathbf {k}_{\mathbf {m}} \cdot (\mathbf {r}_{j0} + \mathbf {v}_{j} t)] + {\mathrm {c.c.}}\,. \end{aligned}$$

(102)

Because of ordering (10), for each beam the corresponding values of $\exp [{\mathrm {i}}\mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{l0}]$ are uniformly distributed on the unit circle, and their global contribution to the 1 factor in the last bracket of Eq. (101) vanishes. Therefore Eq. (102) yields the compact expression

$$\begin{aligned} \varDelta \ddot{\mathbf {r}}_j = - \frac{2 \omega _{{\mathrm {p}}}^2}{N k_{\mathbf {m}}^2} \mathbf {k}_{\mathbf {m}} \ \sum _{l = 1}^N \cos [\mathbf {k}_{\mathbf {m}} \cdot [\mathbf {r}_{j0} - \mathbf {r}_{l0} + (\mathbf {v}_{j} - \mathbf {v}_{l}) \, t \, ] \, \mathbf {k}_{\mathbf {m}} \cdot \varDelta {\mathbf {r}}_l(t). \end{aligned}$$

(103)

This defines a system of N linear differential equations whose coefficients are entire functions of t. Therefore, the $\varDelta {\mathbf {r}}_l(t)$’s are entire functions. Through Eq. (101), this property is transferred to $\tilde{\varphi } (\mathbf {m},t)$ and to the amplitude $A(t) = \tilde{\varphi } (\mathbf {m},t) \exp ({\mathrm {i}}\omega t)$ for $\omega$ real: A(t) is an entire function.

If the $\mathbf {k}_{\mathbf {m}} \cdot \mathbf {v}_{j}$’s are multiples of a given number, Eq. (103) has coefficients with some period T, and belongs to the Floquet class of differential equations. Then its solutions are of the type

$$\begin{aligned} U(t) = V(t) {\mathrm {e}}^{\alpha t} , \end{aligned}$$

(104)

where V(t) is a vector of period T, and $\alpha$ a complex number. The corresponding Floquet exponents are the $\beta _\sigma$’s introduced at the end of Sect. 5.3. A similar equation was met in the self-consistent wave–particle approach introduced in Sect. 9.1 (section 3 of Elskens and Escande 2002). 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.

Solution to all orders

By expanding $A(t-\tau )$ in Taylor series, the second term of Eq. (16) becomes (with $k = k_{\mathbf {m}}$ and $v_\phi = \omega /k)$)

$$\begin{aligned} S= & {} - {\mathrm {i}}\omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \int \frac{\partial }{\partial \varOmega } \int _0^t (-1)^n\frac{\tau ^n}{n!} \exp ( - {\mathrm {i}}\varOmega \tau ) \, {\mathrm {d}}\tau \, g(v) {\mathrm {d}}v \nonumber \\= & {} - {\mathrm {i}}\omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \int \frac{\partial }{\partial \varOmega } \int _0^t \frac{{\mathrm {i}}^{3n}}{n!} \frac{\partial ^n}{\partial \varOmega ^n} \exp ( - {\mathrm {i}}\varOmega \tau ) \, {\mathrm {d}}\tau \, g(v) {\mathrm {d}}v \nonumber \\= & {} \omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {i}}^{3n}}{n!} \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \int \frac{\partial ^{n+1}}{\partial \varOmega ^{n+1}} \frac{\exp ( - {\mathrm {i}}\varOmega t)-1}{\varOmega } \, g(v) {\mathrm {d}}v \nonumber \\= & {} \omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {i}}^{3n}}{n! k^{n+1}} \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \int (-1)^{n+1} g^{(n+1)}(v) \, \frac{\exp ( - {\mathrm {i}}\varOmega t)-1}{\varOmega } {\mathrm {d}}v \nonumber \\= & {} - \omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {i}}^{n}}{n! k^{n+1}} \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \int g^{(n+1)}(v) \frac{\exp ( - {\mathrm {i}}\varOmega t)-1}{\varOmega } {\mathrm {d}}v \nonumber \\&\quad \sim _{t\rightarrow +\infty } \omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty } \frac{{\mathrm {i}}^{n}}{n! k^{n+1}} \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \left[ {\mathrm{P}} \! \! \! \int \frac{g^{(n+1)}(v)}{\varOmega } {\mathrm {d}}v + {\mathrm {i}}\frac{\pi }{k} g^{(n+1)}(v_\phi ) \right] . \end{aligned}$$

(105)

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

$$\begin{aligned} 1 = \frac{\omega _{\mathrm {p}}^2}{A} \sum _{n=0}^{+\infty } \frac{{\mathrm {i}}^{n}}{n! k^{n+1}} \frac{{\mathrm {d}}^n A(t)}{{\mathrm {d}}t^n} \left[ {\mathrm{P}} \! \! \! \int \frac{g^{(n+1)}(v)}{\varOmega } {\mathrm {d}}v + {\mathrm {i}}\frac{\pi }{k} g^{(n+1)}(v_\phi ) \right] . \end{aligned}$$

(106)

Looking for a solution $A(t)=A(0)e^{\nu t}$ with $\nu$ real, Eq. (106) reads

$$\begin{aligned} 1 = \omega _{\mathrm {p}}^2 \sum _{n=0}^{+\infty }\frac{{\mathrm {i}}^{n}}{n! k^{n+1}} \nu ^n \left[ {\mathrm{P}} \! \! \! \int \frac{g^{(n+1)}(v)}{\varOmega } {\mathrm {d}}v + {\mathrm {i}}\frac{\pi }{k} g^{(n+1)}(v_\phi ) \right] , \end{aligned}$$

(107)

whose real part is

$$\begin{aligned} 1 = \omega _{\mathrm {p}}^2 \sum _{p=0}^{+\infty } \frac{(-1)^p}{(2p)! k^{2p+1}} \nu ^{2p} \left[ {\mathrm{P}} \! \! \! \int \frac{g^{(2p+1)}(v)}{\varOmega } {\mathrm {d}}v - \frac{\pi \nu }{(2p+1)k^2} g^{(2p+2)}(v_\phi )\right] , \end{aligned}$$

(108)

and imaginary part is

$$\begin{aligned} 0 = \sum _{p=0}^{+\infty } \frac{(-1)^p\nu ^{2p}}{(2p+1)!k^{2p}} \left[ \nu {\mathrm{P}} \! \! \! \int \frac{g^{(2p+2)}(v)}{\varOmega } {\mathrm {d}}v +\pi (2p+1)g^{(2p+1)}(v_\phi ) \right] . \end{aligned}$$

(109)

Since $g^{(n+1)}(v) \sim v_{{\mathrm {T}}}^{-n}g'(v)$, Eqs. (108) and (109) may be seen as expansions in the small parameter $\varepsilon \equiv \nu / (k \lambda _{\mathrm {D}}\omega _{{\mathrm {p}}})$. Then, these equations may be solved order by order. The first two orders were used in Sect. 3.1.

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

$$\begin{aligned} 1 = \frac{\omega _{\mathrm {p}}^2}{k} \left[ {\mathrm{P}} \! \! \! \int \frac{g'(v + {\mathrm {i}}\nu / k)}{\varOmega } {\mathrm {d}}v + {\mathrm {i}}\frac{\pi }{k} g'(v_\phi +i\nu /k) \right] , \end{aligned}$$

(110)

which is exactly the formula obtained by Landau.

Appendix 2: Relevance of transients before Landau damping

In this Appendix, we discuss the relevance of transients that would occur before the self-consistent 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–18). Such a generalization has already been performed in Bénisti and Gremillet (2007) when the self-consistent field was slowly driven. An important point of Bénisti and Gremillet (2007) was to prove that a wave may be considered as slowly varying, provided that its complex amplitude did not change much during a time interval of the order of $\tau _{\mathrm{mix}}=(k v_{\mathrm {T}})^{-1}$. More precisely, the typical wave growth rate, $\gamma$, had to be such that $\vert \gamma \vert \tau _{\mathrm{mix}}<0.1$.

We now further discuss the importance of the product $\vert \gamma \vert \tau _{\mathrm{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<\tau _{\mathrm {d}}$ (see Bénisti and Gremillet 2007 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 $\varphi =A(t) \exp [{\mathrm {i}}(\mathbf {k_m}.\mathbf {r}-\omega t)] + {\mathrm {c.c.}}$ (which has been proved to be correct in Bénisti and Gremillet 2007 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, $\mathbf {r}_j=\mathbf {r}_{j0}$ when $t=0$. Then, Eq. (14) is changed into

$$\begin{aligned} \varDelta \mathbf {r}_{j1}(t) = \alpha \mathbf {k}_{\mathbf {m}} \int _{0}^{t} \tau A(t-\tau ) \exp [{\mathrm {i}}(\varOmega _j (t-\tau ) + \mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{j0})] {\mathrm {d}}\tau + {\mathrm {c.c.}} \end{aligned}$$

(111)

If the wave is slowly driven, A(t) is a slowly varying function, and one may stop the Taylor expansion of $A(t-\tau )$ at first order, which yields

$$\begin{aligned} A(t)= & {} \int \frac{\omega _{{\mathrm {p}}}^2}{k_{\mathbf {m}}} A(t) g'(v) \frac{1-\cos (\varOmega t)+ {\mathrm {i}}\sin (\varOmega t)}{\varOmega } {\mathrm {d}}v \nonumber \\&- \, {\mathrm {i}}\frac{\omega _{{\mathrm {p}}}^2}{k_{\mathbf {m}}^2} \dot{A}(t) \int g^{\prime \prime }(v) \frac{\cos (\varOmega t)- {\mathrm {i}}\sin (\varOmega t) -1}{\varOmega } {\mathrm {d}}v, \end{aligned}$$

(112)

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 $\delta \mathbf {r}_{ j}=\mathbf {0}$. Hence, unlike in Eq. (18), no transient is expected from this term. Now, when $t>\tau _{\mathrm {d}}$ it is clear that A(t) is nothing but the amplitude of the self-consistent potential, since the electrons no longer feel the effect of the drive. Moreover, because it is slowly driven, only when $\tau _{\mathrm {d}}\gg \tau _{\mathrm{mix}}$ may the self-consistent field reach a significant amplitude, and may effectively be Landau damped. Hence, when $t>\tau _{\mathrm {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 slowly varying when $t < \tau _{\mathrm {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\le \tau _{\mathrm {d}})$. To do so, we still assume that the electrons are at equilibrium when $t=0$. Moreover, we use Eq. (111), which is exact, to calculate the shift in their positions at $t=\tau _{\mathrm {d}}$, induced by the drive and the self-consistent electric fields. This yields

$$\begin{aligned} \varDelta \mathbf {r}_{j1}(t)= & {} \alpha \mathbf {k}_{\mathbf {m}} \int _{0}^{\tau _{\mathrm {d}}} \tau A(\tau _{\mathrm {d}}-\tau ) \exp [{\mathrm {i}}(\varOmega _j (\tau _{\mathrm {d}}-\tau ) + \mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{j0})] {\mathrm {d}}\tau + {\mathrm {c.c.}}\nonumber \\\equiv & {} \delta \mathbf {r}_{ j} \sin (\mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{j0}+\psi _0). \end{aligned}$$

(113)

As noted in Sect. 3.1, the phase $\psi _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>\tau _{\mathrm {d}}$, which will let us conclude about the evolution of the wave amplitude. Using again Eq. (111), we find

$$\begin{aligned} \varDelta \mathbf {r}_{j_1}(t)= & {} \delta \mathbf {r}_{ j} \sin (\mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{j0}) + \left\{ \alpha \mathbf {k_m}\int _{\tau _{\mathrm {d}}}^t \tau \mathcal {A}(t-\tau ) {\mathrm {d}}\tau \right. \nonumber \\&+ \left. \alpha \mathbf {k_m} \int _0^{\tau _{\mathrm {d}}} \tau \left[ \mathcal {A}(t-\tau )-\mathcal {A}(\tau _{\mathrm {d}}-\tau ) \right] {\mathrm {d}}\tau + {\mathrm {c.c.}}\right\} \ , \end{aligned}$$

(114)

where we introduced

$$\begin{aligned} \mathcal {A}(t) \equiv A(t)\exp [{\mathrm {i}}(\varOmega _j t + \mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_{j0})]. \end{aligned}$$

(115)

In Sect. 3.1, the limit $\tau _{\mathrm {d}}\rightarrow 0$ was considered, which yields Eq. (14). In this Appendix, we specify how small $\tau _{\mathrm {d}}$ has to be for the results of Sect. 3.1 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\gg \tau _{\mathrm {d}}$ and if A(t) does not abruptly vanish within a time interval smaller than $\tau _{\mathrm {d}}$. The first condition implies that the results of Sect. 3.1 are only valid when $t\gg \tau _{\mathrm {d}}$. The second condition is only true provided that $\gamma _\mathrm{L} \tau _{\mathrm {d}}\ll 1$.

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

$$\begin{aligned} A(t) =&\int \left[ \frac{Ne}{2\varepsilon _0 k_{\mathbf {m}}^2 }\, h(v) \exp (- {\mathrm {i}}\varOmega t) \right. \nonumber \\&\left. + \frac{\omega _{{\mathrm {p}}}^2}{k_{\mathbf {m}}} A(t) g'(v) \frac{{\mathrm {e}}^{-{\mathrm {i}}\varOmega \tau _{\mathrm {d}}}-\cos (\varOmega t)+ {\mathrm {i}}\sin (\varOmega t)}{\varOmega } \right] {\mathrm {d}}v \nonumber \\&- \, {\mathrm {i}}\frac{\omega _{{\mathrm {p}}}^2}{k_{\mathbf {m}}^2} \dot{A}(t) \int g^{\prime \prime }(v) \frac{\cos (\varOmega t)- {\mathrm {i}}\sin (\varOmega t) - {\mathrm {e}}^{-{\mathrm {i}}\varOmega \tau _{\mathrm {d}}}}{\varOmega } {\mathrm {d}}v. \end{aligned}$$

(116)

Then, Eq. (18) is recovered only in the limit $\varOmega \tau _{\mathrm {d}}\rightarrow 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 $\tau _{\mathrm{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 $\tau _{\mathrm{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, $\tau _{\mathrm {d}}$, was much smaller than $\tau _{\mathrm{mix}}$ and the plasma period, or when it was much smaller than $\tau _{\mathrm{mix}}$. The situation when $\tau _{\mathrm {d}}$ is either of the order of the plasma period or of $\tau _{\mathrm{mix}}$ would deserve further investigation.

Appendix 3: Infinite number of beams

This appendix provides calculations similar to those in Dawson (1960), but reformulated in a form suitable for the derivation of Sect. 5, and with the correction of several errors. We first focus on $\epsilon _{{\mathrm {d}}1}(\mathbf {m},\omega )$ defined by Eq. (51). To separate the regular and the singular parts of this quantity, we add and subtract to the right-hand side of Eq. (51) the quantity

$$\begin{aligned} \omega _{{\mathrm {p}}}^2\left[ \frac{\pi ^2 g(\omega /k)}{k^2 \delta \sin ^2(\pi \omega /k \delta )}-\frac{2 \pi g'(\omega /k)}{k^2}\cot (\pi \omega /k \delta )\right] , \end{aligned}$$

(117)

where $k=k_{\mathbf {m}}$. This yields

$$\begin{aligned} \epsilon _{{\mathrm {d}}1}(\mathbf {m},\omega )=&1 - \omega _{{\mathrm {p}}}^2 \left[ \frac{\pi ^2 g(\omega /k)}{k^2 \delta \sin ^2(\pi \omega /k \delta )}-\frac{2 \pi g'(\omega /k)}{k^2} \cot (\pi \omega /k \delta ) \right. \nonumber \\&\left. + \varSigma _{\sigma = - \infty }^{\infty }\left( \frac{[g(\sigma \delta ) -g(\omega /k)]\delta }{ (\omega - \sigma k \delta )^2}+\frac{2 g'(\omega /k) \omega \delta }{k (\omega ^2 - (\sigma k \delta )^2)} \right) \right] , \end{aligned}$$

(118)

where use has been made of the relations (see Abramowitz and Stegun 1967)⁴⁴

$$\begin{aligned} \frac{\pi ^2 }{\sin ^2(\pi x)}= & {} \varSigma _{\sigma = - \infty }^{\infty }\frac{1}{ (x - \sigma )^2} \end{aligned}$$

(119)

$$\begin{aligned} \pi \cot (\pi x)= & {} \varSigma _{\sigma = - \infty }^{\infty }\frac{x}{ (x^2 - \sigma ^2)}. \end{aligned}$$

(120)

We notice that the function of $\sigma \delta$ inside the summation in Eq. (118) has no real poles. Therefore, the sum passes smoothly to an integral as $\delta$ goes to zero⁴⁵. In this limit, Eq. (118) becomes

$$\begin{aligned} \epsilon _{{\mathrm {d}}1}(\mathbf {m},\omega )=&1 - \omega _{{\mathrm {p}}}^2 \left[ \frac{\pi ^2 g(\omega /k)}{k^2 \delta \sin ^2(\pi \omega /k \delta )}-\frac{2 \pi g'(\omega /k)}{k^2} \cot (\pi \omega /k \delta ) \right. \nonumber \\&\qquad \left. + \int _{- \infty }^{\infty }\left( \frac{g'(v)}{ k(\omega - k v)}+\frac{2 g'(\omega /k) \omega }{k (\omega ^2 - k^2 v^2)} \right) {\mathrm {d}}v \right] . \end{aligned}$$

(121)

We now compute the zeros of $\epsilon _{{\mathrm {d}}1}(\mathbf {m},\omega )$ and write $\omega = \alpha + {\mathrm {i}}\beta$. We first consider those with $\beta$ vanishing when $\delta$ goes to zero. If $\beta$ vanished like or faster than $\delta$, the first term in the bracket of Eq. (121) would diverge, while the second and third one would remain finite, which is impossible. Therefore, $\beta$ vanishes slower than $\delta$, which forces the cotangent to converge toward $- \mu \, {\mathrm {i}}$, where $\mu = \pm 1$ is the sign of $\beta$. Then, to stay finite, the first term requires $\beta$ to scale like $\delta \, |\ln (\delta /v_{{\mathrm {T}}})\,|$. With this in mind, and looking for solutions in the vicinity of $n k \delta$, Eq. (121) requires

$$\begin{aligned}&\omega _{{\mathrm {p}}}^2 \left[ - \frac{4 \pi ^2 g(n \delta ) \exp [2 \pi ({\mathrm {i}}\mu \alpha _1 - |\beta |)/k \delta ]}{k^2 \delta }+\frac{{\mathrm {i}}\mu \pi g'(n \delta )}{k^2}\right] \nonumber \\&\quad = 1 + {\mathrm{P}} \! \! \! \int _{- \infty }^{\infty } \frac{\omega _{{\mathrm {p}}}^2 g'(v)}{ k^2(n \delta - v)} \, {\mathrm {d}}v, \end{aligned}$$

(122)

for $n \delta$ in the support of g, with $\alpha _1 = \alpha - n k \delta$. Equation (121) provides two contributions to the term in $g'(n \delta )$: 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'(\omega /k)$ bring contributions cancelling each other. Solving for $\alpha _1$ and $\beta$ yields

$$\begin{aligned}&\tan \frac{2 \pi \alpha _1}{k \delta } = - \frac{\pi \omega _{{\mathrm {p}}}^2}{k^2} g'(n \delta )/\left[ 1 + {\mathrm{P}} \! \! \! \int _{- \infty }^{\infty } {\mathrm {d}}v \frac{\omega _{{\mathrm {p}}}^2 g'(v)}{ k^2(n \delta - v)}\right] \, , \end{aligned}$$

(123)

$$\begin{aligned}& \beta = \mu \frac{k \delta }{2 \pi } \ln \left\{ \left[ \frac{k^2 \delta }{4 \pi ^2 \omega _{{\mathrm {p}}}^2 g(n \delta )}\right] \right. \times \nonumber \\& \qquad \qquad \qquad \left. \left[ \left( 1 \! \! + {\mathrm {P}}\! \! \! \int _{- \infty }^{\infty } \! \! \! {\mathrm {d}}v \frac{\omega _{{\mathrm {p}}}^2 g'(v)}{ k^2(n \delta - v)} \right) ^2 \! \! + \! \left( \frac{\pi \omega _{{\mathrm {p}}}^2 g'(n \delta )}{k^2} \right) ^2 \right] \right\}. \end{aligned}$$

(124)

Equation (123) yields $2 \pi \alpha _1/k \delta$ modulo $\pi$, and the right solution is obtained by requiring $\cos (2 \pi \mu \alpha _1/k \delta )$ to have the opposite sign to the denominator of this equation.

Like the natural frequencies of the beams, the roots are spaced $k \delta$ apart in $\alpha$. Therefore, the above zeros have real parts between these natural frequencies. There are two roots for each beam, since $\mu$ 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 $\frac{\partial \epsilon }{\partial \omega } (\mathbf {m}, \omega )$ defined in Eq. (54) for the case of a vanishing imaginary part of $\omega$ when $\delta$ goes to zero. Here again, we handle the singularity by adding and subtracting to the right-hand side of Eq. (54) the quantity

$$\begin{aligned} 2 \omega _{{\mathrm {p}}}^2\left[ - \frac{\pi ^3 \cos (\pi \omega /k \delta ) g(\omega /k)}{k^3 \delta ^2 \sin ^3(\pi \omega /k \delta )} - \frac{\pi ^2 g'(\omega /k)}{k^2 \delta \sin ^2(\pi \omega /k \delta )} + \frac{\pi g^{\prime \prime }(\omega /k)}{2 k^3}\cot (\pi \omega /k \delta )\right] . \end{aligned}$$

(125)

These terms can also be written in the form of sums using again Eqs. (119 and 120), and (see Abramowitz and Stegun 1967)

$$\begin{aligned} \pi ^3 \frac{\cos (\pi x) }{\sin ^3(\pi x)}= & {} \varSigma _{\sigma = - \infty }^{\infty }\frac{1}{ (x - \sigma )^3}. \end{aligned}$$

(126)

Using the latter expression, for $\delta$ small we find that $\epsilon '_{\sigma ,\mu }$ is given by Eq. (56).

Appendix 4: 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 $\varDelta \mathbf {v}$ distributed with the continuous distribution $P(\varDelta \mathbf {v})$, instead of a discrete distribution in the case of the pulverization⁴⁶. 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 Eq. (47) can be performed again, but the summation over particles now involves an integral over the nearby velocities of the coarse-grained system. Then $\varphi ^{({\mathrm{bal}})}(\mathbf {m},\omega )$ is substituted with

$$\begin{aligned} \varphi _{{\mathrm{cg}},j}^{({\mathrm{bal}})}(\mathbf {m},\omega ) = - \frac{{\mathrm {i}}e}{\varepsilon _0 k_{\mathbf {m}}^2} \int \frac{\exp [- {\mathrm {i}}\mathbf {k}_{\mathbf {m}} \cdot \mathbf {r}_j(0)]}{\omega -\mathbf {k}_{\mathbf {m}} \cdot ({\dot{\mathbf {r}}}_j(0) + \mathbf {u})} P(\mathbf {u}) {\mathrm {d}}^3 \mathbf {u}, \end{aligned}$$

(127)

and

$$\begin{aligned} f_0(\mathbf {v}) = \sum _{\sigma = 1}^{n_{\mathrm {b}}} N_\sigma P(\mathbf {v}-\mathbf {w}_\sigma ) , \end{aligned}$$

(128)

where $n_{\mathrm {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(\mathbf {u})$ may be taken as a Dirac distribution in Eq. (127), which then becomes Eq. (46).

Appendix 5: Shielded Coulomb potential by a singular limit of the many-beam description

If in Eq. (65) we take first the limit $\delta \rightarrow 0$ and then the limit $\varepsilon \rightarrow 0^{+}$, according to Eq. (121), $\epsilon _{\mathrm {d}}(\mathbf {m},\mathbf {k}_{\mathbf {m}} \cdot \mathbf {v}+ {\mathrm {i}}\varepsilon )$ converges toward

$$\begin{aligned} \epsilon _{\mathrm{lim}} (\mathbf {m},\mathbf {k}_{\mathbf {m}} \cdot \mathbf {u}) = 1 + \omega _{{\mathrm {p}}}^2 \, {\mathrm{P}} \! \! \! \int \frac{\frac{\partial g}{\partial v}(v)}{k_{\mathbf {m}}(\mathbf {k}_{\mathbf {m}} \cdot \mathbf {u}- k_{\mathbf {m}} v)} {\mathrm {d}}v - {\mathrm {i}}\mu \frac{\pi \omega _{{\mathrm {p}}}^2}{k_{\mathbf {m}}^2} \frac{\partial g}{\partial v}\left( \frac{\mathbf {k}_{\mathbf {m}} \cdot \mathbf {u}}{k_{\mathbf {m}}}\right) . \end{aligned}$$

(129)

This is nothing but the contribution of $\epsilon (\mathbf {m},\mathbf {k}_{\mathbf {m}} \cdot \mathbf {u}+ {\mathrm {i}}\varepsilon )$ in the same limit (see for instance equation (9.12) of Nicholson 1983). Therefore, for $\delta$ small enough, Eq. (64) becomes Eq. (68), and Eq. (65) becomes Eq. (69).