1 Introduction
Collisionless damping of electrostatic waves is one of the most fundamental phenomena in plasma physics and a starting point for the studies on the wave–particle interaction. From its first report by Landau (Reference Landau1946) until its experimental verification by Malmberg & Wharton (Reference Malmberg and Wharton1964), this phenomenon was considered with some caution by the plasma community, in particular due to the mathematical aspects (specifically the use of complex integral and the large time limit in the Laplace transform) of its derivation and its paradoxical nature.
A very physical and intuitive approach was proposed by Dawson (Reference Dawson1961), who explained this effect as resulting from a near-resonance mechanism of energy (and momentum) transfer between wave and particles. Dawson also pointed out a time threshold for a breakdown (due to trapping effects) of Landau’s linear analysis. Later on, O’Neil (Reference O’Neil1965) extended the analysis for arbitrary times and found that the nonlinearities may lead the wave to a time-asymptotic behaviour with a constant non-zero amplitude. O’Neil’s picture, introduced in the sixties, gave rise to what we call today ‘nonlinear Landau damping’ and induced a reconsideration of the phenomenon from experimental (Malmberg & Wharton Reference Malmberg and Wharton1967; Franklin, Hamberger & Smith Reference Franklin, Hamberger and Smith1972), numerical (Brodin Reference Brodin1997; Manfredi Reference Manfredi1997) and theoretical (Mouhot & Villani Reference Mouhot and Villani2010, Reference Mouhot and Villani2011) viewpoints. In the last decades, the existence of a critical initial amplitude that distinguishes whether the asymptotic field Landau damps to zero or evolves to a steady non-zero value has been proved, by solving the Vlasov–Poisson system, by Brunetti, Califano & Pegoraro (Reference Brunetti, Califano and Pegoraro2000) and Lancellotti & Dorning (Reference Lancellotti and Dorning2003) and from a statistical-physics-like viewpoint by Firpo & Elskens (Reference Firpo and Elskens2000, Reference Firpo and Elskens2001).
Another prominent interpretation for linear Landau damping (beside Dawson’s mechanical picture) explains it through the phase mixing of generalized (singular) eigenmodes of the corresponding linear dynamical system. This formulation, resting on a strictly linear theory and originally proposed by van Kampen (Reference van Kampen1955), shows that for a given perturbation with wavenumber
$k$
there is a continuous spectrum of real frequencies allowed for the plasma oscillations. In this context, Landau damping emerges from a special superposition of stationary (distribution-like) solutions that now go under the name of van Kampen modes. Though such superpositions account for the exponential decay of initial perturbations at Landau’s rate, special initial data can also be constructed which cause a decay at slower rates (Belmont et al.
Reference Belmont, Chust, Mottez and Hess2011).
Nowadays, it is well known that other kinds of systems also admit normal modes analogous to the van Kampen modes for plasmas (see e.g. Pen Reference Pen1994; Vekstein Reference Vekstein1998; Vandervoort Reference Vandervoort2003). Indeed, similar phenomena are found in a wide variety of systems such as e.g. dusty plasmas (Bliokh, Sinitsin & Yaroshenko Reference Bliokh, Sinitsin and Yaroshenko1995), galaxies (Kandrup Reference Kandrup1998), two-dimensional fluid flows (del Castillo-Negrete & Firpo Reference del Castillo-Negrete and Firpo2002) and liquids containing gas bubbles (Smereka Reference Smereka2002). Despite its being an old problem, Landau damping remains intensely investigated experimentally (Doveil, Escande & Macor Reference Doveil, Escande and Macor2005) and theoretically (see e.g. Dougherty Reference Dougherty1998; Ryutov Reference Ryutov1999; del Castillo-Negrete Reference del Castillo-Negrete, Dauxois, Ruffo, Arimondo and Wilkens2002; Elskens Reference Elskens, Passot, Sulem and Sulem2005; Mouhot & Villani Reference Mouhot and Villani2010, Reference Mouhot and Villani2011; Bratanov et al. Reference Bratanov, Jenko, Hatch and Brunner2013; Bénisti Reference Bénisti2016) with new aspects and viewpoints still being explored.
In contrast to the difficult reception of Landau damping by the physics community, the weak warm beam instability, captured by the same formula (with an opposite sign reflecting its physics), was immediately accepted. We show in this paper that its status is somewhat more subtle, and that this Landau instability has more to share with damping than might be thought initially. Indeed, the aim of this paper is to provide further theoretical and numerical insight into how both Landau effects emerge from the phase mixing among the van Kampen-like modes in the Hamiltonian approach.
The key element in plasmas is their being an
$N$
-body systemFootnote
1
. Therefore, a microscopic description involving actual particles, rather than a continuum represented by a smooth distribution function in
$(\boldsymbol{r},\boldsymbol{v})$
space, is their physical fundamental model (Elskens & Escande Reference Elskens and Escande2003; Elskens, Escande & Doveil Reference Elskens, Escande and Doveil2014). As deriving a kinetic model for singular interactions, with the Coulomb divergence at short range, remains a challenge (Chaffi, Casta & Brenig Reference Chaffi, Casta and Brenig2014; Kiessling Reference Kiessling2014; Brenig, Chaffi & Rocha Filho Reference Brenig, Chaffi and Rocha Filho2016), understanding basic plasma phenomena from the
$N$
-body picture is important. Moreover, the
$N$
-body approach yields its own benefits, such as showing unexpected connections between Landau damping and Debye screening (Escande, Doveil & Elskens Reference Escande, Doveil and Elskens2016).
The formalism adopted here is closely related to the fluid model used by Dawson (Reference Dawson1960), where electrons are distributed into (almost) monokinetic beams, and the central issue consists in analysing normal modes that fully describe the evolution of the system. This approach is mathematically elementary, involving no partial differential equations or functional analysis and not even an analytic continuation or
$\text{i}\unicode[STIX]{x1D700}$
prescription enforcing a causality argument. Its concepts belong in basic classical mechanics.
In § 2, we revisit the Hamiltonian model for plasmas and waves (Escande, Zekri & Elskens Reference Escande, Zekri and Elskens1996; Elskens & Escande Reference Elskens and Escande2003) and discuss connections with the van Kampen–Case approach (van Kampen Reference van Kampen1955; Case Reference Case1959) when the continuous limit is taken. In § 3, we obtain the spectrum of the discrete analogue to the van Kampen frequencies by solving accurately the dispersion relation for systems composed of up to 2000 beams. We also investigate differences between the damping and growth regimes and the consistency between the discrete and the continuous systems by monitoring the evolution of the wave intensity.
Our main results show that not only the damping but also the Landau instability emerge as consequences of a phase mixing mechanism among the eigenmodes of the linear system. We highlight that while in the stable case the phase mixing to generate Landau damping involves all van Kampen-like eigenmodes, in the unstable case the pure Landau growth results from a destructive interference effect, leaving a single eigenmode having a dominant (exclusive, in the continuum limit) contribution to the wave amplitude. This behaviour, observed for dense spectrum systems, is described through (3.2) and its asymptotic form given in (3.3) and illustrated in figures 5 and 6.
The outcomes reported in this paper for systems with many degrees of freedom were made possible only by the development of a new technique to compute complex roots. This new root finding method is further detailed in appendix A.
2 Formalism of monokinetic beams
Our model system is composed of
$N$
charged resonant particles interacting with a single electrostatic Langmuir wave with natural frequency
$\unicode[STIX]{x1D714}_{0}$
, in one space dimension (Onishchenko et al.
Reference Onishchenko, Linetskiĭ, Matsiborko, Shapiro and Shevchenko1970; O’Neil, Winfrey & Malmberg Reference O’Neil, Winfrey and Malmberg1971). The evolution of this system is generated by the self-consistent Hamiltonian (Mynick & Kaufman Reference Mynick and Kaufman1978; Tennyson, Meiss & Morrison Reference Tennyson, Meiss and Morrison1994; Escande et al.
Reference Escande, Zekri and Elskens1996; Elskens & Escande Reference Elskens and Escande2003; Escande Reference Escande, Dauxois, Ruffo and Cugliandolo2010)
$$\begin{eqnarray}H_{sc}=\mathop{\sum }_{l=1}^{N}\frac{p_{l}^{2}}{2}+\unicode[STIX]{x1D714}_{0}\frac{X^{2}+Y^{2}}{2}+\unicode[STIX]{x1D700}k_{w}^{-1}\mathop{\sum }_{l=1}^{N}(Y\sin k_{w}x_{l}-X\cos k_{w}x_{l}),\end{eqnarray}$$
where
$X$
and
$Y$
correspond to the Cartesian components of the complex wave amplitude
$Z=X+\text{i}Y$
, also expressed as
$Z=\sqrt{2I}\text{e}^{-\text{i}\unicode[STIX]{x1D703}}$
in terms of the phase
$\unicode[STIX]{x1D703}$
and intensity
$I$
. We assume that particles have periodic boundary conditions in the interval of length
$L$
, and the wavenumber is
$k_{w}=2\unicode[STIX]{x03C0}j/L$
for some integer
$j$
. Parameter
$\unicode[STIX]{x1D700}$
is the wave–particle coupling constant (which may be determined from further consideration of the underlying plasma).
The first and second terms in (2.1) represent the kinetic energy of the resonant particles (generating ballistic motion) and the energy of the free wave (related to the vibratory motion of the non-resonant bulk particles in the plasma), respectively. The latter (also responsible for the nonlinear nature of the dynamics) corresponds to the interaction energy between the wave and resonant particles.
The equations of motion are directly derived from the Hamiltonian (2.1),
$$\begin{eqnarray}\displaystyle & {\dot{x}}_{l}=p_{l},\quad 1\leqslant l\leqslant N, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\dot{p}}_{l}=\unicode[STIX]{x1D700}\,\text{Im}(Z\text{e}^{\text{i}k_{w}x_{l}}),\quad 1\leqslant l\leqslant N, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{Z}}=-\text{i}\unicode[STIX]{x1D714}_{0}Z+\text{i}\unicode[STIX]{x1D700}k_{w}^{-1}\mathop{\sum }_{l^{\prime }=1}^{N}\text{e}^{-\text{i}k_{w}x_{l^{\prime }}}. & \displaystyle\end{eqnarray}$$
This dynamical system admits as equilibrium states the configurations where the electrostatic field has zero amplitude and the particles are distributed in monokinetic beams (as in Dawson’s beams model), labelled
$1\leqslant s\leqslant b$
, characterized by velocities
$v_{s}$
and number of particles
$N_{s}$
. Assuming
$$\begin{eqnarray}x_{ns}^{(0)}(t)=v_{s}t+nL/N_{s}+\unicode[STIX]{x1D719}_{s}\end{eqnarray}$$
as the position of the
$n$
th particle of beam
$s$
with the number of particles per beam satisfying the condition
$N_{s}>j$
, the sum in (2.4) vanishes at any instant in time and, consequently,
$(x_{l}(t)=x_{ns}^{(0)}(t),p_{l}(t)=v_{s},Z(t)=0)$
corresponds to an exact, equilibrium solution with vanishing wave. No similar solution exists with non-zero
$Z$
(Elskens Reference Elskens2001).
The dynamics preserves total energy
$H$
and total momentum
$P=k_{w}I+\sum _{l}p_{l}$
. The reversibility of Hamiltonian dynamics is expressed by the invariance of system (2.2)–(2.4) under the time-reversal map
$(x^{\prime },p^{\prime },X^{\prime },Y^{\prime },t^{\prime },H^{\prime },k_{w}^{\prime })=(x,-p,-X,Y,-t,H,-k_{w})$
. Here, the reversal of
$k_{w}$
accounts for the reversal of total momentum
$P$
under this map, as well as invariance of total energy
$H$
.
Considering the system slightly displaced from this equilibrium configuration, the evolution of the wave amplitude and the perturbations in particle orbits obey the linearized equations of motion
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}{\dot{x}}_{ns}=\unicode[STIX]{x1D6FF}p_{ns}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}{\dot{p}}_{ns}=-\unicode[STIX]{x1D700}\,\text{Im}(Z\text{e}^{-\text{i}k_{w}x_{ns}^{(0)}}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{Z}}=-\text{i}\unicode[STIX]{x1D714}_{0}Z+\unicode[STIX]{x1D700}\mathop{\sum }_{s=1}^{b}\mathop{\sum }_{n=1}^{N_{s}}\unicode[STIX]{x1D6FF}x_{ns}\text{e}^{-\text{i}k_{w}x_{ns}^{(0)}}, & \displaystyle\end{eqnarray}$$
where
$b$
is the number of beams (
$\sum _{s=1}^{b}N_{s}=N$
),
$\unicode[STIX]{x1D6FF}x_{ns}=x_{ns}-x_{ns}^{(0)}$
and
$\unicode[STIX]{x1D6FF}p_{ns}=p_{ns}-v_{s}$
.
The perturbations in positions and momenta can be expressed in Fourier series with coefficients
$C$
and
$A$
(Escande et al.
Reference Escande, Zekri and Elskens1996),
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}x_{ns}(t)=2\,\text{Im}(C_{s}(t)\text{e}^{\text{i}k_{w}x_{ns}^{(0)}(t)})-\text{i}\mathop{\sum }_{m\in \unicode[STIX]{x1D707}_{s}^{\prime }}C_{ms}(t)\text{e}^{\text{i}k_{m}x_{ns}^{(0)}(t)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FF}p_{ns}(t)=2\,\text{Re}(A_{s}(t)\text{e}^{\text{i}k_{w}x_{ns}^{(0)}(t)})+\mathop{\sum }_{m\in \unicode[STIX]{x1D707}_{s}^{\prime }}A_{ms}(t)\text{e}^{\text{i}k_{m}x_{ns}^{(0)}(t)}, & \displaystyle\end{eqnarray}$$
where the set
$\unicode[STIX]{x1D707}_{s}^{\prime }$
is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D707}_{s}^{\prime }=\left\{\begin{array}{@{}ll@{}}\{m\in \mathbb{Z}:|m|\leqslant (N_{s}-1)/2~\text{and}~m\neq \pm j\}, & \text{if}~N_{s}~\text{is odd,}\\ \{m\in \mathbb{Z}:1-N_{s}/2\leqslant m\leqslant N_{s}/2~\text{and}~m\neq \pm j\}, & \text{if}~N_{s}~\text{is even.}\end{array}\right.\end{eqnarray}$$
The first terms in (2.9) and (2.10), which appear outside the sum, correspond to what we call the wave-like part of the solution. We omit the subscript
$j$
in the coefficients since they are the only components of the wave-like kind. The second terms, involving the sums over
$\unicode[STIX]{x1D707}_{s}^{\prime }$
, corresponds to the ballistic parts whose evolutions are readily expressed in terms of the initial conditions,
$$\begin{eqnarray}\displaystyle & C_{ms}(t)=(C_{ms}(0)+\text{i}A_{ms}(0)t)\text{e}^{-\text{i}k_{m}v_{s}t} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & A_{ms}(t)=A_{ms}(0)\text{e}^{-\text{i}k_{m}v_{s}t}, & \displaystyle\end{eqnarray}$$
because these terms do not couple with
$Z$
as
$\sum _{n}\text{e}^{-\text{i}k_{m}x_{ns}^{(0)}}=0$
for each
$m$
,
$s$
. On the contrary, due to their dependence on the wave amplitude
$Z$
, the wave-like components of the Fourier coefficients have a more intricate form that requires calculating the eigenvalues of the linearized system and expanding initial data in terms of the corresponding eigenvectors.
2.1 Wave-like solution and dispersion relation
The ‘lattice’ Fourier coefficients of the wave-like part satisfy the system of differential equations
$$\begin{eqnarray}\displaystyle & {\dot{C}}_{s}=-\text{i}k_{w}v_{s}C_{s}+\text{i}A_{s},\quad 1\leqslant s\leqslant b, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{A}}_{s}=-\text{i}k_{w}v_{s}A_{s}+\text{i}\frac{\unicode[STIX]{x1D700}}{2}Z,\quad 1\leqslant s\leqslant b, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{Z}}=\text{i}\unicode[STIX]{x1D714}_{0}Z-\text{i}\unicode[STIX]{x1D700}\mathop{\sum }_{s=1}^{b}C_{s}N_{s}. & \displaystyle\end{eqnarray}$$
Assuming solutions of the form
${\sim}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}t}$
, the problem of solving the
$2b+1$
equations of motion reduces to solving the algebraic linear system
$$\begin{eqnarray}\unicode[STIX]{x1D70E}{\mathcal{C}}=\unicode[STIX]{x1D648}\boldsymbol{\cdot }{\mathcal{C}},\end{eqnarray}$$
where
${\mathcal{C}}\in \mathbb{C}^{2b+1}$
denotes
${\mathcal{C}}=[c_{1},\ldots ,c_{b},a_{1},\ldots ,a_{b},z]^{\text{T}}$
, whose values of
$c_{s}$
and
$a_{s}$
represent, respectively, the time-independent part of the Fourier coefficients for the positions and momenta. The real matrix
$\unicode[STIX]{x1D648}\in \mathbb{R}^{(2b+1)\times (2b+1)}$
reads
$$\begin{eqnarray}\unicode[STIX]{x1D648}=\left(\begin{array}{@{}ccccccc@{}}k_{w}v_{1} & & 0 & -1 & & 0 & 0\\ & \ddots & & & \ddots & & \vdots \\ 0 & & k_{w}v_{b} & 0 & & -1 & 0\\ 0 & & 0 & k_{w}v_{1} & & 0 & -\unicode[STIX]{x1D700}/2\\ & \ddots & & & \ddots & & \vdots \\ 0 & & 0 & 0 & & k_{w}v_{b} & -\unicode[STIX]{x1D700}/2\\ \unicode[STIX]{x1D700}N_{1} & \ldots \, & \unicode[STIX]{x1D700}N_{b} & 0 & \ldots \, & 0 & \unicode[STIX]{x1D714}_{0}\end{array}\right).\end{eqnarray}$$
Denoting by
${\mathcal{C}}_{r}=[c_{r1},\ldots ,c_{rb},a_{r1},\ldots ,a_{rb},z_{r}]^{\text{T}}$
the
$r$
th eigenvector of
$\unicode[STIX]{x1D648}$
, its associated eigenvalue
$\unicode[STIX]{x1D70E}_{r}$
is obtained through the characteristic equation
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{r}=\unicode[STIX]{x1D714}_{0}+\frac{\unicode[STIX]{x1D700}^{2}}{2}\mathop{\sum }_{s=1}^{b}\frac{N_{s}}{(\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s})^{2}}.\end{eqnarray}$$
For a given wavenumber
$k_{w}$
, (2.19) states the condition to be satisfied by the eigenfrequencies
$\unicode[STIX]{x1D70E}_{r}$
for the system (2.14)–(2.16) to admit non-trivial solutions. This equation is thus a dispersion relation for the eigenmodes of the linearized system. It can be transformed into a polynomial equation of degree
$2b+1$
with real coefficients, and therefore admits generically
$2b+1$
complex roots, with at least one of them purely real (see figure 1
b).
Figure 1. Illustration of a possible configuration for a system composed of one wave with natural frequency
$\unicode[STIX]{x1D714}_{0}$
and ten monokinetic beams. The beams velocities
$v_{s}=v_{1}+(s-1)\unicode[STIX]{x0394}p$
are set around the wave phase velocity
$\unicode[STIX]{x1D714}_{0}/k$
. (a) Distribution of particles velocities with a constant slope, given according to
$f(v_{s})=f(v_{1})+(s-1)\unicode[STIX]{x0394}pf^{\prime }$
(yielding
$N_{s}=N\unicode[STIX]{x0394}pf(v_{s})$
as the number of particles) in
$I_{v}=[v_{1}-\unicode[STIX]{x0394}p/2,v_{b}+\unicode[STIX]{x0394}p/2]$
and set to zero outside this interval. (b) Graphical representation of the dispersion relation (2.19) through the line
$y=\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{0}$
and the many branches curve
$y=\unicode[STIX]{x1D712}(\unicode[STIX]{x1D714})\equiv (\unicode[STIX]{x1D700}^{2}/2)\sum _{s}(\unicode[STIX]{x1D714}-kv_{s})^{-2}N_{s}$
. The intercept point locates the single purely real solution to this equation.
The continuous limit corresponds to letting
$N\rightarrow \infty$
while keeping
$\unicode[STIX]{x1D700}^{2}N$
fixed. Particles are described with the distribution function (Firpo & Elskens Reference Firpo and Elskens1998)
$$\begin{eqnarray}{\mathcal{F}}(x,p)=\lim _{N\rightarrow \infty }\frac{1}{N}\mathop{\sum }_{l}\unicode[STIX]{x1D6FF}(x-x_{l})\unicode[STIX]{x1D6FF}(p-p_{l})\end{eqnarray}$$
and the wave is rescaled to
$$\begin{eqnarray}{\mathcal{Z}}=\unicode[STIX]{x1D700}Z\end{eqnarray}$$
so that the dynamics (2.2)–(2.3) generates the characteristics of the Vlasov equation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}{\mathcal{F}}+p\unicode[STIX]{x2202}_{x}{\mathcal{F}}-\text{Im}({\mathcal{Z}}\text{e}^{-\text{i}k_{w}x})\unicode[STIX]{x2202}_{p}{\mathcal{F}}=0\end{eqnarray}$$
coupled with the wave evolution
$$\begin{eqnarray}\dot{{\mathcal{Z}}}=-\text{i}\unicode[STIX]{x1D714}_{0}{\mathcal{Z}}+\text{i}k_{w}^{-1}\int _{0}^{L}\int _{-\infty }^{\infty }\text{e}^{-\text{i}k_{w}x}{\mathcal{F}}(x,p)\,\text{d}p\,\text{d}x.\end{eqnarray}$$
The equilibrium reference state is
$({\mathcal{Z}}=0,{\mathcal{F}}(x,v)=L^{-1}f(v))$
, with the velocity distribution function such that
$N_{s}=N\int _{v_{s}-\unicode[STIX]{x0394}p/2}^{v_{s}+\unicode[STIX]{x0394}p/2}f(p)\,\text{d}p$
: in the continuous limit, beam velocities are continuously distributed.
The dispersion relation in this limit follows from (2.19) by replacing the sum over the number of particles with an integral over the interval containing the beams velocities weighted by the normalized distribution
$f(v_{s})$
. Formally, this leads to a singular integral equation if
$\unicode[STIX]{x1D70E}$
is real and
$f(\unicode[STIX]{x1D70E}/k_{w})$
does not vanish; with an
$\text{i}\unicode[STIX]{x1D716}$
prescription, Landau (Reference Landau1946) obtains a solution
$\unicode[STIX]{x1D70E}_{L}$
with imaginary part
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{L}=\frac{\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}^{2}N}{2k_{w}^{2}}f^{\prime }(\unicode[STIX]{x1D714}_{0}/k_{w}),\end{eqnarray}$$
implying that the equilibrium is unstable for positive
$f^{\prime }(\unicode[STIX]{x1D714}_{0}/k_{w})$
whereas perturbations are damped for negative
$f^{\prime }(\unicode[STIX]{x1D714}_{0}/k_{w})$
.
For equally spaced beams,
$v_{s}-v_{s^{\prime }}=(s-s^{\prime })\unicode[STIX]{x0394}p$
, analytical estimates (Dawson Reference Dawson1960; Escande et al.
Reference Escande, Zekri and Elskens1996; Elskens & Escande Reference Elskens and Escande2003) show that a pair of eigenvalues
$\unicode[STIX]{x1D714}_{r}\pm \text{i}\unicode[STIX]{x1D6FE}_{r}$
is located near each beam (with
$|v_{s}-\unicode[STIX]{x1D714}_{r}/k_{w}|<\unicode[STIX]{x0394}p/2$
), with imaginary part
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{r}\simeq -\frac{1}{\unicode[STIX]{x1D70F}_{rec}}\ln \left|\frac{\unicode[STIX]{x1D6FF}\sqrt{1+\unicode[STIX]{x1D703}_{r}^{2}}}{4\unicode[STIX]{x03C0}}\right|+\cdots ,\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}:=f^{\prime }(v_{s})\unicode[STIX]{x0394}p/f(v_{s})$
and
$\unicode[STIX]{x1D703}_{r}=O(1)$
. Their real part
$\unicode[STIX]{x1D714}_{r}$
is near
$(v_{s}+(\unicode[STIX]{x0394}p/4)\,\text{sign}\,\unicode[STIX]{x1D6FF})k_{w}$
. The recurrence time
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{rec}=\frac{2\unicode[STIX]{x03C0}}{k_{j}\unicode[STIX]{x0394}p}\end{eqnarray}$$
characterizes the time scale on which the approximation of the many beams by a smooth distribution function loses its validity.
In the continuum limit, the imaginary parts
$\pm \unicode[STIX]{x1D6FE}_{r}$
of these eigenfrequencies as well as their spacing go to zero, so they approach a continuum spectrum of real eigenfrequencies which corresponds to the analogue of van Kampen modes found in the Vlasovian approach.
When the distribution function has a negative slope, there are
$2b$
such van Kampen-like eigenvalues, and just one real eigenvalue beyond the fastest beam (see figure 1
b). In contrast, for a positive
$f^{\prime }(\unicode[STIX]{x1D714}_{0}/k_{w})$
, besides these solutions condensing to the real axis in the limit, (2.19) is also satisfied by a particular pair of eigenvalues
$\unicode[STIX]{x1D714}_{L}\pm \text{i}\unicode[STIX]{x1D6FE}_{L}$
whose imaginary part is given (in the limit
$\unicode[STIX]{x0394}p\rightarrow 0$
) by (2.24).
A first benefit of the discretized model with monokinetic beams is its avoiding singular integrals and prescriptions of proper integration contour inherent to the kinetic approach. It also stresses a deep difference between the damping and growth cases: despite their being encompassed by a single formula (2.24), they must involve different physics, because the instability is associated with a single eigenvalue solving (2.19), while in the damping case (2.19) admits no solution with negative imaginary part near
$\unicode[STIX]{x1D6FE}_{L}$
.
This paradox reflects the reversibility of Hamiltonian dynamics: if damping resulted from a genuine eigenmode of the system, the dispersion relation would have a conjugate eigenvalue with positive real part, leading to an instability. This paradox of the Vlasov equation was solved by van Kampen (Reference van Kampen1955, Reference van Kampen1957) and Case (Reference Case1959).
A second indication that Landau damping does not result from the ‘dominant eigenmode’ in kinetic theory is that some initial data lead to a damping with smaller decay rate than Landau’s: when present, they dominate over the Landau behaviour for long times (Belmont et al. Reference Belmont, Chust, Mottez and Hess2011).
Actually, we shall see that both Landau damping and growth are related to a phase mixing effect among the van Kampen-like modes when the general solution of the linearized system is written in terms of its eigenfunctions.
2.2 Normal modes expansion
Denote by
${\mathcal{C}}_{r}$
and
$\unicode[STIX]{x1D70E}_{r}$
the
$r$
th eigenvector of
$\unicode[STIX]{x1D648}$
and its eigenvalue. Then, by the linearity of (2.17), the general solution
${\mathcal{G}}(t)=[C_{1}(t),\ldots ,C_{b}(t),A_{1}(t),\ldots ,A_{b}(t),Z(t)]^{\text{T}}$
is a superposition of the eigensolutions
${\mathcal{C}}_{r}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t}$
, i.e.
$$\begin{eqnarray}{\mathcal{G}}(t)=\mathop{\sum }_{r=1}^{2b+1}\unicode[STIX]{x1D709}_{r}{\mathcal{C}}_{r}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t},\end{eqnarray}$$
with components
$$\begin{eqnarray}\displaystyle & \displaystyle C_{s}(t)=\frac{\unicode[STIX]{x1D700}}{2}\mathop{\sum }_{r=1}^{2b+1}\frac{\unicode[STIX]{x1D709}_{r}}{(\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s})^{2}}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle A_{s}(t)=\frac{\unicode[STIX]{x1D700}}{2}\mathop{\sum }_{r=1}^{2b+1}\frac{\unicode[STIX]{x1D709}_{r}}{\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s}}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Z(t)=\mathop{\sum }_{r=1}^{2b+1}\unicode[STIX]{x1D709}_{r}\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t}. & \displaystyle\end{eqnarray}$$
The coefficients
$\unicode[STIX]{x1D709}_{r}$
of the linear combination are obtained from the initial condition and the left eigenvectors of matrix
$\unicode[STIX]{x1D648}$
,
$\unicode[STIX]{x1D70E}_{r}{\mathcal{C}}_{r}^{\prime }=\unicode[STIX]{x1D648}^{\text{T}}\cdot {\mathcal{C}}_{r}^{\prime }$
, through the inner product
$$\begin{eqnarray}\unicode[STIX]{x1D709}_{r}={\mathcal{C}}_{r}^{\prime \text{T}}\cdot {\mathcal{G}}(0)=z_{r}^{\prime }\left[Z(0)+\unicode[STIX]{x1D700}\mathop{\sum }_{s=1}^{b}N_{s}\left(\frac{C_{s}(0)}{\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s}}-\frac{A_{s}(0)}{(\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s})^{2}}\right)\right],\end{eqnarray}$$
where the scaling factor
$z_{r}^{\prime }$
, given by
$$\begin{eqnarray}z_{r}^{\prime }=\left[1+\unicode[STIX]{x1D700}^{2}\mathop{\sum }_{s=1}^{b}\frac{N_{s}}{(\unicode[STIX]{x1D70E}_{r}-k_{w}v_{s})^{3}}\right]^{-1},\end{eqnarray}$$
is chosen so that the normalization condition
${\mathcal{C}}_{r^{\prime }}^{\prime \text{T}}\cdot {\mathcal{C}}_{r}=\unicode[STIX]{x1D6FF}_{r^{\prime },r}$
is satisfied.
In the linear regime, the system of beams interacting with a single Langmuir wave is therefore analytically solvable in terms of a normal modes expansion of the linearized equations of motion. Within the present approach, specifically through (2.9)–(2.13) and (2.27)–(2.32), note that once the full spectrum of eigenfrequencies
$\{\unicode[STIX]{x1D70E}_{r}\}$
and the initial condition of the dynamics (the coefficients
$\unicode[STIX]{x1D709}_{r}$
) are known, the evolution of the wave amplitude and the perturbations on the particles orbits are determined.
In particular, the scaling factor
$z_{r}^{\prime }$
coincides with the weight of eigenvector
${\mathcal{C}}_{r}$
in the decomposition of a ‘quiet start’ initial condition (
$\unicode[STIX]{x1D6FF}x_{ns}(0)=0$
,
$\unicode[STIX]{x1D6FF}p_{ns}(0)=0$
, viz.
$C_{s}(0)=0,A_{s}(0)=0$
).
3 Numerical results and discussion
In this section, we discuss the spectrum of the van Kampen-like eigenfrequencies obtained by solving numerically (2.19) for many-beam systems and stress some microscopic aspect of the dynamics. As we shall see, both wave damping and growth arise as consequence of an interference among the
$2b+1$
eigenmodes. However, in respect to this interference (or phase mixing) effect, these two cases exhibit a remarkable distinctive behaviour as the number of beams becomes large.
The monokinetic beams approach presented in § 2 breaks down after times of the order of the bouncing time
$\unicode[STIX]{x1D714}_{b}^{-1}=(\unicode[STIX]{x1D700}|Z|k_{w})^{-1/2}$
, where the particle orbit may no longer be considered approximately ballistic. To ensure that trapping effects play a negligible role and the structure of beams is preserved, all the following results were obtained in the weak amplitude regime (easily implemented in taking
$\unicode[STIX]{x1D700}\ll 1$
,
$N\sim \unicode[STIX]{x1D700}^{-2}\rightarrow \infty$
). To simplify our numerical calculations, we take
$\unicode[STIX]{x1D714}_{0}=0.0$
which, from the physical viewpoint, amounts to performing a Galileo transformation to a reference frame moving with the wave (nominal) phase velocity. We set the equilibrium state with particle velocities equally spaced (
$v_{s}=v_{1}+(s-1)\unicode[STIX]{x0394}p$
) inside an interval of length
$I_{v}$
centred at the origin and the distribution function
$f(v_{s})=f(v_{1})+(s-1)f^{\prime }\unicode[STIX]{x0394}p$
, where
$f^{\prime }$
is constant and
$\unicode[STIX]{x0394}p=I_{v}/b$
is the velocity interval between two adjacent beams. The data
$f(v_{1})=1/I_{v}-(b-1)\unicode[STIX]{x0394}pf^{\prime }/2$
ensure that the distribution of velocities is normalized to unity. Thus, the number of particlesFootnote
2
in beam
$s$
is given by
$N_{s}=N\unicode[STIX]{x0394}pf(v_{s})$
.
Figure 2. Damping case: discrete analogue of van Kampen eigenfrequencies in the complex
$\unicode[STIX]{x1D70E}$
– plane (c), the phase (b) and the modulus (a) of the scaling factor. The curves in the two upper panels are guides to the eye. The red losangle locates the Landau value
$\unicode[STIX]{x1D70E}_{L}$
(which corresponds to no eigenmode in the finite-
$N$
model). The results were obtained for systems composed of 50 and 200 beams with parameters set to
$I_{v}=1.6$
,
$\unicode[STIX]{x1D714}_{0}=0.0$
,
$k_{w}=1.0$
,
$\unicode[STIX]{x1D700}Nf^{\prime }=-4.8$
and
$\unicode[STIX]{x1D6FE}_{L}=-2.26\times 10^{-2}$
.
In the weak-coupling regime and for many-beam systems, the spectrum of eigenfrequencies is very close to the real axis, and therefore the numerical solutions of (2.19) must be computed with high accuracy. To obtain the spectrum of eigenfrequencies for such systems with an accuracy down to
$10^{-7}$
(
$|\unicode[STIX]{x1D70E}_{r}^{approx}-\unicode[STIX]{x1D70E}_{r}|\leqslant 10^{-7}$
), we developed a root-finding scheme based on Cauchy’s residue theorem with triangular contours that enabled us to implement a bisection method in the complex
$\unicode[STIX]{x1D70E}$
-plane. In this paper, we shall not dwell on the technical details of the method; its key ideas are discussed in appendix A.
3.1 Spectra and scaling factors
In the lowest graph in figure 2, we plot the spectrum of van Kampen-like eigenfrequencies for the damping case (
$f^{\prime }<0$
) for systems composed of 50 and 200 beams. The upper two plots display the logarithm of the modulus of the scaling factor
$z_{r}^{\prime }$
and its phase
$\text{Arg}\,z_{r}^{\prime }$
, obtained from the spectra and (2.32). These curves correspond only to the eigenmodes with
$\unicode[STIX]{x1D6FE}_{r}\geqslant 0$
: thanks to the time reversibility of Hamiltonian dynamics implying real coefficients in the algebraic equation (2.32), the stable eigenmodes differ only by the opposite sign of their phase. Regardless of the number of beams, the phase velocities
$\unicode[STIX]{x1D714}_{r}/k_{w}$
of the damped and growing eigenmodes are confined to the interval
$[-0.8,0.8]$
that contains the resonant particles. The more beams we consider, denser the spectrum becomes and it approaches as a whole the real axis. The dependence of the mean value of the imaginary part
$\bar{\unicode[STIX]{x1D6FE}_{r}}=b^{-1}\sum _{r,\unicode[STIX]{x1D6FE}_{r}>0}\unicode[STIX]{x1D6FE}_{r}$
(which amounts to the
$\ell ^{1}$
norm of the sequence of these imaginary parts) on the beam spacing is illustrated in figure 3, indicating how, in the continuous limit (
$\unicode[STIX]{x0394}p\rightarrow 0$
), these eigenfrequencies accumulate toward the real axis forming what corresponds to the continuous spectrum of van Kampen frequencies derived from a Vlasovian description. The spacing between the real parts of successive eigenfrequencies also goes to zero like
$\unicode[STIX]{x0394}p$
because there is always a real
$\unicode[STIX]{x1D714}_{r}/k_{w}$
between two successive beam velocities, so that the condensation of eigenvalues onto the real axis tends to a cut along the interval
$[\inf (v_{s}),\sup (v_{s})]$
: this cut is the locus of van Kampen’s singular spectrum. This behaviour
$\unicode[STIX]{x1D6FE}_{r}\sim \unicode[STIX]{x0394}p\,|\ln \unicode[STIX]{x0394}p|$
, (2.25) was computed by Dawson (Reference Dawson1960).
Figure 3. Damping case: behaviour of the average distance to real axis
$\bar{\unicode[STIX]{x1D6FE}_{r}}\equiv b^{-1}\sum _{r,\unicode[STIX]{x1D6FE}_{r}>0}\unicode[STIX]{x1D6FE}_{r}$
as a function of the beam spacing. The linear regression
$y=(0.177\pm 0.001)x\pm 0.001$
shows that
$\unicode[STIX]{x1D6FE}_{r}$
goes to zero like
$x\equiv \unicode[STIX]{x0394}p\,|\log \unicode[STIX]{x0394}p|$
as computed by Dawson (Reference Dawson1960). The error bars represent the tolerance
$10^{-3}$
used in the root-finding routine.
Landau’s value
$\unicode[STIX]{x1D70E}_{L}=\unicode[STIX]{x1D714}_{L}+\text{i}\unicode[STIX]{x1D6FE}_{L}$
, represented by the red diamond point in figure 2, is obtained after taking the continuous limit of (2.19). The interaction between the wave and resonant particles is responsible for the small deviation
$|\unicode[STIX]{x1D714}_{L}-\unicode[STIX]{x1D714}_{0}|=2.30\times 10^{-2}$
from the free wave frequency, which can also be estimated analytically (Elskens & Escande Reference Elskens and Escande2003).
We also observe, in the vicinity of
$\unicode[STIX]{x1D714}_{L}$
, a striking behaviour for the modulus and for the phase of the scaling factor. The top plot evidences a peak slightly shifted from the origin, and the centre plot shows a change in the phase’s sign. These two behaviours together indicate that the modes with greater contribution to the initial data
$(\unicode[STIX]{x1D6FF}x_{ns}(0),\unicode[STIX]{x1D6FF}p_{ns}(0),Z(0))$
are those with a frequency close to the Landau value, in agreement with the estimate (Elskens & Escande Reference Elskens and Escande2003)
$$\begin{eqnarray}z_{r}^{\prime }\simeq -\frac{1}{2\unicode[STIX]{x03C0}}\frac{k_{w}\unicode[STIX]{x0394}p}{\unicode[STIX]{x1D6FE}_{L}\pm \text{i}(\unicode[STIX]{x1D714}_{r}-\unicode[STIX]{x1D714}_{L})}\end{eqnarray}$$
from (2.32). The plus-or-minus sign in this equation take care of the two
$r$
eigenmodes with same phase velocity.
Figure 4. Growth case: (a) modulus and phase of the scaling factor (the lines are guides to the eye) and (b) spectrum of the van Kampen-like eigenfrequencies in the complex
$\unicode[STIX]{x1D70E}$
plane for systems composed of 100, 500 and 2000 beams. The red losangle locates the Landau mode, which is slightly shifted from the origin by
$\unicode[STIX]{x1D714}_{L}=-4.87\times 10^{-3}$
. The parameters were set to
$I_{v}=1.0$
,
$\unicode[STIX]{x1D714}_{0}=0.0$
,
$k_{w}=1.0$
,
$\unicode[STIX]{x1D700}Nf^{\prime }=16.0$
and
$\unicode[STIX]{x1D6FE}_{L}=2.51\times 10^{-3}$
.
For the growth case (
$f^{\prime }>0$
), the spectrum and the scaling factor shown in figure 4 exhibit two distinctive features (compared to the damping regime) for systems with many beams. The first one is the presence of two specific eigenfrequencies, complex conjugate to each other, that do not approach the real axis as the number of beams increase. We denote the eigenfrequency close to the Landau frequency by
$\unicode[STIX]{x1D70E}_{bL}$
(Landau-like eigenmode, with subscript
$b$
recalling the finite number of beams) and its complex conjugate by
$\unicode[STIX]{x1D70E}_{bL^{\ast }}$
(anti-Landau eigenmode), both highlighted in the figure. The second striking difference is the spike (note the log scale) in the curve for the modulus of the scaling factor. It shows, along with
$\text{Arg}(z_{bL}^{\prime })\approx 0$
, that the Landau and the anti-Landau eigenmodes have a dominant contribution to the initial data.
3.2 Quiet start evolution
In order to monitor the contributions of the
$bL$
and
$bL^{\ast }$
eigenmodes to the wave amplitude, we consider a single realization of the system where the wave is launched with initial amplitude
$Z(0)=1$
and particles start from the equilibrium configuration of monokinetic arrays. For this specific quiet start realization, the evolution of the wave amplitude is given according to (2.30) by
$$\begin{eqnarray}Z(t)=z_{bL}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{bL}t}+z_{bL^{\ast }}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{bL^{\ast }}t}+\mathop{\sum }_{\substack{ r=1 \\ r\neq bL,bL^{\ast }}}^{2b+1}z_{r}^{\prime }\text{e}^{-\text{i}\unicode[STIX]{x1D70E}_{r}t},\end{eqnarray}$$
recalling that the scaling factor corresponds to a weight factor indicating how much each eigenmode contributes to construct the initial amplitude.
Figure 5. Growth case: evolution of real (a) and imaginary (b) parts of the Landau (
$bL$
), anti-Landau (
$bL^{\ast }$
) and van Kampen-like components in the wave complex amplitude
$Z=X+\text{i}Y$
for a system of 2000 beams in the unstable case. The symmetry between the red and green curves shows the compensation between the anti-Landau eigenmode and the van Kampen spectrum. The Landau (blue) line is growing exponentially, exceeding largely the range of our ordinate axis.
Due to the fact that for many-beam systems
$z_{bL}^{\prime }\approx 1$
and
$\unicode[STIX]{x1D70E}_{bL}\approx \unicode[STIX]{x1D70E}_{L}$
, as seen from the black graphs in figure 4, one would expect the Landau-like (
$bL$
) eigenmode to provide, by itself, the correct growth of the wave. However, the contribution of the other
$2b$
eigenmodes still have to be taken into account. Thanks to a destructive interference, these
$2b$
eigenmodes superposed contribute only with a small oscillatory part that does not compromise the exponential growth
${\sim}\text{e}^{\unicode[STIX]{x1D6FE}_{L}t}$
of the wave. Figure 5 illustrate this destructive interference by monitoring the evolution of the Cartesian components of each term of (3.2) for a system of 2000 beams. This graph exhibits a clear symmetry between the red and green curves, showing that the
$O(1)$
contribution of the anti-Landau eigenmode is cancelled out by the superposition of the (dense, individually small) van Kampen-like eigenmodes. The dashed curves that appear as an envelope in figure 5 are given by
$\pm \text{e}^{-\unicode[STIX]{x1D6FE}_{L}t}$
and indicate that, in the same way as for the anti-Landau eigenmode, the superposition of the van Kampen-like eigenmodes also damps with the anti-Landau rate.
Analytically, indeed, the asymptotic form (3.1) for the scaling factor implies that the third term in (3.2) behaves like the Fourier transform of a Lorentz function, namely as
$-Z(0)\text{e}^{-\unicode[STIX]{x1D6FE}_{L}\,|t|}$
. Therefore, (3.2) reduces to
$$\begin{eqnarray}Z(t)\simeq Z(0)\text{e}^{\unicode[STIX]{x1D6FE}_{L}t}+Z(0)\text{e}^{-\unicode[STIX]{x1D6FE}_{L}t}-Z(0)\text{e}^{-\unicode[STIX]{x1D6FE}_{L}|t|}\end{eqnarray}$$
which reduces to
$Z(0)\text{e}^{\unicode[STIX]{x1D6FE}_{L}|t|}$
for both positive and negative
$t$
.
This form rescues the time reversibility in the solution to the initial value problem. If the response to the initial perturbation were described with only the Landau-like eigenmode, it would lead to a decreasing
$|Z(t)|$
for
$t\rightarrow -\infty$
. Yet time-reversal symmetry (unbroken by the quiet start initial condition) requires the wave to increase as
$t\rightarrow -\infty$
just as it does for
$t\rightarrow \infty$
. Therefore, the system evolution for all times must involve both the Landau-like and anti-Landau modes, the sum of which reads
$2Z(0)\cosh \unicode[STIX]{x1D6FE}_{L}t$
. But this sum does not start exactly exponentially at
$t=0$
, so that the van Kampen-like eigenmodes are needed to cancel the anti-Landau mode for
$t>0$
, and to cancel the Landau-like mode for
$t<0$
.
This peculiarity of the growth regime can be seen numerically only when considering a larger number of beams. The necessity of working with such many-beam systems motivated us to developed the Cauchy root-finding scheme once that traditional computer algebra systems failed in providing the full spectrum.
Note that the damping case is analytically simpler: it involves only van Kampen-like eigenmodes, for which the scaling factors given by (3.1) again yield the Fourier transform of a Lorentz function. However,
$\unicode[STIX]{x1D6FE}_{L}$
is now negative, so that
$z_{r}^{\prime }>0$
for
$\unicode[STIX]{x1D714}_{r}\approx \unicode[STIX]{x1D714}_{L}$
and (3.2) reduces to the well known
$$\begin{eqnarray}Z(t)\simeq Z(0)\text{e}^{\unicode[STIX]{x1D6FE}_{L}|t|}.\end{eqnarray}$$
Figure 6. Evolution of the wave intensity for systems composed of 30 (red), 50 (blue) and 200 (orange) beams for the stable (a) and the unstable (b) cases with the dashed lines corresponding to the Landau prediction. The same parameters that generated the spectra in figures 2 and 4 were used here. As initial condition we consider
$Z(0)=1.0$
and unperturbed trajectories
$\unicode[STIX]{x1D6FF}x_{ns}(0)=\unicode[STIX]{x1D6FF}p_{ns}(0)=0$
with particles starting in the configuration of equally spaced monokinetic beams.
Figure 6 show the evolution of the wave intensity
$I=Z^{\ast }Z/2$
, obtained through the superposition of the van Kampen-like eigenmodes for systems composed of 30, 50 and 200 beams. We observe that the Langmuir wave damps (or grows) initially with the expected Landau prediction
$\log I_{L}(t)=\log I(0)\pm 2|\unicode[STIX]{x1D6FE}_{L}|t$
, represented by the dashed black lines. For the unstable case, even for a small number of beams, where
$\max \{\unicode[STIX]{x1D6FE}_{r}\}>\unicode[STIX]{x1D6FE}_{L}$
and the
$bL$
and
$bL^{\ast }$
modes are not prominent, the discrete and continuous systems agree.
The divergences observed in the graphs occur for times close to the recurrence time
$\unicode[STIX]{x1D70F}_{rec}\sim 2\unicode[STIX]{x03C0}/\unicode[STIX]{x0394}\unicode[STIX]{x1D714}_{r}\sim 2\unicode[STIX]{x03C0}/(k_{w}\unicode[STIX]{x0394}p)$
where the lack of an effective phase mixing of the van Kampen-like eigenmodes along with the large values of the modulus of the complex amplitude
$z_{r}^{\prime }\text{e}^{\unicode[STIX]{x1D6FE}_{r}t}$
of the unstable eigenmodes make the wave intensity depart from the Landau line. After this characteristic time, the approximation of a discrete system composed of monokinetic beams by a continuous one is no longer valid (Firpo & Elskens Reference Firpo and Elskens1998). It is worth noting that the departure of the orange curve from the Landau line in the top part of figure 6 is not related to the breakdown of the phase mixing but to the small values of
$I_{L}(t)$
compared to the amplitude of the oscillations.
4 Summary and perspectives
Our results in this paper provide an accurate numerical support for studying the phase mixing mechanism in Hamiltonian models. Success in investigating this process for many-beams systems rests on the development of a new root finding scheme based on the Cauchy’s integral theorem enabling one to compute all the roots of the dispersion relation (2.19).
Analysing the spectrum of eigenfrequencies of the van Kampen-like modes and the evolution of the Cartesian components of the wave amplitude highlights a remarkable aspect in which the stable and unstable regimes differ. In the stable case, the normal modes simply contribute collectively in providing the Landau damping; however, in the unstable case, the pure exponential growth of the wave amplitude is the result of a destructive interference effect (between its dual exponentially decaying mode, with same initial amplitude, and a flurry of van Kampen modes with small individual amplitudes) that enables a single eigenmode to provide a dominant contribution on the dynamics from the very beginning. We conclude the paper pointing out the consistence between the discrete and the continuous systems (Firpo & Elskens Reference Firpo and Elskens1998), and verifying that the superposition of the van Kampen-like eigenmodes indeed provides the expected exponential Landau damping/growth rate for the wave intensity.
The quiet start as well as the equally spaced beams setting were chosen in order to simplify the calculations. However, one expects similar results for other discretizations and initial conditions of the system in the many beams regime.
Acknowledgements
The first author thanks Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) for financing his stay at Aix-Marseille Université under the Programa de Doutorado Sanduíche no Exterior (PDSE), process no. 99999.004512/2014-06. Valuable discussions with B. V. Ribeiro and comments from D. F. Escande, M. A. Amato, F. Doveil and D. F. G. Minenna are gratefully acknowledged. The authors also acknowledge the reviewers for their useful comments.
Appendix A. Bisection-like root finding method in the complex plane
The Cauchy integral theorem states that if
$F$
is an analytic function on a simply connected domain
${\mathcal{D}}$
bounded by a simple closed contour
$\unicode[STIX]{x2202}{\mathcal{D}}$
, then
$\oint _{\unicode[STIX]{x2202}D}F(z)\,\text{d}z=0$
, and if
$F$
has a finite number of simple poles
$z_{j}$
in
${\mathcal{D}}$
with residues
$R_{j}$
, then
$\oint _{\unicode[STIX]{x2202}D}F(z)\,\text{d}z=2\unicode[STIX]{x03C0}\text{i}\sum _{j}R_{j}$
. In this appendix, we expose a direct scheme based on this elementary tool of complex analysis with the aim of computing numerically the full set of van Kampen-like eigenfrequencies.
Figure 7. Illustration of the bisection-like method in the complex
$\unicode[STIX]{x1D70E}$
-plane by supposing that the root is located to the right side of the line corresponding to the
$s$
th beam. The shaded triangles highlight the regions for which the contour integral is non-zero and that, consequently, contain a root inside.
With
$\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E})$
being the right-hand side of (2.19), the method consists in obtaining the root of
$\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E})$
, inside a given region, by searching for the pole of
$(\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E}))^{-1}$
. To find a root in the vicinity of the velocity of the
$s$
th beam, we define two rectangles: the left one with vertices
$A_{1}$
to
$A_{4}$
and the right one with vertices
$A_{1}^{\prime }$
to
$A_{4}^{\prime }$
, both illustrated in figure 7(a). Firstly, we estimate the Cauchy integral along these rectangular contours in order to identify to which side from the
$s$
th beam the root is located. If one of these rectangles yields
$$\begin{eqnarray}\left|\oint _{\unicode[STIX]{x2202}{\mathcal{D}}}(\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E}))^{-1}\,\text{d}\unicode[STIX]{x1D70E}\right|\leqslant C_{crit},\end{eqnarray}$$
with
$C_{crit}$
representing the Cauchy criterion to assume the integral numerically zero, it is discarded as a possible region containing the root.
After identifying the presence of a root inside a rectangle, the next step consists in implementing the bisection-like method which is based on successive divisions of the region of search. In figure 7, we illustrate this procedure with the shaded triangles in figures 7(a) and 7(b) showing, respectively, the configuration at the beginning and after 6 iterations of the method. Iterations are stopped when the largest side of the shaded triangle is smaller than the tolerance of the method (which is an input data). Following so, the approximate root of
$\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E})$
is given by the geometrical centre of the triangle.
The bisection-like method is displayed through the pseudo-code in Algorithm 1, where we called the function ‘Larger_side (A,B,C)’, that yields the length of the largest side of triangle ABC, and the function ‘exist_pole_inside (A,B,C)’ whose Boolean output informs, by means of the Cauchy criterion (A 1), whether there is a pole of
$(\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D712}(\unicode[STIX]{x1D70E}))^{-1}$
inside the triangle ABC. The code generating the results in this paper was implemented in C language with the ‘complex.h’ library.
Although the method was built specifically to tackle the problem of monokinetic beams presented in this paper, it can a priori be used in an efficient way to find the complex roots of any function
$F(\unicode[STIX]{x1D70E})$
which is well behaved in the integration domains and such that the poles of
$F^{-1}(\unicode[STIX]{x1D70E})$
have non-vanishing residues.
