Notes for Cours Ecole Doctorale

Advanced Theory of Plasmas

Ecole Polytechnique Fédérale de Lausanne

Stephan Brunner
stephan.brunner@epfl.ch

Centre de Recherches en Physique des Plasmas

Association Euratom-Confédération Suisse
Ecole Polytechnique Fédérale de Lausanne
CRPP-PPB, CH-1015 Lausanne, Switzerland

1
Contents

1 Non-Linear Effects in Plasmas 3

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Non-Linear Evolution of an Electron Plasma Wave . . . . . . . . . . . . . 6
1.2.1 Motivation/Illustration . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Re-Deriving Linear Landau Damping Invoking Energy Conservation 8
1.2.3 Limit of Linear Landau Damping . . . . . . . . . . . . . . . . . . . 13
1.2.4 “Non-Linear Landau Damping” . . . . . . . . . . . . . . . . . . . . 15
1.2.5 BGK waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.2.6 Non-Linear Frequency Shift . . . . . . . . . . . . . . . . . . . . . . 26
1.2.7 Stability of BGK Mode: The Trapped Particle Instability . . . . . . 32
1.2.8 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
1.3 Three Wave Interactions and Parametric Instabilities . . . . . . . . . . . . 33
1.3.1 System of Three Coupled Oscillators . . . . . . . . . . . . . . . . . 34
1.3.2 Illustration of Three Wave Coupling: Stimulated Raman Scattering 37
1.3.3 Matching Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 45
1.3.4 Manley-Rowe Relations . . . . . . . . . . . . . . . . . . . . . . . . . 48
1.3.5 Non-Linear Analytic Solution to the Space Independent Coupling
Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
1.3.6 Saturation Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . 58
1.3.7 Illustrations from Simulations . . . . . . . . . . . . . . . . . . . . . 58
1.3.8 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

A Particles in a Sinusoidal Potential 61

A.1 Trapped and Untrapped Particles . . . . . . . . . . . . . . . . . . . . . . . 61
A.2 Deeply Trapped Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
A.3 General Case of Untrapped Particles . . . . . . . . . . . . . . . . . . . . . 63
A.4 General Case of Trapped Particles . . . . . . . . . . . . . . . . . . . . . . . 63

B Solving the (Linearized) Vlasov Equation by Integrating Along Trajec-

tories 65
B.1 Brief Review of the Vlasov Equation . . . . . . . . . . . . . . . . . . . . . 65
B.2 Solving the Vlasov Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 67
B.3 Solving the Linearized Vlasov Equation . . . . . . . . . . . . . . . . . . . . 67

2
Chapter 1

Non-Linear Effects in Plasmas

1.1 Introduction
For many plasmas of interest, the Vlasov-Maxwell system of equations can be considered
as the “fundamental” kinetic description:
• The Vlasov equation describes how the evolution of the phase space distribution
fα (~x, ~v , t) for each plasma species α (= electrons “e”, ions “i”) is subject to macro-
scopic electric and magnetic fields (E, ~ B):
~

∂fα ∂fα qα  ~ ~ · ∂fα = 0.


+ ~v · + E + ~v × B (1.1)
∂t ∂~x mα ∂~v
This equation is justified in the limit of weakly coupled plasmas, characterized by
a small value p  1 of the plasma parameter p = 1/(N λ3D ) (N is the density,
and λD the Debye length). One recalls, that the parameter p is a measure of the
relative fluctuation level of interaction energy due to particle discreteness compared
to the kinetic energy. The weakly coupled approximation is justified for handling
most plasmas of interest in magnetic fusion.
~ x, t), B(~
• In turn, Maxwell’s equations describe the evolution of the fields [E(~ ~ x, t)]:

~ ~
~ = − ∂B ,
∇×E ∇×B~ = µ0 ~j + 1 ∂ E , (1.2)
∂t c2 ∂t
~ ρ
~ = ,
∇ · B = 0, ∇·E (1.3)
0

where the charge density ρ and current density ~j are generated by the plasma itself:
X Z X Z
ρint = qα d~v fα , ~j int = qα d~v ~v fα , (1.4)
α α

as well as by possible external sources (ρext , ~j ext ), so that in general:

ρ = ρint + ρext , ~j = ~j int + ~j ext . (1.5)

3
The system of equations (1.1)-(1.5) clearly show how the plasma is both subject to, and
~ B)
the source of, electromagnetic fields. The distributions fα and the fields (E, ~ must
therefore be solved for self-consistently.

Note furthermore, that in terms of the unknown quantities fα and (E, ~ B),
~ equations
(1.1)-(1.5) form a non-linear system of integro-differential equations. Indeed, the last
term in Vlasov’s equations (1.1) represents a quadratic non-linearity. This is in fact the
only non-linear term in this system, as Maxwell’s equation are themselves fully linear
[including the evaluation of the sources (ρint , ~j int ) in terms of fα ].

Although collisional effects will not be discussed in any detail in this chapter, let us
nonetheless briefly comment here on the additional non-linearity found in the more gen-
eral Fokker-Planck equation. Indeed, in cases for which the fluctuations due to binary
scattering effects can not be fully neglected, the Vlasov equation is replaced by the more
general Fokker-Planck equation:
∂fα ∂fα qα  ~ ~ · ∂fα =
 X
+ ~v · + E + ~v × B C[fβ , fα ],
∂t ∂~x mα ∂~v β

where the collision operator C[fβ , fα ] represents the scattering of species α off of species
β. This collision operator can in many cases be modeled by the Landau-type operator:
Z !
∂ 0 0 1 ∂ 1 ∂
C[fβ , fα ] = Γα,β · d~v U(~v − ~v ) · 0
− fβ (~v 0 )fα (~v ), (1.6)
∂~v mβ ∂~v mα ∂~v

where Γα,β = qα2 qβ2 ln Λ/(8π20 mα ), ln Λ is the Coulomb logarithm, and one has defined
the tensor U(~u) = (u2 1 − ~u : ~u)/u3 . The collision operator scales as C ∼ νc f , where
the collision frequency νc is of the order νc /ωp ∼ O(p ), ωp being the plasma frequency.
The collision operator (1.6) obviously provides an additional non-linearity to the Fokker-
Planck equation. The origin of this non-linearity is similar to the one in the Vlasov part
of the equation, i.e. the effect of self-consistent electromagnetic fields on the distribution.
The non-linearity in C however results from the Coulomb forces relative to random bi-
nary interactions, while the non-linearity in the Vlasov part results from the macroscopic
electro-magnetic fields from collective phenomena.

For understanding certain basic mechanisms, one can justify linearizing the Vlasov-Maxwell
system of equations. It is assumed here that the reader is familiar with this linear ap-
proximation, in particular as a first approach to studying the dispersion and dissipation
of small amplitude waves. However, plasmas are fundamentally non-linear in nature, and
the non-linearities pointed out above are thus essential for describing a whole set of im-
portant plasma phenomena.

Non-linear plasma theory is a vast topic, and this chapter only provides an introduc-
tion to the subject through a couple of specific examples. The first part of this chapter
considers the non-linear evolution of a single, finite amplitude Langmuir wave, and in

4
particular points out the break down of linear Landau damping. To model the non-
linear evolution of resonant wave-particle interaction naturally requires the framework of
a kinetic description. The second part of this chapter addresses the issue of non-linear
wave-wave interaction, and in particular the mechanism of parametric instabilities. These
wave-wave interactions do not involve resonant particles in their basic form, and thus can
be derived from fluid equations.

The examples of non-linear phenomena considered in this chapter are all cases of weak
non-linearity, for which the wave energy remains small compared to the total plasma en-
ergy. Under these conditions, the dispersion of the considered waves remains near the
dispersion predicted by linear theory (some non-linear corrections may nonetheless occur,
see the section on non-linear frequency shift), so that non-linear effects mainly affect the
evolution of the amplitudes of the waves, by either altering damping in the case of non-
linearly interacting resonant particles, or by providing coupling between waves in the case
of wave-wave interactions.

General Reference: Sagdeev and Galeev Ref. [1].

5
1.2 Non-Linear Evolution of an Electron Plasma Wave
1.2.1 Motivation/Illustration
In the following, the terms “Langmuir wave” and “electron plasma wave” (EPW) are used
interchangeably to describe the high frequency (near ωp ), electrostatic plasma oscillations.

Figure 1.1 presents results from non-linear Vlasov-Poisson simulations of Langmuir waves
using the SAPRISTI code1 . In this case, the code simply solves the 1-Dim Vlasov-Poisson
system of equations, where only the electron distribution fe (x, v, t) is evolved:
 Z 
∂f ∂f e ∂f ∂E 1
+v − E = 0, = −e dvf + qi Ni,0 . (1.7)
∂t ∂x m ∂v ∂x 0
For studying Langmuir waves, it is a good approximation to assume that the ions form a
fixed, homogeneous, neutralizing background, with density Ni,0 . Here, and in the follow-
ing, one drops the subscript “e” for electronic quantities, unless required for clarity.

For practical reasons, it is simpler to initiate a standing wave than a propagating wave.
This is done by choosing the initial electron distribution as a Maxwellian with a sinusoidal
density perturbation δN :
" #
δN N v2
f0 (x, v) = 1 + cos(k0 x) fM (v), fM (v) = exp(− 2
).
N (2π)1/2 vth 2vth

Figures 1.1 a) and b) present the evolution of the potential energy, the variation of kinetic
energy, and the total energy for initial conditions with perturbation levels δN/N = 0.01
and 0.1 respectively. For these illustrations, the wavelength is chosen such that k0 λD =
0.3, where λD = vth /ωp is the Debye length, and ωp2 = N e2 /m0 the squared plasma
frequency. By numerically solving the linear dispersion relation for Langmuir waves, one
obtains for k0 λD = 0.3: Frequency ω0 /ωp = 1.1598, and linear Landau damping rate
γL /ωp = 1.2620 · 10−2 .

Note, that throughout the simulation the total energy Etot (= Kinetic energy Kin of
electrons + Electrostatic potential energy Pot):
Z Z Z
m λ0
2 0 λ0
Etot = Kin + Pot = dx dv v f + dx E 2 ,
2 0 2 0

remains essentially constant. This is naturally expected as Etot is conserved by the system
of Eqs. (1.7) (check it!).

As shown in Fig. 1.1.a, the very first stage of the wave’s evolution in the non-linear
1
The acronym stands for (S)emi-Lagrangian (A)dvection code for (P)a(R)ametric (I)n(ST)ab(I)lity
studies. You are welcome to use this Fortran 90 code to familiarize with the different physical mechanisms
addressed in this chapter. It is available on the SVN web server of the CRPP.

6
 −4
x 10 δ N / N = 0.01 δ N / N = 0.1
Kinetic 0.03
Electrostatic
Total
Linear Damping

Energy / N T
2
Energy / N T

0.02

Linear Landau
damping rate γL
1 0.01

(a) (b)
0 0
0 50 100 150 200 250 0 50 100 150 200 250
t ωp tω
p

t ω = 250
th

t ωp = 250 p
v/v

5
v / vth

4.2 2∆v
trap
2∆v
trap
4
v vφ 4
φ
3.8

3.6
3

3.4
(c) (d)
0 5 10 15 20 0 5 10 15 20
x/λ x/λ
D D

Figure 1.1: Results from Vlasov-Poisson simulations. Non-linear evolution of a standing Lang-
muir wave. Potential energy, and variation of kinetic energy for initial density perturbation (a)
δN/N = 0.01 and (b) δN/N = 0.1 respectively. Distribution in phase space at tω p = 250 for (c)
δN/N = 0.01 and (d) δN/N = 0.1 respectively.

7
simulation is clearly exponential, and in agreement with the estimate γL for the Landau
damping from the linear theory. At later times however, the wave’s amplitude starts to
oscillate. This oscillation ultimately dies out, and the wave amplitude settles at a finite
value. The asymptotic state of the system is thus an undamped mode. These are obvi-
ously features not predicted by the linear theory of Landau damping. Figures 1.1 c) and
d) show the electron distribution at time t = 250 ωp−1 in phase space (x, v) near the phase
velocity vφ = ω0 /k0 for both cases δN/N = 0.01 and 0.1 respectively. These phase-space
plots clearly illustrate trapping for particles with velocity v in an interval of width 2∆vtrap
around vφ .

The purpose of the following sections is to study the non-linear effects illustrated by
these simulation results. For experimental confirmation of these effects, see the work by
Danielson [3] and the references therein.

In the following sections, the evolution of a traveling Langmuir wave will be consid-
ered, while the simulation results just presented here involve standing waves. So, will the
following discussion apply to these results? Yes, but why? A standing wave can naturally
be considered as the superposition of a forward traveling wave with phase velocity v φ ,
and a backward traveling wave with phase velocity −vφ . However, as one considers here
non-linear phenomena, one must be careful before invoking a superposition principal. As
will be discussed in detail, the non-linearities affecting each of these waves involve the
corresponding resonant particles, i.e. particles with velocities within an interval of order
∆vtrap from the corresponding phase velocity. As in general ∆vtrap /vφ  1, the reso-
nant particles relative to the forward and backward traveling waves thus usually form
two distinct groups with velocities in the vicinity of vφ and −vφ respectively. For this
reason the two traveling components forming a standing Langmuir wave can be assumed
non-interacting.

1.2.2 Re-Deriving Linear Landau Damping Invoking Energy Con-

servation
One considers an initially Maxwellian, homogeneous plasma. The study is limited here to
the evolution of an electron plasma wave, so that ions may be assumed fixed, providing a
neutralizing background for the mobile electrons.

One considers a slab-like system, so that the electron distribution f (x, v, t) verifies the
1-Dim Vlasov equation:
∂f ∂f e ∂f
+v − E = 0.
∂t ∂x m ∂v
One furthermore assumes the electrostatic wave to be essentially monochromatic of the
form:
E(x, t) = E0 (t) sin(k0 x − ω0 t),
where the amplitude E0 (t) evolves at a slow time scale compared to the frequency ω0 . i.e.
|(dE0 /dt)/E0 |  ω0 .

8
One is interested here in the rate of change of the kinetic energy Kin of the electrons
in the field E(x, t). The spatially averaged kinetic energy is defined by
Z
m
Kin = dv v 2 hf ix ,
2
where the brackets hix stand for the spatial average:
Z
1 λ0
hix = dx,
λ0 0

λ0 = 2π/k0 being one wavelength of the electrostatic field, which can be considered here
as the length of the periodic system.

The rate of change is then given by

Z
d Kin m ∂f
= dv v 2 h ix . (1.8)
dt 2 ∂t
One starts by considering a perturbative approach, so that the distribution is expanded
as
f = f0 (v) + f1 (x, v, t) + f2 (x, v, t) + . . . ,
where fn ∼ O(E n ), and f0 is the initial unperturbed state. Expanding the Vlasov Eq.
for the first two perturbation orders leads to
∂f1 ∂f1 e ∂f0
+v − E = 0, (1.9)
∂t ∂x m ∂v
∂f2 ∂f2 e ∂f1
+v − E = 0. (1.10)
∂t ∂x m ∂v
Obviously h∂f0 /∂tix = 0 as f0 6= f0 (t), and also h∂f1 /∂tix = 0 as f1 ∼ E ∼ cos(k0 x − ω0 t)
from Eq. (1.9). Thus, to lowest order in the perturbation, equation (1.8) is evaluated by
Z Z
d Kin m ∂f2 e ∂f1
= dv v 2 h ix = dv v 2 hE ix
dt 2 Z ∂t 2 ∂v
= −e dv vhE f1 ix , (1.11)

having made use of Eq. (1.10), h∂f2 /∂xix = 0, and performed an integration by parts in
the last step. Note that the last equality in Eq. (1.11) simply corresponds to dKin/dt =
~ x , where ~j is the electronic charge current.
h~j · Ei

One now addresses the problem of deriving hE f1 ix . The linear perturbation f1 can
be obtained from Eq.(1.9) by integrating along the unperturbed trajectories. Indeed, Eq.
(1.9) can be written as
df1 e ∂f0
= E , (1.12)
dt u.t. m ∂v

9
where d/dt|u.t. stands for the total time derivative along the unperturbed trajectories (for
more details on solving Vlasov-type equations by integrating along trajectories, see Ap-
pendix B). For the here considered homogeneous, unmagnetized plasma, the unperturbed
trajectories are simply given by (free streaming):
dx0 dv 0
= v0, = 0,
dt0 dt0
with initial conditions x0 (t) = x, v 0 (t) = v,

whose solution is x0 = x + v(t0 − t) and v 0 = v. Equation(1.12) can thus be solved by

integrating along these characteristics:
Z t df1 (x0 , v 0 , t0 ) e Zt 0 ∂f0 (v 0 )
f1 (x, v, t) − f1 (x − vt, v, 0) = dt0 = dt E(x0 , t0 )
0 dt0 u.t.
m 0 ∂v
Z
e ∂f0 t
= dt0 E0 (t0 ) sin [k0 x + k0 v(t0 − t) − ω0 t0 ] .(1.13)
m ∂v 0

As E(x, t) is assumed a self-consistent field, the initial perturbation f1 (x, v, 0) must be

such that Poisson’s equation is verified at t = 0:
Z
∂E(x, 0) −e
= k0 E0 (0) cos(k0 x) = dv f1 (x, v, 0).
∂x 0
This equation is verified for

f1 (x, v, 0) = f1,0 (v) cos(k0 x), (1.14)

with f1,0 (v) such that

Z
0 k 0
dv f1,0 (v) = −
E0 (0).
e
Combining (1.13) and (1.14), one thus obtains for the linear perturbation:
e ∂f0 Z t 0
f1 (x, v, t) = f10 (v) cos[k0 (x − vt)] + dt E0 (t0 ) sin [k0 x + k0 v(t0 − t) − ω0 t0 ] .
m ∂v 0
(1.15)
The first term on the right hand side of Eq.(1.15) is a free streaming term, and therefore
is a transient, as will appear clearly further on. The second term is related to the actual
coherent wave.

Equation (1.15) can now be used for computing hE f1 ix :

hE f1 ix = E0 (t)f10 (v) hsin(k0 x − ω0 t) cos[k0 (x − vt)]ix

Z t
e ∂f0
+ E0 (t) dt0 E0 (t0 ) hsin(k0 x − ω0 t) sin [k0 x + k0 v(t0 − t) − ω0 t0 ]ix
m ∂v 0
1
= E0 (t)f10 (v) sin[(k0 v − ω0 )t]
2
Z t
1 e ∂f0
+ E0 (t) dt0 E0 (t0 ) cos [(k0 v − ω0 )(t0 − t)] , (1.16)
2 m ∂v 0

10
having applied the relation hsin(k0 x + α) sin(k0 x + β)ix = (1/2) cos(α − β).

The time integral in the second term of this last relation can be carried out after Taylor
expanding the electrostatic field amplitude E(x, t0 ) = E(x, t) + (dE(t)/dt)(t0 − t) + . . .,
which is justified by invoking the assumption |(dE0 /dt)/E0 |  ω0 , so that
Z t
dt0 E0 (t0 ) cos [(k0 v − ω0 )(t0 − t)]
0
Z t
0 dE0 (t) Z t 0 0
0
= E0 (t) dt cos [(k0 v − ω0 )(t − t)] + dt (t − t) cos [(k0 v − ω0 )(t0 − t)]
0 dt 0
sin[(k0 v − ω0 )t] dE0 (t) ∂ 1 − cos[(k0 v − ω0 )t]
= E0 (t) + , (1.17)
k0 v − ω 0 dt ∂ω0 k0 v − ω 0
R R
having used dτ τ cos[(k0 v − ω0 )τ ] = −(∂/∂ω0 ) dτ sin[(k0 v − ω0 )τ ].

Combining (1.11), (1.16) and (1.17) finally provides:

Z
d Kin e
= − E0 (t) dv v f10 (v) sin[(k0 v − ω0 )t]
dt 2
e2 2
Z
∂f0 sin[(k0 v − ω0 )t]
− E0 (t) dv v
2m ∂v k0 v − ω 0
2 2 Z
e dE0 (t) ∂ ∂f0 1 − cos[(k0 v − ω0 )t]
− dv v .
4m dt ∂ω0 ∂v k0 v − ω 0
In this last relation, one notes that the first term phase mixes to zero for times t >
1/(k0 vth ), where the thermal velocity vth is the typical variation scale of the distribution
in velocity. In the same time limit, one also has sin(Ωt)/Ω → πδ(Ω) and [1−cos(Ωt)]/Ω →
P/Ω, where P stands for principal value, so that

d Kin e2 2
Z
∂f0 e2 dE02 (t) ∂
Z
v ∂f∂v
0

= − E0 (t) dv v πδ(k0 v − ω0 ) − P dv
dt 2m ∂v 4m dt ∂ω0 k0 v − ω 0
 
2 2 ∂f 0 /N
 0 ωp 2 0 dE02 (t) ∂  ωp
Z
∂f0 /N ∂v
= −π E (t) v − ω0 2 P dv , (1.18)
2 k0 0 ∂v v=ω/k0 4 dt ∂ω0 k0 v − ω0 /k0
| {z }| {z }
resonant bulk
where ωp2 = N e2 /(m0 ) is the squared plasma frequency, and having used δ(k0 v − ω0 ) =
(1/k0 )δ(v − ω0 /k0 ). The first term on the right hand side of Eq.(1.18) is clearly the contri-
bution from the resonant particles with velocities matching the phase velocity vφ = ω0 /ko ,
while the second term corresponds to the contribution from the bulk of the distribution.
In a sinusoidal wave with fixed amplitude, the bulk particles simply oscillate back and
forth in the electrostatic field, and experience no secular gain or loss in energy. The
resonant particles however experience a nearly constant field, and so can be efficiently
accelerated or decelerated.

11
 f 0(v)

Non−Resonant

Resonant

v
v =ω /k
φ 0 0

Figure 1.2: Separation of the distribution into a non-resonant and resonant part.

Invoking conservation of the total energy Etot , the kinetic energy gained/lost by the
particles must be compensated by the loss/gain of potential energy:
dEtot d Kin dPot
= + = 0, (1.19)
dt dt dt
where the space averaged potential energy Pot, i.e. the electrostatic energy, is given by
0 0 0
Pot = h E 2 (x, t)ix = E02 (t) hsin2 (k0 x − ω0 t)ix = E02 (t). (1.20)
2 2 4
Combining (1.18), (1.19) and (1.20) leads to
  
∂f0 /N
0 dE02 (t) ∂   ω2 Z   ω2 ∂f0 /N
ω0 1 − 2p P dv ∂v  = π 0 p E 2 (t) v . (1.21)
0
4 dt ∂ω0  k0 v − ω0 /k0  2 k0 ∂v v=vφ

On the right hand side of Eq.(1.21) one finds again the variation of kinetic energy of the
resonant particles Kinres , while on the left hand side one identifies the time variation of
the wave energy Ewave (= bulk kinetic energy + Pot), given by the general relation:
1 ∂
Ewave = Kinbulk + Pot = 0 [ω R (ω)] E02 (t), (1.22)
4 ∂ω ω0

where R is the real, i.e. non-resonant, part of the media’s dielectric function. Thus
Eq.(1.21) reads
dEwave dKinres
=− . (1.23)
dt dt

For the Langmuir wave model considered here, with mobile electrons and fixed ions,
one indeed has:
∂f0 /N
ωp2 Z ∂v
R (k, ω) = 1 − 2 P dv . (1.24)
k v − ω/k

12
 In the limit vφ /vth ' 1/k0 λD  1, one can consider the cold fluid approximation
R = 1 − ωp2 /ω 2 of Eq.(1.24), so that ∂(ω R )/∂ω|ω0 ' 2, as ω0 ' ωp .

Assuming an exponential decay of the wave:

d E02 (t)
E0 (t) = E0 e−γL t =⇒ = −2γL E02 (t),
dt
one then recovers from Eq.(1.21) the well-known relation for the linear Landau damping
(in the resonant approximation):
π ωp2 ∂f0 /N
γL = − v . (1.25)
2 k0 ∂v v=vφ
√
In the case of a Maxwellian distribution f0 (v) = N/( 2πvth ) exp −v 2 /(2vth
2
), Eq.(1.25)
becomes r  
γL π 1 ω0 1 ω0 2
= exp − . (1.26)
ω0 8 (k0 λD )3 ωp 2 k0 vth
Naturally, Eq.(1.25) leads to growth, i.e. instability, in case ∂f0 /∂v|v=vφ > 0. This is the
bump on tail instability.

Note that to re-derive the linear Landau damping relation (1.25) invoking energy con-
servation, one had to consider perturbation terms of the distribution up to second order,
i.e. f2 . This is due to the fact that energy is intrinsically a non-linear quantity, as in par-
ticular the potential energy Pot is quadratic in the perturbation field, as appears clearly
in Eq.(1.20).

Exercises:
1.2.2.1 Re-derive relation (1.25) for the linear Landau damping, starting from the linearized
Vlasov-Poisson equations, computing the dispersion function (k, ω), establishing
the appropriate dispersion relation, and solving using the resonant approximation.
1.2.2.2 Derive the general relation (1.22) for the wave energy of an electrostatic wave. In
fact, Eq. (1.22) can be further generalized to the case of an electromagnetic wave,
for which:
( )
1 ~? + 1 H
~ 0 · ∂ [ω H (ω)] · E ~ 0 · ∂ [ω µH (ω)] · H
~? ,
Ewave = 0 E 0 0
4 ∂ω ω0
µ0 ∂ω ω0

1.2.3 Limit of Linear Landau Damping

The derivation of the linear Landau damping presented in the previous section naturally
breaks under conditions for which one reaches the limits of the considered perturbative
approach. This is the case when ∂f1 /∂v becomes of the same order as ∂f0 /∂v:
∂f1 ∂f0
Linear Landau damping derivation breaks down when ∼ .
∂v ∂v
13
To estimate under which conditions this linear limit is met, one considers equation (1.15)
for f1 . For our purpose here, one can neglect the time dependence of the field envelope,
and thus directly carry out the remaining time integral in (1.15) to obtain:
e ∂f0 cos[k0 (x − vt)] − cos(k0 x − ω0 t)
f1 (x, v, t) = f10 (v) cos[k0 (x − vt)] + E0 . (1.27)
m ∂v k0 v − ω 0
Away from resonance (Ω = k0 v − ω0 6= 0), the derivative ∂f1 /∂v obviously produces
secular terms, i.e. which grow linearly in time t. Near resonance, both the numerator and
denominator in the second term on the right hand side of (1.27) go to zero, and so one
must expand the numerator for small Ω = k0 v − ω0 to address the variation of f1 in this
region. Noting that near resonance one has
Ω
z }| {
cos[k0 (x − vt)] = cos[k0 x − ω0 t − (k0 v − ω0 ) t]
= cos(k0 x − ω0 t) + sin(k0 x − ω0 t) (k0 v − ω0 )t
1
− cos(k0 x − ω0 t) (k0 v − ω0 )2 t2 + . . . ,
2
equation (1.27) can be written in this region as
 
e ∂f0 1
f1 (x, v, t) = f10 (v) cos[k0 (x−vt)]+ E0 t sin(k0 x − ω0 t) − (k0 v − ω0 )t2 cos(k0 x − ω0 t) + . . . .
m ∂v 2
The term in t2 rapidly becomes dominant when estimating ∂f1 /∂v:
∂f1 1 e ∂f0
'− k0 E0 t2 cos(k0 x − ω0 t),
∂v 2 m ∂v
so that  1/2
∂f1 ∂f0 m 1
At resonance: ∼ ⇐⇒ t> = , (1.28)
∂v ∂v ek0 E0 ωb
where ωb2 = ek0 E0 /m. The frequency ωb is identified as the bounce frequency of electrons
deeply trapped in the potential wells of the wave. It is indeed important to note, that the
electrons moving in the field E = E0 sin(k0 x − ω0 t) are separated into two groups: Passing
and trapped. This is discussed in more detail, in appendix A, where useful quantities such
as ωb are derived.

In the bulk of the distribution, i.e. away from resonance, one has k0 v − ω0 ' −ω0
and from Eq.(1.27):
∂f1 e ∂f0 sin[k0 (x − vt)]
'− k 0 E0 t ,
∂v m ∂v ω0
so that
∂f1 ∂f0 mω0 ω0
In bulk: ∼ ⇐⇒ t > = 2, (1.29)
∂v ∂v ek0 E0 ωb
Assuming ω0  ωb , one sees from Eqs. (1.28) and (1.29) that the linear approximation
breaks down much later in the bulk of the distribution than in the resonant region.

14
 For the linear result, and in particular the handling of resonant particles, to be valid
for describing the full evolution of the damping of the waves thus requires

γL > ωb . (1.30)

Under this condition, the wave has indeed damped out before the linear approximation
breaks down. For a given value of k0 λD , this last condition is equivalent to an upper
limit on the amplitude of the wave, which can be written in the case of a Maxwellian
distribution as:
s r " #
ωb δN γL π 1 1 1
= < ' exp − , (1.31)
ωp N ωp 8 (k0 λD )3 2 (k0 λD )2

where δN is the density perturbation amplitude of the Langmuir wave, and having used
relation (1.26) for γL in the case of a Maxwellian plasma.

Exercises:

1.2.3.1 Show that for Langmuir wavesq one has the bounce frequency of deeply trapped
electrons verifying ωb /ωp = δN/N . Show also that the trapping width ∆vtrap is
q
such that ∆vtrap = 2 ωb /k0 ' 2 δN/N vφ .

1.2.3.2 Verify the derivations in appendix A.

1.2.3.3 Draw the parallel between charged particles trapped in a sinusoidal electrostatic
field, and particles trapped in the magnetic well of a large aspect ratio tokamak.

1.2.3.4 Show that the condition described by Eq. (1.30) is equivalent to imposing

∆Kin > Ewave ,

where ∆Kin is the variation in kinetic energy which would result from flattening the
electron distribution in the resonant region, and Ewave is the average wave energy
of the Langmuir wave.

1.2.4 “Non-Linear Landau Damping”

In the previous section it was shown that the linear derivation of Landau damping is valid
for describing the full collisionless attenuation of the wave in the limit ωb  γL  ω0 .
We shall now consider the case ωb2 /ω0  γL  ωb  ω0 . In this limit the trapped
electrons have time to bounce back and forth many times in the potential wells of the
electrostatic field E(x, t) before the amplitude of the field is damped significantly. In this
case, at least for the resonant particles, one must correctly account for the full non-linear
trajectories of the electrons in the sinusoidal field E = E0 sin(k0 x − ω0 t). Note however,
that according to Eq.(1.29), the scaling ωb2 /ω0  γL ensures that the linear response of

15
the bulk distribution remains valid over the characteristic time scale of damping. Thus,
the decomposition illustrated by Fig. 1.2 is still preserved, i.e. of a bulk, non-resonant
distribution supporting, through a linear response, the oscillatory motion of the plasma
wave (= dispersion), and of a relatively small fraction of resonant particles leading to
damping/growth. As will appear clearly in the following, let us point out already that
the effect of resonant particles include significant contributions from both the trapped
electrons, as well as passing electrons near the separatrix shown in Fig. A.1.

The basic procedure for handling the non-linear evolution of the wave’s amplitude is
essentially the same as for the linear regime addressed in Sec. 1.2.2. Let us summarize.
One again invokes total energy conservation by equating the rate of change of the wave
energy Ewave to the variation of the kinetic energy of resonant particles Kinres , as already
written in Eq.(1.23):
dEwave dKinres
=− .
dt dt
As just pointed out, the bulk of the distribution may still be assumed to respond linearly
under the considered scaling, so that relation (1.22) for the wave energy still holds:

1 ∂
Ewave = 0 [ω R (ω)] E02 (t).
4 ∂ω ω0

Recalling that for Langmuir waves one has ∂(ω R )/∂ω|ω0 ' 2, and allowing for a time
dependant damping/growth rate γ(t) for the wave amplitude:
Rt
dt0 γ(t0 )
E0 (t) = E0 e− 0 ,

the rate of change of wave energy is thus again given by

dEwave
= −γ(t)0 E0 (t)2 . (1.32)
dt
Equation (1.8) for the variation of the kinetic energy of resonant particles is naturally still
valid here:
d Kinres
Z
m ∂f
= dv v 2 h ix . (1.33)
dt 2 res ∂t
Thus from the above relations, the time dependant rate γ(t) is computed from
1 d Kinres
Z
m ∂f
γ(t) = 2
= dv v 2 h ix . (1.34)
0 E0 dt 20 E02 res ∂t
The non-linear calculation differs however from the linear case in the way the distribution
f res of resonant particles is computed. Indeed, here the distribution cannot be derived
applying a perturbative approach as in Sec. 1.2.2, but must be calculated directly from
the non-linear Vlasov equation:
df ∂f ∂f e ∂f
= +v − E = 0, (1.35)
dt n.l.t.
∂t ∂x m ∂v

16
Figure 1.3: Evolution in resonant region of a Maxwellian distribution interacting with a sinu-
soidal wave. Results as seen from the wave frame. Black lines represent orbits of particles. Color
coding reflects amplitude of distribution. Note how the density in phase space is preserved along
the trajectories, but how the difference in the bounce/transit period between neighboring orbits
leads to a filamentation of the distribution.

Damping rate of a finite amplitude wave in limit ω >> γ

B L
1

0.8
0.6
0.4

0.2
γ (t) / γL

−0.2
a b c d e f
−0.4
−0.6
−0.8
Results from Matlab script
~/TEX/COURS_3EME_CYCLE/NonLinear/MODEL/NonLinLandau_ONeil_V2.m
−1
0 0.5 1 1.5 2 2.5 3 3.5 4
Time t ωB/ 2*π
)

Figure 1.4: Instantaneous rate of change γ(t) of wave amplitude as a function of time, as
computed from Eq. (1.34). The rate γ(t) is normalized with respect to the linear Landau
damping rate γL , while time is given in units of the deeply trapped bounce period τ b = 2π/ωb .
The rate γ(t) reflects the change in kinetic energy of the resonant particles, as clearly illustrated
by identifying the states of the distribution shown in Fig. 1.3 to the times a-f pointed out in
this graph: At times when a majority of the resonant electrons are being accelerated (resp.
decelerated), which corresponds to an increase (resp. decrease) in kinetic energy, one indeed
observes positive damping (resp. negative damping = growth) of the wave.

17
Figure 1.5: State of the distribution after tω b = 50, clearly illustrating the increased filamen-
tation in time of the distribution. At this stage, on average, the resonant particles are neither
accelerated nor decelerated, and the rate of change γ(t) of the wave amplitude becomes zero.

where d/dt|n.l.t. stands for the total time derivative along the non-linear trajectories, i.e.
along the full characteristics in the electrostatic field E.

The Vlasov equation (1.35) can always be solved formally by integrating along the char-
acteristics (see Appendix B):

f (x, v, t) = f [x0 (0; x, v, t), v 0 (0; x, v, t), 0], (1.36)

where f (x, v, 0) is the initial distribution, and [x0 (t0 ; x, v, t), v 0 (t0 ; x, v, t)] the non-linear
trajectories verifying:
dx0 dv 0 e
0
= v0, 0
= − E0 sin(k0 x0 − ω0 t0 ), (1.37)
dt dt m
with initial conditions x0 (t) = x, v 0 (t) = v. (1.38)

Note that one makes use of the time scale separation ωb  γL by considering the ampli-
tude E0 of the field fixed when integrating the non-linear trajectories (1.37)-(1.38), which
are then used for computing (1.36) and (1.34). The time dependence of E0 (t) is then
taken account for iteratively when computing dEwave /dt through Eq. (1.32).

The above system of equations can easily be solved numerically, which provides a useful
illustration of the mechanism underlying the non-linear evolution of the wave. This has
been done to obtain the results presented in Figs. 1.3, 1.4 and 1.5.

Figures 1.3 and 1.4 show how the wave amplitude decreases [γ(t) > 0] as a majority
of the resonant particles are accelerated, which corresponds to an increase in kinetic en-
ergy. Inversely, the wave amplitude increases [γ(t) < 0] as a majority of the resonant

18
particles are decelerated, which corresponds to a decrease in kinetic energy. This is obvi-
ously the origin of the oscillations in kinetic and potential energy already observed in the
full simulation results shown in Fig. 1.1.

According to the Vlasov equation, the density in phase space remains invariant along
the trajectories of the particles. But as a result of the fact that neighboring orbits of
passing particles (resp. trapped particles) have different transit times (resp. bounce pe-
riods), as shown in Fig. A.2, one observes a filamentation of the distribution over time.
This is clearly illustrated in Fig. 1.5, which shows the distribution as computed from
Eq.(1.36) at time tωb = 50. At this stage, the resonant particles are neither accelerated
nor decelerated on average, and the rate of change γ(t) ofR the wave amplitude therefore
tends to zero. As a result, the attenuation factor exp[− 0∞ dtγ(t)] is non-zero, so that
asymptotically in time one has an undamped mode. This resulting wave is a BGK mode.
BGK modes are discussed in some detail in Sec. 1.2.5. As filamentation happens over the
time scale of a few bounce periods ωb , this phenomena is obviously only observed under
the assumed scaling ωb  γ. Note the similarity between Figs. 1.3 & 1.5 and the phase
space plots from the full simulations shown in Fig. 1.1.

The above system of equations (1.34), (1.36)-(1.38) can in fact also be solved analytically.
The corresponding derivation, which is a somewhat lengthy exercise involving Jacobian
elliptic functions [5], is described to some detail in Ref. [6]. From this derivation, the
final result for the damping rate is
h i h i
nπ ωb t (2n+1)π ωb t
∞
X 64
Z 1 2nπ 2 sin κF
(2n + 1)π 2 κ sin 2F
γ(t) = γL dκ { + },
n=0 π 0 κ5 F 2 (1 + q 2n )(1 + q −2n ) F 2 (1 + q 2n+1 )(1 + q −2n−1 )
| {z } | {z }
passing trapped
(1.39)
where q = exp(πF 0 /F ), √F = F (κ2 ) is the complete elliptic integral of the Rfirst kind [5],
and F 0 = F 0 (K 2 ) = F ( 1 − K 2 ). Equation (1.39) includes the integral dκ over the
energy variable κ, defined by (A.3), characterizing the different orbits of the resonant
particles in phase space. The first term in the integrand corresponds to the contribution
from the resonant passing particles, while the second term is related to trapped parti-
P
cles. Notice also the sum n over harmonics of the transit period τt (κ) = 2κF/ωb (resp.
bounce period τb (κ) = 4F/ωb ) of passing (resp. trapped) particles. These relations for τt
and τb for arbitrary energy levels κ are derived in Appendix A.

In the first stage of the evolution, that is in the limit t → 0, one can show that Eq.
(1.39) indeed recovers the linear damping rate γ(t) → γL , as illustrated in Fig. 1.4.
One needs to be careful however in taking this limit, as this is a typical case where
R R
limt→0 dκ 6= dκ limt→0 .

Due to the dependence in κ of τt and τb , as shown in Fig. A.2, one can see that the
integrals in κ over the terms sin[nπ ωb t/κF ] and sin[(2n + 1)π ωb t/2F ] phase mix to zero
as t becomes large, and thus one can show that limt→∞ γ(t) = 0. This can also be seen in

19
 2.5
untrapped
2

1.5
δv
1

0.5

0
trapped

V in units ω / k
W
W+δW
b
0

−0.5

−1

−1.5

−2

−2.5
−0.5 0 0.5
X in units λ0

Figure 1.6: Initial sub-volumes in phase space within the region limited by orbits with energy
levels W and W = δW for both passing and trapped particles. In time, due to the difference in
periods of the neighboring trajectories, these sub-volumes are stretched into ever finer filaments
that ultimately uniformly fill the full volume between the two orbits.

Fig. 1.4.

The full derivation of Eq.(1.39) is not presented

R∞
here. However, a direct calculation
of the time integrated damping exponent 0 dt γ(t), which is perhaps one of the most
useful results, is now carried out. From Eq.(1.34) one can write:
Z ∞ 1
dt γ(t) = ∆Kinres , (1.40)
0 0 E02
where
m Z λ0 /2 res
Z
∆Kin = dX dV (V + vφ )2 (f∞ − f0 ), (1.41)
2λ0 −λ0 /2 res

f0 = f (X, V, t = 0) and f∞ = f (X, V, t = ∞) being respectively the initial and time

asymptotic distribution in wave frame variables (X = x − vφ t, V = v − vφ ). For conve-
nience, the following derivation is indeed carried out in wave frame variables. However, to
lighten notations one reverts to the notation (x, v) for position and velocity in the wave
frame.

The main issue here is to evaluate f∞ (x, v). For this purpose, let us analyze in some-
what more detail the mechanism of filamentation. Figure 1.6 shows in the wave frame
a small compact sub-volume in phase space within the region limited by orbits with en-
ergy levels W and W + δW , both for the case of trapped, as well as forward/backward
passing particles. In all cases, the difference in periods of the trajectories with energies
between W and W + δW lead to a stretching over time of the sub-volume into an ever
finer filament that ultimately evenly fills the whole region between the two orbits. As a
result of the incompressibility of phase space, the density f remains invariant between
the initial and final stretched state of the sub-volume. This stretching obviously happens

20
for all the other initially compact sub-volumes within the orbits W and W + δW , so
that asymptotically in time all these sub-volumes are finely intertwined with each other.
This effect is clearly illustrated in Fig. 1.5. The coarse grain distribution fˆ obtained by
averaging the actual distribution f over phase elements sufficiently large to be crossed by
many of these filaments thus becomes uniform within the orbits W and W + δW and cor-
responds to the average of the initial distribution f0 within these two orbits. Obviously,
the size of the coarse-graining for computing fˆ can be taken to zero asymptotically in time.

Note, that in a real physical system, even with very low collisionality, these ever finer
phase space structures in fact always end up getting smeared out in a finite time through
collisional diffusion. Indeed, for a given collisionality νc , a structure in the distribution
with scale λv in velocity diffuses in a characteristic time tc ∼ νc (vth /λv)2 .

On the basis of the above arguments, the distribution for untrapped particles [W >
max(−eφ)] is obtained asymptotically in time from
R λ0 /2
−λ0 /2 dx δv(x, W ) f0 [σv(x, W )]
u
f∞ (W, σ) = fˆ∞ = R λ0 /2 ,
−λ0 /2 dx δv(x, W )

where the initial distribution f0 is in fact only function of velocity v. In the above
relation, σ = sign(v), v(x, W ) is the velocity at point x of the trajectory with energy W ,
by convention chosen positive, and δv(x, W ) the separation in velocity between orbits W
and W + δW . Note that for untrapped particles each energy level W corresponds to both
a forward passing (σ = +1) as well as a backward passing (σ = −1) trajectory, which
must in general be treated separately. From energy conservation for a single particle, one
has
 1/2 
1 2 2
mv − eφ(x) = W =⇒ v(x, W ) = (W + eφ)
2 m
∂v δW
=⇒ δv(x, W ) = δW = .
∂W m v(x, W )
The asymptotic untrapped distribution thus can be written
R λ0 /2
u −λ0 /2 f0 [σv(x, W )] dx/v(x, W )
f∞ (W, σ) = R λ0 /2 . (1.42)
−λ0 /2 dx/v(x, W )

For trapped particles [min(−eφ) < W < max(−eφ)] each energy level W represents a
single orbit, which includes however both a forward as well as a backward going segment,
so that for this group of particles the asymptotic distribution is computed from
P R x2
f0 [σv(x, W )] dx/v(x, W )
= fˆ∞ =
t σ=±1 x1
f∞ (W ) R x2 . (1.43)
2 x1 dx/v(x, W )
where x1,2 are the turning points, i.e. where v(x1,2 , W ) = 0.

21
As the time asymptotic distributions (1.42) and (1.43) are to be used for computing the
variation in kinetic energy (1.41) of the resonant particles, the distribution f0 in (1.42)
and (1.43) can be Taylor-expanded around the wave phase velocity, which in the wave
frame is v = 0. This naturally assumes that the resonant region, of the order ∆vtrap , is
such that ∆vtrap /vth  1. For the purpose of computing ∆Kinres , the leading order effect
is obtained by considering a Taylor expansion to first order:

df0 (0)
f0 (v) = f0 (0) + v + .... (1.44)
dv
Inserting (1.44) in (1.42) thus leads to
Rλ /2
u df0 (0) σ −λ0 0 /2 dx df0 (0) σλ0
f∞ (W, σ) = f0 (0) + R λ0 /2 + . . . = f0 (0) + + ...,
dv −λ /2 dx/v(x, W ) dv τt
0

where τt stands for the transit time. Restricting the derivation to sinusoidal waves, one
makes use of Eq. (A.4) to obtain:

u df0 (0) π ∆vtrap

f∞ (κ, σ) = f0 (0) + σ + ...,
dv 2 κ F (κ2 )

so that together with (1.44):

!
df0 (0) π ∆vtrap
∆f u = f∞
u
− f0u = σ − v + .... (1.45)
dv 2 κ F (κ2 )

For trapped particles, the term in df0 /dv obviously cancels out:
R
df0 (0) σ=±1 σ xx12 dx
P
t
f∞ (W ) = f0 (0) + R + . . . = f0 (0) + . . . ,
dv 2 xx12 dx/v(x, W )

so that
df0 (0)
∆f t = f∞
t
− f0t = − v + .... (1.46)
dv
To the considered order, ∆f = f∞ −f0 is clearly odd in v, so that Eq.(1.41) can be written
Z Z
res 2m vφ λ0 /2 ∞
∆Kin = dx dv v ∆f. (1.47)
λ0 −λ0 /2 0

Note, that although one is only computing the resonant particle contribution to dKin/dt,
the integration in velocity has been prolonged here to infinity. This is justified by proving
that the relation ∆f u for passing particles derived in Eq. (1.45) actually goes to zero for
v → ∞ (κ → 0), and in this way the velocity integral in (1.47) finds its own cutoff (check
it!).

22
 2
By making the change of variables v ↔ κ, so that from Eq.(A.2) dv = ∆vtrap /(vκ3 ) dκ,
one can write:
 
Z Z Z Z 
 1 dκ 1 λ0 /2 ∞ dκ 1 x2 
∆Kinres = 2m vφ ∆vtrap
2  dx ∆f u + dx ∆f t 
, (1.48)

 0 κ3 λ 0 −λ0 /2 1 κ3 λ 0 x1
| {z } | {z }
passing trapped

where the first term in (1.48) corresponds to the contribution from untrapped electrons,
and the second term to trapped. Expliciting the space integral for passing particles leads
to
Z " #
1 λ0 /2
u df0 (0) π ∆vtrap
dx ∆f = − < v >x
λ0 −λ0 /2 dv 2 κ F (κ2 )
" #
df0 (0) ∆vtrap π 2 E(κ2 )
= − , (1.49)
dv κ 2 F (κ2 ) π

having made use of Eq. (1.45) for ∆f u and Eq. (A.5) for the spatial averaged velocity
< v >x . The corresponding integral for trapped particles becomes:
Z
1 x2 df0 (0)
dx ∆f t = − < v >x
λ0 x1 dv
 
df0 (0) 2∆vtrap 1 1 1
= − E( 2 ) + ( 2 − 1)F ( 2 ) , (1.50)
dv π κ κ κ
having made use of Eq. (1.46) for ∆f t and Eq. (A.7) for < v >x .

Inserting Eqs. (1.49) and (1.50) into (1.48) then provides

(Z " # )
E(κ2 )
Z 
df0 (0) 1 dκ π ∞ dκ 1 1 1 1
∆Kinres 3
= 4m vφ ∆vtrap 4 2
− − 3
E( 2 ) + ( 2 − 1)F ( 2 )
dv 0 κ 4 F (κ ) π 1 κ π κ κ κ
( " # )
Z 1 2
γL 64 1 E(κ ) π κ h i
= 0 E02 dκ 4 − 2
+ E(κ2 ) + (κ2 − 1)F (κ2 ) , (1.51)
ωb π 0 κ π 4 F (κ ) π

having made the change of variable κ → 1/κ for the trapped contribution, and made use
of Eq. (1.25) for the linear Landau damping γL .

Finally, inserting (1.51) in (1.40) leads to:

 

 

 " # 
1 E(κ2 )
Z Z  i
∞ γL 64 1  π κh 
dt γ(t) = dκ 4
− 2)
+ E(κ2 ) + (κ2 − 1)F (κ2 ) .
0 ωb π 0 

 κ π 4 F (κ π
| {z

}


| {z } 

passing trapped

(1.52)
Note how this result scales as γL /ωb , which is assumed small in the considered ordering.
Recall that for a Maxwellian, γL /ωp is essentially a function of kλD [see Eq.(1.26)], while

23
 q
ωb /ωp = δN/N .

In the constant factor on the right hand side of Eq. (1.52), one can still distinguish the
contribution of untrapped (first term) from the one of trapped particles
R
(second term).
By numerical integration, this constant can be estimated: (64/π) 01 dκ. . . ' 1.96, with
the contribution from untrapped particles being 0.52 and from trapped 1.44.

Exercises:
1.2.4.1 Recover (1.52) from (1.39) making use of the following relations:
∞
2π X 1 E(κ2 ) π
2 2n
= − ,
−2n
F (κ ) n=1 (1 + q )(1 + q ) π 4 F (κ2 )
∞
2π 2 X
 
1 1 1
2 2 2n−1
= 1− F (κ2 ) + E(κ2 ).
κ F (κ ) n=1 (1 + q )(1 + q −2n+1 ) κ2 κ2

1.2.4.2 Compare the “exact” non-linear simulation results shown in Fig. 1.1 to the analytic
relations (1.39) and (1.52). Is there qualitative agreement? Quantitative agreement?
Recall that the simulation results in figure 1.1 are for standing Langmuir waves.

1.2.5 BGK waves

From the previous section, where the regime γL  ωb has been studied, it appears that
the non-linear evolution of a Langmuir wave leads, asymptotically in time, towards a finite
amplitude, undamped mode. Let us now see in more detail how such a state is charac-
terized. In fact, one wants to further convince oneself that such an undamped mode can
actually exist as an exact self-consistent solution of the non-linear Vlasov-Poisson system.

In the wave frame, an undamped mode corresponds to a stationary solution of the Vlasov-
Poisson system. Such a state is thus characterized by a time independent distribution
f = f (x, v) for each species, and a time independent potential field φ = φ(x). The dis-
tribution f for each species, with charge q and mass m, is solution to the corresponding
stationary Vlasov equation:
" #
∂ q dφ(x) ∂
v − f (x, v) = 0.
∂x m dx ∂v
As a result, f must be a function of the invariants of motion. For particles trapped
in the troughs of the potential φ, the only conserved quantity is the total energy W =
mv 2 /2 + qφ(x), while for untrapped particles the sign σ of the velocity is an additional
invariant. The stationary distribution must therefore be of the form:

f (x, v) = f (W, σ) = f t (W ) + f u (W, σ), (1.53)

where f t (W ) is the distribution of trapped particles, non-zero for min(qφ) < W <
max(qφ), and f u (W, σ) is the distribution of untrapped particles, non-zero for W >

24
max(qφ).

The potential φ(x) is solution of Poisson’s equation:

d2 φ(x) 1 X Z +∞
= − q dv f (x, v)
dx2 0 species −∞
1 X Z +∞
= − q dv f (W = mv 2 /2 + q φ(x), σ)
0 species −∞
Z ∞ P
1 X dW σ=±1 f (W, σ)
= − q . (1.54)
0 species qφ [2m(W − qφ)]1/2

Thus, in principal, given any distribution of the form (1.53), Eq. (1.54) defines a second
order ordinary differential equation (ODE) for φ(x). The so-obtained set [{f }species , φ]
corresponds to a non-linear, undamped state. Such undamped states are the so-called
Bernstein-Greene-Kruskal (BGK) modes [7].
P
Note from Eq. (1.54), that the potential φ(x) in fact depends only on σ=±1 f (W, σ). As
a result, φ(x) remains invariant if one varies the partition between forward and backward
P
passing particles while keeping σ=±1 f (W, σ) unchanged. Thus, for a given φ one can
arbitrarily vary the pase velocity in the lab frame (defined as the frame in which the
system has zero average momentum).

Equation (1.54) can be written as

d2 φ(x)
= F (φ), (1.55)
dx2
where F (φ) is the right hand side of (1.54), and can be viewed as the “force” acting on
a “particle” with “position” φ and evolving in “time” x. By multiplying Eq.(1.55) by
dφ/dx and integrating in x, one thus obtains

1 dφ 2
( ) + V (φ) = const., (1.56)
2 dx
where the “potential” V (φ) is such that dV (φ)/dφ = −F (φ), and is given by
Z  1/2 X
1 X ∞ 2
V (φ) = − dW (W − qφ) f (W, σ).
0 species qφ m σ=±1

Obviously, if V (φ) has a minimum, one can obtain periodic solutions for φ(x).

Equation (1.56) can be solved by quadrature:

Z φ dφ
± = x − x0
φ0 [2 (V (φ0 ) − V (φ))]1/2

25
where the “initial conditions” have been chosen such that at position x = x0 one has
φ(x0 ) = φ0 and dφ(x0 )/dx = 0.

Inversely, one can demonstrate [7] that from any given potential field φ(x) and given dis-
tributions for the various species (electrons and ions) and groups (passing and trapped),
except one (e.g. the distribution fet (W ) of trapped electrons), one can build a BGK-type
mode by solving for this undetermined distribution.

Thus, it obviously appears that one can build BGK modes of quite arbitrary shape of
the potential φ(x) and in particular wavelength. Furthermore, as already pointed out,
for a given φ the wave velocity in the lab frame can also be chosen arbitrarily. Hence,
in general, the wave numbers and frequencies of BGK modes do not need to obey the
dispersion relations derived in the framework of linear theory. The larger this deviation
from the linear dispersion relation, the stronger the distribution of the BGK mode must
be deformed from the initial, unperturbed equilibrium distribution (usually a Maxwellian)
considered in linear theory. To reach a state which is far from equilibrium naturally re-
quires the system to be very strongly driven.

1.2.6 Non-Linear Frequency Shift

In Sec. 1.2.4, the non-linear evolution of the amplitude of a Langmuir wave was studied
in the regime ωb2 /ω0  γL  ωb  ω0 . It was shown how such a wave ultimately evolves
towards a finite amplitude, undamped mode, a so-called BGK mode. The characteristics
of these BGK modes were discussed in Sec. 1.2.5, and it was pointed out how in general
they do not necessarily obey the dispersion relations from linear theory. In the case of
interest in section 1.2.4, let us recall however that the assumed scaling ensures that the
main part of the distribution, the so-called bulk, still responds linearly. One therefore
expects in this case, that the frequency and wavelength of the final BGK state verify a
dispersion relation which corresponds to the linear one with at most a small correction
due to the minority fraction of non-linearly behaving resonant particles. The purpose of
this section is to derive this non-linear correction. Note that the evolution of the electron
plasma wave considered in 1.2.4 is an initial value problem, for which the fundamental
wavenumber k0 remains fixed. A non-linear correction to the dispersion relation will thus
affect the frequency, i.e. potentially leading to a non-linear frequency shift.

One starts from Poisson’s equation for the potential φ of the final, BGK-like state of
the electron plasma wave considered in Sec. 1.2.4, i.e. a case with dynamic electrons and
a fixed, neutralizing ion background:

d2 φ
 Z 
1 +∞
2
= e dv f∞ − qi Ni,0 , (1.57)
dx 0 −∞

26
where f∞ is the time asymptotic distribution given by Eqs. (1.42) and (1.43) for passing
and trapped particles respectively, qi is the charge of ions and Ni,0 their uniform density.

To explicit the purely non-linear response, one adds and subtracts in (1.57) the charge
density related to the linear response of the system:
d2 φ
 Z Z 
1 +∞ +∞
= e dv f L − q N
i i,0 + e dv (f∞ − fL )
dx2 0 −∞ −∞
Z +∞ Z +∞ 
e
= dv δf + dv (f∞ − fL ) , (1.58)
0 −∞ −∞

where fL = f0 + δf , fR0 being the initial unperturbed electron distribution verifying the
neutrality condition e dvf0 = qi Ni,0 , and δf the actual linear perturbation. For a wave
with wave number k and frequency ω in the lab frame, one has
e kφ ∂f0
δf = − , (1.59)
m kv − ω ∂v
which is obtained from the linearized Vlasov equation:
∂δf ∂δf e ∂φ ∂f0
+v + = 0.
∂t ∂x m ∂x ∂v
Inserting (1.59) in (1.58), and making use of the fact that the wave remains mainly
sinusoidal with wave number k0 , enables to write:
Z
e +∞
−k02 L (k0 , ω)φ = dv ∆fNL , (1.60)
0 −∞

having defined the non-linear deviation of the distribution ∆fNL = f∞ − fL . In equation

(1.60) L (k, ω) is the linear dielectric function
∂f0 /N
ωp2 Z ∂v
L (k, ω) = 1 − dv .
k2 v − ω/k

The non-linear charge density on the right hand side of Eq. (1.60) appears as an external
source term to a system responding linearly with dielectric L .
R λ0 /2
One now projects Eq. (1.60) onto φ, i.e. carries out hφ . . .ix = (1/λ0 ) −λ0 /2 dx φ . . .
on each side of the equation, to obtain:
Z Z
2 1 λ0 /2 +∞
L (k0 , ω) = − dx eφ dv ∆fNL
0 E02 λ0 −λ0 /2 −∞
Z Z +∞
2 1 λ0 /2 dW X
= − dx eφ ∆fNL , (1.61)
0 E02 λ0 −λ0 /2 −eφ mv(x, W ) σ=±1

having used hφ2 ix = E02 /(2k02 ). To obtain the last equality in (1.61) one has changed to
wave frame variables, and made a variables transformation from velocity v to wave frame

27
energy W = mv 2 /2 − eφ, so that v(x, W ) = +[(2/m)(W + eφ(x))]1/2 . As will clearly
appear in the following, this change of frame and variables is convenient for expliciting
∆fNL = f∞ − fL .

The distribution f∞ is given in wave frame variables by Eqs. (1.42) and (1.43) for passing
P
and trapped particles respectively. As appears from (1.61), it is in fact σ=±1 f∞ which is
required. This quantity can actually be represented for both groups of particles, passing
and trapped, by the single, compact relation:
P R λ0 /2
X σ −λ0 /2 dxf0 [σv(x, W )] H(W +eφ)
v(x,W )
f∞ (W, σ) = R λ0 /2
σ=±1 −λ0 /2 dx H(W +eφ)
v(x,W )

hf0 [σv(x, W )] H(W +eφ)

P
σ v(x,W ) x
i
= , (1.62)
h H(W +eφ)
v(x,W ) x
i

where H(w) is the Heaviside step function, defined as H(w) = 0 for w < 0, and H(w) = 1
for w > 0.

One expects the non-linear deviation ∆fNL to be essentially non-zero in the resonant
region, i.e. around v = 0 in the wave frame. Thus, as in section 1.2.4 for addressing the
non-linear damping, one can here again Taylor-expand f0 around v = 0. However, as a
result of the summing over σ = ±1, which cancels out the first order terms, the expansion
must be considered in this case to second order:
df0 (0) 1 d2 f0 (0) 2
f0 (v) = f0 (0) + v+ v + .... (1.63)
dv 2 dv 2
Inserting (1.63) in (1.62) provides:
X d2 f0 (0) hHv(x, W )ix d2 f0 (0) v̄
f∞ (W, σ) = 2f0 (0) + H + . . . ' 2f 0 (0) + , (1.64)
σ=±1 dv 2 h v(x,W )
i x dv 2 mv̄ 0

having used the notation v̄ = hHvix for the average velocity and v̄ 0 = dv̄/dW =
(1/m)hH/v(x, W )ix . Naturally, H is the shorter notation for H(W + eφ).

The linear distribution fL = f0 + δf is also expanded up to order d2 f0 /dv 2 , so that

using (1.59) in wave frame variables one obtains:
" #
X X e φ df0 (v)
fL = f0 (v) −
σ=±1 σ m v dv
( " #)
X df0 (0) 1 d2 f0 (0) 2 e φ df0 (0) d2 f0 (0)
' f0 (0) + v+ v − + v
σ dv 2 dv 2 mv dv dv 2
!
d2 f0 (0) 2 eφ d2 f0 (0) 2
= 2f0 (0) + v − 2 = 2f 0 (0) + W. (1.65)
dv 2 m dv 2 m

28
Combining (1.64) and (1.65) thus provides

d2 f0 (0) 2 d2 f0 (0)
 
X X v̄ 2 v̄
∆fNL = (f∞ − fL ) ' 2 0
− W =− 0 2
(W v̄ 0 − ). (1.66)
σ=±1 σ dv mv̄ m mv̄ dv 2
P P
Figure 1.7 shows the distributions σ f∞ and σ fL in the resonant region.
P
The above relation derived for σ ∆fNL is clearly only function of W , so that in Eq.(1.61)
one can invert the integrals over x and W as follows:
Z Z
2 +∞ X 1 λ0 /2 eφ H(W + eφ)
L (k0 , ω) = − 2
dW ∆fNL dx
0 E0 min(−eφ) σ λ0 −λ0 /2 m v(x, W )
Z +∞
2 v̄ X
= 2
dW (W v̄ 0 − ) ∆fNL , (1.67)
0 E0 min(−eφ) 2 σ

having used

−(W − mv 2 /2) H(W + eφ)

Z Z
1 λ0 /2 eφ H(W + eφ) 1 λ0 /2
dx = dx
λ0 −λ0 /2 mv(x, W ) λ0 −λ0 /2 mv
W H 1 v̄
= − h ix + hHvix = −(W v̄ 0 − ).
m v 2 2

Inserting (1.66) into (1.67) thus gives

4 d2 f0 (0) Z +∞ 1 v̄
L (k0 , ω) = − 2 2
dW 0 (W v̄ 0 − )2 . (1.68)
m0 E0 dv min(−eφ) v̄ 2

Assuming that the potential field φ(x) is still essentially sinusoidal, the average velocity
v̄ is given by Eqs. (A.5) and (A.7) for passing and trapped particles respectively. Fur-
thermore, for passing particles one has hH/v(x, W )ix = τt /λ0 where τt is the transit time
given by relation (A.4), while for trapped particles hH/v(x, W )ix = τb /(2λ0 ), where τb is
the bounce period given by (A.6). In terms of the energy variable κ, defined in Eq.(A.3),
these terms thus read for passing particles (0 < κ < 1):
2 ∆vtrap
v̄ = E(κ2 ),
π κ
2 κF (κ2 )
v̄ 0 = ,
π m∆vtrap

and for trapped particles (κ > 1):

 
2∆vtrap 1 1 1
v̄ = E( 2 ) + ( 2 − 1)F ( 2 ) ,
π κ κ κ
1
2 F ( κ2 )
v̄ 0 = .
π m∆vtrap

29
 By inserting these last relations into Eq. (1.68), one then finally obtains:
5
(Z
1 m∆vtrap d2 f0 1 dκ 1 h 2 2 2 2
i
L (k0 , ω) = − (2 − κ )F (κ ) − 2E(κ )
2π 0 E02 dv 2 vφ 0 κ6 F (κ2 )
Z  2 )
∞ dκ 1 1 1
+ 1 F ( ) − 2E( ) ,
1 κ3 F ( κ 2 ) κ2 κ2
 

 

ωp2 d2 (f0 /N )
Z  
8 1  1 h 2
i2 κ 
= − ∆v trap dκ 2(F − E) − κ F + (F − 2E)2 , (1.69)
k02 dv 2 vφ
π 0 


6
κ F } |F

}

| {z {z 
passing trapped

having used the shortened notations F = F (κ2 ) and E = E(κ2 ) for the final equal-
ity. To obtain (1.69) one has transformed integration variable W to κ, so that W =
2
(m∆vtrap /4)(2/κ2 − 1) and dW = (m∆vtrap 2
/κ3 )dκ. For the trapped particle contri-
bution, oneR has furthermore carried out the transformation κ → 1/κ. The integral
α = (8/π) dκ . . . in (1.69) is a constant, and can be integrated numerically, providing
the value α = 0.823, which is composed of the contribution αu = 0.117 from untrapped
particles and αt = 0.705 from trapped particles.

For a mode with given wave number k0 , the non-linear correction derived in Eq.(1.69)
to the linear dispersion relation L (k0 , ωL ) = 0, obviously leads to a shift δω of the lin-
ear frequency ωL . By assuming the non-linear effect small, so that one may expand
L (ωL + δω) = L (ωL ) + δω ∂L (ωL )/∂ω + . . . ' δω ∂L (ωL )/∂ω, one obtains:
α ωp2 d2 (f0 /N ) α ωp3 d2 (f0 /N )
δω = − ∂L (ωL ) 2 ∆vtrap =− ∆v trap ,
∂ω
k0 dv 2 vφ
2 k02 dv 2 vφ

having again made use of the cold fluid approximation L ' 1 − ωp2 /ω 2 for the dispersion
function of electron plasma waves. This frequency shift is shown for Langmuir waves in
Fig. 1.8 for both the initial value problem considered here (=”sudden”), as well as for
the case of a wave turned on adiabatically (=”adiabatic”. See exercise 1.2.6.1).

Exercises:
1.2.6.1 Derive the non-linear frequency shift in the case of a Langmuir wave turned on
adiabatically instead of “suddenly”, as in the initial value problem which has been
considered in this section.
R
Start by deriving the distribution f adiab. by invoking the
invariance of the action dxv as the amplitude of the wave is slowly turned on. This
distribution is to be compared to the asymptotic distribution f∞ considered in the
initial value case. This exercise illustrates the non-uniqueness of the trapped particle
distribution for a given wave amplitude, and thus reflects the potential diversity of
BGK modes.
1.2.6.2 Adapt the derivation in this section to handle the non-linear frequency shift of an
ion acoustic wave. Point out the contribution to the frequency shift from both the
ions and the electrons.

30
 Distributions in resonant region
1

[ Σ f(W,σ) − 2 f (v ) ] / ∆ v2 f ’’(v )
p
Linear
t

0.5 Sudden
Adiabatic
p
0

0
σ

−0.5
−1 −0.5 0 0.5 1 1.5 2
Energy W / (e φ )
0

Figure 1.7: Distributions in resonant region for the linear, “sudden”, and “adiabatic” cases.
For the non-linear distributions, a sinusoidal wave is assumed. The “sudden” distribution cor-
responds to the initial value case considered in the main text, while the “adiabatic” case corre-
sponds to the distribution obtained by adiabatically turning on the wave, and is addressed in
exercise 1.2.6.1.
Comparison with theoretical result by G. Morales & T. O’Neil, PRL 28, p.417 (1972)
and R. Dewar, Phys. Fluids 15, p.712 (1972)
0
Theory, Sudden
−0.005 Theory, Adiabatic
Numerical, Sudden
−0.01 Numerical, Driven

−0.015
p
Frequency Shift δω / ω

−0.02

−0.025

−0.03

−0.035

−0.04

−0.045

0 0.01 0.02 0.03 0.04 0.05 0.06

Density Amplitude δN / N

Figure 1.8: Non-linear frequency shift of a Langmuir wave as a function of the density pertur-
bation amplitude δN/N . The considered wave number is k 0 λD = 1/3. The numerical simulation
results from the SAPRISTI code are compared to the analytic predictions for both the “sudden”
and “adiabatic” case.

31
 Large Amplitude Wave λ0 = 2 π / k0

ωb Trapped
macro−particle

∆v
Sideband

λs = 2 π / ks

Figure 1.9: Resonant interaction of trapped particles with sidebands.

1.2.7 Stability of BGK Mode: The Trapped Particle Instability

In general, the undamped BGK modes are not stable equilibrium states. Indeed, they are
subject to be destabilized by the growth of sidebands, which may be resonantly driven
by the particles trapped in the principal wave. The resonant mechanism of this trapped
particle instability (also called modulational instability) is schematized in Fig. 1.9.

1.2.8 Further Reading

• Non-linear Landau damping and frequency shift: See O’Neil [6], Morales and
O’Neil [8], and Dewar [9].

• BGK waves: Bernstein, Greene, and Kruskal [7].

• Trapped Particle Instability: Kruer et. al [10], Goldman [11], and Dewar et. al
[12].

• Experimental evidence of non-linear evolution: Danielson et. al [3] and the

references therein.

32
Figure 1.10: Basic mechanism of a parametric instability that may affect a laser beam in a
plasma: The incident light (= wave #1), scatters off a plasma wave (= wave # 2), i.e. either
an electron plasma wave (EPW) or an ion acoustic wave (IAW). If the matching conditions are
met, the scattered light ( = wave # 3) and the incident light may beat together in such a way
as to reinforce the plasma wave.

1.3 Three Wave Interactions and Parametric Insta-

bilities
Parametric instabilities result from the resonant interaction between three non-linearly
coupled waves. The basic mechanism is schematized in figure 1.10 in the case of a laser
beam, i.e. a transverse electromagnetic wave, propagating through an unmagnetized,
under-dense plasma (frequency ω0 of incident light > plasma frequency ωp ). This in-
cident light will reflect off any electron density fluctuation, in particular perturbations
related to electron plasma waves (EPWs) or ion acoustic waves (IAWs). Under condi-
tions of appropriate phase matching, the scattered and incident light may beat together
in such a way as to reinforce the plasma wave ( = EPW or IAW). This reinforced plasma
wave will in turn lead to a higher level of scattering, and this increased scattered light will
lead to a stronger beating with the incident light, which will further amplify the plasma
wave, and so on, obviously initiating an instability. The case of an electromagnetic wave
scattering off an EPW is called Stimulated Raman Scattering (SRS), while the case of
an electromagnetic wave scattering off an IAW is called Stimulated Brillouin Scattering
(SBS).

The cases of SRS and SBS, as parametric instabilities affecting laser light propagating
through a plasma, are of particular concern in the context of inertial fusion, as they lead
to a loss of control of the laser energy deposition. Due to the relative simplicity in mod-
eling these particular parametric instabilities (mainly the fact that they involve waves in
a non-magnetized plasma), they will also be considered as illustrations in the following.
However, it is important to emphasize the generality of the underlying instability mech-
anism which may affect any set of three linearly coupled waves verifying the appropriate
matching conditions.

33
1.3.1 System of Three Coupled Oscillators
The problem of three non-linearly interacting waves is closely related to the dynamics of
a system of three coupled oscillators. One therefore starts by considering this somewhat
simpler system. This exercise will in particular point out the necessary condition of phase
matching in time.

A system of three non-linearly coupled oscillators can be characterized by a Hamilto-

nian of the form: !
X3
p2j 2
2 xj
H(x, p) = + ωj + V x 1 x2 x3 , (1.70)
j=1 2 2
where xj , pj , and ωj are respectively the position, momentum, and eigenfrequency of the
jth oscillator. The strength of the non-linear coupling is defined by the constant V .

The Hamilton-Jacobi equations:

dxj ∂H dpj ∂H
= , =− ,
dt ∂pj dt ∂xj
lead in the case of H given by (1.70) to the following set of coupled equations:

ẍ1 + ω12 x1 = −V x2 x3 , (1.71)

ẍ2 + ω22 x2 = −V x1 x3 , (1.72)
ẍ3 + ω32 x3 = −V x1 x2 . (1.73)

The left hand sides of these equations clearly correspond to the linear equation of motion
for each independent harmonic oscillator, while the right hand sides model the non-linear
coupling. The doted quantities correspond to the standard notation for time differentia-
tion.

If the non-linear coupling remains relatively small (|V xj |  |ωj20 |), a logical represen-
tation for the solution to the system (1.71)-(1.73) is of the form:
1h i
xj (t) = Aj (t)eiωj t + A?j (t)e−iωj t , (1.74)
2
where the (possibly complex) amplitude Aj (t) of the jth oscillator is expected to vary
slowly compared to the frequency ωj , i.e. |(dAj /dt)/Aj |  |ωj |.

Note, that if one wants to make use of a complex representation, one must be some-
what more careful for the non-linear physical system considered here than in the case
of a linear one. Indeed, if C is a complex solution to a linear, physical (implying with
real coefficients) set of equations, than the real part Re(C) and imaginary part Im(C)
are obviously solutions as well. This property is not verified in the case of a non-linear
system as a result of Re(z1 z2 ) which in general is not equal to Re(z1 )Re(z2 ), where z1,2
are complex values. When dealing with a non-linear physical system, one must therefore

34
ensure reality from the start, which explains the presence of the complex conjugate term
in Eq. (1.74).

Differentiating relation (1.74) with respect to time leads to

1h i
ẋj (t) = (Ȧj + iωj Aj )eiωj t + c.c. , (1.75)
2
1h i
ẍj (t) = (Äj + 2iωj Ȧj − ωj2 Aj )eiωj t + c.c.
2
1h i
= (Äj + 2iωj Ȧj )eiωj t + c.c. − ωj2 xj , (1.76)
2
where “c.c.” stands for complex conjugate. Inserting (1.76) into (1.71) provides an
equation for oscillator #1:
V h i
(Ä1 + 2iω1 Ȧ1 )eiω1 t + c.c. = − A2 A3 ei(ω2 +ω3 )t + A2 A?3 ei(ω2 −ω3 )t + c.c. ,
2
which, after multiplying by exp(−iω1 t), becomes:

V h
(Ä1 + 2iω1 Ȧ1 ) + (Ä?1 − 2iω1 Ȧ?1 )e−2iω1 t = − A2 A3 e−i(ω1 −ω2 −ω3 )t + A2 A?3 e−i(ω1 −ω2 +ω3 )t
2 i
+A?2 A3 e−i(ω1 +ω2 −ω3 )t + A?2 A?3 e−i(ω1 +ω2 +ω3 )t .(1.77)

To obtain a slow time scale equation for the amplitude A1 (t), one averages this last rela-
tion over the fast time scale of the eigenfrequencies ωj . Note, that the assumption of time
scale separation |Ȧj /Aj |  |ωj | also enables to neglect the terms Äj compared to 2iωj Ȧj
on the right hand side of Eq.(1.77).

In the case of frequency mismatch, such that all the phase factors exp −i(ω1 ± ω2 ± ω3 )t
on the right hand side of Eq.(1.77) vary on the fast time scale, i.e.

|ω1 ± ω2 ± ω3 | ∼ |ωj |,

the averaging process applied to Eq.(1.77) simply leads to

2iω1 Ȧ1 = 0 =⇒ A1 = const.,

Equivalent equations for A2 and A3 are also obtained by starting from Eqs. (1.72) and
(1.73) respectively. This result clearly shows that in the case of frequency mismatch the
oscillators are essentially decoupled, so that in the absence of damping there amplitude
remains constant.

However, if the condition of resonant coupling is verified, i.e.

ω1 = ω2 + ω3 + δω, (1.78)

35
the fast time scale averaging of Eq. (1.77) gives:
V
2iω1 Ȧ1 = − A2 A3 e−i δω t .
2
In Eq. (1.78) one has nonetheless allowed for the possibility of a small mismatch δω,
but this mismatch is assumed small such that |δω|  |ωj |. Similar equations for A2 and
A3 are obtained from Eqs. (1.72)-(1.73), providing the following system of non-linearly
coupled equations for the slow time scale variation of the oscillator amplitudes Aj :
V
2iω1 Ȧ1 = − A2 A3 e−i δω t , (1.79)
2
V
2iω2 Ȧ2 = − A1 A?3 e+i δω t , (1.80)
2
V
2iω3 Ȧ3 = − A1 A?2 e+i δω t . (1.81)
2
From the set of equations (1.79)-(1.81) one can now start by carrying out a stability
analysis of the state of the system in which one assumes that the oscillator #1 has been
initially excited with a much larger amplitude than the two other ones. Thus, considering
A2 and A3 as small perturbations compared to A1 (|A2,3 |  |A1 |), one can linearize the
system (1.79)-(1.81), which leads to

2iω1 Ȧ1 = 0 =⇒ A1 = A1,0 = const,

V
2iω2 Ȧ2 = − A1,0 A?3 e+i δω t ,
2
V
2iω3 Ȧ3 = − A1,0 A?2 e+i δω t .
2
By considering the Ansatz A2,3 (t) = a2,3 exp(γ + i δω/2)t, one then obtains:
δω V
2iω2 (γ + i ) a2 = − A1,0 a?3 ,
2 2
δω V
2iω3 (γ + i ) a3 = − A1,0 a?2 ,
2 2
To obtain a solution with non-zero amplitudes a2,3 , the determinant of this last linear
system must be zero, which provides the following relation for γ:
2 !2
|A1,0 |2

2 V δω
γ = − . (1.82)
4 ω2 ω3 2

To have instability requires γ 2 > 0. From Eq. (1.82), a necessary condition for instability
is thus
ω2 ω3 > 0, (1.83)
which is also a sufficient condition if there is no mismatch (δω = 0). The two conditions
(1.78) and (1.83) thus lead to
|ω1 | > |ω2 |, |ω3 |.

36
 Here the conditions for instability must be interpreted as the conditions for effective
energy transfer from oscillator #1 to the oscillators #2 and #3. We thus conclude that
for this transfer to be effective, the frequency ω1 of oscillator #1 must be larger than the
frequencies of the other two oscillators. In the frame of a quantum description, conditions
(1.78) and (1.83) correspond to the conditions of energy conservation between the energy
quanta of the three oscillators:

h̄|ω1 | = h̄|ω2 | + h̄|ω3 |.

In the presence of a frequency mismatch, condition (1.83) is obviously not sufficient for
instability. From Eq. (1.82), one clearly sees that there is an additional condition on the
amplitude of oscillator #1:
2|δω|(ω2ω3 )1/2
|A1,0 | > .
V
Exercises:
1.3.1.1 Re-derive the system of equations (1.79)-(1.81) but furthermore assuming that the
j’th oscillator undergoes damping with rate γj . Repeat the linear stability analysis
against parametric instabilities (for this, neglect the damping of oscillator #1), and
show that the relation (1.82) for the rate γ now becomes:
 !2 1/2
2
|A1,0 |2 

γ2 + γ 3  γ2 − γ 3 δω V
γ=− ± +i +
4 4 2 4 ω2 ω3

Re-asses the conditions for instability in this case. In particular, show that there is
now an amplitude threshold even in the case of perfect frequency match (δω = 0).

1.3.2 Illustration of Three Wave Coupling: Stimulated Raman

Scattering
As has been done in the previous section for the amplitudes of three coupled harmonic
oscillators, one now derives the equations governing the evolution of the envelopes of three
non-linearly interacting waves. This derivation is carried out here for the particular case
of stimulated Raman scattering (SRS), which, as we recall, involves the interaction of
1. An incident electromagnetic wave, the so-called pump.

2. A scattered electromagnetic wave.

3. An electron plasma wave (EPW).

The specific case of SRS enables to illustrate the derivation of a system of coupled equa-
tions for the wave amplitudes whose form is generic for any set of three non-linearly
interacting waves.

SRS involves two types of waves: Transverse electromagnetic waves, and electron plasma

37
 z

A
B y

x
Ees
~ = Ez ~ez , B
Figure 1.11: Fields (E ~ = By ~ey ) of linearly polarized electromagnetic waves repre-
~ and electrostatic field E
sented by the vector potential A ~ es = Ex~ex of plasma wave, for slab
model of SRS and SBS.

waves. One thus needs to derive the corresponding wave equations, and in particular
identify the dominant non-linear coupling terms. As the coupling involves no essential
wave-particle resonance, these equations are derived, for simplicity, in the frame of a fluid
description. One assumes furthermore the system to be one-dimensional (slab), so that
all fields depend on a single spatial variable x.

Equation for Transverse EM Waves

One starts by deriving the wave equation for the transverse EM waves. The corresponding
~ B)
field components (E, ~ verify Maxwell’s equations, and in particular

~
∇×E ~ = − ∂B ,
∂t
~
∇ · B = 0.

The general solution to these two equations can be expressed in terms of the vector
~ and scalar potential φ:
potential A
~
~ = − ∂ A − ∇φ,
E ~ = ∇ × A.
B ~
∂t
~ = 0, one has for a transverse wave
In Coulomb gauge, for which ∇ · A
~ = −4φ = 0
∇·E =⇒ φ = 0.

One shall furthermore assume here that the transverse wave is linearly polarized along
~ = Az (x, t) ~ez , and
the direction Oz (see Fig. 1.11), so that A

~ = Ez ~ez = − ∂Az ~ez ,

E ~ = By ~ey = − ∂Az ~ey .
B (1.84)
∂t ∂x
38
Inserting these relations into Ampere’s law

1 ∂E~
~ = µ0 ~j +
∇×B ,
c2 ∂t
provides
1 ∂ 2 Az ∂ 2 Az
= + µ 0 jz , (1.85)
c2 ∂t2 ∂x2
where c is the speed of light, and jz is the current along ~ez induced by the EM field in
the plasma. This plasma current is derived from the equation of motion for the particles:
d~v ~ + ~v × B).
~
m = q (E
dt
Projecting onto ~ez and inserting relations (1.84) provides:

dvz ∂Az ∂Az dAz d

m = q (Ez + vx By ) = −q ( + vx ) = −q =⇒ (mvz + qAz ) = 0.
dt ∂t ∂x dt dt

This last relation expresses the invariance of momentum p~ = m~v + q A ~ along Oz, which
results from the fact that the system is translationally invariant in this direction. If the
particles are essentially immobile before the passage of the wave, one obtains
qAz
vz = − , (1.86)
m
so that the transverse current can be written
X X N q2 N e e2
jz = N qvz = −Az '− Az , (1.87)
species species m me

having only kept the electron contribution, as the ion contribution is smaller by the ratio
me /mi . To lighten notations, the subscript “e” for electron values will be dropped from
here on.

By considering the electron density N = N0 + δN as the superposition of an initial,

homogeneous background N0 and of fluctuations δN , one finally obtains by inserting
(1.87) in (1.85):
∂ 2 Az 2
2 ∂ Az δN
2
− c 2
+ ωp2 Az = −ωp2 Az , (1.88)
∂t ∂x N0
where ωp2 = N0 e2 /m0 is the plasma frequency squared. The left hand side of this last
equation is simply the linear wave equation for transverse electromagnetic waves, giving
rise to the dispersion relation ω 2 = ωp2 + (kc)2 , and in particular to the well-known
condition for propagation |ω| > ωp . Longitudinal plasma waves (either EPWs or IAWs)
are naturally a source of density fluctuations δN , so that the term on the right hand side
of Eq. (1.88) provides the non-linear coupling term between the transverse EM waves and
the plasma waves.

39
Equation for Electron Plasma Waves
For modeling the high frequency EPWs, one can treat the massive ions as an immobile,
uniform, neutralizing background. The electrons are modeled here with a warm fluid
description. For this, one considers the continuity equation and momentum equation
along Ox:
∂N ∂
+ (N vx ) = 0, (1.89)
∂t ∂x
∂vx ∂vx ∂p
mN ( + vx ) = −eN (Ex − vz By ) − , (1.90)
∂t ∂x ∂x
as well as the closure from the equation of state:
p
= const., (1.91)
Nγ
where p is the electronic pressure. For EPWs, the phase velocity vφ = ω/k is such that
|vφ |  vth . One thus considers the adiabatic equation of state for which γ = (D + 2)/D,
where D is the number of degrees of freedom. For wave propagation the appropriate value
is D = 1, so that γ = 3.

The electric component Ex on the right hand side of Eq. (1.90) corresponds to the
longitudinal, electrostatic field related to the EPWs, and verifies Poisson’s equation:
∂Ex e δN
=− . (1.92)
∂x 0
where δN is the electron density perturbation related to the EPWs.

The next term on the right hand side of Eq. (1.90) is the Lorentz force Fp = e vz By ,
which results from the transverse oscillatory motion represented by Eq. (1.86), and gives
rise to the so-called ponderomotive force:

eAz ∂Az 1 ∂ e2 A2z

Fp = e v z B y = e (− )=− ( ).
m ∂x 2 ∂x m
The ponderomotive force provides the non-linear coupling of the EPWs with the transverse
EM waves. This is therefore the only non-linearity which is retained, and all other terms
in Eqs. (1.89)-(1.91) are linearized with respect to the EPW perturbation terms δN , v x ,
Ex and δp:
∂ δN ∂ vx
+ N0 = 0, (1.93)
∂t ∂x
∂ vx 1 ∂ e2 A2z 1 ∂δp
m = −e Ex − ( )− , (1.94)
∂t 2 ∂x m N0 ∂x
p0
δp = γ δN = 3 T0 δN, (1.95)
N0

40
having used p0 = N0 T0 , where T0 is the background electron temperature.

Inserting Eq. (1.92) into Eqs. (1.93) and (1.95) enables to express the longitudinal
velocity and pressure perturbations in terms of Ex :
0 ∂Ex
vx = ,
e N0 ∂t
T0 ∂Ex
δp = −3 0 .
e ∂x
These relations are in turn inserted into Eq. (1.94), finally providing:

∂ 2 Ex 2
2 ∂ Ex 2
2
2 1 ∂ eAz
− 3v th + ω p E x = −ω p . (1.96)
∂t2 ∂x2 2 ∂x m
The left hand side of Eq. (1.96) is clearly the linear wave equation for EPWs, leading to
the Bohm-Gross dispersion relation ω 2 = ωp2 + 3(kvth )2 . The right hand side represents
the non-linear coupling with the transverse EM waves.

Coupled Three Wave Equations

To summarize the above results, one rewrites here the wave equations for the transverse
EM waves with vector potential component Az and EPWs with electrostatic component
Ex :
∂ 2 Az 2
2 ∂ Az e ∂Ex
2
− c 2
+ ωp2 Az = Az , (1.97)
∂t ∂x m ∂x
∂ 2 Ex 2
2 ∂ Ex 2 2 1 ∂ eAz
2
− 3v th + ω E
p x = −ω p , (1.98)
∂t2 ∂x2 2 ∂x m
having used (1.92) to replace δN by Ex in Eq. (1.88).

To obtain a model for SRS, one now intends to derive from Eqs. (1.97)–(1.98) a system
of coupled equations for the amplitudes A0 , As and E of the incident EM wave, scattered
EM wave, and the EPW respectively. The method for deriving these equations is similar
to the one used for obtaining the system (1.79)–(1.81) for the amplitudes of three cou-
pled harmonic oscillators, with the additional complication however of spatial dependence.

Under the assumption of small non-linear coupling, it is logical to assume that each
wave involved in the SRS mechanism has a wave number k and frequency ω still verifying
the linear dispersion relation, but with an amplitude that may vary slowly both in space
and time. The vector potential field Az in Eqs. (1.97)–(1.98), which is the superposition
of the incident and scattered EM waves, is thus written:
1h i 1h i
Az (x, t) = A0 (x, t) ei(k0 x−ω0 t) + c.c. + As (x, t) ei(ks x−ωs t) + c.c. , (1.99)
|2 {z } |2 {z }
Incident EM Scattered EM

41
while a similar Ansatz is considered for the electrostatic field:
1h i
Ex (x, t) = E(x, t) ei(ke x−ωe t) + c.c. , (1.100)
|2 {z }
EPW
where the wave number-frequency pairs (k, ω) of all three waves verify their respective
linear dispersion relations:

ω02 = ωp2 + (k0 c)2 , (1.101)

ωs2 = ωp2 2
+ (ks c) , (1.102)
ωe2 = ωp2 2
+ 3(ke vth ) . (1.103)

As an example, the condition of slowly varying envelopes reads in the case of the incident
EM wave:
1 ∂A0 1 ∂A0
| |  |k0 |, and | |  |ω0 |.
A0 ∂x A0 ∂t
Similar scalings hold for As and E.

To obtain the equations for the amplitudes A0 and As for the incident and scattered
EM waves, one starts by differentiating Eq. (1.99) with respect to x and t:
X 1h i
∂x Az = (∂x A + ikA) ei(kx−ωt) + c.c. ,
0,s 2
X1h i
∂xx Az = (∂xx A + 2ik ∂x A − k 2 A) ei(kx−ωt) + c.c. ,
0,s 2
X1h i
∂t Az = (∂t A − iωA) ei(kx−ωt) + c.c. ,
0,s 2
X1h i
∂tt Az = (∂tt A − 2iω ∂t A − ω 2 A) ei(kx−ωt) + c.c. .
0,s 2

Differentiating Eq.(1.100) provides similar relations for Ex , which can then all be inserted
into Eq.(1.97), leading to:
1h i 1h i
(∂tt A0 − 2iω0 ∂t A0 ) ei(k0 x−ω0 t) + c.c. − c2 (∂xx A0 + 2ik0 ∂x A0 ) ei(k0 x−ω0 t) + c.c.
2 2
1h i 1h i
+ (∂tt As − 2iωs ∂t As ) ei(ks x−ωs t) + c.c. − c2 (∂xx As + 2iks ∂x As ) ei(ks x−ωs t) + c.c.
2 2
e 1h i
= (∂x E + ike E) ei(ke x−ωe t) + c.c. ×
m 2
1h i 1h i
i(k0 x−ω0 t) i(ks x−ωs t)
A0 e + c.c. + As e + c.c. (1.104)
2 2
having made use of Eqs. (1.101)-(1.102). Invoking the assumption of slow variation of the
envelope, one has |∂tt A|  |ω ∂t A|, |∂xx A|  |k ∂x A| and |∂x E|  |ke E|, which justifies

42
neglecting these smaller terms in Eq.(1.104).

By multiplying Eq. (1.104) by exp −i(k0 x − ω0 t) and averaging over the fast space and
time scales of the wave numbers k and wave frequencies ω, one obtains a slow variation
scale equation for the incident wave amplitude A0 (x, t). Similarly as for the system of
coupled harmonic oscillators, one sees that the non-linear coupling terms on the right
hand side of Eq. (1.104) will in general average out to zero unless certain resonant con-
ditions are met. These conditions correspond to phase matching of the three waves both
in space and time, and can be written in terms of the wave numbers and frequencies as:

k0 = k s + k e , (1.105)
ω0 = ωs + ωe + δω, (1.106)

having allowed in (1.106) for a possible frequency mismatch δω, such that |δω|  |ω0,s,e|.
Under these matching conditions, one then obtains after averaging:
e
−2iω0 ∂t A0 − c2 2ik0 ∂x A0 = ike EAs eiδωt .
2m
In the same way, by multiplying Eq. (1.104) by exp −i(ks x−ωs t) and performing the same
averaging leads to the corresponding equation for the scattered wave amplitude A s (x, t):
e
−2iωs ∂t As − c2 2iks ∂x As = − ike E ? A0 e−iδωt .
2m
The equation for the EPW envelope E(x, t) is obtained through a similar derivation by
inserting (1.99) and (1.100) into (1.98), multiplying by exp −i(ke x − ωe t), and averag-
ing. Again invoking the assumption of phase matching, one then obtains (check it as an
exercise):
2 e 2
−2iωe ∂t E − 3vth 2ike ∂x E = − ω ike A0 A?s e−iδωt .
2m p
The system of non-linearly coupled equations for the three wave amplitudes can thus be
summarized as follows:
e ke
∂t A0 + vg,0 ∂x A0 = − E As eiδωt , (1.107)
4m ω0
e ke
∂t As + vg,s ∂x As = A0 E ? e−iδωt , (1.108)
4m ωs
e ke 2
∂t E + vg,e ∂x E = ω A0 A?s e−iδωt , (1.109)
4m ωe p

having used the notations vg,0,s = dω0,s /dk0,s = c2 k0,s /ω0,s and vg,e = dωe /dke = 3vth
2
ke /ωe
for the group velocities of the three waves.

In case the phase matching conditions (1.105)–(1.106) are not verified, the non-linear
coupling terms on the right hand sides of Eqs. (1.107)–(1.109) are absent. The resulting

43
 equations then simply represent the independent advection of the three wave envelopes
at their respective group velocity. For example, for the incident wave:

∂t A0 + vg,0 ∂x A0 = 0 =⇒ A0 (x, t) = A0 (x − vg,0 t, t = 0). (1.110)

In the case of space independent wave amplitudes, one can easily show the equivalence
of the system of coupled equations (1.107)-(1.109) with the system (1.79)-(1.81) for the
amplitudes of the coupled harmonic oscillators. For this, one simply needs to identify the
corresponding terms (harmonic oscillators ↔ waves): (A1 , ω1 ) ↔ (−iA0 , −ω0 ), (A2 , ω2 ) ↔
(−iAs , −ωs ), (A3 , ω3 ) ↔ (−iE/ωp , −ωe ), V ↔ e ke ωp /m, and δω ↔ −δω.

Linear Analysis of the Parametric Instability

As in section 1.3.1, one can now start by carrying out a linear stability analysis of the
state of the system in which the incident EM wave has an amplitude A0 significantly
larger than the amplitudes As and E of the scattered wave and EPW respectively. In
a real physical system, A0 could represent a high intensity laser beam, while the initial
values of As and E could be at the level of thermal fluctuations. One furthermore assumes
here that there is no spatial variations of the envelopes, so that only temporal variations
are considered. Thus, considering As and E as small perturbations, one can linearize the
system (1.107)-(1.109) with respect to these terms, which leads to

∂t A0 = 0 =⇒ A0 = const, (1.111)
e ke
∂t As = A0 E ? e−iδωt , (1.112)
4m ωs
e ke 2
∂t E = ω A0 A?s e−iδωt . (1.113)
4m ωe p
Thus, considering A0 as a constant, Eqs. (1.112)-(1.113) become linear in As and E. The
linear solution for these two fields can be easily found by expressing As , E ∼ exp(γ −
iδω/2)t, which leads to the rate:
!2 !2
2 ke vos ωp2 δω
γ = − , (1.114)
4 ωs ωe 2

where one has defined vos = e|A0 |/m the velocity oscillation amplitude of electrons in the
incident EM wave. Note again the analogy between (1.114) and (1.82).

From Eq. (1.114), one sees that a necessary condition for instability is thus again ωs ωe > 0,
which together with (1.106) implies that all three frequencies ω0,s,e must have same sign,
and that |ω0 | > |ωs |, |ωe |. From here on, as a convention, all three frequencies ω0,s,e can
thus be assumed positive.

Exercises:
1.3.2.1 Carry out the derivation of equation (1.109).

44
1.3.2.2 How do Eqs. (1.107)-(1.109) generalize in the presence of possible damping mecha-
nisms for the various waves?

1.3.2.3 Derive the non-linear coupling equations for the envelopes of the three waves in-
volved in stimulated Brillouin scattering (SBS). Compute the corresponding linear
growth rate of the parametric instability.

1.3.3 Matching Conditions

Let us make some additional comments here relative to the matching conditions for non-
linear three wave coupling. For illustration purposes, one pursues this discussion in the
particular case of SRS.

For the one-dimensional slab model of SRS that was considered in Sec. 1.3.2, it ap-
pears clearly, that having fixed one of the 6 real wave number-frequency values (k 0 , ω0 ),
(ks , ωs ) and (ke , ωe ) the other values are in general determined by the system of 5 equa-
tions formed by the 3 dispersion relations (1.101)-(1.103) and the 2 matching conditions
(1.105) and (1.106).

Note however, that the dispersion relations contain quadratic terms, and therefore, for
a given incident wave with wave number-frequency pair (k0 , ω0 ) their may be zero, one
or two solutions for (ks , ωs ) and (ke , ωe ) verifying the matching conditions. Obviously,
for SRS, as a result of the properties ωs > ωp and ωe > ωp , one must have from (1.106)
ω0 > 2ωp . By defining the critical density Nc such that ω02 = Nc e2 /m0 = (Nc /N )ωp2 , one
thus obtains a necessary condition for SRS to be able to develop (density below quarter
critical):
Nc
N< .
4
The geometrical solution to the matching conditions (1.106)-(1.105) in a one-dimensional
system appears as a sum of vectors in the (k, ω) plane: (k0 , ω0 ) = (ks , ωs ) + (ke , ωe ),
where the three vectors (k0,s,e, ω0,s,e ) must lie on the curves of their respective dispersion
relation. This is shown in the case of SRS in Fig. 1.12. From this figure, one clearly sees
∼
that for ω0 > 2ωp there are in general two solutions to the matching conditions for a given
incident wave (k0 , ω0 ). One solution is such that k0 ks < 0, and is called Backward SRS
(BSRS), as it corresponds to the scattered EM wave propagating opposite to the incident
wave. The other solution with k0 ks > 0 is called Forward SRS (FSRS), and corresponds
to the scattered wave propagating in the same direction as the incident one. According
to (1.114), the SRS growth rate related to backward scattering tends to be larger than
for forward scattering, as |keBSRS | > |keFSRS |. This is indeed true as long as damping of
the EPW is not important. However, as Landau damping increases with ke λD , FSRS can
become competitive when keBSRS λD becomes large.

Given the frequency ω0 of the incident laser light, the matching conditions for SRS can
thus be solved together with the dispersion relations to obtain the wave number-frequency

45
 BSRS
5
Incident
Dispersion rel.
Laser Light
for EM waves (k0, ω0)
4

Scattered
Laser Light
3
(ks, ωs)
ω / ωp

2
Dispersion rel.
for EPW

1 EPW
(k , ω )
e e

(a)
0
−0.5 0 0.5
k λD
FSRS
5
Incident
Dispersion rel.
Laser Light
for EM waves (k0, ω0)
4
Scattered
Laser Light

3 (k , ω )
s s
ω / ωp

2
Dispersion rel.
for EPW

1
EPW
(k , ω )
e e
(b)
0
−0.5 0 0.5
k λD

Figure 1.12: Geometrical solution to the matching conditions for (a) backward and (b) forward
SRS. The blue curve corresponds to the dispersion relation ω 2 = ωp2 + (kc)2 for electromagnetic
waves, and the red curve shows the Bohm-Gross dispersion relation ω 2 = ωp2 + 3(kvth )2 for
EPWs. Space and time phase matching: The purple vector, representing the wave number and
frequency (ω0 , k0 ) of the incident laser light, must be the sum of the vectors (k s , ωs ) and (ke , ωe )
representing the scattered light and EPW respectively.

46
pairs of all waves involved in the SRS mechanism. This is achieved by first inserting (1.105)
and (1.106) into (1.102), which (neglecting possible mismatch δω) leads to:

(ω0 − ωe )2 = ωp2 + (k0 − ke )2 c2 .

Making use of Eq. (1.101) for (k0 , ω0 ), this relation can then be reduced to
ωe (ωe − 2ω0 )
ke2 − 2k0 ke − = 0,
c2
which in turn is solved for ke :
1 2 2
ke = k 0 ± (k c + ωe2 − 2ω0 ωe )1/2
c 0
ω0 ωp
' k0 ± (1 − 2 )1/2 , (1.115)
c ω0
having approximated ωe ' ωp in the last step, and having again invoked Eq. (1.101).
From the matching condition on the wave numbers, one has ks = k0 − ke , so that
ω0 ωp
ks = ∓ (1 − 2 )1/2 .
c ω0
For k0 > 0, the solution ke with the + sign in (1.115) is thus related to BSRS, while the
solution with the − sign corresponds to F SRS. From Eq.(1.115) it also appears clearly
that keBSRS varies from 2k0 in the limit of low density (N  Nc /c ⇔ ω0  2ωp ), down to
k0 in the limit N ' Nc /4 ⇔ ω0 ' 2ωp .

If the waves are allowed to propagate in more than one dimensions, the wave number
matching condition (1.105) is replaced by the vector equation on the wave vectors:
~k0 = ~ks + ~ke . (1.116)

In this case, having fixed for instance the frequency ω0 and direction ~n0 = ~k0 /k0 of the
incident wave, there remain 9 parameters to be determined for the wave vectors and fre-
quencies of the three waves involved in the parametric instability: k0 , (~ks , ωs ), and (~ke , ωe ).
These 9 parameters must verify the system of 7 equations formed by the 3 dispersion re-
lations (1.101)-(1.103) and the 4 matching conditions (1.106) and (1.116). Thus, there
remain 2 degrees of freedom. These can be chosen as the scattering direction ~ns = ~ks /ks .
Therefore, in general, scattering can occur in all spatial directions. Note however that
each scattering direction has its own growth rate and interaction length.

In the context of a quantum description, the matching condition (1.116) on the wave
vectors and the matching condition (1.106) on the frequencies correspond respectively to
the conservation of momentum and energy of an incident photon scattering off a plasmon:

h̄~k0 = h̄~ks + h̄~ke , (1.117)

h̄ω0 = h̄ωs + h̄ωe , (1.118)

47
 where (h̄~k0 , h̄ω0 ), (h̄~ks , h̄ωs ), and (h̄~ke , h̄ωe ) are respectively the (momentum, energy) of
the incident photon, scattered photon, and plasmon involved in the process.

Exercises:

1.3.3.1 Draw the geometrical solution to the matching conditions for SBS.

1.3.3.2 From the appropriate plot, convince yourself that matching conditions can be ver-
ified for the potential decay of an EPW into another EPW and an IAW. This
parametric instability is the so-called Langmuir Decay Instability (LDI).

1.3.4 Manley-Rowe Relations

One can show that the non-linear system of coupled equations (1.107)–(1.109) for the
amplitudes of the three waves involved in the parametric instability verifies certain con-
servation laws. These laws are conveniently derived after an appropriate normalization of
the wave amplitudes.

Inspired by the equivalence between the matching conditions on the wave vectors/frequencies
of the waves in the classical description and the conservation laws (1.117)-(1.118) of the
corresponding quantum process, one defines for each wave the complex action amplitude
a, such that the corresponding action density n = |a2 | = a a? is given from the wave
energy density Ewave by
Ewave = a a? ω = n ω. (1.119)
To a factor h̄, the action density n can thus be interpreted as the density of wave quanta.

The wave energy density Ewave for each wave is derived from the general relation for
the wave energy in a dispersive media:
( )
1 ~? + 1 H
~ 0 · ∂ [ω H (ω)] · E ~ 0 · ∂ [ω µH (ω)] · H
~? ,
Ewave = 0 E 0 0 (1.120)
4 ∂ω µ0 ∂ω

where H and µH are the hermitian parts of the dielectric and magnetic permeability
~ 0 and B
tensors respectively, E ~ 0 = µH
~ 0 are the complex amplitudes of the electric and
magnetic field components of the wave respectively, and ω is the frequency of the wave.
Note, that in a plasma, the magnetic permeability tensor µ = 1. Proof of this relation
was the goal of exercise 1.2.2.1.

In the case of transverse electromagnetic waves, the dielectric tensor is given by  =

(1 − ωp2 /ω 2 )1, and, according to (1.84), for a given wave with wave number k and fre-
quency ω the amplitude of the electric and magnetic fields are related to the amplitude
A of the vector potential by

E0 = iωA, B0 = −ikA,

48
so that by applying Eq. (1.120) one obtains:
" ! # " #
1 ∂ ωp2 1 2 2 1 ωp2 + (kc)2
Ewave = 0 ω 2 |A|2 ω− + k |A| = 0 ω 2 |A|2 1 +
4 ∂ω ω µ0 4 ω2
1
= 0 ω 2 |A|2 , (1.121)
2
having made use of the dispersion relation ω 2 = ωp2 + (kc)2 for transverse waves. Com-
paring Eq. (1.121) with (1.119), the action amplitudes for the incident and scattered
electromagnetic waves participating in the SRS mechanism are thus given respectively
by:
   
0 ω0 1/2 0 ωs 1/2
a0 = A0 , as = As . (1.122)
2 2
For EPWs, the dielectric function in the warm fluid limit considered here (see model for
EPW in Sec. 1.3.2) is given by  = 1 − ωp2 /[ω 2 − 3(kvth )2 ], so that for a given wave with
frequency ω and amplitude E of the electric field, one obtains from Eq. (1.120):

1 ∂ 1 ∂ 1 ω2
Ewave = 0 (ω) |E|2 = 0 ω |E|2 = 0 2 |E|2 , (1.123)
4 ∂ω 4 ∂ω 2 ωp

having made use of the dispersion relation  = 0, i.e. ω 2 = ωp2 + 3(kvth )2 , for the EPWs.
Note that in this case of a longitudinal wave, there is no magnetic field contribution to
the energy. Comparing Eq. (1.123) with (1.119), the action amplitude for the EPW in
the SRS mechanism is thus given by:
 1/2
 0 ωe E
ae = . (1.124)
2 ωp

Inserting relations (1.122) and (1.124) into Eqs. (1.107)–(1.109), one obtains the normal-
ized system of coupled equations for the action amplitudes of the three waves participating
in SRS:

∂t a0 + vg,0 ∂x a0 = −Γ ae as eiδωt , (1.125)

∂t as + vg,s ∂x as = +Γ a0 a?e e−iδωt , (1.126)
∂t ae + vg,e ∂x ae = +Γ a0 a?s e−iδωt , (1.127)

where one has defined the normalized coupling parameter

1 e ωp ke
Γ= √ √ .
2 2 m 0 (ω0 ωs ωe )1/2

Note the symmetry in equations (1.125)–(1.127). Symmetry which was lacking in Eqs.
(1.107)–(1.109).

The form of Eqs. (1.125)–(1.127) is now convenient for deriving conservation relations, in

49
particular the conservation of action. Indeed, by multiplying Eqs. (1.125) by a?0 , and by
multiplying the complex conjugate of Eq. (1.126) by as , one obtains:

a?0 ∂t a0 + vg,0 a?0 ∂x a0 = −Γ a?0 as ae eiδωt ,

as ∂t a?s + vg,s as ∂x a?s = +Γ a?0 as ae eiδωt .

Then, by adding these two equations:

a?0 ∂t a0 + vg,0 a?0 ∂x a0 + as ∂t a?s + vg,s as ∂x a?s = 0. (1.128)

and also considering the complex conjugate of this last relation:

a0 ∂t a?0 + vg,0 a0 ∂x a?0 + a?s ∂t as + vg,s a?s ∂x as = 0, (1.129)

one then finally obtains by adding Eqs. (1.128) and (1.129):

∂t |a0 |2 + vg,0 ∂x |a0 |2 + ∂t |as |2 + vg,s ∂x |as |2 = 0,

which can also be written in terms of the action densities as:

∂t n0 + ∂x (vg,0 n0 ) = − [∂t ns + ∂x (vg,s ns )] . (1.130)

Noting, that vg,0 n0 and vg,s ns are the action fluxes for the incident and scattered waves
respectively, one identifies the left hand side of Eq. (1.130) as the continuity equation for
the action in the incident wave, while the right hand side corresponds to the continuity
equation for action in the scattered wave. Equation (1.130) clearly states that the sink of
action in the incident wave is the source of action in the scattered wave. The integral form
of this local conservation law is obtained by integrating Eq. (1.130) over space, which
leads to Z
(n0 + ns ) dx = const., (1.131)

having assumed that waves do not propagate in or out of the system (true in particular
for a periodic system).

Starting from Eqs. (1.125) and (1.127), one obtains by a similar derivation the equa-
tions for action transfer between the incident wave and the EPW, both in local form:

∂t n0 + ∂x (vg,0 n0 ) = − [∂t ne + ∂x (vg,e ne )] . (1.132)

and in global form: Z

(n0 + ne ) dx = const. (1.133)

Equations (1.130)-(1.133) are the so-called Manley-Rowe relations (in local and global
form). Equations (1.131) and (1.133) state that for each quantum disappearing in the
incident wave, a quantum appears both in the scattered as well as in the electron plasma
wave.

50
 Together with the frequency matching condition (1.106), the Manley-Rowe equations also
imply energy conservation. Indeed, multiplying Eq. (1.130) by ωs and Eq. (1.132) by ωe
and then adding these two relations, leads to

(ωs + ωe ) [∂t n0 + ∂x (vg,0 n0 )] + ωs [∂t ns + ∂x (vg,s ns )] + ωe [∂t ne + ∂x (vg,e ne )] = 0. (1.134)

Invoking the matching condition ω0 = ωs + ωe , Eq. (1.134) can finally be written

X X
∂t ( Ewave ) + ∂x ( vg Ewave ) = 0, (1.135)
0,s,e 0,s,e

having made use of the relation Ewave = nω between the action density and wave density.
P
Identifying 0,s,e vg Ewave as the total energy flux from the three waves, Eq. (1.135)
obviously corresponds to the local energy conservation law. The corresponding global
conservation law is again obtained by integrating Eq. (1.135) over space, leading to
XZ
Ewave dx = const. (1.136)
0,s,e

Starting from Eqs.(1.79)–(1.81), relations analog to the global (i.e. space independant)
conservation laws (1.131), (1.133) and (1.136) can be derived for the system of three non-
linearly coupled harmonic oscillators (see following exercise).

Exercise:
1.3.4.1 For the system of three non-linearly coupled harmonic oscillators, define the appro-
priate action amplitude aj for each oscillator. Starting from Eqs. (1.79)–(1.81) for
the oscillator amplitudes Aj , obtain the more symmetric set of equations for the
action amplitudes aj . Finally, derive the conservation laws for action and energy.

1.3.5 Non-Linear Analytic Solution to the Space Independent

Coupling Equations
By carrying out a linear stability analysis of the three wave system (1.107)–(1.109) at
the end of Sec. 1.3.2 [which is naturally equivalent to the stability analysis of the system
(1.125)–(1.127)], one identified under which conditions a parametric instability may de-
velop. One shall now consider the non-linear evolution of the three wave interaction. To
facilitate the analytic derivation, one will assume that the envelopes of the waves are space
independant (valid in an infinitely long or periodic system), and furthermore assume that
there is no frequency mismatch. In this case, Eqs. (1.125)–(1.127) become

a˙1 = −Γ a2 a3 , (1.137)
a˙2 = +Γ a1 a?3 , (1.138)
a˙3 = +Γ a1 a?2 , (1.139)

having replaced the subscripts 0, s, and e used for labeling the waves involved in the SRS
mechanism by the subscripts 1, 2, and 3, so as to emphasize the generality of the following

51
results to any set of three non-linearly interacting waves. Note that one again uses here
the doted notation for time differentiation. Recall as well that Eqs. (1.137)–(1.139) have
been derived under the assumption of the wave frequencies ω1,2,3 being all positive, which,
together with the frequency matching condition ω1 = ω2 + ω3 , ensures |ω1 | > |ω2 |, |ω3 |.
This is the necessary condition for possible decay of quanta from wave #1 into quanta of
wave #2 and #3.

In order to derive the non-linear time evolution of the three wave amplitudes governed
by Eqs. (1.137)–(1.139), one starts by expliciting the modulus αj and phase φj for each
action amplitude aj :
aj (t) = αj (t) eiφj (t) , j = 1, 2, 3. (1.140)
Inserting (1.140) into (1.137)–(1.139), and defining the phase difference θ = φ1 − φ2 − φ3 ,
one obtains:
α˙1 + i α1 φ̇1 = −Γ α2 α3 e−iθ ,
α˙2 + i α2 φ̇2 = +Γ α1 α3 e+iθ ,
α˙3 + i α3 φ̇3 = +Γ α1 α2 e+iθ .
Taking the real parts of these relations provides:
α˙1 = −Γ α2 α3 cos θ, (1.141)
α˙2 = +Γ α1 α3 cos θ, (1.142)
α˙3 = +Γ α1 α2 cos θ, (1.143)
and the imaginary parts:
α1 φ̇1 = Γ α2 α3 sin θ, (1.144)
α2 φ̇2 = Γ α1 α3 sin θ, (1.145)
α3 φ̇3 = Γ α1 α2 sin θ. (1.146)
From Eqs. (1.141)–(1.143) one can naturally recover the Manley-Rowe relations:
d
(n1 + n2 ) = 2 α1 α̇1 + 2 α2 α̇2 = 0,
dt
d
(n1 + n3 ) = 2 α1 α̇1 + 2 α3 α̇3 = 0,
dt
having used the relation n = |a|2 = α2 between the action density n and the modulus α
of the action amplitude a. From Eqs. (1.144)–(1.146), one can obtain an equation for the
time variation of the phase difference θ:
 
α2 α3 α1 α3 α1 α2
θ̇ = φ̇1 − φ̇2 − φ̇3 = Γ − − sin θ,
α1 α2 α3
and making again use of Eqs. (1.141)–(1.143) one can further write:
 
α̇1 α̇2 α̇3 d
θ̇ = − + + tan θ = − tan θ ln(α1 α2 α3 ),
α1 α2 α3 dt

52
which leads to another invariant quantity:

θ̇ d d
= − ln(α1 α2 α3 ) =⇒ ln(α1 α2 α3 sin θ) = 0.
tan θ dt dt
The invariants for the system (1.137)–(1.139) can thus be summarized as follows:
.
n1 + n2 = const. = m2 , (1.147)
.
n1 + n3 = const. = m3 , (1.148)
.
α1 α2 α3 sin θ = const. = C. (1.149)

These constants are now used in deriving the non-linear time evolution of the spatially
uniform action densities nj . Starting from Eq. (1.141), one computes the time variation
of n1 :
ṅ1 = 2 α1 α̇1 = −2 Γ α1 α2 α3 cos θ.
Writing cos θ = ±(1 − sin2 θ)1/2 and making use of Eq. (1.149) then leads to:
!1/2
C2
ṅ1 = ± 2 Γ α1 α2 α3 1− 2 2 2 = ± 2 Γ (n1 n2 n3 − C 2 )1/2 .
α1 α2 α3

Finally, one invokes the invariants given by Eqs. (1.147) and (1.148) to obtain:
h i1/2
ṅ1 = ± 2 Γ n1 (m2 − n1 ) (m3 − n1 ) − C 2 .

This last equation depends only on n1 and can thus be solved by quadrature:
Z
1 n1 (t) dn1
= ± Γ (t − t0 ), (1.150)
2 n1 (t0 ) [n1 (m2 − n1 ) (m3 − n1 ) − C 2 ]1/2
| {z }
I

where for now, t0 represents an arbitrary reference time. The integral I on the left hand
side of Eq. (1.150) is an elliptic integral [5], and can in fact be expressed in terms of the
elliptic integral of the first kind, as is now shown. For this one notes na , nb and nc the
roots of the cubic equation

n1 (m2 − n1 ) (m3 − n1 ) − C 2 = 0. (1.151)

These roots can be ordered such that

0 ≤ n a ≤ nb ≤ nc . (1.152)

The integral I can thus be written:

Z
1 n1 (t) dn1
I= .
2 n1 (t0 ) [(n1 − na ) (n1 − nb ) (n1 − nc )]1/2

53
 Elliptic Integral of First Kind F(φ, m) Complete Elliptic Integral of First Kind F(m)
3
m=0
F(m) = F(π/2, m)
2.5 m = 0.4
u = F(φ , m) units of π

1.5
m = 0.8
2 m = 0.95

F(m) / π
1.5
1

0.5
(a) 0.5 (b)
0 0 0.2 0.4 0.6 0.8 1
0 0.5 1 1.5 2
φ/π m
Jacobian elliptic function sn(u, m) Jacobian elliptic function cn(u, m)
m=0
1 1 m = 0.4
0.8 m = 0.8
m = 0.95
0.6
0.5
0.4
sn(u, m)

0.2
cn(u, m)

0 0
−0.2
−0.4
−0.5 m=0
m = 0.4 −0.6
m = 0.8 −0.8
m = 0.95
−1 −1
(c) (d)
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
u/π u/π

Figure 1.13: (a) Elliptic integral of the first kind F (ϕ, m). (b) Complete elliptic integral of
the first kind F (m) = F (π/2, m). (c) Jacobian elliptic function sn(u, m). (d) Jacobian elliptic
function cn(u, m).

54
By further making the change of variables
 1/2
n1 − n a
y= , (1.153)
nb − n a
one obtains (check it):
Z
1 y(t) dy 1
I= = F [arcsin y(t), µ2 ], (1.154)
(nc − na )1/2 0 [(1 − y 2 )(1 − µ2 y 2 )]1/2 (nc − na ) 1/2

having defined the parameter

nb − n a
µ2 = , (1.155)
nc − n a
which, as a result from the ordering (1.152), is such that 0 ≤ µ2 ≤ 1. To obtain the first
equality in (1.154), one has now defined the time t0 such that

n1 (t0 ) = na ,

which led to the lower boundary of the integral in y being equal to zero. The second
equality in (1.154) was obtained using the definition of the elliptic integral of the first
kind [5]:
Z ϕ Z sin ϕ
dθ dx
F (ϕ, m) = 2 = ,
0 (1 − m sin θ) 1/2 0 [(1 − x2 )(1 − mx2 )]1/2
with x = sin θ, and 0 ≤ m ≤ 1.

Inserting (1.154) into (1.150) gives then:

F [arcsin y(t), µ2 ] = ± (nc − na )1/2 Γ (t − t0 ).

This last relation can be inverted for y(t):

h i
y(t) = ± sn (nc − na )1/2 Γ (t − t0 ), µ2 . (1.156)

Here sn(u, m), 0 ≤ m ≤ 1, is one of the Jacobian elliptic functions defined by

u = F (ϕ, m) ⇐⇒ sn(u, m) = sin ϕ. (1.157)

Note that one has sn(u, m = 0) = sin u, and for 0 < m < 1 the elliptic function sn(u, m)
appears as a “flattened” sin function (see Fig. 1.13). Let us immediately introduce here
cn(u, m = 0), which is another Jacobian elliptic function, defined by

u = F (ϕ, m) ⇐⇒ cn(u, m) = cos ϕ. (1.158)

Comparing definitions (1.157) and (1.158) for sn(u, m) and cn(u, m) respectively, one ob-
viously has the relation sn2 (u, m) + cn2 (u, m) = 1 for any argument u and parameter
0 ≤ m ≤ 1.

55
Now, inserting (1.156) into (1.153) finally provides the solution for the action density
of wave #1: h i
n1 (t) = na + (nb − na ) sn2 (nc − na )1/2 Γ (t − t0 ), µ2 . (1.159)
Note that n1 (t) oscillates between the values na and nb with period T = 2F (µ2 )/[Γ(nc −
na )1/2 ], where F (m) = F (π/2, m) is the complete elliptic integral of the first kind, and
that one indeed has n1 (t0 ) = na .

The time evolution of the action densities for wave #2 and #3 are then simply obtained
from (1.147) and (1.148):
h i
n2 (t) = m2 − n1 (t) = (m2 − nb ) + (nb − na ) cn2 (nc − na )1/2 Γ (t − t0 ), µ2 , (1.160)
h i
n3 (t) = m3 − n1 (t) = (m3 − nb ) + (nb − na ) cn2 (nc − na )1/2 Γ (t − t0 ), µ2 . (1.161)

Equations (1.159)–(1.161) clearly provide solutions for the action densities of the three
waves in terms of the initial conditions. Note in particular, that the values na , nb and nc ,
being solution of Eq. (1.151), are function of the invariants of motion m2 , m3 and C, and
thus function of the initial conditions of the system.

One now considers the solutions (1.159)–(1.161) for the action densities of the three non-
linearly coupled waves in two particular cases of initial conditions.

Case 1. n3 (0)  n2 (0) > n1 (0) = 0.

For these initial conditions, the invariants (1.147)–(1.149) become

m2 = n1 (0) + n2 (0) = n2 (0),

m3 = n1 (0) + n3 (0) = n3 (0),
C = [n1 (0) n2 (0) n3 (0)]1/2 sin[θ(0)] = 0.

The properly ordered roots na,b,c to the cubic equation (1.151) thus are given by

[na = 0] < [nb = m2 = n2 (0)] < [nc = m3 = n3 (0)],

and the parameter µ2 = n2 (0)/n3 (0)  1, so that sn(u, µ2 ) ' sin(u).

In this case, the solutions (1.159)-(1.161) can be written

q
n1 (t) ' n2 (0) sin2 [ n3 (0) Γ t], (1.162)
q
n2 (t) ' n2 (0) cos2 [ n3 (0) Γ t], (1.163)
q
n3 (t) ' n3 (0) − n2 (0) sin2 [ n3 (0) Γ t]. (1.164)

This scenario is plotted in figure 1.14.a. It clearly illustrates stability in the case for which
the frequency of the large amplitude wave (here wave # 3) is lower than at least one of
the frequencies of the two other waves.

56
 n (t)
3
n (t)
Action density n(t)

Action density n(t)

n (t)
2
n (t) n (t)
1 2

0 0
n (t)
3
(a) (b)
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5
Time t, units π / [ Γ n (0)
1/2
] Time t, units 2 F(µ2) / [ Γ (n (0) + n (0))1/2]
3 1 2

Figure 1.14: Time evolution of the action densities for the three non-linearly coupled waves in
the case where the initially finite amplitude wave has a frequency (a) smaller than at least one
of the frequencies of the two other waves, (b) larger than the two other frequencies.

Case 2. n1 (0)  n2 (0) > n3 (0) = 0.

For these initial conditions, the invariants (1.147)–(1.149) are given by

m2 = n1 (0) + n2 (0),
m3 = n1 (0),
C = 0,

and the properly ordered roots na,b,c to Eq. (1.151) by

[na = 0] < [nb = m3 = n1 (0)] < [nc = m2 = n1 (0) + n2 (0)],

so that the parameter µ2 = n1 (0)/[n1 (0) + n2 (0)] ' 1.

In this case, the solutions (1.159)–(1.161) can thus be written

q 
n1 (t) = n1 (0) sn2 n1 (0) + n2 (0) Γ (t − t0 ), µ2 , (1.165)
q 
n2 (t) = n2 (0) + n1 (0) cn2 n1 (0) + n2 (0) Γ (t − t0 ), µ2 , (1.166)
q 
n3 (t) = n1 (0) cn2 n1 (0) + n2 (0) Γ (t − t0 ), µ2 , (1.167)

where t0 = F (µ2 )/Γ [n1 (0) + n2 (0)]1/2 , so that one indeed has n1 (t = 0) = n1 (0).

This scenario is plotted in figure 1.14.b. It corresponds to the unstable case where the
frequency of the large amplitude wave (here wave # 1) is larger than the frequencies of
the two other waves. This result illustrates the most obvious saturation mechanism for
the parametric instability: Depletion of the pump. Indeed, once there is no more energy

57
in the pump wave, the instability clearly stops to grow. The saturation time is given by
t0 ∼ 1/γ = 1/[Γ n1 (0)1/2 ], i.e. scales as the inverse of the growth rate γ obtained from
the linear stability analysis of system (1.137)–(1.139) in the case |a1 |  |a2 |, |a3 |.

The oscillatory nature of solution (1.165)–(1.167) furthermore points out that the de-
cay of quanta from wave #1 into quanta of wave #2 and #3 is in fact reversible.

1.3.6 Saturation Mechanisms

Pump depletion turns out to be the only saturation mechanism included in the simple
fluid model for the three wave interaction considered in the previous sections. In reality
many other saturation mechanisms, involving additional waves and/or kinetic effects,
may potentially play a role. In the case of SRS for example, let us point out the following
saturation processes:
• As the SRS instability develops and the amplitude of the EPW increases, a non-
linear frequency shift of this wave may be induced (see Sec. 1.2.6). This frequency
shift thus leads to a frequency mismatch of the three wave interaction, and, according
to Eq. (1.114), to a less efficient growth of the parametric instability.
• SRS may saturate as a result of the EPW, generated by parametric instability,
undergoing a secondary instability once it reached a certain amplitude threshold:
The so-called Langmuir Decay Instability (LDI, see exercise 1.3.3.2), in which the
primary EPW decays into another EPW and an IAW.
• The EPW may also be subject to another type of secondary instability once its
amplitude is such that the wave has trapped a significant fraction of electrons: The
Trapped Particle Instability (TPI, see Sec. 1.2.7) [17].

1.3.7 Illustrations from Simulations

Bursting behavior of SRS
Raman Amplifier

1.3.8 Further Reading

• Parametric instabilities in laser plasma interaction (LPI): See the excellent
introductory book by Kruer [13].
• Parametric instability in magnetized plasma: See illustration in Sec. 2 of
Chap. I of Ref. [1] treating the case of coupling between Alvén and sound waves.
• Parametric instabilities in inhomogeneous plasmas: Article by Rosenbluth
[14].
• Systematic study of parametric instabilities affecting EM waves in plas-
mas: See the much referenced article by Drake et. al [15].

58
 /RamAmp/RUN3/RES_RamAmp1_full_000
t ω = 2250 t ω = 2700 t ω = 3000
p Incident p p
10 Scattered 10 10
Amplitude

Amplitude

Amplitude
EPW

5 5 5

0 0 0
0 5 10 0 5 10 0 5 10
x/λ 4 x/λ 4 x/λ 4
De x 10 De x 10 De x 10
t ω = 3450 t ω = 3750 t ω = 4200
p p p
10 10 10
Amplitude

Amplitude

Amplitude

5 5 5

0 0 0
0 5 10 0 5 10 0 5 10
x / λDe 4 x / λDe 4 x / λDe 4
x 10 x 10 x 10
t ωp = 4500 t ωp = 4950 t ωp = 5250
10 10 10
Amplitude

Amplitude

Amplitude

5 5 5

0 0 0
0 5 10 0 5 10 0 5 10
x/λ 4 x/λ 4 x/λ 4
De x 10 De x 10 De x 10

59
• Simulations: See for example Vu et. al [16] as an example for PIC-type simula-
tions, and Johnston et. al and Brunner & Valeo [17] as examples of Eulerian-type
calculations.

• Raman amplifiers: For theory see work by Shvets, Fish, Malkin et. al [18, 19],
for simulation results the thesis by Clark [20], and for experimental results see the
work by Ping [21].

60
Appendix A

Particles in a Sinusoidal Potential

A.1 Trapped and Untrapped Particles

One considers here the trajectories of electrons in a sinusoidal wave E(x, t) = E0 sin(k0 x−
ω0 t). One notes that the electrons moving in the field E are separated into two groups:
Passing and trapped. This is clearly seen by working in the frame of reference moving
with the wave, i.e. at velocity vφ = ω0 /k0 with respect to the lab frame. In this wave
frame, the electrostatic field becomes E = E0 sin(k0 x) = −∂φ/∂x, with the potential
φ = (E0 /k0 ) cos(k0 x), so that the total energy of an electron is given by

1 eE0
W = m v2 − cos(k0 x). (A.1)
2 k0
If one assumes the amplitude E0 to be time independent, the energy W is conserved for
each particle. As illustrated in Fig. A.1, electrons with energy levels −eE0 /k0 < W <
eE0 /k0 are trapped, while particles with energy levels W > eE0 /k0 are untrapped. The
terms “untrapped” and “passing” are used interchangeably.

A.2 Deeply Trapped Particles

To start, one considers the case of electrons with energies near Wmin = −eE0 /k0 , which
remain at the bottom of the potential wells, i.e. around the positions xmin = nλ0 =
n 2π/k0 , n integer. These are the so-called deeply trapped particles. By expanding the
potential to second order at the bottom of the well, the total energy (A.1) can be written
1 1
W = m v 2 + ek0 E0 (x − xmin )2 + const,
2 2
which is simply the energy of a harmonic oscillator with frequency ωb = (ek0 E0 /m)1/2 ,
the so-called deeply trapped bounce frequency.

61
 2

Energy W in units e E 0 / k 0
W
untrap
1.5
−e φ(x)
1 W
lim

0.5

0
Wtrap
−0.5

−1
−1 −0.5 0 0.5 1

x1 x2

2
V in units ω b / k 0

untrapped
1
separatrix
trapped 2 ∆ vtrap
0

−1

−2

−3
−1 −0.5 0 0.5 1
X in units λ
0

Figure A.1: a) Sinusoidal potential −eφ(x) in wave frame. b) Untrapped and trapped orbits of
particles.

Period T as a function of energy variable κ

2.5

2
B

Trapped
Period T units 2π /ω

1.5

Untrapped

0.5

0
0 0.5 1 1.5 2 2.5 3 3.5 4
κ

Figure A.2: Transit time and bounce period of particles in a sinusoidal potential as a function
of the energy variable κ.

62
A.3 General Case of Untrapped Particles
One now considers the case of untrapped electrons in a sinusoidal wave with no further
assumption. From Eq. (A.1), one obtains
  1/2
2 eE0 ∆vtrap
v(x, W ) = W+ cos(k0 x) = (1 − κ2 sin2 ξ)1/2 , (A.2)
m k0 κ
having defined ξ = k0 x/2, as well as the transformed energy variable
2eE0
κ2 = . (A.3)
k0 W + eE0
Here one has also made use of the notation ∆vtrap = 2ωb /k0 for the trapping width in
velocity (see Fig. A.1). In terms of κ, the passing condition becomes 0 < κ < 1.

The transit time τt , i.e. the time required for an untrapped particle to cover one pe-
riod λ0 = 2π/k0 , is derived as follows
Z Z
λ0 /2 dx κ 2 π/2 dξ 2κ
τt (κ) = = 2 = F (κ2 ), (A.4)
−λ0 /2 v(x, κ) ∆vtrap k0 −π/2 (1 − κ2 sin ξ) 1/2 ωb
R π/2
where F (m) = 0 dθ/(1 − m sin2 θ)1/2 , 0 < m < 1, is the complete elliptic integral of
the first kind [5]. The transit time τt is plotted as a function of κ in Fig. A.2. Note how
τt becomes infinite in the limit κ → 1−, i.e. for marginally passing particles.

The space averaged velocity is another useful quantity:

Z Z
1 λ0 /2 ∆vtrap 2 π/2
v̄(κ) = < v >x = dx v(x, κ) = dξ (1 − κ2 sin2 ξ)1/2
λ0 −λ0 /2 κ λ0 k 0 −π/2
2 ∆vtrap
= E(κ2 ), (A.5)
π κ
R π/2
where E(m) = 0 dθ(1 − m sin2 θ)1/2 , 0 < m < 1, is the complete elliptic integral of the
second kind [5].

A.4 General Case of Trapped Particles

Here one treats the case of trapped electrons which are not necessarily deeply trapped.
Noting that equation (A.2) is still valid here, and that in terms of κ the trapping condition
becomes κ > 1, the bounce period τb is derived as follows:
Z Z
dx x2 κ 2 ξ2 dξ
τb (κ) = 2 =2
x1 v(x, κ) ∆vtrap k0 ξ1 (1 − κ2 sin2 ξ)1/2
Z
2 π/2 dη 4 1
= 1 2 = F ( ), (A.6)
ωb −π/2 (1 − κ2 sin η)1/2 ω b κ2

63
where the turning points x1,2 = 2ξ1,2 /k0 verify cos(k0 x1,2 ) = −k0 W/eE0 ⇐⇒ κ sin ξ1,2 =
±1. One has also made use here of the change of variables ξ ↔ η defined by κ sin ξ = sin η,
so that κdξ/(1 − κ2 sin2 ξ)1/2 = dη/(1 − sin2 η/κ2 )1/2 . Note that in deriving (A.6) one took
account of both the time for the forward and backward segment of the trapped particle
orbit [factor 2 in first equality of relation(A.6)]. The bounce period τb is also plotted as
a function of κ in Fig. A.2. In the limit κ → 1+, i.e. for marginally trapped particles, τb
becomes infinite.

Deeply trapped particles correspond to κ → ∞, so that, as expected, one recovers

from Eq.(A.6) the bounce period for deeply trapped particles derived previously: τbdeep =
limκ→∞ τb = 2π/ωb , having used F (0) = π/2.

The space averaged velocity is again obtained in a similar way:

1 Z x2 ∆vtrap 2 Z ξ2
v̄(κ) = < v >x = dx v(x, κ) = dξ (1 − κ2 sin2 ξ)1/2
λ0 x1 κ λ0 k 0 ξ 1
cos2 η
Z
∆vtrap π/2
= dη
πκ2 −π/2 (1 − κ12 sin2 η)1/2
 
2∆vtrap 1 1 1
= E( 2 ) + ( 2 − 1)F ( 2 ) , (A.7)
π κ κ κ
having used cos2 η = (1 − κ2 ) + κ2 (1 − sin2 η/κ2 ) for obtaining the final relation.

64
Appendix B

Solving the (Linearized) Vlasov

Equation by Integrating Along
Trajectories

One discusses here the general approach to solving the non-linear Vlasov equation or its
linearized form. One thus considers the Vlasov equation for the distribution f (~x, ~v , t) in
phase space (~x, ~v ) of a given species with charge q and mass m evolving in the electro-
~ x, t), B(~
magnetic fields [E(~ ~ x, t)]:
∂f ∂f q ~ ~ · ∂f = 0.

+ ~v · + E + ~v × B (B.1)
∂t ∂~x m ∂~v
To lighten notations in the following, one writes F~ = q(E
~ + ~v × B)
~ the total force on the
particles.

B.1 Brief Review of the Vlasov Equation

Let us recall, that the Vlasov equation results from the incompressibility of phase space
~ (~z, t) = [~v , F~ (~x, ~v , t)/m], which in general is a 6-dimensional vector in phase space
flux V
~z = (~x, ~v ). This is briefly reviewed here.

Indeed, if the total number of particles remains constant (no sources or sinks for the
considered species, such as ionization or recombination processes), one can write a conti-
nuity equation for the phase space density f (~x, ~v , t) of particles:
∂f ∂ ~
+ ~ f ) = ∂f + ∂ · (~v f ) + ∂ · ( F f ) = 0.
· (V (B.2)
∂t ∂~z ∂t ∂~x ∂~v m
One notes that for the considered electromagnetic forces F~ , related to the macroscopic
~ B),
fields (E, ~ the phase space flux V ~ is such that:
∂ ~ ∂ ∂ q ~ ~ = 0,
·V = · ~v + · (E + ~v × B) (B.3)
∂~z ∂~x
| {z } | ∂~
v m {z }
=0 =0

65
 ~ x, t) is velocity independent, and
having in particular used the fact that E(~

∂ ~ · ( ∂ × ~v ) −~v · ( ∂ × B)
~ =B ~ = 0.
· (~v × B)
∂~v ∂~v
| {z } | ∂~
v {z }
=0 =0

~.
Equation (B.3) describes the incompressibility of phase space flux V

As a result of relation (B.3), the continuity equation (B.2) becomes:

∂f ∂ ~ f ) = ∂f + ( ∂ · V ~ · ∂f = 0.
~ )f + V
+ · (V
∂t ∂~z ∂t | ∂~z{z } ∂~z
=0

This last relation clearly provides Vlasov’s equation (B.1).

Vlasov’s equation in fact states that the distribution f remains invariant along the particle
trajectories in phase space, which naturally again reflects phase space incompressibility.
One can indeed write Eq. (B.1) as

df
= 0, (B.4)
dt traj.

where d/dt|traj is the total time derivative along the trajectories. To convince oneself that
Eq. (B.4) is equivalent to the Vlasov equation, one considers the trajectory [~x 0 (t0 ), ~v 0 (t0 )]
in phase space of any given particle. This trajectory must verify the equations of motion:

d~x 0 (t0 )
0
= ~v 0 (t0 ), (B.5)
dt
d~v 0 (t0 ) q ~ 0 0
0
= F [~x (t ), ~v 0 (t0 ), t0 ]. (B.6)
dt m
The value of the distribution along the particle trajectory, i.e. f [~x 0 (t0 ), ~v 0 (t0 ), t0 ], is func-
tion of the single variable t0 . The variation in time of this function is thus:

d 0 0 0 0 0 ∂f (~x 0 , ~v 0 , t0 ) d~x 0 ∂f (~x 0 , ~v 0 , t0 ) d~v 0 ∂f (~x 0 , ~v 0 , t0 )

f [~
x (t ), ~
v (t ), t ] = + 0 · + 0 ·
dt0  ∂t0 dt ∂~x 0  dt ∂~v 0
∂f ∂f F~ (~x, ~v , t) ∂f 
=  + ~v · + · = 0, (B.7)
∂t ∂~x m ∂~v
(~
x ,~v ,t )
0 0 0

having used Eqs.(B.5)-(B.6) and the Vlasov equation (B.1). This last result clearly proves
Eq. (B.4), and validates the notation for the Vlasov operator:

d ∂ ∂ F~ ∂
= + ~v · + · .
dt traj.
∂t ∂~x m ∂~v

66
B.2 Solving the Vlasov Equation
One considers solving the full Vlasov equation

df ∂f ∂f F~ ∂f
= + ~v · + · = 0,
dt traj.
∂t ∂~x m ∂~v

for the initial condition

f (~x, ~v , 0) = f0 (~x, ~v ).
Formally at least, the solution f (~x, ~v , t) can be written in terms of the particle trajectories,
by invoking the invariance of f along these characteristics. Indeed, by integrating (B.7)
between time t0 = 0 and t0 = t for the trajectory [~x 0 (t0 ; ~x, ~v , t), ~v 0 (t0 ; ~x, ~v , t)] verifying the
particular initial conditions:

~x 0 (t0 = t) = ~x, and ~v 0 (t0 = t) = ~v ,

one simply obtains:

Z t d
f (~x, ~v , t) − f0 [~x 0 (0; ~x, ~v , t), ~v 0 (0; ~x, ~v , t)] = dt0 f [~x 0 (t0 ), ~v 0 (t0 ), t0 ] = 0,
0 dt0
=⇒ f (~x, ~v , t) = f0 [~x 0 (0; ~x, ~v , t), ~v 0 (0; ~x, ~v , t)].
One must however point out here, that computing the trajectories (~x 0 , ~v 0 ) may not be
straightforward, as the electromagnetic forces F~ = q(E+~
~ v ×B)~ determining these trajecto-
ries are themselves function of f [the plasma itself provides sources to the electromagnetic
~ B)].
fields (E, ~ This issue reflects the non-linear nature of the plasma.

B.3 Solving the Linearized Vlasov Equation

The method of integrating along particle trajectories can also be applied when solving
~ 0, B
the linearized Vlasov equation. Let us assume that (f0 , E ~ 0 ) is a known unperturbed
state of the Vlasov-Maxwell system. Note that this state need not be necessarily time
independent. In particular, Vlasov’s equation reads
 
df0 ∂ ∂ F~0 ∂ 
=  + ~v · + · f0 = 0,
dt u.traj.
∂t ∂~x m ∂~v

where (d/dt)|u.traj. stands for the total time derivative along the unperturbed trajectories
in the force field F~0 = q(E
~ 0 + ~v × B
~ 0 ).

~ δ B)
One now considers a perturbation (δf, δ E, ~ to this state. The Vlasov equation for
the full distribution f = f0 + δf reads
 
df ∂ ∂ F~0 + δ F~ ∂ 
=  + ~v · + · (f0 + δf ) = 0, (B.8)
dt traj.
∂t ∂~x m ∂~v

67
where (d/dt)|traj. stands for the total time derivative along the trajectories in the full (i.e.
including perturbations) force field F~ = F~ + δ F~ , with δ F~ = q(δ E
~ + ~v × δ B).
~ Assuming
small amplitude perturbation, one can justify linearizing Eq. (B.8), which leads to the
linearized Vlasov equation for δf :
 
dδf ∂ ∂ F~0 ∂  δ F~ ∂f0
=  + ~v · + · δf = − · , (B.9)
dt u.traj.
∂t ∂~x m ∂~v m ∂~v

noting that on the left hand side of Eq. (B.9) one finds again the total time derivative
along the unperturbed trajectories. Note, that as a result of the non-zero right hand side
term in Eq. (B.9), δf is not invariant along the unperturbed trajectories. Only the full
distribution f = f0 + δf is invariant along the full trajectories as stated by Eq. (B.8).

Equation (B.9) can nonetheless be solved for δf by integrating along trajectories. The
non-zero right hand side provides however an additional contribution. Indeed, by evaluat-
ing Eq. (B.9) along an unperturbed trajectory [~x 0 (t0 ; ~x, ~v , t), ~v 0 (t0 ; ~x, ~v , t)], and integrating
again from time t0 = 0 to t0 = t, one obtains:
Z
0 0 dt
δf (~x, ~v , t) − δf0 [~x (0; ~x, ~v , t), ~v (0; ~x, ~v , t)] = dt0
δf [~x 0 (t0 ), ~v 0 (t0 ), t0 ]
0 dt0
Z t
δ F~ ∂f0
= − dt0 · ,
0 m ∂~v 0 0 0
(~
x ,~v ,t )

δ F~ ∂f0
Z t
=⇒ δf (~x, ~v , t) = δf0 [~x 0 (0; ~x, ~v , t), ~v 0 (0; ~x, ~v , t)] − dt0 · ,
0 m ∂~v
(~
x 0 ,~v 0 ,t0 )

where δf0 (~x, ~v ) = δf (~x, ~v , t = 0) is the initial distribution perturbation, and the unper-
turbed trajectory now verifies

d~x 0 (t0 )
= ~v 0 (t0 ),
dt0
d~v 0 (t0 ) q ~ 0 0
= F0 [~x (t ), ~v 0 (t0 ), t0 ].
dt 0 m
with initial conditions

~x 0 (t0 = t) = ~x, and ~v 0 (t0 = t) = ~v .

68
Bibliography

[1] R. Z. Sagdeev and A. A. Galeev, Nonlinear Plasma Theory (W. A. Benjamin, Inc.,
New York, 1969).

[2] K. Appert, Théorie des Plasmas Chauds (EPFL-Repro, EPFL, 2003).

[3] J. R. Danielson, F. Anderegg, and C. F. Driscoll, Physical Review Letters 92,

#245003 (2004).

[4] L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media, Vol. 8 of

Course of Theoretical Physics (Pergamon Press, Oxford, 1960).

[5] M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions with Formulas,

Graphs, and Mathematical Tables (John Wiley & Sons, New York, 1964).

[6] T. O’Neil, Physics of Fluids 8, 2255 (1965).

[7] I. B. Bernstein, J. M. Greene, and M. D. Kruskal, Physical Review 108, 546 (1957).

[8] G. J. Morales and T. O’Neil, Physical Review Letters 28, 417 (1972).

[9] R. L. Dewar, Physics of Fluids 15, 712 (1972).

[10] W. L. Kruer, J. M. Dawson, and R. N. Sudan, Physical Review Letters 23, 838
(1969).

[11] M. Goldman, Physics of Fluids 13, 1281 (1970).

[12] R. L. Dewar, W. L. Kruer, and W. M. Manheimer, Physical Review Letters 28, 215
(1972).

[13] W. L. Kruer, The Physics of Laser Plasma Interactions, Frontiers in Physics

(Addison-Wesley Publishing Co., Redwood City, California, 1988).

[14] M. N. Rosenbluth, Physical Review Letters 29, 565 (1972).

[15] J. F. Drake et al., Physics of Fluids 17, 778 (1974).

[16] H. X. Vu, D. F. DuBois, and B. Bezzerides, Physical Review Letters 86, 4306 (2001).

[17] S. Brunner and E. J. Valeo, Physical Review Letters 93, 145003 (2004).

69
[18] G. Shvets, N. J. Fisch, A. Pukhov, and J. M. ter Vehn, Physical Review Letters 81,
4879 (1998).

[19] V. M. Malkin, G. Shvets, and N. J. Fisch, Physical Review Letters 82, 4448 (1999).

[20] D. S. Clark, Ph.D. thesis, Princeton University, 2003.

[21] Y. Ping et al., Physical Review Letters 92, #175007 (2004).

70