Schekochihin Houches 2017
Schekochihin Houches 2017
Schekochihin Houches 2017
Alexander A. Schekochihin†
The Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3NP, UK
Merton College, Oxford OX1 4JD, UK
These are the notes for my lectures on Kinetic Theory of Plasmas and on Magnetohy-
drodynamics, taught since 2014 as part of the MMathPhys programme at Oxford. Part
I contains the lectures on plasma kinetics that formed part of the course on Kinetic
Theory, taught jointly with Paul Dellar and James Binney. Part II is an introduction
to magnetohydrodynamics, which was part of the course on Advanced Fluid Dynam-
ics, taught jointly with Paul Dellar. These notes have evolved from two earlier courses:
“Advanced Plasma Theory,” taught as a graduate course at Imperial College in 2008,
and “Magnetohydrodynamics and Turbulence,” taught as a Mathematics Part III course
at Cambridge in 2005-06. I will be grateful for any feedback from students, tutors or
sympathisers.
PART I
Kinetic Theory of Plasmas
1. Kinetic Description of a Plasma
We would like to consider a gas consisting of charged particles—ions and electrons. In
general, there may be many different species of ions, with different masses and charges,
and, of course, only one type of electrons.
We shall index particle species by α (α = e for electrons, α = i for ions). Each is
characterised by its mass mα and charge qα = Zα e, where e is the magnitude of the
electron charge and Zα is a positive or negative integer (e.g., Ze = −1).
1.1. Quasineutrality
We shall always assume that plasma is neutral overall:
X X
qα Nα = eV Zα n̄α = 0, (1.1)
α α
where Nα is the number of the particles of species α, n̄α = Nα /V is their mean number
density and V the volume of the plasma. This condition is known as quasineutrality.
(α)
where by r i we mean the position of the ith particle of species α. We can safely
anticipate that we will only be able to have a nice closed kinetic description if the gas is
approximately ideal, i.e., if particles interact weakly, viz.,
e2
kB T Φ ∼ ∼ e2 n1/3 , (1.3)
∆r
where kB is the Boltzmann constant, which will henceforth be absobed into the tem-
perature T , and ∆r ∼ n−1/3 is the typical interparticle distance. Let us see what this
condition means and implies physically.
Using also the obvious fact that the solution of Eq. (1.5) must be spherically symmetric,
we recast this equation as follows
1 ∂ 2 ∂ϕ 1
2
r − 2 ϕ = −4πqα δ(r). (1.7)
r ∂r ∂r λD
The solution to this that asymptotes to the Coulomb potential ϕ → qα /r as r → 0 and
to zero as r → ∞ is
qα −r/λD
ϕ= e . (1.8)
r
Thus, in a quasineutral plasma, charges are shielded on typical distances ∼ λD .
Obviously, this calculation only makes sense if the “Debye sphere” has many particles
in it, viz., if
nλ3D 1. (1.9)
Let us check that this is the case: indeed,
3/2 3/2
3 T T
nλD ∼ n = 1, (1.10)
ne2 n1/3 e2
provided the weak-interaction condition (1.3) is satisfied. The quantity nλ3D is called the
plasma parameter.
Before we use this approach to construct a description of the plasma as a continuum (on scales
& λD ), let us check that particles travel sufficiently long distances between collisions in order to
4 A. A. Schekochihin
feel the macroscopic fields, viz., that their mean free paths λmfp λD . The mean free path can
be estimated in terms of the collision cross-section σ:
1 T2
λmfp ∼ ∼ (1.13)
nσ ne4
because σ ∼ d2 and the effective distance d by which the particles have to approach each other
in order to have significant Coulomb interaction is inferred by balancing the Coulomb potential
energy with the particle temperature, e2 /d ∼ T . Using Eqs. (1.13) and (1.6), we find
1/2
T 2 ne2
λmfp
∼ ∼ nλ3D 1, q.e.d. (1.14)
λD ne4 T
Thus, it makes sense to talk about a particle travelling long distances experiencing the macro-
scopic fields exerted by the rest of the plasma collectively before being deflected by a much
larger, but also much shorter-range, microscopic field of another individual particle.
(α) (α)
where r i (t) and v i (t) are the exact phase-space coordinates of particle i of species α
at time t, i.e., these are solutions of the exact equations of motion for all these particles
moving in microscopic fields E (micro) (t, r) and B (micro) (t, r). The function Fα is called
the Klimontovich distribution function and it is a random object (i.e., it fluctuates on
scales λD ) because it depends on the exact particle trajectories, which depend on the
exact microscopic fields. In terms of this distribution function,
X Z
σ (micro) (r, t) = qα d3 v Fα (r, v, t), (1.20)
α
X Z
(micro)
j (r, t) = qα d3 v vFα (r, v, t). (1.21)
α
We now need to average these quantities for use in Eqs. (1.15) and (1.18). We shall
assume that the average over microscales (1.12) and the ensemble average (i.e., average
over many different initial conditions) are the same. The ensemble average of Fα is an
object familiar from the kinetic theory of gases, the so-called one-particle distribution
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 5
function:
hFα i = f1α (r, v, t) (1.22)
(we shall henceforth omit the subscript 1). If we learn how to compute fα , then we
can average Eqs. (1.20) and (1.21), substitute into Eqs. (1.15) and (1.18), and have the
following set of macroscopic Maxwell’s equations:
X Z
∇ · E = 4π qα d3 v fα (r, v, t), (1.23)
α
∇ · B = 0, (1.24)
1 ∂B
∇×E+ = 0, (1.25)
c ∂t Z
1 ∂E 4π X
∇×B− = qα d3 v vfα (r, v, t). (1.26)
c ∂t c α
Any other macroscopic force that the plasma might be subject to (e.g., gravity) can
be added to the Lorentz force in the third term on the left-hand side, as long as its
divergence in velocity space is (∂/∂v) · force = 0. Eq. (1.30) is closed by Maxwell’s
equations (1.23–1.26).
mα v 2 ∂fα
XZ
d3 v = 0; (1.33)
α
2 ∂t c
cannot decrease, and, as S is conserved by all the collisionless terms in Eq. (1.30), the
collision operator must have the property that
XZ Z
dS ∂fα
=− d3 r d3 v ln fα > 0, (1.35)
dt α
∂t c
(αα0 ) (αα0 )
where the drag Ai [fα0 ] and diffusion Dij [fα0 ] coefficients are integral (in v space)
functionals of fα0 . The Fokker–Planck form (1.36) of the Landau operator means that it
describes diffusion in velocity space and so will erase sharp gradients in fα with respect
to v—a property that we will find very important in §4.
(α) (α)
First, because r i (t) and v i (t) obviously do not depend on the phase-space variables r and
v, the derivatives ∂/∂r and ∂/∂v can be pulled outside, so the right-hand side of Eq. (1.37) can
be written as a divergence in phase space. Secondly, the particle equations of motion give us
(α)
dr i (t) (α)
= v i (t), (1.38)
dt " #
(α) (α) (α)
B (micro) r i (t), t
dv i (t) qα (α) v i (t) ×
E (micro) r i (t), t +
= , (1.39)
dt mα c
which we substitute into the right-hand side of Eq. (1.37)—after it is written in the divergence
(α) (α)
form. As the time derivatives of r i (t) and v i (t) inside the divergence multiply delta functions
(α) (α) (α) (α)
identifying r i (t) with r and v i (t) with v, we may replace r i (t) with r and v i (t) with
v in the right-hand sides of Eqs. (1.38) and (1.39) when they go into Eq. (1.37). This gives
(wrapping up all the sums of delta functions back into Fα )
v × B (micro) (r, t)
∂Fα ∂ qα
= −∇ · (vFα ) − · E (micro) (r, t) + Fα . (1.40)
∂t ∂v mα c
Finally, because r and v are independent variables and the Lorentz force has zero divergence in
v space, we find that Fα satisfies exactly
v × B (micro)
∂Fα qα ∂Fα
+ v · ∇Fα + E (micro) + · = 0. (1.41)
∂t mα c ∂v
† The simplest example that I can think of in which the collision operator is a velocity-space
diffusion opeartor of this kind is the gas of Brownian particles [each with velocity described by
Langevin’s equation (10.20)]. This is treated in detail in §6.9 of Schekochihin (2016).
8 A. A. Schekochihin
This is the Klimontovich equation. There is no collision integral here because microscopic fields
are explicitly present. The equation is closed by the microscopic Maxwell’s equations:
X Z 3
∇ · E (micro) = 4π qα d v Fα (r, v, t), (1.42)
α
∇ · B (micro) = 0, (1.43)
(micro)
1 ∂B
∇ × E (micro) + = 0, (1.44)
c ∂t
1 ∂E (micro)
Z
(micro) 4π X
∇×B − = qα d3 v vFα (r, v, t). (1.45)
c ∂t c α
Now we separate the microscopic fields into mean (macroscopic) and fluctuating parts ac-
cording to Eq. (1.11) and also
Fα = hFα i + δFα . (1.46)
| {z }
≡ fα
Maxwell’s equations are linear, so averaging them gives the same equations for E and B in
terms of fα [Eqs. (1.23–1.26)] and for δE and δB in terms of δFα . Averaging the Klimontovich
equation (1.45) gives the Vlasov–Landau equation:
∂fα qα v×B ∂fα
+ v · ∇fα + E+ ·
∂t mα c ∂v
qα v × δB ∂δFα ∂fα
=− δE + · ≡ . (1.47)
mα c ∂v ∂t c
The macroscopic fields in the left-hand side satisfy the macroscopic Maxwell’s equations (1.23–
1.26). The microscopic fluctuating fields δE and δB inside the average in the right-hand side
satisfy microscopic Maxwell’s equations with fluctuating charge and current densities expressed
in terms of δFα . Thus, the right-hand side is quadratic in δFα . In order to close this equation, we
need an expression for the correlation function hδFα δFα0 i in terms of fα and fα0 . This is basically
what the BBGKY procedure plus truncation of velocity integrals based on an expansion in 1/nλ3D
achieve. The result is the Landau collision operator (or the more precise Lenard–Balescu one;
see Balescu 1963).
Further details are complicated, but my aim here was just to show how the fields are split into
macroscopic and microscopic ones, with the former appearing explicitly in the kinetic equation
and the latter wrapped up inside the collision operator. The presence of the macroscopic fields
and the consequent necessity for coupling the kinetic equation with Maxwell’s equations for
these fields is the main mathematical difference between the kinetics of neutral gases and the
kinetics of plasmas.
• First, particles are charged, so they interact via Coulomb potential. The collision
operator is, therefore, different: the cross-section is the Rutherford cross-section, most
collisions are glancing (with interaction on distances up to the Debye length), leading to
diffusion of the particle distribution function in velocity space. Mathematically, this is
manifested in the collision operator in (1.30) having the Fokker–Planck structure (1.36).
One can spin out of the Vlasov–Landau equatipn (1.30) a theory that is analogous to
what is done with Boltzmann’s equation in gas kinetics (Dellar 2016): derive fluid equa-
tions, calculate viscosity, thermal conductivity, Ohmic resistivity, etc., of a collisionally
dominated plasma, i.e., of a plasma in which the collision frequency of the particles is
much greater than all other relevant time scales. This is done in the same way as in gas
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 9
kinetics, but now applying the Chapman–Enskog procedure to the Landau collision oper-
ator. This is quite a lot of work—and constitutes core textbook plasma-physics material
(see Parra 2017a). In magnetised plasmas especially, the resulting fluid dynamics of the
plasma are quite interesting and quite different from neutral fluids—we shall see some of
this in Part II of these Lectures, while the classic treatment of the transport theory can
be found in Braginskii (1965); a great textbook on this is Helander & Sigmar (2005).
• Secondly, Coulomb potential is long-range, so the electric and magnetic fields have a
macroscopic (mean) part on scales longer than the Debye length—a particle experiencing
these fields is not undergoing a collision in the sense of bouncing off another particle, but is
rather interacting, via the fields, with the collective of all other particles. Mathematically,
this manifests itself as a Lorentz-force term appearing in the right-hand side of the
Vlasov–Landau kinetic equation (1.30). The macroscopic E and B fields that figure in
it are determined by the particles via their mean charge and current densities that enter
the macroscopic Maxwell’s equations (1.23–1.26).
In the case of neutral gas, all the interesting kinetic physics is in the collision operator,
hence the focus on transport theory in gas-kinetic literature (see, e.g., the classic mono-
graph by Chapman & Cowling 1991 if you want an overdose of this). In the collisionless
limit, the kinetic equation for a neutral gas,
∂f
+ v · ∇f = 0, (1.48)
∂t
simply describes particles with some initial distribution ballistically flying in straight
lines in their initial directions of travel. In contrast, for a plasma, even the collisionless
kinetics (and, indeed, especially the collisionless—or weakly collisional —kinetics) are in-
teresting and nontrivial because, as the initial distribution starts to evolve, it gives rise
to charge densities and currents, which modify E and B, which modify fα , etc. This
opens up a whole new conceptual world and it is on these effects involving interactions
between particles and fields that we shall focus here, in pursuit of maximum novelty.†
We shall also be in pursuit of maximum simplicity (well, “as simple as possible, but
not simpler”!) and so will mostly restrict our consideration to the “electrostatic approx-
imation”:
B = 0, E = −∇ϕ. (1.49)
This, of course, eliminates a huge number of interesting and important phenomena with-
out which plasma physics would not be the voluminous subject that it is, but we cannot
do them justice in just a few lectures (most of Parra 2017b is devoted to collisionless
magnetised plasmas).
Thus, we shall henceforth focus on a simplified kinetic system, called the Vlasov–
Poisson system:
∂fα qα ∂fα
+ v · ∇fα − (∇ϕ) · = 0, (1.50)
∂t mα ∂v
X Z
−∇2 ϕ = 4π qα d3 v fα . (1.51)
α
† Similarly interesting things happen when the field tying the particles together is gravi-
ty—an even more complicated situation because, while the potential is long-range, rather like
the Coulomb potential, gravity is not shielded and so all particles feel each other at all distances.
This gives rise to remarkably interesting theory (Binney 2016).
10 A. A. Schekochihin
and the Vlasov equation (1.50) written for k = 0 (i.e., the spatial average of the equa-
tion) is
∂f0 ∂δfk=0 q X ∂δfk
+ =− ϕ−k ik · . (2.10)
∂t ∂t m ∂v
k
Averaging also over time according to Eq. (2.7) eliminates fast variation and gives us
∂f0 q X ∗ ∂δfk
=− ϕk ik · (2.11)
∂t m ∂v
k
The three terms that control the evolution of the perturbed distribution function in
Eq. (2.12) represent the three physical effects that we shall focus on in these Lectures. The
second term on the left-hand side represents free ballistic motion of particles (“stream-
ing”). It will give rise to the phenomenon of phase mixing (§4) and, in its interplay with
plasma waves, to Landau damping and kinetic instabilities (§3). The first term on the
right-hand side contains the interaction of the electric-field perturbations (waves) with
the equilibrium particle distribution. The second term on the right-hand side of Eq. (2.12)
has nonlinear interactions between fluctuating fields and the perturbed distribution—it
is negligible when fluctuation amplitudes are small enough (which, sadly, they rarely are)
and responsible for plasma turbulence when they are not.
The programme for determining the slow evolution of the equilibrium is “simple”: solve
Eq. (2.12) together with Eq. (2.9), calculate the correlation function of the fluctuations,
hϕ∗k δfk i, as a functional of f0 , and use it to close Eq. (2.11); then proceed to solve the
latter. Obviously, this is impossible to do in most cases. But it is possible to construct
a hierarchy of approximations to the answer and learn much interesting physics in the
process.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 13
2.4. Hierarchy of Approximations
2.4.1. Linear Theory
Consider first infinitesimal perturbations of the equilibrium. All nonlinear terms can
then be ignored, Eq. (2.11) turns into f0 = const and Eq. (2.12) becomes
∂δfk q ∂f0
+ ik · v δfk = ϕk ik · , (2.13)
∂t m ∂v
the linearised kinetic equation. Solving this together with Eq. (2.9) allows one to find
oscillating and/or growing/decaying† modes of the plasma perturbing a particular equi-
librium. The theory for doing this is very well developed and contains some of the core
ideas that give plasma physics its intellectual shape (§3).
Physically, the linear solutions will describe what happens over short term, viz., on
times t such that
ω −1 t teq or tnl , (2.14)
where ω is the characteristic frequency of the perturbations, teq is the time over which the
equilibrium starts getting modified by the perturbations (which depends on the amplitude
to which they can grow: see the right-hand side of Eq. (2.11); if perturbations grow, they
can modify the equilibrium by this mechanism so as to render it stable), and tnl is the
time at which perturbation amplitudes become large enough for nonlinear interactions
between individual modes to matter (second term on the right-hand side of Eq. (2.12);
if perturbations grow, they can saturate by this mechanism).
For compactness of notation, we will drop both the species index α and the wave number
k in the subscripts, unless they are necessary for understanding.
We will discover that that electrostatic perturbations in a plasma described by Eqs. (3.1)
and (3.2) oscillate, can pass their energy to particles (damp) or even grow, sucking energy
from the particles. We will also discover that it is useful to know some complex analysis.
Figure 3. Layout of the complex-p plane: δfˆ(p) is analytic for Re p > σ. At Re, p < σ, δfˆ(p)
may have singularities (poles).
It is a mathematical certainty that if there exists a real number σ > 0 such that
then the integral (3.4) exists (i.e., is finite) for all values of p such that Re p > σ. The
inverse Laplace transform, giving us back our distribution function as a function of time,
is then
Z i∞+σ
1
δf (t) = dp ept δfˆ(p), (3.6)
2πi −i∞+σ
where the integral is along a straight line in the complex plane parallel to the imaginary
axis and intersecting the real axis at Re p = σ (Fig. 3).
Since we expect to be able to recover our desired time-dependent function δf (v, t) from
its Laplace transform, it is worth knowing the latter. To find it, we Laplace-transform
Eq. (3.1):
Z ∞ Z ∞
∂δf ∞
dt e−pt = e−pt δf 0 + p dt e−pt δf = −g + p δfˆ,
l.h.s. =
0 ∂t 0
q ∂f0
r.h.s. = −ik · v δfˆ + ϕ̂ ik · . (3.7)
m ∂v
Equating these two expressions, we find the solution:
1 q ∂f0
δfˆ(p) = i ϕ̂(p) k · +g . (3.8)
p + ik · v m ∂v
16 A. A. Schekochihin
The Laplace transform of the potential, ϕ̂(p), itself depends on δfˆ via Eq. (3.2):
Z ∞ Z
4π X
ϕ̂(p) = dt e−pt ϕ(t) = 2 qα d3 v δfˆα (p)
0 k α
Z
4π X 3 1 qα ∂f0α
= 2 qα d v i ϕ̂(p) k · + gα . (3.9)
k α p + ik · v mα ∂v
This is an algebraic equation for ϕ̂(p). Collecting terms, we get
" #
X 4πq 2 Z 1 ∂f0α 4π X
Z
gα
α 3
1− 2m
i d v k · ϕ̂(p) = 2
q α d3 v . (3.10)
α
k α p + ik · v ∂v k α
p + ik · v
| {z }
≡ (p, k)
The prefactor in the left-hand side, which we denote (p, k), is called the dielectric func-
tion, because it encodes all the self-consistent charge-density perturbations that plasma
sets up in response to an electric field. This is going to be an important function, so let
us write it out beautifully:
X ωpα 2 Z
i 1 ∂f0α
(p, k) = 1 − 2 n
d3 v k· , (3.11)
α
k α p + ik · v ∂v
2 4πqα2 nα
ωpα = . (3.12)
mα
The solution of Eq. (3.10) for the potential is
4π X Z gα
ϕ̂(p) = 2 qα d3 v . (3.13)
k (p, k) α p + ik · v
To calculate ϕ(t), we need to inverse-Laplace-transform ϕ̂: similarly to Eq. (3.6),
Z i∞+σ
1
ϕ(t) = dp ept ϕ̂(p). (3.14)
2πi −i∞+σ
How do we do this integral? Recall that δfˆ and, therefore, ϕ̂ only exists (i.e., is finite)
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 17
for Re p > σ, whereas at Re p < σ, it can have singularities, i.e., poles—let us call them
pi , indexed by i. If we analytically continue ϕ̂(p) everywhere to Re p < σ except those
poles, the result must have the form
X ci
ϕ̂(p) = + A(p), (3.15)
i
p − pi
where ci are some coefficients (residues) and A(p) is the analytic part of the solution. The
integration contour in Eq. (3.14) can be shifted to Re p → −∞ but with the proviso that
it cannot cross the poles, as shown in Fig. 4).† Then the contributions to the integral
from the vertical segments of the contour are exponentially small, the contributions from
the segments leading towards and away from the poles cancel, and the contributions from
the circles around the poles can, by Cauchy’s formula, be expressed in terms of the poles
and residues:
X
ϕ(t) = ci epi t . (3.16)
i
Thus, initial perturbations of the potential will evolve ∝ epi t , where pi are poles of ϕ̂(p).
In general, pi = −iωi + γi , where ωi is a real frequency (giving wave-like behaviour
of perturbations), γi < 0 represents damping and γi > 0 growth of the perturbations
(instability).
Going back to Eq. (3.13), we realise that the poles of ϕ̂(p) are zeros of the dielectric
function:
(pi , k) = 0 ⇒ pi = pi (k) = −iωi (k) + γi (k). (3.17)
To find the wave frequencies ωi and the damping/growth rates γi , we must solve this
equation, which is called the plasma dispersion relation.
We are not particularly interested in what ci ’s are, but they can be computed from the initial
conditions, also via Eq. (3.13). The reason we are not particularly interested is that if we set
up an initial perturbation with a given k and then wait long enough, only the fastest-growing
or, failing growth, the slowest-damped mode will survive, with all others having exponentially
small amplitudes. Thus, a typical outcome of the linear theory is a mode oscillating at some
frequency and growing or decaying at some rate. Since this is a linear solution, the prefactor in
front of the exponential can be scaled arbitrarily and so does not matter.
† This is proven by making a closed loop out of the old and the new contours, joining them
at ±i∞, and noting that this loop contains no poles.
18 A. A. Schekochihin
Let us prove that this is indeed an analytic continuation, i.e., that the integral (3.18), adjusted
to be along CL , is analytic for all values of p. Let us cut the Landau contour at vz = ±R and
close it in the upper half-plane with a semicircle CR of radius R such that R > σ/k (Fig. 6).
Then, with integration running along the truncated CL and counterclockwise along CR , we get,
‡ NB: This means that in what follows, k > 0 by definition.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 19
by Cauchy’s formula,
F 0 (vz ) F 0 (vz )
Z Z
ip
dvz + dvz = 2πi F 0 . (3.19)
CL vz − ip/k CR vz − ip/k k
Since analyticity is guaranteed for Re p > σ, the integral along CR is analytic. The right-hand
side is also analytic, by assumption. Therefore, the integral along CL is analytic—this is the
integral along the Landau contour if we take R → ∞.
where the integrals are again over the real axis and the imaginary bits come from the
contour making a half (when Re p = 0) or a full (when Re p < 0) circle around the pole.
In the case of Re p = 0, or ip = ω, the integral along the real axis is formally divergent
and so we take its principal value, defined as
Z +∞ "Z Z +∞ #
ω/k−ε
F 0 (vz ) F 0 (vz )
P dvz = lim + dvz . (3.21)
−∞ vz − ω/k ε→0 −∞ ω/k+ε vz − ω/k
The difference between Eq. (3.21) and the usual Lebesgue definition of an integral is that the
latter would be
" Z +∞ #
F 0 (vz ) F 0 (vz )
Z +∞ Z ω/k−ε1
dvz = lim + lim dvz , (3.22)
−∞ vz − ω/k ε1 →0 −∞ ε2 →0 ω/k+ε
2
vz − ω/k
and this, with, in general, ε1 6= ε2 , diverges logarithmically, whereas in Eq. (3.21), the divergences
neatly cancel.
The Re p = 0 case in Eq. (3.20),
F 0 (vz ) F 0 (vz )
Z Z +∞ ω
dvz =P dvz + iπF 0 , (3.23)
CL vz − ω/k −∞ vz − ω/k k
which tends to be of most use in analytical theory, is a particular instance of Plemelj’s formula:
for a real ζ and a well-behaved function f (no poles on or near the real axis),
Z +∞ Z +∞
f (x) f (x)
lim dx =P dx ± iπf (ζ), (3.24)
ε→+0 −∞ x − ζ ∓ iε −∞ x−ζ
also sometimes written as
1 1
lim =P ± iπδ(x − ζ), (3.25)
ε→+0 x − ζ ∓ iε x−ζ
Finally, armed with Landau’s prescription, we are ready to calculate. Eq. (3.11) becomes
2
X ωpα Fα0 (vz )
Z
1
(p, k) = 1 − dv z (3.26)
α
k 2 nα CL vz − ip/k
20 A. A. Schekochihin
and, analogously, Eq. (3.13) becomes
Z
4πi X Gα (vz )
ϕ̂(p) = − 3 qα dvz , (3.27)
k (p, k) α CL vz − ip/k
R R
where Gα (vz ) = dvx dvy gα (vx , vy , vz ).
Setting the real part of Eq. (3.28) to zero gives the equation for the real frequency:
Re (−iω, k) = 0 . (3.30)
Thus, we now only need (p, k) with p = −iω. Using Eq. (3.23), we get
2
X ωpα F 0 (vz )
Z
1
Re = 1 − 2
P dvz α , (3.31)
α
k nα vz − ω/k
X ωpα 2
π 0 ω
Im = − F . (3.32)
α
k 2 nα α k
Let us now consider a two-species plasma, consisting of electrons and a single species
of ions. There will be two interesting limits: p
• “High-frequency” electron waves: ω kvthe , where vthe = 2Te /me is the “thermal
speed” of the electrons†; this limit will give us Langmuir waves (§3.4), slowly damped or
growing (§3.5).
• “Low-frequency” ion waves: a particularly tractable limit will
p be that of “hot” elec-
trons and “cold” ions, viz., kvthe ω kvthi , where vthi = 2Ti /mi is the “thermal
speed” of the ions; this limit will give us the sound (“ion-acoustic waves”; §3.8), which
also can be damped or growing (§3.9).
√ p
where λDe = vthe / 2 ωpe = Te /4πe2 ne is the “electron Debye length” [cf. Eq. (1.6)].
Eq. (3.39) is the Bohm & Gross (1949a) dispersion relation, describing an upgrade of
p
† For hydrogen plasma, mi /me ≈ 42, more or less, the answer to the Ultimate Question of
Life, Universe and Everything (Adams 1979).
22 A. A. Schekochihin
the Langmuir oscillations to dispersive Langmuir waves, which have a non-zero group
velocity (this effect is due to electron pressure: see Exercise 3.1).
Note that all this is only valid for ω kvthe , which we now see is equivalent to
kλDe 1. (3.40)
Exercise 3.1. Langmuir hydrodynamics. Starting from the linearised kinetic equation for
electrons and ignoring perturbations of the ion distribution functions completely, work out the
fluid equations for electrons (i.e., evolution equations for the electron density ne and velocity
ue ) and show that you can recover the Langmuir waves (3.39) if you assume that electrons
behave as a 1D adiabatic fluid (i.e., have the equation of state pe n−γ
e = const with γ = 3). You
can prove that they indeed do this by calculating their density and pressure directly from the
Landau solution for the perturbed distribution function (see §§4.3 and 4.6), ignoring resonant
particles. The “hydrodynamic” description of Langmuir waves will reappear in §9.
where ω is given by Eq. (3.39). Provided ωF 0 (ω/k) < 0 (as would be the case, e.g., for
any distribution monotonically decreasing with |vz |; see Fig. 7a), γ < 0 and so this is
indeed a damping rate, the celebrated Landau damping (Landau 1946; it was confirmed
experimentally two decades later, by Malmberg & Wharton 1964).
The same theory also describes a class of kinetic instabilities: if ωF 0 (ω/k) > 0, γ > 0,
so perturbations grow exponentially with time. An iconic example is the bump-on-tail
instability (Fig. 7b), which arises when a high-energy (vz vthe ) electron beam is
injected into a plasma and which we will study in great detail in §7.
We see that the damping or growth of plasma waves occurs via their interaction with
the particles whose velocities coincide with the phase velocity of the wave (“Landau
resonance”). Because such particles are moving in phase with the wave, its electric field
is stationary in their reference frame and so can do work on these particles, giving its
energy to them (damping) or receiving energy from them (instability). In contrast, other,
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 23
out-of-phase, particles experience no mean energy change over time because the field that
they “see” is oscillating. It turns out (§3.6) that the process works in the spirit of socialist
redistrbution: the particles that are slightly lagging behind the wave will, on average,
receive energy from it, damping the wave, whereas those overtaking the wave will give
up their energy, amplifying the wave. The condition ωF 0 (ω/k) < 0 corresponds to the
stragglers being more numerous than the strivers, leading to net damping; ωF 0 (ω/k) > 0
implies the opposite, leading to an instability (which then leads to a flattening of the
distribution; see §7).
Let us note again that these results are quantitatively valid only in the limit (3.33),
or, equivalently, (3.40). It makes sense that damping should be slow (γ ω) in the
limit where the waves propagate much faster than the majority of the electrons (ω/k
vthe ) and so can interact only with a small number of particularly fast particles (for a
2 2 2
Maxwellian equilibrium distribution, it is an exponentially small number ∼ e−ω /k vthe ).
If, on the other hand, ω/k ∼ vthe , the waves interact with the majority population and
the damping should be strong: a priori, we might expect γ ∼ kvthe .
A few years ago, Landau damping became a cause célèbre in the mathematics community and
in the wider science world with the award of the Fields Medal in 2010 to Cédric Villani, who
proved (with C. Mouhot) that, basically, Landau’s solution of the linearised Vlasov equation
survived as a solution of the full nonlinear Vlasov equation for small enough and regular enough
initial perturbations: see a “popular” account of this by Villani (2014). The regularity restriction
is apparently important and the result can break down in interesting ways: see Bedrossian
(2016). The culprit is plasma echo, of which more will be said in §§6.2 and 10 (without claim to
mathematical rigour; see also Schekochihin et al. 2016 and Adkins & Schekochihin 2017).
Exercise 3.2. Stability of isotropic distributions. Prove that if f0e (vx , vy , vz ) = f0e (v),
i.e., is 3D-isotropic distribution, monotonic or otherwise, the Langmuir waves at kλDe 1 are
always damped (this is solved in Lifshitz & Pitaevskii 1981; the statement of stability of isotropic
distributions is in fact valid much more generally than just for long-wavelength Langmuir waves,
as we will see in §5).
The work done by the field on the electron per unit time, averaged over time, is the power gained
by the electron:
δP (v0 ) = −e hE(z(t), t)vz (t)i
∂E
≈ −e E(v0 t, t) + δz(t) (v0 t, t) [v0 + δvz (t)]
∂z
* +
εt
= −eE0 e v0 cos[(ω − kv0 )t] + δvz (t) cos[(ω − kv0 )t] + v0 δz(t)k sin[(ω − kv0 )t]
| {z } | {z } | {z }
vanishes only cos term only sin term
under from from
averaging Eq. (3.47) Eq. (3.48)
survives survives
averaging averaging
e2 E02 2εt
ε 2kv0 ε(ω − kv0 )
= e +
2me (ω − kv0 )2 + ε2 [(ω − kv0 )2 + ε2 ]2
e2 E02 2εt d εv0
= e . (3.49)
2me dv0 (ω − kv0 )2 + ε2
| {z }
≡ χ(v0 )
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 25
† The fact that we are working in 1D means that we are restricting our consideration to
perturbations whose wave numbers k are parallel to the beam’s velocity. In general, allowing
transverse wave numbers brings into play the transverse (electromagnetic) part of the dielectric
tensor (see Q2). However, for non-relativistic beams, the fastest-growing modes will still be the
longitudinal, electrostatic ones (see, e.g., Alexandrov et al. 1984, §32).
26 A. A. Schekochihin
Figure 10. Sketch of the growth rate of the hydrodynamic and kinetic beam instabilities:
Eq. (3.57) for k < ωpe /ub , Eq. (3.58) for k = ωpe /ub , and Eq. (3.41) for ωpe /ub < k < ωpe /v0 ,
where v0 is the point of the minimum of the distribution in Fig. 7(b) and ub is the point of the
maximum of the bump.
of a thermal spread around ub . It turns out that if the width of the beam is sufficiently small,
another instability appears, whose origin is hydrodynamic rather than kinetic. In the interest of
having a full picture, let us work it out.
Consider a very simple limiting case of the distribution (3.53): let vb → 0 and nb ne . Then
(Fig. 9)
Fe (vz ) = FM (vz ) + nb δ(vz − ub ), (3.54)
where FM is the bulk Maxwellian from (3.53) (with density ≈ ne , neglecting nb in comparison).
Let us substitute the distribution (3.54) into the dielectric function (3.26), seek solutions with
p/k vthe , expand the part containing FM in the same way as we did in §3.4†, neglect the
ion contribution for the same reason as we did there, and deal with the term in the dielectric
function containing δ 0 (vz − ub ) by integrating by parts. The resulting dispersion relation is
2 2
ωpe nb ωpe
≈1+ − = 0. (3.55)
p2 ne (kub − ip)2
Since nb ne , the last term can only matter for those perturbations that are close to resonance
with the beam (this is called Cherenkov’s resonance):
p = −ikub + γ, γ kub . (3.56)
This turns Eq. (3.55) into
2 2 r −1/2
ωpe nb ωpe nb 1 1
1− + =0 ⇒ γ=± − . (3.57)
k2 u2b ne γ 2 ne k2 u2b 2
ωpe
The expression under the square root is positive and so there is indeed a growing mode only if
k < ωpe /ub . This is contrast to the case of a hot (or warm) beam in §3.5: there, having a kinetic
instability required ωFe0 (ω/k) > 0, which was only possible at k > ωpe /ub (the phase speed
of the perturbations had to be to the left of the bump’s maximum). The new instability that
we have found—the hydrodynamic beam instability—has the largest growth rate at kub = ωpe ,
i.e., when the beam and the plasma oscillations are in resonance, in which case, to resolve the
singularity, we need to retain γ in the second term in Eq. (3.55). Doing so and expanding in γ,
we get
2 2 2 √ 1/3
ωpe nb ωpe 2iγ nb ωpe ± 3+i nb
≈ 1− 2 2
+ 2
≈ + =0 ⇒ γ= , −i ωpe .
(ωpe + iγ) ne γ ωpe ne γ 2 2 2ne
(3.58)
† We can treat the Landau contour as simply running along the real axis because we are
expecting to find a solution with Re p > 0 [see Eq. (3.20)], for reasons independent of the
Landau resonance.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 27
(a) Cold streams, Eq. (3.59) (b) Hot streams, Eq. (5.16)
Figure 11. Two streams.
The unstable root (Re γ > 0) is the interesting one. The growth rate of the combined beam
instability, hydrodynamic and kinetic, is sketched in Fig. 10.
Exercise 3.3. This instability is called “hydrodynamic” because it can be derived from fluid
equations (cf. Exercise 3.1) describing cold electrons (vthe = 0) and a cold beam (vb = 0).
Convince yourself that this is the case.
Exercise 3.4. Using the model distribution (3.53), work out the conditions on vb and nb that
must be satisfied in order for our derivation of the hydrodynamic beam instability to be valid,
i.e., for Eq. (3.55) to be a good approximation to the true dispersion relation. What is the effect
of finite vb on the hydrodynamic instability? Sketch the growth rate of unstable perturbations
as a function of k, taking into account both the hydrodynamic instability and the kinetic one,
as well as the Landau damping.
Exercise 3.5. Two-stream instability. This is a popular instability† that arises, e.g., in a
situation where the plasma consists of two cold counterstreams of electrons propagating against
a quasineutrality-enforcing background of effectively immobile ions (Fig. 11a). Model the corre-
sponding electron distribution by
ne
Fe (vz ) = [δ(vz − ub ) + δ(vz + ub )] (3.59)
2
and solve the resulting dispersion relation (where the ion terms can be neglected for the same
reason as in §3.4). Find the wave number at which perturbations grow fastest and the corre-
sponding growth rate. What happens, qualitatively, if the two streams are warm (have a finite
thermal spread vb ; Fig. 11b)? A nice fully tractable quantitative model of the latter situation is
the double-Lorentzian distribution (5.16). Solve the dispersion relation for it to verify that your
qualitative expectations are correct.
where
r
ZTe
cs = ωpi λDe = (3.68)
mi
is the sound speed, called that because, if we take kλDe 1, Eq. (3.67) describes a wave
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 29
Figure 12. Ion-acoustic resonance: damping (cs > ue ) or instability (cs < ue ). Ion Landau
damping is weak because cs vthi , so in the tail of Fi (vz ); electron damping/instability are
also weak because ue , cs vthe , so close to the peak Fe (vz ).
Exercise 3.7. Find the ion contribution to the damping of ion-acoustic waves. Under what
conditions does it become comparable to, or larger than, the electron contribution?
Exercise 3.8. Derive the dispersion relation for ion Langmuir waves. Investigate their damp-
ing/instability.
plasma. These are summarised in Fig. 13. In the limit of short wavelengths, kλDe 1 and
kλDi 1, the electron and ion branches, respectively, becomes dispersive, their damping
rates increase and eventually stop being small. This corresponds to waves having phase
speeds that are comparable to the speeds of the particles in the thermal bulk of their
distributions, so a great number of particles are available to have Landau resonance with
the waves and absorb their energy—the damping becomes strong.
Note that if the cold-ion condition Ti Te is not satisfied, the sound speed is com-
parable to the ion thermal speed cs ∼ vthi , and so the ion-acoustic waves are strongly
damped at all wave numbers—it is well-nigh impossible to propagate sound through a
collisionless hot plasma (no one will hear you scream)!
and similar dispersion relations arising from, e.g., considering electromagnetic perturbations (see
Q2), magnetised plasmas (see Parra 2017b), different equilibria Fα (see Q3), etc.
A Maxwellian equilibrium is obviously an extremely important special case because that is,
after all, the distribution towards which plasma is pushed by collisions on long time scales:
nα 2 2 nα 2 2
f0α (v) = 2
e−v /vthα ⇒ Fα (vz ) = √ e−vz /vthα . (3.77)
(πvthα )3/2 πvthα
For this case, we would like to introduce a new “special function” that would incorporate the
Landau prescription for calculating the velocity integral in Eq. (3.76) and that we could in some
sense “tabulate” once and for all.†
† In the olden days, one would literally tabulate it (Fried & Conte 1961). In the 21st century,
we could just teach a computer to compute it [see Eq. (3.86)] and have an app.
32 A. A. Schekochihin
Taking Fα to be a Maxwellian and letting u = vz /vthα and ζα = ip/kvthα , we can rewrite the
velocity integral in Eq. (3.76) as follows
2
Fα0 (vz ) u e−u
Z Z
1 2 2
dvz = −√ 2 du = − 2 [1 + ζα Z(ζα )] , (3.78)
nα CL vz − ip/k π vthα CL u − ζα vthα
where the plasma dispersion function is defined to be
2
e−u
Z
1
Z(ζ) = √ du . (3.79)
π u−ζ
Note that if the Maxwellian distribution (3.77) has a mean flow, as it did, e.g., in Eq. (3.61),
this amounts to a shift by some mean velocity uα and all one needs to do to adjust the above
results is to shift the argument of Z accordingly:
uα
ζα → ζα − . (3.81)
vthα
This turns out to be a uniformly valid expression for Z(ζ): our function is simply a complex erf!
Here is a Mathematica script for calculating it:
Z[zeta ] := I Sqrt[Pi] Exp[−zeta2 ](1 + I Erfi[zeta]). (3.86)
You can use this to code up Eq. (3.80) and explore, e.g., the strongly damped solutions (ζ ∼ 1,
γ ∼ ω).
√ 2ζ 2 4ζ 4 8ζ 6
2
Z(ζ) = i π e−ζ − 2ζ 1 − + − + ... . (3.87)
3 15 105
√
2 1 1 3 15
Z(ζ) = i π e−ζ − 1 + 2 + 4 + 6 + ... . (3.88)
ζ 2ζ 4ζ 8ζ
All the results (for a Maxwellian equilibrium) that I derived in §§3.4–3.10 can be readily obtained
from Eq. (3.80) by using the above limiting cases (see Q1). It is, indeed, a general practical
strategy for studying this and similar plasma dispersion relations to look for solutions in the
limits ζα → 0 or ζα → ∞, then check under what physical conditions the solutions thus obtained
are valid (i.e., they indeed satisfy |ζα | 1 or |ζα | 1, |Re ζ| |Im ζ|), and then fill in the
non-asymptotic blanks in the same way that an experienced hunter espying antlers sticking out
above the shrubbery can reconstruct, in contour outline, the rest of the hiding deer.
Exercise 3.9. Work out the Taylor series (3.87). A useful step might be to prove this interesting
formula (which also turns out to be handy in other calculations; see, e.g., Q7):
2
dm Z (−1)m Hm (u) e−u
Z
m
= √ du , (3.89)
dζ π CL u−ζ
Exercise 3.10. Work out the asymptotic series (3.88) using the Landau prescription (3.20)
and expanding the principal-value integral similarly to the way it was done in Eq. (3.35). Work
out also (or look up; e.g., Fried & Conte 1961) other asymptotic forms of Z(ζ), relaxing the
condition |Re ζ| |Im ζ|.
where j is the current density. So the rate of change of the electric field is minus the rate
at which electric field does work on the charges, a.k.a. Joule heating—not a surprising
result. The energy goes into accelerating particles, of course: the rate of change of their
kinetic energy is
dK X ZZ mα v 2 ∂fα
= d3 r d3 v
dt α | 2 {z ∂t}
use
Eq. (4.1)
m2α
X ZZ
3 3 qα ∂fα ∂fα
= d rd v −v · ∇fα + (∇ϕ) · +
2 mα ∂v ∂t
α | {z } | {z } | {z }c
vanishes by parts in v vanishes
because because
full energy is
divergence conserved
X ZZ Z
3 3
=− qα d r d v fα v · ∇ϕ = d3 r E · j. (4.4)
α
Combining Eqs. (4.3) and (4.4) gives us the law of energy conservation:
2
Z
d 3 E
K+ d r =0 . (4.5)
dt 8π
Exercise 4.1. Demonstrate energy conservation for the more general case in which magnetic-
field perturbations are also allowed.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 35
Thus, if the perturbations are damped, the energy of the particles must increase—
this is usually called heating. Strictly speaking, heating refers to a slow increase in the
mean temperature of the thermal equilibrium. Let us make this statement quantitative.
Consider a Maxwellian plasma, homogeneous in space but possibly with some slow de-
pendence on time (cf. §2):
3/2
nα −v 2 /vthα
2 mα 2
f0α = 2 3/2
e = n α e−mα v /2Tα . (4.6)
(πvthα ) 2πTα
In a homogeneous system with a fixed volume, the density nα is constant in time be-
cause the number of particles is constant: d(V nα )/dt = 0. We allow, however, that the
temperature is Tα = Tα (t). The total kinetic energy of the particles is
XZ mα v 2 X ZZ mα v 2
K=V d3 v f0α + d3 r d3 v δfα . (4.7)
α |
2 α
2
{z }
3
= nα Tα
2
Let us average this over time, as per Eq. (2.7): the perturbed part vanishes and we have
X3
hKi = V nα Tα . (4.8)
α
2
The Landau damping rate of the electric-field perturbation is the heating rate of the
equilibrium.†
This result, while at first glance utterly obvious, might, on reflection, appear to be
paradoxical: surely, the heating of the equilibrium implies increasing entropy—but the
damping that is leading to the heating is collisionless and, in a collisionless system, in
view of the H theorem, how can the entropy change?
4.2. Entropy and Free Energy
The kinetic entropy for each species of particles is defined to be
ZZ
Sα = − d3 r d3 v fα ln fα . (4.11)
† Obviously, the damping of waves on particles of species α increases only the temperature
of that species.
36 A. A. Schekochihin
This quantity (or, indeed, the full-phase-space integral of any quantity that is a function
only of fα ) can only be changed by collisions and, furthermore, the plasma-physics version
of Boltzmann’s H theorem establishes that
X ZZ
d X 3 3 ∂fα
Sα = − d rd v ln fα > 0, (4.12)
dt α α
∂t c
where equality is achieved iff all fα are Maxwellian with the same temperature Tα = T .
Thus, if collisions are ignored, the total entropy stays constant and everything that
happens is, in principle, reversible. So how can there be net damping of waves and,
worse still, net heating of the equilibrium particle distribution?! Presumably, any damping
solution can be turned into a growing solution by reversing all particle trajectories—so
should the overall perturbation level stay constant?
As we already noted in §4.1, strictly speaking, heating is the change in the equilibrium
temperature—and, therefore, equilibrium entropy. Indeed, for each species, the equilib-
rium entropy is
mv 2
ZZ ZZ
3 3 3 3 m 3/2 3
S0 = − d r d v f0 ln f0 = − d r d v f0 ln n − ln T −
2π 2 2T
m 3/2 3
3
= V −n ln n + n ln T + n , (4.13)
2π 2 2
where we have used d3 v (mv 2 /2)f0 = (3/2)nT . Since n = const,
R
dS0 3 dT
T =V n , (4.14)
dt 2 dt
so heating is indeed associated with the increase of S0 .
Since, according to Eq. (4.10), this can be successfully accomplished by collisionless
damping and since entropy overall can only increase due to collisions, we must search for
the “missing entropy” (or, rather, for the missing decrease in entropy) in the perturbed
part of the distribution. The mean entropy associated with the perturbed distribution is
ZZ
hδSi = hS − S0 i = − d3 r d3 v h(f0 + δf ) ln(f0 + δf ) − f0 ln f0 i
δf 2
ZZ
3 3 δf
=− d r d v (f0 + δf ) ln f0 + − + . . . − f0 ln f0
f0 2f0
hδf 2 i
ZZ
≈− d3 r d3 v , (4.15)
2f0
after expanding to second order in small δf /f0 and using hδf i = 0. The total entropy
of each species, S = S0 + δS, can only by changed by collisions, so, if collisions are
ignored, any heating of a given species, i.e., an increase in its S0 [see Eq. (4.14)] must be
compensated by a decrease in its δS. The latter can only be achieved by increasing the
mean square amplitude of the perturbations of the distribution function: indeed, using
Eqs. (4.14) and (4.15), we find†
T hδf 2 i
ZZ
dS0 dhδSi 3 dT d
T + =V n − d3 r d3 v
dt dt 2 dt dt 2f0
ZZ
∂f
=− d3 r d3 v T ln f . (4.16)
∂t c
† T can be brought inside the time derivative in the second term because hδf 2 i/f0 f0 .
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 37
If the right-hand side is ignored, T can only increase if hδf 2 i increases too.
It is useful to work out the collision term in Eq. (4.16) in terms of f0 and δf : using the fact that
hδf i = 0 by definition an d the conservation of the particle number by the collision operator,
we get
ZZ ZZ
∂f ∂f0 T δf ∂δf
d3 r d3 v T ln f ≈ d3 r d3 v T ln f0 +
∂t c ∂t c f0 ∂t c
2
Z ZZ
mv ∂f 0 T δf ∂δf
= V d3 v + d3 r d3 v . (4.17)
2 ∂t c f0 ∂t c
The second term is the collisional damping of δf , of which more will be said soon. The first
term is the collisional energy exchange between the equilibrium distributions of different species
(intra-species collisions conserve energy, but there can be friction between
P species). If the species
under consideration is α, this energy exchange can be represented as α0 ναα0 (Tα − Tα0 ) (e.g.,
Helander & Sigmar 2005) and will act to equilibrate temperatures between species as the system
strives toward thermal equilibrium. If the collision frequencies ναα0 are small, this will be a slow
effect. Due to overall energy conservation, the energy-exchange terms vanish exactly if Eq. (4.17)
is summed over species.
Finally, let us sum Eq. (4.16) over species and use Eq. (4.9) to relate the total heating
to the rate of change of the electric-perturbation energy:
" #
2 2
3 Tα hδfα i 3 hE i
ZZ Z X ZZ
d X 3 Tα δfα ∂δfα
d rd v + d r = d3 r d3 v 60 ,
dt α 2f0α 8π α
f0α ∂t c
| {z }
≡W
(4.18)
where we used Eq. (4.17) in the right-hand side (with the total equilibrium collisional
energy-exchange terms vanishing upon summation over species). The right-hand side
must be non-positive-definite because collisions cannot decrease entropy [Eq. (4.12)].
Eq. (4.18) is a way to express the idea that, except for the effect of collisions, the
change in the electric-perturbation energy (= −heating) must be compensated by the
change in hδf 2 i, in terms of a conservation law of a quadratic positive-definite quantity,
W , that measures the total amount of perturbation in the system (a quadratic norm of
the perturbed solution).† Examining the quantity under the time derivative in Eq. (4.18),
we realise that it is the free energy of the perturbed state, comprising the entropy of the
perturbed distribution and the energy of the electric field:
hE 2 i
X Z
W =E− Tα hδSα i, E = d3 r . (4.19)
α
8π
Figure 14. Shifting the integration contour in Eq. (4.23). This is analogous to Fig. 4 but note
the additional “kinetic” pole.
ρu2
Z Z
d
d3 r = −µ d3 r |∇u|2 6 0. (4.21)
dt 2
The conserved quadratic quantity is kinetic energy and the negative-definite dissipation (leading
to net entropy production) is due to viscosity.
Thus, as the electric perturbations decay via Landau damping, the perturbed distri-
bution function must grow. This calls for going back to our solution for it (§3.1) and
analysing carefully the behaviour of δf .
To compute the inverse Laplace transform, Eq. (3.6), we adopt the same method as in
§3.1 (Fig. 4), viz., shift the contour to large negative Re p as shown in Fig. 14 and use
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 39
Cauchy’s formula:
Z i∞+σ
1
δf (t) = dp ept δfˆ(p)
2πi −i∞+σ
" #
q X ci epi t ∂f0 q X c i ∂f0
=i k· + e−ik·vt g − i k· . (4.23)
m i pi + ik · v ∂v m i pi + ik · v ∂v
| {z } | {z }
eigenmode-like ballistic response,
solution, comes from comes from
the poles of ϕ(p) p = −ik · v
Note that the analytic part of ϕ̂(p) must be A(−ik ·v) = 0 in order for the above solution
to satisfy δf (t = 0) = g.
The solution (4.23) teaches us two important things.
1) First, the Landau-damped solution is not an eigenmode. Even though the evolution
of the potential, Eq. (3.16), does look like a sum of damped eigenmodes of the form
ϕ ∝ epi t , Re pi < 0, the full solution of the Vlasov–Poisson system does not decay: there
is a part of δf (t), the “ballistic response” ∝ e−ik·vt , that oscillates without decaying—in
fact, we shall see in §4.6 that δf even has a growing part! It is this part that is responsible
for keeping free energy conserved, as per Eq. (4.18) without collisions. Thus, you may
think of Landau damping as a process of transferring (free) energy from the electric-field
perturbations to the perturbations of the distribution function.
In contrast to the case of damping, a growing solution (Re pi > 0) can be viewed as an eigenmode
because, after a few growth times, the first term in Eq. (4.23) will be exponentially larger than
the ballistic term. This will allow us to ignore the latter in our treatment of QLT (§7.1)—a
handy, although not necessary (see Q8) simplification. Note that reversibility is not an issue for
the growing solutions: so, there may be (and often are) damped solutions as well, so what? We
only care about the growing modes because they will be all that is left if we wait long enough.
The velocity integral over the fine structure increasingly cancels as time goes on—a
perturbation initially “visible” as ϕ phase-mixes away, disappearing into the negative
entropy associated with the fine velocity dependence of δf [Eq. (4.15)].
More generally speaking, one can similarly argue that the refinement of velocity de-
pendence of δf causes lower velocity moments of δf (density, flow velocity, temperature,
heat flux, and so on) to decrease with time, transferring free energy to higher moments
(ever higher as time goes on). One way to formalise this statement neatly is in terms
of Hermite moments: since Hermite polynomials are orthogonal, the free energy of the
perturbed distribution can be written as a sum of “energies” of the Hermite moments
[see Eq. (10.32)]. It is then possible to represent the Landau-damped perturbations as
having a broad spectrum in Hermite space, with the majority of the free energy resid-
ing in high-order moments—infinitely high in the formal limit of zero collisionality and
infinite time (see Q7 and Kanekar et al. 2015).
Since the mth-order Hermite moment can, for m √ 1, be asymptotically represented as a cos
function in v space oscillating with the “frequency” 2m/vth [Eq. (10.33)], Eq. (4.24) implies
that the typical moment in which the free energy resides grows with time as m ∼ (kvth t)2 .
Note that tc ν −1 provided ν kvth , i.e., tc is within the range of times over which
our “collisionless” theory is valid. After time tc , “collisionless” damping becomes irre-
versible because the part of δf that is fast-varying in velocity space is lost (entropy has
grown) and so it is no longer possible, even in principle, to invert all particle trajectories,
have the system retrace back its steps, “phase-unmix” and thus “undamp” the damped
perturbation.
Some rather purist theoreticians sometimes choose to replace collisional estimates of the type
discussed above by a stipulation that δf (v) must be “coarse-grained” beyond some suitably
chosen scale in v (Fig. 16)—this is equivalent to saying that the formation of the fine-structured
phase-space part of δf constitutes a loss of information and so leads to growth of entropy (i.e., the
loss of negative entropy associated with hδf 2 i). Somewhat non-rigorously, this means that we can
just consider the ballistic term in Eq. (4.23) to have been wiped out and use the coarse-grained
(i.e., velocity-space-averaged) version of δf :†
q X ci epi t ∂f0
δf = i k· . (4.30)
m i pi + ik · v ∂v
We can check that the correct solution for the potential, Eq. (3.16), can be recovered from this:
Z
4π X
ϕ= 2 qα d3 v δf α
k α
X p t X 4πqα 2 Z
X
3 1 ∂f0α
= ci e i i d v k · − 1 +1 = ci epi t . (4.31)
i α
mα k2 pi + ik · v ∂v i
| {z }
= −(pi , k) = 0 by definition of pi
† With an understanding that any integral involving the resonant denominator must be taken
along the Landau contour (see Q8).
42 A. A. Schekochihin
If you are wondering how this works without the coarse-graining kludge, read on.
The second term is the ballistic evolution of the initial perturbation (particles flying
apart in straight lines at their initial velocities)—the homogeneous solution of the kinetic
equation (3.1). This develops a lot of fine-scale velocity-space structure, but obviously
does not grow. The first term, the particular solution arising from the (linear) wave-
particle interaction, is more interesting, especially around the resonances Re pi +k·v = 0.
Consider one of the modes, pi = −iω + γ, and assume γ k · v ∼ ω. This allows us
to introduce “intermediate” times:
1 1
t . (4.33)
k·v γ
This means that the wave has had time to oscillate, phase mixing has got underway,
but the perturbation has not yet been damped away completely. We have then, for the
relevant piece of the perturbed distribution (4.32),
epi t − e−ik·vt eγt − e−i(k·v−ω)t 1 − e−i(k·v−ω)t
δf ∝ = −ie−iωt ≈ −ie−iωt , (4.34)
pi + ik · v k · v − ω − iγ k·v−ω
with the last, approximate, expression valid at the intermediate times (4.33), assuming
also that, even though we might be close to the resonance, we shall not come closer than
γ, viz., |k · v − ω| γ. Respecting this ordering, but taking |k · v − ω| 1/t, we find
δf ∝ t e−iωt . (4.35)
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 43
Thus, δf has a peak that grows with time, emerging from the sea of fine-scale but constant-
amplitude structures (Fig. 17). The width of this peak is obviously |k · v − ω| ∼ 1/t and
so δf around the resonance develops a sharp structure, which, in the formal limit t → ∞
(but respecting γt 1, i.e., with infinitesimal damping), tends to a delta function:
1 − e−i(k·v−ω)t
δf ∝ −ie−iωt → e−iωt πδ(k · v − ω) as t → ∞. (4.36)
k·v−ω
T |δfk |2 d |E k |2 |E k |2
Z
d
d3 v =− = −2γk . (4.38)
dt 2f0 dt 8π 8π
Ignoring the term in Eq. (4.32) involving g as it obviously cannot give us a growing amplitude,
letting the relevant root of the dispersion relation be pi = −iωpe + γk , where γk is given
by Eq. (3.41), and assuming a Maxwellian f0 , we may write the solution (4.32) for electrons
(q = −e) as
† Here we come into passing contact with an alternative (to Landau’s) formalism, due to
van Kampen, for treating collisionless plasmas. The objective is more mathematical rigour—but
even if this is of limited appeal to you, the book by van Kampen & Felderhof (1967) is still a
good read and a good chance to question and reexamine your understanding of how it all works.
44 A. A. Schekochihin
We are going to have to compute |δfk |2 and squaring delta functions is a dangerous game
belonging to the class of games that one must play veeery carefully. Here is how:
2
1 − e−ixt 1 − eixt 1 − e−ixt 1 − eixt
= ≈ iπδ(x) = tπδ(x). (4.40)
x x x x }
| {z
= −it at
x=0
Using this prescription,
T |δfk |2 e2 |ϕk |2
Z Z
d3 v = d3 v 2
(k · v)2 tπδ(k · v − ωpe )f0
2f0 me vthe
4
k2 |ϕk |2 ωpe π ω
pe |E k |2
= 2t 2
F = −2γ k t . (4.41)
8π k3 ne vthe k 8π
| {z }
= −γk , Eq. (3.41)
Thus, the entropic part of the free energy grows secularly with time and its time derivative
satisfies Eq. (4.38), q.e.d.
the question of whether the dispersion relation (3.17) has any unstable solutions: roots
with growth rates γi (k) > 0.
It is going to be useful to write the dielectric function (3.26) as follows
2 Z
ωpe F̄ 0 (vz )
(p, k) = 1 − 2 dvz , (5.1)
k CL vz − ip/k
1 X 2 me Fe Zme Fi
F̄ = Zα Fα = + , (5.2)
ne α mα ne mi ni
where the last expression in (5.2) is for the case of a two-species plasma. Thus, the dis-
tribution functions of different species come into the linear problem additively, weighted
by their species’ charges and (inverse) masses.
Let us develop a method (due to Nyquist 1932) for counting zeros of (p) (we will
henceforth suppress k in the argument) in the half-plane Re p > 0 (the unstable roots
of the dispersion relation). We observe that (p) is analytical (be virtue of our efforts in
§3.2 to make it so) and that if p = pi is its zero of order Ni , then in its vicinity,
∂ ln (p) Ni
(p) = const (p − pi )Ni + . . . ⇒ = + ..., (5.3)
∂p p − pi
so zeros of are poles of ∂ ln (p)/∂p; the latter function has no other poles because is
analytical. If we now integrate this function over a closed contour CR running along the
imaginary axis (and just to the right of it: p = −iω + 0) in the complex p plane from iR
to −iR and then along a semicircle of radius R back to iR (Fig. 18), we will, in the limit
R → ∞, capture all these poles:
Z
∂ ln (p) X
lim dp = 2πi Ni = 2πiN, (5.4)
R→∞ C
R
∂p i
where N is the total number of zeros of (p) in the half-plane Re p > 0. It turns out
(as we shall prove in a moment) that the contribution to the integral over CR from the
semicircle vanishes at R → ∞ and so we need only integrate along the imaginary axis:
Z −i∞+0
∂ ln (p) (−i∞)
2πiN = dp = ln . (5.5)
+i∞+0 ∂p (+i∞)
46 A. A. Schekochihin
Proof. All we need to show is that
∂ ln (p)
|p| →0 as |p| → ∞, Re p > 0. (5.6)
∂p
Indeed, using Eq. (5.1) and the Landau integration rule (3.20), we have in this limit:
2 Z +∞
ωpe ik ikvz 1 X 2
(p) = 1 − 2 dvz F̄ 0 (vz ) 1− + ... ≈ 1 + 2 ωpα , (5.7)
k −∞ p p p α
R
where we have integrated by parts and used dvz Fα = nα . Manifestly, the condition (5.6) is
satisfied.
Note that, along the imaginary axis p = −iω, by the same expansion and using also the
Plemelj formula (3.23), we have
2
1 X 2 ωpe 0 ω
(−iω) ≈ 1 − ωpα − iπ F̄ → 1 ∓ i0 as ω → ∓∞. (5.8)
ω2 α k2 k
In view of Eq. (5.8) and of our newly proven formula (5.5), as the function (−iω) runs
along the real line in ω, it changes from
Here P (v∗ ) is (minus) the principal-value integral in Eq. (5.11), which can be manipulated
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 47
as follows:
+∞ Z +∞
F̄ 0 (vz )
Z
1 ∂
P (v∗ ) = −P dvz = −P dvz F̄ (vz ) − F̄ (v∗ )
−∞ vz − v∗ −∞ vz − v∗ ∂vz
+∞
F̄ (v∗ ) − F̄ (vz )
Z
= dvz , (5.13)
−∞ (vz − v∗ )2
where we have integrated by parts; the additional term F̄ (v∗ ) was inserted under the
derivative in order to eliminate the boundary terms arising in this integration by parts
around the pole vz = v∗ .†
Now we are ready to analyse particular (and, as we shall see, also generic) equilibrium
distributions F̄ (vz ).
5.1.1. Stability of Monotonically Decreasing Distributions
Consider first a distribution function that has a single maximum at vz = v0 and
monotonically decays in both directions away from it (Fig. 20a): F̄ 0 (v0 ) = 0, F̄ 00 (v0 ) < 0.
This means that, besides at ω = ∓∞, Im (−iω) ∝ F̄ 0 (ω/k) also vanishes at ω = kv0 . It
is then clear that
2
ωpe
(−ikv0 ) = 1 + 2 P (v0 ) > 1 (5.14)
k
because F̄ (v0 ) > F̄ (v) for all vz and so P (v0 ) > 0. Thus, the Nyquist curve departs from
1 − i0 at ω = −∞, intersects the real line once at ω = kv0 and then comes back to
1 + i0 without circling zero; the corresponding Nyquist digram is sketched in Fig. 19(a).
Conclusion:
Monotonically decreasing distributions are stable against electrostatic perturbations.
We do not in fact need all this mathematical machinery just to prove the stability of monotoni-
cally decreasing distributions (in §5.2.1, we will see that this is a very robust result)—but it will
come handy when dealing with less simple cases. Parenthetically, let me give two alternative
proofs of stability, the first of which (Exercise 5.1) is useful technically and the second (§5.1.2)
conceptually.
† Note that in the final expression in Eq. (5.13), there is no longer a need for principal-value
integration because, v∗ being a point of extremum of F̄ , the numerator of the integrand is
quadratic in vz − v∗ in the vicinity of v∗ .
48 A. A. Schekochihin
(a) (b)
(c) (d)
Figure 21. Various possible forms of the Nyquist diagram for a single-minimum distribution
sketched in Fig. 20b: (a) (−ikv0 ) > 1, stable; (b) (−ikv0 ) < 0, (−ikv2 ) > 1, unstable;
(c) (−ikv0 ) < (−ikv2 ) < 0, stable; (d) (−ikv0 ) < 0 < (−ikv2 ) < 1, unstable.
of F̄ and so P (v∗ ) > P (v0 ) > 0 for all other extrema. In this situation, illustrated in
Fig. 21(a), the Nyquist curve never circumnavigates zero and, therefore, P (v0 ) > 0 is
a sufficient condition of stability. It is also the necessary one, which is proved in the
following way.
Suppose P (v0 ) < 0. Then, in Eq. (5.12), we can always find a range of k that are
small enough that (−ikv0 ) < 0, so the downward crossing at v0 happens on the negative
side of zero in the plane. After this downward crossing, the Nyquist curve will make
more crossings, until it finally comes to rest at 1 + i0 as ω = +∞. Let us denote by
v2 > v0 the point of extremum for which the corresponding crossing occurs at a point
on the Re axis that is closest to (but always will be to the right of) (−ikv0 ) < 0. If
(−ikv2 ) > 0, then there is no way back, zero has been fully circumnavigated and so
there must be at least one unstable root (see Fig. 21b,d). If (−ikv2 ) < 0, there is in
principle some wiggle room for the Nyquist curve to avoid circling zero (see Fig. 21c for a
single-minimum distribution of Fig. 20b—or Fig. 19b for some serious wiggles). However,
since P (v2 ) > P (v0 ) for any v2 (because v0 is the global minimum of F̄ ), we can always
increase k in Eq. (5.12) just enough so (−ikv2 ) > 0 even though (−ikv0 ) < 0 still (this
corresponds to turning Fig. 21c into Fig. 21d). Thus, if P (v0 ) < 0, there will always be
some range of k inside which there is an instability.
We have obtained a sufficient and necessary condition of instability of an equilibrium
50 A. A. Schekochihin
F̄ (vz ) against electrostatic perturbations: if v0 is the point of global minimum of F̄ ,†
+∞
F̄ (v0 ) − F̄ (vz )
Z
P (v0 ) = dvz <0 ⇔ F̄ is unstable . (5.15)
−∞ (vz − v0 )2
This is the famous Penrose’s instability criterion (famous criterion, not the famous Pen-
rose; it was proved by Oliver Penrose 1960, in a stylistically somewhat different way
than I did it here). Note that considerations of the kind presented above can be used
to work out the wave-number intervals, corresponding to various troughs in F̄ , in which
instabilities exist.
Intuitively, the criterion (5.15) says that, in order for a distribution to be unstable, it
needs to have a trough and this trough must be deep enough. Thus, if F (v0 ) = 0, i.e.,
if the distribution has a “hole”, it is always unstable (an extreme example of this is the
two-stream instability; see Exercise 3.5). Another corollary is that you cannot stablise a
distribution by just adding some particles in a narrow interval around v0 , as this would
create two minima nearby, which, the filled interval being narrow, are still going to
render the system unstable. To change that, you must fill the trough substantially with
particles—hence the tendency to flatten bumps into plateaux, which we will discover
in §7 (this answers, albeit in very broad strokes, the question posed at the beginning of
§5 about the types of stable distributions towards which the unstable ones will be pushed
as the instabilities saturate).
Exercise 5.2. Consider a single-minimum distribution like one in Fig. 20(b), but with the
global maximum on the right and the lesser maximum on the left of the minimum. Draw various
possible Nyquist diagrams and convince yourself that Penrose’s criterion works. If you enjoy
this, think of a distribution that would give rise to the Nyquist diagram in Fig. 19(b).
Exercise 5.3. What happens if the distribution function F̄ has an inflection point, i.e., F̄ (v0 ) =
0, F̄ 0 (v0 ) = 0, F̄ 00 (v0 ) = 0?
Exercise 5.4. What happens if the distribution function has a trough with a flat bottom (i.e.,
a flat minimum over some interval of velocities)?
† Another way of putting this is: a distribution F̄ is unstable iff it has a minimum at some
v0 for which P (v0 ) < 0. Obviously, if P (v0 ) < 0 at some minimum, it is also so at the global
minimum.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 51
Figure 22. Combined distribution (5.2) for cold ions and hot electrons (cf. Fig. 12).
Exercise 5.5. Use Nyquist’s method to work out the range of wave numbers at which perturba-
tions will grow for the two-stream instability (you will find the answer in Jackson 1960—yes, that
Jackson). Convince yourself that this is all in accord with the explicit solution of the dispersion
relation that you might have already obtained in Exercise 3.5.
Exercise 5.6. Construct an equilibrium distribution to model your favorite plasma system
with flows and/or beams and investigate its stability: find the growth rate as a function of wave
number, instability conditions, etc.
5.1.5. Anisotropies
So we have found that various holes, bumps, streams, beams, flows, currents and
other such nonmonotonic features in the (combined, multispecies) equilibrium distribu-
tion present an instability risk, unless they are sufficiently small, shallow, wide and/or
close enough to the thermal bulk. All of these are, of course, anisotropic features—indeed,
as we saw in Exercise 5.1(c,d), 3D-isotropic distributions are harmless, instability-wise.
It turns out that anisotropies of the distribution function in velocity space are dangerous
even when the distribution decays monotonically in all directions.† However, the insta-
bilities that occur in such situations are electromagnetic, rather than electrostatic, and
so require an investigation into the properties of the transverse dielectric function of the
kind derived in Q2 or Q3, but for a general equilibrium. The corresponding instability
criterion is derived in Q4, by a somewhat adjusted version of Nyquist’s method. A nice
treatment of anisotropy-driven instabilities can be found in Krall & Trivelpiece (1973)
and an even more thorough one in Davidson (1983). In §5.2.4, we will show in quite a
simple way that, at least in principle, there is energy to be extracted from anisotropic
distributions.
where f∗ is some trial distribution, which will represent our best guess about the prop-
erties of the stable distribution towards which the system will want to evolve and/or in
the general vicinity of which we are interested in investigating stability. The function
A(r, v, f ) is chosen in such way that for any f ,
A[f, f∗ ] > 0. (5.20)
If it is also chosen so that H is conserved by the (collisionless) Vlasov–Maxwell equations,
then H(t) = H(0) and the inequality (5.20) gives us the following bound on the field
† In Q3, you have an opportunity to derive the most famous of the instabilities triggered by
anisotropy.
‡ These ideas appear to have crystallised in the papers by T. K. Fowler in early 1960s (see
his review, Fowler 1968; his reminiscences and speculations on the subject 50 years later can be
found in Fowler 2016), although a number of founding fathers of plasma physics were thinking
along these lines around the same time (references are given in opportune places below).
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 53
energy at time t:
It is obvious that we cannot extract from a distribution more energy than K(0), but
the above tells us that, in fact, one can extract less. Thus, A[f0 , f∗ ] can be thought of
as an upper bound on the available energy of the distribution f0 . The sharper it can be
made, the closer we are to learning something useful. Thus, the idea is to identify some
suitable functional A[f, f∗ ] for which H is conserved and some class of trial distributions
f∗ for which (5.20) holds, then minimise A[f0 , f∗ ] within that class, subject to whatever
physical constraints one can reasonably expect to hold: e.g., conservation of particles,
momentum, any other (possibly approximate) invariants that the system might possess
(e.g., its adiabatic invariants; see Helander 2017).†
To make some steps towards practical implementation of this programme, let us in-
vestigate how to choose A in such way as to ensure conservation of H:
mα v 2 ∂fα
ZZ X ZZ
dH dE X ∂A ∂fα ∂A
= + d3 r d3 v = d3 r d3 v − = 0.
dt dt α
∂fα ∂t α
∂fα 2 ∂t
(5.23)
The second equality was obtained by using the conservation of total energy,
d X ZZ mα v 2
(E + K) = 0, K = d3 r d3 v fα , (5.24)
dt α
2
where K is the kinetic energy of the particles. Eq. (5.23) tells us how to choose A:
X mα v 2
A(r, v, f ) = fα + Gα (fα ) , (5.25)
α
2
where Gα (fα ) are arbitrary functions of fα . These can be added here because Vlasov’s
equation has an infinite number of invariants: for any (sufficiently smooth) Gα (fα ),
ZZ
d
d3 r d3 v Gα (fα ) = 0. (5.26)
dt
This follows from the fact that, in the absence of collisions, the kinetic equation (1.30)
expresses the conservation of phase volume in (r, v) space (the flow in this phase space
is divergence-free).
† Krall & Trivelpiece (1973) comment with a slight air of resignation that, with the rules
of the game much vaguer than in linear theory, the thermodynamical approach to stability is
“more art than science”. In the Russian translation of their textbook, this statement provokes
a stern footnote from the scientific editor (A. M. Dykhne), observing that the right way to put
it would be “more art than craft”.
54 A. A. Schekochihin
Exercise 5.7. Prove the conservation law (5.26), assuming that the system is isolated.
The existence of an infinite number of conservation laws suggests that the evolution of a colli-
sionless system in phase space is much more constrained than that of a collisional one. In the
latter case, the evolution is constrained only by conservation of particles, momentum and energy
and the requirement that entropy must not decrease. We shall return shortly to the question of
how available energy might be related to entropy.
A quick sanity check is to choose Gα (fα ) = 0. The inequality (5.20) is then certainly
satisfied for f∗ ∝ δ(v) and the bound (5.21) becomes
E(t) − E(0) 6 K(0), (5.27)
i.e., one cannot extract any more energy than the total energy contained in the distribution—
indeed, one cannot. Let us now move on to more nontrivial results.
∂f0α
<0 ⇒ stability . (5.28)
∂εα
Proof. For every species (suppressing species indices), let us again take G(f ) = 0 in
Eq. (5.25), but construct a nontrivial f∗ that satisfies (5.20) with f (t) at any time t since
the beginning of its evolution from the initial distribution f0 .
Let us define f∗ to be a monotonically decreasing function of v 2 (i.e., energy) such that
for any Λ > 0, the volume of the region in the phase space (r, v) where f∗ > Λ is the
same as the volume of the phase-space region where f0 > Λ. Then f∗ is the distribution
with the smallest kinetic energy [Eq. (5.24)], denoted here by K∗ , that can be reached
† The stability of Maxwellian equlibria against small perturbations was first proved by
W. Newcomb, whose argument was published as Appendix I of Bernstein (1958) (and followed
by Fowler 1963, who proved stability against large perturbations). Gardner (1963) attributes the
first appearance of the stability condition (5.28) to an obscure 1960 report by M. N. Rosenbluth,
although the same condition was derived also by Kruskal & Oberman (1958), more or less in the
manner described in §5.2.2. Many great minds were clearly thinking alike in those glory days of
plasma physics.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 55
from f0 while preserving phase-space volume:
K(t) > K∗ . (5.29)
Indeed, while the phase-space volume occupied by any given value of the probability
density is the same for f0 and for f∗ , the corresponding energy is always lower for f∗
than for f0 or for any other f that can evolve from it, because in f∗ , the values of the
probability density are rearranged in such a way as to put the largest of them at the
lowest values of v 2 , thus minimising the velocity integral in Eq. (5.24).
A vivid analogy is to think of the evolution of f under the collisionless kinetic equation
(1.28) as the evolution of a mixture of “fluids” of different densities (values of f ) advected
in a 6D phase-space (r, v) by a divergence-free flow (ṙ, v̇). The lowest-energy state is the
one in which these fluids are arranged in layers of density decreasing with increasing v 2 ,
the heaviest at the bottom, the lightest at the top (Fig. 23).
In view of (5.29) and since A is given by Eq. (5.25) with G(f ) = 0,
A[f, f∗ ] = K(t) − K∗ > 0, (5.30)
so (5.20) holds and the bound (5.21) follows. When f0 = f∗ , i.e., the equilibrium distri-
bution satisfies (5.28), the system is stable, q.e.d.
Note that the condition (5.28) is sufficient, but not necessary, as we already know from,
e.g., Exercise 5.1(c,d).
In a fresh-off-the-press paper, Helander (2017) has developed a beautiful scheme for calculating
“ground states” (the states of minimum energy) of Vlasov’s equation, i.e., for determining f∗
and then calculating K∗ to work out specific values of the available energy implied by the upper
bound (5.21), A[f0 , f∗ ] = K(0) − K∗ .
The condition (5.28) makes H positive definite and so no wonder the system is sta-
ble: perturbations around f0 have a conserved norm! For a Maxwellian equilibrium,
−∂f0α /∂εα = f0α /Tα , so this H is none other than W , (the electromagnetic version of)
our free energy (4.19), and so it is tempting to think of Eq. (5.35) as providing a natural
generalisation of free energy to non-Maxwellian plasmas.
In Q5, the results of this section are obtained in a more straightforward way, directly from the
Vlasov–Maxwell equations.
This style of thinking has been having a revival lately: see, e.g., the discussion of firehose and
mirror stability of a magnetised plasma in Kunz et al. (2015). Generalised energy invariants like
H are important not just for stability calculations, but also for theories of turbulence in near-
collisionless environments that are not particularly close to thermal (Maxwellian) equilibrium,
e.g., the solar wind (see, e.g., Schekochihin et al. 2009).
Exercise 5.8. Prove that if Gα and f∗α are given by Eq. (5.36), then
X Z Z 3 3 mα v 2
A[f, f∗ ] = d rd v (fα − f∗α ) + Gα (fα ) − Gα (f∗α ) > 0 (5.37)
α
2
Thus, Eq. (5.21) provides an upper bound on the energy of any electromagnetic fields
that can be extracted from the initial distribution f0α . In order to make this bound as
sharp as possible, one picks the constants Cα and Tα so as to minimise A[f0 , f∗ ] subject
to constraints that cannot change: e.g., the number of particles of each species:
3/2 Z Z
mα 1
Cα = d3 r d3 v f0α . (5.38)
2πTα V
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 57
5.2.4. Anisotropic Equilibria
Let me give an example of the use of this scheme for deriving an upper bound on the
energy of unstable perturabtions for a case of an anisotropic initial distribution—the case
that, at the end of §5.1, I had to relegate to Q4 as it needed substantial extra work if it
were to be handled by the method developed there.
Consider a bi-Maxwellian distribution, a useful and certainly the simplest model for
anisotropic equilibria:
!
m 3/2
α 1 2
m α v⊥ mα vk2
f0α = nα 1/2
exp − − , (5.39)
2π T⊥α T 2T⊥α 2Tkα
kα
where T⊥α and Tkα are the “temperatures” of particle motion perpendicular and parallel
to some special direction. Is this distribution unstable? (Yes: see Q3.) To obtain an
upper bound on the energy available for extraction from it, we substitute the distribution
function (5.39) into Eq. (5.36), use also Eq. (5.38), and find
3/2
X Tα Tkα 3
A[f0 , f∗ ] = V nα ln 1/2
+ T⊥α + − Tα . (5.40)
α T⊥α T 2 2
kα
2/3 1/3
This is minimised by Tα = T⊥α Tkα , resulting in the following estimate of the available
energy:†
3 X 2 1 2/3 1/3
E(t) − E(0) 6 min A[f0 , f∗ ] = V nα T⊥α + Tkα − T⊥α Tkα . (5.41)
Tα 2 α
3 3
The bound is zero when T⊥α = Tkα and is always positive otherwise because it is the
difference between an arithmetic and a geometric mean of the two temperatures. We do
not, of course, have any way of knowing how good an approximation this is to the true
saturated level of whatever instability might exist here in any particular physical regime
(if any),‡ but this does hint rather suggestively that temperature anisotropy is a source
of free energy.
Further examples of such calculations can be found in Krall & Trivelpiece (1973, §9.14) and
Fowler (1968). A certain further development of the methodology discussed above allows one
to derive upper bounds not just on the energy of perturbations but also on their growth rates
(Fowler 1964, 1968).
7. Quasilinear Theory
7.1. General Scheme of QLT
In §§3 and 4, we discussed at length the structure of the linear solution corresponding to
a Landau-damped initial perturbation. This could be adequately done for a Maxwellian
plasma and we have found that, after some interesting transient time-dependent phase-
space dynamics, perturbations damp away and their energy turns into heat, increasing
somewhat the temperature of the equilibrium (see, however, Q8).
We now consider a different problem: an unstable (and so decidedly non-Maxwellian)
equilibrium distribution giving rise to exponentially growing perturbations. The specific
example on which we shall focus is the bump-on-tail instability, which involves generation
of unstable Langmuir waves with phase velocities corresponding to instances of positive
derivative of the equilibrium distribution function (Fig. 24). The energy of the waves
grows exponentially:
3
∂|E k |2 π ωpe 1 0 ωpe
= 2γk |E k |2 , γk = F , (7.1)
∂t 2 k 2 ne k
R R
where F (vz ) = dvx dvy f0 (v) [see Eq. (3.41)]. In the absence of collisions, the only
way for the system to achieve a nontrivial steady state (i.e., such that |E k |2 is not just
zero everywhere) is by adjusting the equilibrium distribution so that
ω
pe
γk = 0 ⇔ F 0 =0 (7.2)
k
at all k where |E k |2 6= 0, say, k ∈ [k2 , k1 ]. If we translate this range into velocities,
v = ωpe /k, we see that the equilibrium must develop a flat spot:
ωpe ωpe
F 0 (v) = 0 for v ∈ [v1 , v2 ] = , . (7.3)
k1 k2
This is called a quasilinear plateau (§7.4). Obviously, the rest of the equilibrium distri-
bution may (and will) also be modified in some, to be determined, way (§§7.6, 7.7).
These modifications of the original (initial) equilibrium distributions can be accom-
plished by the growing fluctuations via the feedback mechanism already discussed in §2.3,
¶ The third kind is asking for general criteria of certain kinds of behaviour, such as stability
or otherwise—we dabbled in this type of nonlinear theory in §5.2.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 59
namely, the equilibrium distribution will evolve slowly according to Eq. (2.11):
∂f0 q X ∗ ∂δfk
=− ϕk ik · . (7.4)
∂t m ∂v
k
−1
The time averaging here [see Eq. (2.7)] is over ωpe ∆t γk−1 .
The general scheme of QLT is:
We keep only the fastest growing mode (all others are exponentially small after a
while), and so the solution (3.16) for the electric perturbations is
In the solution (4.23) for the perturbed distribution function, we may ignore the ballistic
term because the exponentially growing piece (the first term) will eventually leave all
this velocity-space structure behind,† so
† See, however, Q8, where we work out how to avoid having to wait for this to happen: in
fact, the results below are valid for γk t . 1 as well.
60 A. A. Schekochihin
This is a diffusion equation in velocity space, with a velocity-dependent diffusion matrix
q2 X 1
D(v) = − ikk|ϕk |2
m2 k · v − ωk − iγk
k
q2 X
21 1 1
=− 2 ikk|ϕk | +
m 2 k · v − ωk − iγk −k · v − ω−k − iγ−k
k | {z }
here we changed
variables k → −k
q 2 X kk
2i 1 1
=− 2 |E k | −
m k2 2 k · v − ωk − iγk k · v − ωk + iγk
k
q 2 X kk 1
= 2 2
|E k |2 Im
m k k · v − ωk − iγk
k
q 2 X kk γk
= 2 |E k |2 . (7.8)
m k 2 (k · v − ωk )2 + γk2
k
To obtain these expressions, we used the fact that the wave-number sum could just as
well be over −k instead of k and that ω−k = −ωk , γ−k = γk (because ϕ−k = ϕ∗k , where
ϕk is given by Eq. (7.5)). The matrix D is manifestly positive definite—this adds credence
to our a priori expectation that a plateau will form: diffusion will smooth the bump in
the equilibrium distribution function.
The question of validity of the QL approximation is quite nontrivial and rife with subtle issues, all
of which I have swept under the carpet. They mostly have to do with whether coupling between
waves [the last term in Eq. (2.12)] will truly remain unimportant throughout the quasilinear
evolution, especially as the plateau regime is approached and the growth rate of the waves
becomes infinitesimally small. If you wish to investigate further—and in the process gain a finer
appreciation of nonlinear plasma theory,—the article by Besse et al. (2011) (as far as I know,
the most recent substantial contribution to the topic) is a good starting point, from which you
can follow the paper trail backwards in time and decide for yourself whether you trust the QLT.
When we get to the stage of solving a specific problem (§7.3), we shall see that paying attention
to energy and momentum budgets leads one to important discoveries about the QL evolution of
the particle distribution. With this prospect in mind, as well as by way of a consistency check,
let us show that Eqs. (7.7–7.8) conserve energy and momentum.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 61
7.2.1. Energy Conservation
The rate of change of the particle energy associated with the equilibrium distribution is
mα v 2 X Z Z 3 3 mα v 2 ∂
ZZ
dK d X ∂f0α
≡ d3 r d3 v f0α = d rd v · Dα (v) ·
dt dt α 2 α
2 ∂v ∂v
X ZZ 3 3 ∂f0α
=− d r d v mα v · Dα (v) ·
α
∂v
2 X
|E k |2
Z
X qα k·v ∂f0α
= −V 2
d3 v Im k·
α
m α k k · v − ω k − iγk ∂v
k | {z }
add and substract
ωk + iγk in the
numerator
X |E k |2 X ωpα 2 Z
1 3 1 ∂f0α
= −V Im (ωk + iγk ) d v k ·
4π α
k2 nα k · v − ωk − iγk ∂v
k
| {z }
= 1 − (−iωk + γk , k) = 1
because −iωk + iγk is a solution of
dispersion relation = 0
2 Z 2
X |E k | d E
= −V 2γk =− d3 r , q.e.d., (7.9)
8π dt 8π
k
viz., the total energy K + d3 r E 2 /8π = const. This will motivate §7.6.
R
These two papers, published in the same year, are a spectacular example of the “great minds
think alike” principle. They both appeared in the Proceedings of the 1961 IAEA conference
in Salzburg, one of those early international gatherings in which the Soviets (grudgingly al-
lowed out) and the Westerners were telling each other about their achievements in the recently
declassified controlled-nuclear-fusion research. The entire Proceedings are now online (www-
naweb.iaea.org/napc/physics/FEC/1961.pdf)—they are a remarkable historical document and
a great read, containing, besides the papers (in three languages), a record of discussions that
were held. The Vedenov et al. (1962) paper is in Russian, but you will find a very similar expo-
sition in English in the review by Vedenov (1963) published shortly thereafter. Two other lucid
62 A. A. Schekochihin
accounts of quasilinear theory belonging to the same period are in the books by Kadomtsev
(1965) and by Sagdeev & Galeev (1969).
∂F ∂ ∂F
= D(v) , (7.11)
∂t ∂v ∂v
where F (v) is the 1D version of the distribution function, v = vz and the diffusion
coefficient, now a scalar, is given by
e2 X 1
D(v) = |Ek |2 Im . (7.12)
m2e kv − ωpe − iγk
k
As we explained when discussing Eq. (7.1), if the fluctuation field has reached a steady
state, it must be the case that
∂|Ek |2
= 2γk |Ek |2 = 0 ⇔ |Ek |2 = 0 or γk = 0, (7.13)
∂t
i.e., either there are no fluctuations or there is no growth (or damping) rate. The re-
sult is a non-zero spectrum of fluctuations in the interval k ∈ [k2 , k1 ] and a plateau in
the distribution function, Eq. (7.3), in the corresponding velocity interval v ∈ [v1 , v2 ] =
[ωpe /k1 , ωpe /k2 ] (Fig. 25). The particles in this interval are resonant with Langmuir
waves; those in the (“thermal”) bulk of the distribution outside this interval are non-
resonant. We will have solved the problem completely if we find
we get
∂ e2 4π 2 ωpe 2 k 2
∂F ∂ ωpe ω
pe
= W ne γk = 2γ ω /v W . (7.18)
∂t ∂v m2e v v 3
π ωpe k=ω
∂v me v 3 | pe
v
pe /v {z }
= ∂W/∂t
Rearranging, we arrive at
ω
∂ ∂ ωpe pe
F− W = 0. (7.19)
∂t ∂v me v 3 v
Thus, during QL evolution, the expression in the square brackets stays constant in time.
Since at t = 0, there are no waves, W = 0, we find
∂ ωpe ω
pe
F (t, v) = F (0, v) + W t, → F plateau as t → ∞. (7.20)
∂v me v 3 v
In the saturated state (t → ∞), W (ωpe /v) = 0 outside the interval v ∈ [v1 , v2 ].
Therefore, Eq. (7.20) gives us two implicit equations for v1 and v2 :
Finally, integrating Eq. (7.20) with respect to v and using the boundary condition
† Why the prefactor is 1/4π, rather than 1/8π, will become clear in §7.6.
‡ In fact, the wave-number integral must be taken along the Landau contour (i.e., keeping
the contour below the pole) regardless of the sign of γk : see Q8, where we work out the QL
theory for Landau-damped, rather than growing, perturbations.
¶ This is somewhat reminiscent of the “Maxwell construction” in thermodynamics of real
gases: the plateau sits at such a level that the integral under it, i.e., the number of particles
involved, stays the same as it was for the same velocities in the initial state; see Fig. 24.
64 A. A. Schekochihin
Thus, only half of the energy lost by the resonant particles goes into the electric-field
energy of the waves,
Kres (0) − Kres (∞)
E(∞) = . (7.26)
2
Since the energy must be conserved overall [Eq. (7.9)], we must account for the missing
half: this is easy to do physically, as, obviously, the electric energy of the waves is their
potential energy, which is half of their total energy—the other half being the kinetic
energy of the oscillatory plasma motions associated with the wave (Exercise 3.1 will help
you make this explicit). These oscillations are enabled by the non-resonant, “thermal-
bulk” particles, and so we must be able to show that, as a result of QL evolution, these
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 65
particles pick up the total of E(∞) of energy—in formal terms, one might say that the
plasma is heated.†
e2 X 2 γk e2 X γk
D(v) = 2
|E k | 2 2 ≈ 2
|Ek |2 2
me (kv − ωpe ) + γk me ωpe
k k
e2 X 1 ∂|Ek |2 4πe2 d X |Ek |2 1 dE
= 2 2 = 2 2 = , (7.27)
me ωpe 2 ∂t me ωpe dt 8π me ne dt
k k
independent of v. The QL evolution equation (7.11) for the bulk distribution is then‡
∂F 1 dE ∂ 2 F
= . (7.28)
∂t me ne dt ∂v 2
Eq. (7.28) describes slow diffusion of the bulk distribution, i.e., as the wave field grows,
the bulk distribution gets a little broader (which is what heating is). Namely, the “ther-
mal” energy satisfies
me v 2 me v 2 ∂ 2 F
Z Z
dKth d 1 dE dE
= dv F = dv 2
= . (7.29)
dt dt 2 me ne dt 2 ∂v dt
| {z }
= me ne
(by parts twice)
Integrating this with respect to time, we find that the missing half of the energy lost by
the resonant particles indeed goes into the heating of the thermal bulk:
Kth (∞) + Kres (∞) + E(∞) = Kth (0) + Kres (0), (7.31)
Eq. (7.28) can be explicitly solved: changing the time variable to τ = E(t)/me ne turns it into a
simple diffusion equation
∂F ∂2F
= . (7.32)
∂τ ∂v 2
† This is slightly loose language. Technically speaking, since there are no collisions, this is
not really heating, i.e., the exact total entropy does not increase. The “thermal” energy that
increases is the energy of plasma oscillations, which are mean “fluid” motions of the plasma,
whereas “true” heating would involve an increase in the energy of particle motions around the
mean. R
‡ Note that this implies d dvF (v)/dt = 0, so the number of these particles is conserved,
there is no exchange between the non-resonant and resonant populations.
66 A. A. Schekochihin
Letting the initial distribution be a Maxwellian and ignoring the bump on its tail, the solution is
0 2
e−(v−v ) /4τ v 02 (v − v 0 )2
Z Z
ne
F (τ, v) = dv 0 F (0, v 0 ) √ = dv 0 p 2 exp − 2 −
4πτ πvthe 4πτ vthe 4τ
2
ne v
= p exp − 2 . (7.33)
2
π(vthe + 4τ ) vthe + 4τ
Since
2 2Te 4E(t) 2 2E(t)
vthe + 4τ = + = Te + , (7.34)
me me ne me ne
one concludes that an initially Maxwellian bulk stays Maxwellian but its temperature grows as
the wave energy grows, reaching in saturation
2E(∞)
Te (∞) = Te (0) + . (7.35)
ne
This is negative, so momentum is indeed lost. Since it cannot go into electric field
[Eq. (7.10)], it must all get transferred to the thermal particles. Let us confirm this.
Going back to the QL diffusion equation for the non-resonant particles, Eq. (7.28),
at first glance, we have a problem: the diffusion coefficient is independent of v and so
momentum is conserved. However, one should never take zero for an answer when dealing
with asymptotic expansions—indeed, it turns out here that we ought to work to higher
order in our calculation of D(v). Keeping next-order terms in Eq. (7.27), we get
e2 X e2 X
2 γk 2 γk 2kv
D(v) = 2 |Ek | = 2 |Ek | 2 1 + + ...
me (kv − ωpe )2 + γk2 me ωpe ωpe
k k
" #
4πe2 d X |Ek |2 X k|Ek |2 Z
1 d kW (k)
≈ 2 2 +v = E + v dk . (7.37)
me ωpe dt 8π 4πωpe me ne dt ωpe
k k
Thus, there is a wave-induced drag term in the QL diffusion equation (7.11), which
indeed turns out to impart to the thermal particles the small additional momentum
that, according to Eq. (7.36), the resonant particles lose when rearranging themselves to
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 67
Figure 27. The initial distribution and the final outcome of the QL evolution: its bulk hotter
and shifted towards the plateau in the tail.
This means that the thermal bulk of the final distribution is not only slightly broader
(hotter) than that of the initial one (§7.6), but it is also slightly shifted towards the plateau
(Fig. 27).
In a collisionless plasma, this is the steady state. However, as this steady state is
approached, γk → 0, so the QL evolution becomes ever slower and even a very small
collision frequency can become important. Eventually, collisions will erode the plateau
and return the plasma to a global Maxwellian equilibrium—which is the fate of all things.
8. Kinetics of Quasiparticles
Let us reimagine our system of particles and waves as a mixture of two interacting
gases: “true” particles (electrons) and quasiparticles, or plasmons, which will be the
“quantised” version of Langmuir waves. If each of these plasmons has momentum ~k
and energy ~ωk , we can declare
V |E k |2 /4π
Nk = (8.1)
~ωk
to be the mean occupation number of plasmons with wave number k (in a box of volume
V ). The total energy of the plasmons is then
X X |E k |2
~ωk Nk = V , (8.2)
4π
k k
twice the total electric energy in the system (twice because it includes the energy of the
mean oscillatory motion of electrons within a wave; see discussion at the end of §7.5).
68 A. A. Schekochihin
This is indeed in line with our previous calculations [Eq. (7.39)]. Note that the role of ~
here is simply to define a splitting of wave energy into individual plasmons—this can be
done in an arbitrary way, provided ~ is small enough to ensure Nk 1. Since there is
nothing quantum-mechanical about our system, all our results will in the end have to be
independent of ~, so we will use ~ as an arbitrarily small parmeter, in which it will be
convenient to expand, expecting it eventually to cancel out in all physically meaningful
relationships.
where w are the probabilities of absorption and emission and must be equal:
w(p − k, k → p) = w(p → k, p − k) ≡ w(p, k). (8.7)
The first term in the right-hand side of Eq. (8.6) describes the absorbtion of one of
(indistinguishable) Nk plasmons by one of np−k electrons, the second term desribes the
emission by one of np electrons of one of Nk + 1 plasmons. The +1 is, of course, a
small correction to Nk 1 and can be neglected, although sometimes, in analogous but
more complicated calculations, it has to be kept because lowest-order terms cancel. Using
Eqs. (8.7), (8.5) and (8.4), we find
∂Nk X
= w(p, k)δ(εep − εlk − εep−k )(np − np−k )Nk
∂t p
Z m v
3 e ~k
= V d vw , k δ(~(k · v − ωk )) f0 (v) − f0 v − Nk
~ me
Z m v 1
e ~ ∂f0
= V d3 v w ,k δ(k · v − ωk ) k· Nk
~ ~
me ∂v
V me ωpe 0 ωpe
= w ,k F Nk . (8.8)
m ~k k
| e {z }
= 2γk
Note that ~ has disappeared from our equations, after being used as an expansion pa-
rameter.
Since Nk ∝ |E k |2 [Eq. (8.1)], the prefactor in Eq. (8.8) is clearly just the (twice) growth
or damping rate of the waves. Comparing with Eq. (7.1), we read off the expression for
the absorption/emission probability:
m ω πme ω 3
e pe pe
w ,k = . (8.9)
~k V ne k 2
Thus, our calculation of Landau damping could be thought of as a calculation of this
probability. Whether there is damping or an instability is decided by whether it is absorp-
tion or emission of plasmons that occurs more frequently—and that depends on whether,
for any given k, there are more electrons that are slightly slower or slightly faster than
the plasmons with wave number k. Note that getting the correct sign of the damping
rate is automatic in this approach, since the probability w must obviously be positive.
The evolution equation for the occupation number of electrons can be derived in a
similar fashion, if we itemise the processes that lead to an electron ending up in a state
with a given wave number p = me v/~ or moving from this state to one with a different
wave number. The four relevant diagrams are the two in Fig. 28 and the additional two
shown in Fig. 29. The absorption and emission probabilities are the same as before and
70 A. A. Schekochihin
where we have expanded twice in small k (i.e., in ~). This is a diffusion equation in p
(or, equivalently, v = ~p/me ) space. In view of Eq. (8.4), Eq. (8.10) has the same form
as Eq. (7.7), viz.,
∂f0 ∂ ∂f0
= · D(v) · , (8.11)
∂t ∂v ∂v
where the diffusion matrix is
X ~Nk me v e2 X kk
D(v) = kk 2 w , k δ(k · v − ωk ) = 2 |E k |2 πδ(k · v − ωk ). (8.12)
me ~ me k2
k k
The last expression is identical to the resonant form of the QL diffusion matrix (7.8)
[cf. Eqs. (7.16) and (10.45)]. To derive it, we used the definition (8.1) of Nk and the
absorption/emission probability (8.7), already known from linear theory.
Thus, we are able to recover the (resonant part of the) QL theory from our new
electron-plasmon interaction approach. There is more to this approach than a pretty
“field-theoretic” reformulation of already-derived earlier results. The diagram technique
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 71
and the interpretation of the nonlinear state of the plasma as arising from interactions
between particles and quasiparticles can be readily generalised to situations in which the
nonlinear interactions in Eq. (2.12) cannot be neglected and/or more than one type of
waves is present. In this new language, the nonlinear interactions would be manifested
as interactions between plasmons (rather than only between plasmons and electrons)
contributing to the rate of change of Nk . There are many possibilities: four-plasmon in-
teractions, interactions between plasmons and phonons (sound waves), as well as between
the latter and electrons and/or ions, etc. Some of these will be further explored in §8.2
and onwards. A comprehensive monograph on this subject is Tsytovich (1995) (see also
Kingsep 2004).
This is a good place to stop these lectures, although it is not, of course, the end of
plasma kinetics: weak and strong turbulence theory, magnetised plasma waves, “drift”
and “gyro-” kinetics—there are vast expanses of interesting physics and interesting theory
to explore beyond this basic introduction. Some of these topics are covered by Parra
(2017b) and others call for further reading (e.g., on turbulence: Kadomtsev 1965; Sagdeev
& Galeev 1969; Tsytovich 1995; Krommes 2015; on gyrokinetics: Howes et al. 2006; Abel
et al. 2013).
† In Q3, dealing with the Weibel instability, you will have to do essentially the same calcu-
lation, but with a non-Maxwellian equilibrium. To avoid doing the work twice, you could do
that question first and then specialise to a Maxwellian f0α . However, the algebra is a bit hairier
for the non-Maxwellian case, so it may be useful to do the simpler case first, to train your
hand—and also to have a way to cross-check the later calculation.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 73
(b) Hence solve the transverse dispersion relation, TT (p, k) = 0, and show that, in
the high-frequency limit (|ζe | 1), the resulting waves are simply the light waves,
which, at long wave lengths, turn into plasma oscillations. What is the wave length
above which light can “feel” that it is propagating through plasma?—this is called the
plasma (electron) skin depth, de . Are these waves damped?
(c) In the low-frequency limit (|ζe | 1), show that perturbations are aperiodic (have
zero frequency) and damped. Find their damping rate and show that this result is valid for
perturbations with wave lengths longer than the plasma skin depth (kde 1). Explain
physically why these perturbations fail to propagate.
3. Weibel instability. Weibel (1958) realised that transverse plasma perturbations can
go unstable if the equilibrium distribution is anisotropic with respect to some special
direction n̂, namely if f0α = f0α (v⊥ , vk ), where vk = v · n̂, v⊥ = |v ⊥ | and v ⊥ =
v − vk n̂. The anisotropy can be due to some beam or stream of particles injected into the
plasma, it also arises in collisionless shocks or, generically, when plasma is sheared or non-
isotropically compressed by some external force. The simplest model for an anisotropic
distribution of the required type is a bi-Maxwellian:†
!
nα v⊥2 vk2
f0α = 3/2 2 exp − 2 − 2 , (10.5)
π vth⊥α vthkα vth⊥α vthkα
p p
where, formally, vth⊥α = 2T⊥α /mα and vthkα = 2Tkα /mα are the two “thermal
speeds” in a plasma characterised by two effective temperatures T⊥α and Tkα (for each
species).
(a) Using exactly the same method as in Q2, consider electromagnetic perturbations
in a bi-Maxwellian plasma, assuming their wave vectors to be parallel to the direction
of anisotropy, k k n̂. Show that the dielectric tensor again has the form (10.3) and the
longitudial dielectric function is again given by Eq. (3.80), while the transverse dielectric
function is
" #
1 2 2
X
2 T⊥α
TT (p, k) = 1 + 2 k c + ωpα 1 − [1 + ζα Z(ζα )] . (10.6)
p α
Tkα
(b) Show that in one of the tractable asymptotic limits, this dispersion relation has a
zero-frequency, purely growing solution with the growth rate
kvthke Tke
∆e − k 2 d2e ,
γ= √ (10.7)
π T⊥e
where ∆e = T⊥e /Tke − 1 is the fractional temperature anisotropy, which must be positive
in order for the instability to occur. Find the maximum growth rate and the corresponding
wave number. Under what condition(s) is the asymptotic limit in which you worked
indeed a valid approximation for this solution?
(c) Are there any other unstable solutions? (cf. Weibel 1958)
(d) What happens if the electrons are isotropic but ions are not?
† In Q4, you will need the dielectric tensor in terms of a general equilibrium distribution
f0α (vx , vy , vz ). If you are planning to do that question, it may save time (at the price of a
very slight increase in algebra) to do the derivation with a general f0α and then specialise to the
bi-Maxwellian (10.5). You can check your algebra by looking up the result in Krall & Trivelpiece
(1973) or in Davidson (1983).
74 A. A. Schekochihin
(e∗∗ ) If you want a challenge and a test of stamina, work out the case of perturbations
whose wave number is not necessarily in the direction of the anisotropy (k × n̂ 6= 0). Are
the k k n̂ or some oblique perturbations the fastest growing? This is a lot of algebra,
so only do it if you enjoy this sort of thing. The dispersion relation for this general case
appears to be in the Appendix of Ruyer et al. (2015), but they only solve it numerically;
no one seems to have looked at asymptotic limits. This could be the start of a dissertation.
4.∗ Criterion of instability of anisotropic distributions.
Coming soon. See Krall & Trivelpiece (1973, §9.10).
5. Free energy and stability. (a) Starting from the linearised Vlasov–Poisson system
and assuming a Maxwellian equilibrium, show by direct calculation from the equations,
rather then via expansion of the entropy function and the use of energy conservation (as
was done in §4.2), that free energy is conserved:
" #
2
|∇ϕ|2
Z XZ
d 3 3 Tα δfα
d r d v + = 0. (10.8)
dt α
2f0α 8π
This is an exercise in integrating by parts.
(b) Now consider the full Vlasov–Maxwell equations and prove, again for a Maxwellian
plasma plus small perturbations,
" #
2
|E|2 + |B|2
Z XZ
d 3 3 Tα δfα
d r d v + =0 . (10.9)
dt α
2f0α 8π
(c) Consider the same problem, this time with an equilibrium that is not Maxwellian,
but merely isotropic, i.e., f0α = f0α (v), or, in what will prove to be a more convenient
form,
f0α = f0α (εα ), (10.10)
2
where εα = mα v /2 is the particle energy. Find an integral quantity quadratic in per-
turbed fields and distributions that is conserved by the Vlasov–Maxwell system under
these circumstances and that turns into the free energy (10.9) in the case of a Maxwellian
equilibrium (if in difficulty, you will find the answer in, e.g., Davidson 1983 or in Kruskal
& Oberman 1958, which appears to be the original source). Argue that
∂f0α
<0 (10.11)
∂εα
is a sufficient condition for stability of small (δfα f0α , but not necessarily infinitesimal)
perturbations in such a plasma.
6.∗ Fluctuation-dissipation relation for a collisionless plasma. Let us consider a
linear kinetic system in which perturbations are stirred up by an external force, which we
can think of as an imposed (time-dependent) electric field E ext = −∇χ. The perturbed
distribution function then satisfies
∂δfα qα ∂f0α
+ v · ∇δfα − (∇ϕtot ) · = 0, (10.12)
∂t mα ∂v
where ϕtot = ϕ + χ is the total potential, whose self-consistent part, ϕ, obeys the usual
Poisson equation
X Z
2
−∇ ϕ = 4π qα d3 v δfα (10.13)
α
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 75
and the equilibrium f0α is assumed to be Maxwellian.
(a) By considering an initial-value problem for Eqs. (10.12) and (10.13) with zero initial
perturbation, show that the Laplace transforms of ϕtot and χ are related by
χ̂(p)
ϕ̂tot (p) = , (10.14)
(p)
where (p) is the dielectric function given by (3.80).
(b) Consider a time-periodic external force,
χ(t) = χ0 e−iω0 t . (10.15)
Working out the relevant Laplace transforms and their inverses [see Eq. (3.14)], show
that, after transients have decayed, the total electric field in the system will oscillate at
the same frequency as the external force and be given by
χ0 e−iω0 t
ϕtot (t) = . (10.16)
(−iω0 )
(c) Now consider the plasma-kinetic Langevin problem: assume the external force to
be a white noise, i.e., a random process with the time-correlation function
hχ(t)χ∗ (t0 )i = 2D δ(t − t0 ). (10.17)
Show that the resulting steady-state mean-square fluctuation level in the plasma will be
Z +∞
2 D dω
h|ϕtot (t)| i = . (10.18)
π −∞ |(−iω)|2
Here is a reminder of how the standard Langevin problem can be solved using Laplace transforms.
The Langevin equation is
∂ϕ
+ γϕ = χ(t) , (10.20)
∂t
where ϕ describes some quantity, e.g., the velocity of a Brownian particle, subject to a damping
rate γ and an external force χ. In the case of a Brownian particle, χ is assumed to be a white
noise, as per Eq. (10.17). Assuming ϕ(t = 0) = 0, the Laplace-trasform solution of Eq. (10.20) is
χ̂(p)
ϕ̂(p) = . (10.21)
p+γ
Considering first a non-random oscillatory force (10.15), we have
Z ∞
χ0 χ0
χ̂(p) = dt e−pt χ(t) = ⇒ ϕ̂(p) = . (10.22)
0 p + iω 0 (p + γ)(p + iω0 )
76 A. A. Schekochihin
The inverse Laplace transform of ϕ̂ is calculated by shifting the integration contour to large
negative Re p while not allowing it to cross the two poles, p = −γ and p = −iω0 , in a manner
analogous to that explained in §3.1 (Fig. 4) and shown in Fig. 30. The integral is then dominated
by the contributions from the poles:
−iω0 t
e−γt χ0 e−iω0 t
Z i∞+σ
1 e
ϕ(t) = dp ept ϕ̂(p) = χ0 + → as t → ∞,
2πi −i∞+σ −iω0 + γ −γ + iω0 −iω0 + γ
(10.23)
which is quite obviously the right solution of Eq. (10.20) with a periodic force (the second term
in the brackets is the decaying transient needed to enforce the zero initial condition).
In the more complicated case of a white-noise force [Eq. (10.17)],
*Z 2
+
i∞+σ
2 1 pt χ̂(p)
h|ϕ(t)| i = dp e
(2π)2 −i∞+σ p+γ
Z +∞ Z +∞ ∗ 0
1 0 [−i(ω−ω 0 )+2σ]t hχ̂(−iω + σ)χ̂ (−iω + σ)i
= dω dω e , (10.24)
(2π)2 −∞ −∞ (−iω + σ + γ)(iω 0 + σ + γ)
where we have changed variables p = −iω + σ and similarly for the second integral. The corre-
lation function of the Laplace-transformed force is, using Eq. (10.17),
Z ∞ Z ∞ Z ∞
0∗ 0 0∗ 2D
χ̂(p)χ̂∗ (p0 ) = dt dt0 e−pt−p t χ(t)χ∗ (t0 ) = 2D dt e−(p+p )t = ,
0 0 0 p + p0∗
(10.25)
provided Re p > 0 and Re p0 > 0. Then Eq. (10.24) becomes
Z +∞ Z +∞ 0
D 0 e[−i(ω−ω )+2σ]t
h|ϕ(t)|2 i = dω dω
2π 2 −∞ −∞ [−i(ω − ω 0 ) + 2σ] (−iω + σ + γ)(iω 0 + σ + γ)
Z +∞ 0 Z i∞+σ
D e(iω +σ)t 1 ept
= dω 0 0
dp 0
π −∞ (iω + σ + γ) 2πi −i∞+σ (p + iω + σ)(p + γ)
(iω 0 +σ)t −(iω 0 +σ)t
" #
e−γt
Z +∞
D 0 e e
= dω + , (10.26)
π −∞ (iω 0 + σ + γ) −iω 0 − σ + γ γ + iω 0 + σ
where we have reverted to the p variable in one of the integrals and then performed the integra-
tion by the same manipulation of the contour as in Eq. (10.23). We now note that, since there
are no exponentially growing solutions in this system, σ > 0 can be chosen arbitrarily small.
Taking σ → +0 and neglecting the decaying transient in Eq. (10.26), we get, in the limit t → ∞,
D +∞ dω 0 D +∞ dω 0
Z Z
D
h|ϕ(t)|2 i = = = . (10.27)
π −∞ | − iω 0 + γ|2 π −∞ ω 02 + γ 2 γ
Note that, while the integral in Eq. (10.27) is doable exactly, it can, for the case of weak damping,
also be computed via Plemelj’s formula.
Eq. (10.27) is the standard Langevin fluctuation-dissipation relation. It can also be obtained
without Laplace transforms, either by directly integrating Eq. (10.20) and correlating ϕ(t) with
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 77
itself or by noticing that
t
∂ hϕ2 i
Z
+ γhϕ2 i = hχ(t)ϕ(t)i = dt0 −γϕ(t0 ) + χ(t0 )
χ(t) = D, (10.28)
∂t 2 0
where we have used Eq. (10.17) and the fact that hχ(t)ϕ(t0 )i = 0 for t0 6 t, by causality.
Eq. (10.27) is the steady-state solution to the above, but Eq. (10.28) also teaches us that, if we
interpret hϕ2 i/2 as energy, D is the power that is injected into the system by the external force.
Thus, fluctuation-dissipation relations such as Eq. (10.27) tells us what fluctuation energy will
persist in a dissipative system if a certain amount of power is pumped in.
7.∗ Phase-mixing spectrum. Here we study the velocity-space structure of the per-
turbed distribution function δf derived in Q6.
In a plasma where perturbations are constantly stirred up by a force, Landau damping must
be operating all the time, removing energy from ϕ to provide “dissipation” of the injected
power. The process of phase mixing that accompanies Landau damping must then lead to
a certain fluctuation level h|δfm |2 i in the Hermite moments of δf . Lower m’s correspond to
“fluid” quantities: density (m = 0), flow velocity (m = 1), temperature (m = 2). Higher m’s
correspond to finer structure in velocity space: indeed, for m 1, the Hermite polynomials can
be approximated by trigonometric functions,
m/2
√ √
2m πm u2 /2
Hm (u) ≈ 2 cos 2m u − e , (10.33)
e 2
and so the Hermite
√ transform is somewhat analogous to a Fourier transform in velocity space
with “frequency” 2m/vth .
(a) Show that in the kinetic Langevin problem described in Q6(c), the mean square
fluctuation level of the m-th Hermite moment of the perturbed distribution function is
given by
Z +∞ 2
q2 D ζZ (m) (ζ) ω
|δfm (t)|2 = 2 m dω , ζ= , (10.34)
T π2 m! −∞ (−iω) kvth
where Z (m) (ζ) is the m-th derivative of the plasma dispersion function [note Eq. (3.89)].
78 A. A. Schekochihin
√
(b∗∗ ) Show that, assuming m 1 and ζ 2m,
m/2 √
√
2m 2
Z (m) (ζ) ≈ 2π im+1 eiζ 2m−ζ /2 (10.35)
e
and, therefore, that
const
|δfm (t)|2 ≈ √ . (10.36)
m
Thus, the Hermite spectrum of the free energy is shallow and, in particular, the total
free energy diverges—it has to be regularised by collisions. This is a manifestation of a
copious amount of fine-scale structure in velocity space (note also how this shows that
Landau-damped perturbations involve all Hermite moments, not just the “fluid” ones).
Deriving Eq. (10.35) is a (reasonably hard) mathematical exercise: it involves using Eqs. (3.89)
and (10.33) and manipulating contours in the complex plane. This is a treat for those who like
this sort of thing. Getting to Eq. (10.36) will also require the use of Stirling’s formula.
The Hermite order at which the spectrum (10.36) must be cut off due to collisions can be quickly
deduced as follows. We saw in §4.5 that the typical velocity derivative of δf can be estimated
according to Eq. (4.24) and the time it takes for this perturbation to be wiped out by collisions
is given by Eq. (4.29). But,
√ in view of Eq. (10.33), the velocity gradients probed by the Hermite
moment m are of order 2m/vth . The collisional cut off mc in Hermite space can then be
estimated so:
2
2/3
2 ∂ 2 kvth
mc ∼ vth ∼ (kv t
th c ) ∼ . (10.37)
∂v 2 ν
Therefore, the total free energy stored in phase space diverges: using Eqs. (10.32) and (10.36),
Z mc
δf 2
Z
1 1X const
d3 v = |δfm |2 ∼ dm √ ∝ ν −1/3 → ∞ as ν → +0. (10.38)
n 2f0 2 m m
In contrast, the total free-energy dissipation rate is finite, however small is the collision frequency:
estimating the right-hand side of Eq. (4.18), we get
Z mc
√
Z
1 3 δf ∂δf X 2
d v ∼ −ν m |δfm | ∝ ν dm m ∼ kvth . (10.39)
n f0 ∂t c m
Thus, the kinetic system can collisionally produce entropy at a rate that is entirely independent
of the collision frequency.
If you find phase-space turbulence and generally life in Hermite space as fascinating as I
do, you can learn more from Kanekar et al. (2015) (on fluctuation-dissipation relations and
Hermite spectra) and from Schekochihin et al. (2016) and Adkins & Schekochihin (2017) (on
what happens when a nonlinearity strikes).
Figure 31. Initial stable equilibrium distribution and the final outcome of the QL evolution of
a system with Landau-damped electric perturbations.
distribution F (0, v). What is the energy of the waves in this steady state? What is the
lower bound on initial electric energy E(0) below which the perturbations would just
decay without forming a fully-fledged plateau?
(b) Derive the evolution equation for the thermal (nonresonant) bulk of the distribution
and show that it cools during the QL evolution, with the total thermal energy declining
by the same amount as the electric energy of the waves:
Kth (t) − Kth (0) = − [E(0) − E(t)] . (10.40)
Identify where all the energy lost by the thermal particles and the waves goes and thus
confirm that the total energy in the system is conserved. Why, physically, do thermal
particles lose energy?
(c) Show that we must have
E(0) γk δv
(10.41)
ne Te ωpe v
in order for the wave energy to change only by a small fraction before saturating and
3
E(0) γk δv 1
(10.42)
ne Te ωpe v (kλDe )2
in order for the QL evolution to be faster than the damping. Here δv = v2 − v1 and
v ∼ v1 ∼ v2 .
This question requires some nuance in handling the calculation of the QL diffusion coefficient.
In §7.1, we used the expression for δfk , Eq. (7.6), in which only the eigenmode-like part was
retained, while the phase-mixing terms were dropped on the grounds that we could always just
wait long enough for them to be eclipsed by the term containing an exponentially growing factor
eγk t . When we are dealing with damped perturbations, there is no point in waiting because the
exponential term is getting smaller, while the phase-mixing terms do not decay (except by
collisions, see §§4.3 and 4.5, but we are not prepared to wait for that).
Let us, therefore, bite the bullet and use the full expression for the perturbed distribution
function, Eq. (4.23), where we single out the slowest-damped mode and assume that all others,
if any, will be damped fast enough never to produce significant QL effects:
q 1 − e−i(k·v−ωk )t−γk t ∂f0
δfk = i ϕk k· + e−ik·vt (gk + . . .) , (10.43)
m k · v − ωk − iγk ∂v
where “. . . ” stand for any possible undamped, phase-mixing remnants of other modes. When
the solution (10.43) is substituted into Eq. (7.4), where it is multiplied by ϕ∗k and time averaged
80 A. A. Schekochihin
[according to Eq. (2.7)], the second term vanishes because, for resonant particles (k · v ≈ ωk ),
it contains no resonant denominators and so is smaller than the first term, whereas for the
nonresonant particles, it is removed by time averaging (check that this works at least for |γk |t . 1
and indeed beyond that). Keeping only the first term in the expression (10.43), substituting it
into Eq. (7.4) and going through a calculation analogous to that given in Eq. (7.8), we find that
the diffusion matrix is (check this)
q 2 X kk 1 − e−i(k·v−ωk )t−γk t
2
D(v) = 2 |E k | Im , (10.44)
m k2 k · v − ωk − iγk
k
which is a generalisation of the penultimate line of Eq. (7.8). For nonresonant particles, the
phase-mixing term is eliminated by time averaging and we and up with the old result: the last
line of Eq. (7.8). For resonant particles, assuming |γk | |k · v − ωk | ωk ∼ k · k and |γk |t 1,
we may adopt the approximation (4.37), which we have previously used to analyse the structure
of δf . This gives us
q 2 X kk
D(v) = 2 |E k |2 πδ(k · v − ωk ), (10.45)
m k2
k
which is the same result as Eq. (7.16)—including, importantly, the sign, which we would have
gotten wrong had we just mechanically applied Plemelj’s formula to Eq. (7.12) with γk < 0.
This is equivalent to saying that the k integral in Eq. (7.16) should be taken along the Landau
contour, rather than simply along the real line.
Note that the above construction was done assuming |γk |t 1, i.e., all the QL action has to
occur before the initial perturbations decay away (which is reasonable). Note also that there is
nothing above that would not apply to the case of unstable perturbations (γk > 0) and so we
conclude that results of §7, derived formally for γk t 1, in fact also hold on shorter time scales
(γk t 1, but, obviously, still ωk t 1).
9. Quasilinear theory of Weibel instability. (a) Starting from the Vlasov equa-
tions including magnetic perturbations, show that the slow evolution of the equilibrium
distribution function is described by the following diffusion equation:
∂f0 ∂ ∂f0
= · D(v) · , (10.46)
∂t ∂v ∂v
where the QL diffusion tensor is
q2 X v × B ∗k
1 ∗ v × Bk
D(v) = 2 Ek + Ek + (10.47)
m i(k · v − ωk ) + γk c c
k
and ωk and γk are the frequency and the growth rate, respectively, of the fastest-growing
mode.
(b) Consider the example of the low-frequency electron Weibel instability with wave
numbers k parallel to the anisotropy direction [Eq. (10.7)]. Take k = kẑ and B k = Bk ŷ
and, denoting Ωk = eBk /me c (the Larmor frequency associated with the perturbed
magnetic field), show that Eq. (10.46) becomes
∂f0 ∂ ∂f0 ∂f0 ∂ ∂f0
= Dxx + Dxz + Dzz , (10.48)
∂t ∂vx ∂vx ∂vz ∂vz ∂vz
where the coefficients of the QL diffusion tensor are
X γk X 2γk vx vz X γk vx2
Dxx = |Ωk |2 , Dxz = − |Ωk |2 , Dzz = |Ωk |2 .
k2 k 2 vz2 + γk2 k vz2 + γk2
2
k k k
(10.49)
(c) Suppose the electron distribution function f0 is initially bi-Maxwellian, Eq. (10.5),
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 81
with 0 < T⊥ /Tk − 1 1 (as should be the case for this instability to work). As QL
evolution starts, we may define the temperatures of the evolving distribution according to
m(vx2 + vy2 )
Z Z
1 3 1
T⊥ = d v f0 , Tk = d3 v mvz2 f0 . (10.50)
n 2 n
Show that initially, viz., before f0 has time to change shape significantly so as no longer
to be representable as a bi-Maxwellian, the two temperatures will evolve approximately
(using γk kvth ) according to
Thus, QL evolution will lead, at least initially, to the reduction of the temperature
anisotropy, thus weaking the instability (these equations should not be used to trace
T⊥ /Tk − 1 all the way to zero because there is no reason why the QL evolution should
preserve the bi-Maxwellian shape of f0 ).
Note that, even modulo the caveat about the bi-Maxwellian not being a long-term solution,
this does not give us a way to estimate (or even guess) what the saturated fluctuation level
will be. The standard Weibel lore is that saturation occurs when the approximations that were
used to derive the linear theory (Q3) break down, namely, when magnetic field becomes strong
enough to magnetise the plasma, rendering the Larmor scale ρe = vthe /Ωk associated with the
fluctuations small enough to be comparable to the latter’s wavelengths ∼ k−1 . Using the typical
values of k from Eq. (10.7), we can write this condition as follows
√ vthe 1 B 2 /8π
Ωk ∼ kvthe ∼ ∆e ⇔ ≡ ∼ ∆e . (10.52)
de βe ne Te
Thus, Weibel instability will produce fluctuations the ratio of whose magnetic-energy density
to the electron-thermal-energy density (customarily referred to as the inverse of “plasma beta,”
1/βe ) is comparable to the electron pressure anisotropy ∆e . Because at that point the fluctu-
ations will be relaxing this pressure anisotropy at the same rate as they can grow in the first
place [in Eq. (10.51), λ ∼ γk ], the QL approach is not valid anymore.
These considerations are, however, usually assumed to be qualitatively sound and lead people
to believe that, even in collisionless plasmas, the anisotropy of the electron distribution must be
largely self-regulating, with unstable Weibel fluctuations engendered by the anisotropy quickly
acting to isotropise the plasma (or at least the electrons).
This is all currently very topical in the part of the plasma-astrophysics world preoccupied
with collisionless shocks, origin of the cosmic magnetism, hot weakly collisional environments
such as the intergalactic medium (in galaxy clusters) or accretion flows around black holes and
many other interesting subjects.
(d) Eqs. (10.51) say that the total mean kinetic energy,
2
Tk
Z
3 mv
d v f0 = n T⊥ + , (10.53)
2 2
does not change. But fluctuations are generated and grow at the rate γk ! Without much
further algebra, can you tell whether you should therefore doubt the result (10.51)?
ρ—mass density,
u—flow velocity,
p—pressure,
σ—charge density,
j—current density,
E—electric field,
B—magnetic field.
Our immediate objective is to find a set of closed equations that would allow us to
determine all of these quantities as functions of time and space within the fluid.
As this equation holds for any V , however small, it can be converted into a differential
relation
∂ρ
+ ∇ · (ρu) = 0 . (11.2)
∂t
One part of this equation does have to be calculated from some knowledge of the microscopic
properties of the constituent fluid or gas—the viscous stress. For a gas, it is done in kinetic
theory (e.g., Lifshitz & Pitaevskii 1981; Dellar 2016; Schekochihin 2016):
2
Π = −ρν ∇u + (∇u)T − ∇ · u I , (11.6)
3
where ν is the kinematic (Newtonian) viscosity. In what follows, we will never require the explicit
form of Π (except perhaps in §11.10).
In a magnetised plasma (i.e., such that its collision frequency Larmor frequency of the
gyrating charges), the viscous stress is much more complicated and anisotropic with respect to
the direction of the magnetic field: because of their Larmor motion, charged particles diffuse
84 A. A. Schekochihin
differently across and along the field. This gives rise to the so-called Braginskii (1965) stress
(see, e.g., Helander & Sigmar 2005; Parra 2017a).
and if we sum this over all particles (or, to be precise, average over their distribution and
sum over species), we will get
j×B
F = σE + . (11.8)
c
This body force (force density) goes into Eq. (11.5) and so we must know E, B, σ and
j in order to compute the fluid flow u.
Clearly it is a good idea to bring in Maxwell’s equations:
Normally, the resistivity, like viscosity, has to be computed from kinetic theory (see, e.g., He-
lander & Sigmar 2005; Parra 2017a) or tabulated by assidiuous experimentalists. In a magnetised
plasma, the simple form (11.13) of Ohm’s law is only valid at spatial scales longer than the Lar-
mor radii and time scales longer than the Larmor periods of the particles (see, e.g., Goedbloed
& Poedts 2004; Parra 2017b).
Eqs. (11.9–11.13) can be reduced somewhat if we assume (quite reasonably for most
applications) that our fluid flow is non-relativistic. Let us stipulate that all fields evolve
on time scales ∼ τ , have spatial scales ∼ l and that the flow velocity is
l
u∼ c. (11.14)
τ
Then, from Ohm’s law (11.13),
u
E∼ B B, (11.15)
c
so electric fields are small compared to magnetic fields.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 85
In Ampère–Maxwell’s law (11.12),
1 ∂E 11u
c ∂t B u2
∼ c τ c ∼ 2 1, (11.16)
|∇ × B| 1 c
B
l
so the displacement current is negligible (note that at this point we have ordered out
light waves; see Q2 in Kinetic Theory). This allows us to revert to the pre-Maxwell form
of Ampère’s law:
c
j= ∇×B . (11.17)
4π
Thus, the current is no longer an independent field, there is a one-to-one correspondence
j ↔ B.
Finally, comparing the electric and magnetic parts of the Lorentz force (11.8), and
using Gauss’s law (11.9) to estimate σ ∼ E/l, we get
1 2
|σE| E E2 u2
∼ l ∼ 2 ∼ 2 1. (11.18)
1 1c 2 B c
j×B B
c cl
Thus, the MHD body force is
j×B (∇ × B) × B
F = = . (11.19)
c 4π
This goes into Eq. (11.5) and we note with relief that σ, j and E have all fallen out of
the momentum equation—we only need to know B.
Thus, we learn that the Lorentz force consists of two distinct parts (Fig. 33):
• curvature force, so called because b · ∇b is the vector curvature of the magnetic field
line—the implication being that field lines, if bent, will want to straighten up;
• magnetic pressure, whose presence implies that field lines will resist compression or
rarefication (the field wants to be uniform in strength).
Note that both forces act perpendicularly to B, as they must, since magnetic field never
exerts a force along itself on a charged particle [Eq. (11.7)].
So this is the effect of the field on the fluid. What is the effect of the fluid on the field?
∂B
= ∇ × (u × B) + η∇2 B . (11.24)
∂t | {z } | {z }
advection diffusion
Note that if ∇ · B = 0 is satisfied initially, any solution of Eq. (11.24) will remain
divergence-free at all times.
Generally speaking, when flow velocities are large/distances are large/resistivities are
low, Rm 1 and it makes sense to discuss “ideal MHD,” i.e., the limit η → 0. In fact,
η often needs to be brought back in to deal with instances of large ∇B, which arise
naturally from solutions of ideal MHD equations (see §11.10, §15.2 and Parra 2017a),
but let us consider the ideal case for now to understand what the advective part of the
induction equation does to B.
we arrive at
d B B
= · ∇u . (11.32)
dt ρ ρ
Let us compare the evolution of the vector B/ρ with the evolution of an infinitesimal La-
grangian separation vector in a moving fluid: the convective derivative is the Lagrangian
time derivative, so
d
δr(t) = u(r + δr) − u(r) ≈ δr · ∇u. (11.33)
dt
Thus, δr and B/ρ satisfy the same equation. This means that if two fluid elements are
initially on the same field line,
B
δr = const , (11.34)
ρ
then they will stay on the same field line, q.e.d.
This means that in MHD, the fluid flow will be entraining the magnetic-field lines with
it—and, as saw in §11.4, the field lines will react back on the fluid:
—when the fluid tries to bend the field, the field will want to spring back,
—when the fluid tries to compress or rarefy the field, the field will resist as if it possessed
(perpendicular) pressure.
This is the sense in which MHD fluid is “elastic”: it is threaded by magnetic-field lines,
which move with it and act as elastic bands.
(dS ≡ n̂ dS, where n̂ is a unit normal pointing out of the surface). The flux Φ depends on
the loop ∂S, but not on the choice of the surface spanning it. Indeed, if we consider
R two
surfaces, S1 and S2 , spanning the same loop ∂S (Fig. 35b) and define Φ1,2 = S1,2 B · dS,
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 89
(a) (b)
Figure 35. Magnetic flux through a loop ∂S (a) is independent of the surface spanning the
loop (b).
Alfvén’s Theorem. Flux through any loop moving with the fluid is conserved.
Proof. Let S(t) be a surface spanning the loop at time t. If the loop moves with the
fluid (Fig. 36), at the slightly later time t + dt it is spanned (for example) by the surface
∂B
∂t | ·{z
+ u ∇B} = B
| ·{z∇u} − | {z· u} ,
B∇ (11.44)
advection stretching compression
the three terms in it are responsible for, in order, advection of the field by the flow (i.e.,
the flow carrying the field around with it), “stretching” (amplification) of the field by
velocity gradients that make fluid elements longer and, finally, compression or rarefication
of the field by convergent or divergent flows (unless ∇·u = 0, as it is in an incompressible
fluid).
Hence arises the famous problem of MHD dynamo: are there fluid flows that lead to
sustained amplification of the magnetic fields? The answer is yes—but the flow must be
3D (the absence of dynamo action in 2D is a theorem, the simplest version of which is due
to Zeldovich 1956; see Q4). Magnetic fields of planets, stars, galaxies, etc. are all believed
to owe their origin and sustainement to this effect. This topic requires (and merits) a
more detailed treatment (§11.10), but for now let us flag two important aspects:
• resistivity, however small, turns out to be impossible to neglect because large gradi-
ents of B appear as the field is advected by the flow;
• the amplification of the field is checked by the Lorentz force once the field is strong
enough that it can act back on the flow, viz., when their energy densities become com-
parable:
B2 ρu2
∼ . (11.45)
8π 2
Let us summarise the equations that we have derived so far, namely Eqs. (11.2), (11.5)
and (11.24), expressing conservation of
∂ρ
mass + ∇ · (ρu) = 0, (11.46)
∂t
∂u (∇ × B) × B
momentum ρ + u · ∇u = −∇p − ∇ · Π +
∂t 4π
B2
BB
= −∇ · p + I− +Π , (11.47)
8π 4π
| {z }
total stress
∂B
and flux = ∇ × (u × B) + η∇2 B. (11.48)
∂t
To complete the system, we need an equation for p, which has to come from the one
conservation law that we have not yet utilised: conservation of energy.
The total energy density is
ρu2 p 2
E B2
ε= + + + , (11.49)
2
|{z} γ−1 8π
|{z} 8π
|{z}
| {z }
kinetic internal electric magnetic
where the electric energy can (and, for consistency with §11.3, must) be neglected because
E 2 /B 2 ∼ u2 /c2 1. We follow the same logic as we did in §§11.1 and 11.2:
ρu2
Z Z Z
d p
d3 r ε = − + u · dS − [(p I + Π) · u] · dS
dt 2 γ−1
| V {z } | ∂V {z } | ∂V
{z }
energy inside fluid flow carrying its work done on boundary
a volume of kinetic and internal energy by pressure and viscous
fluid through boundary stress
Z Z
c
− q · dS − (E × B) · dS . (11.50)
4π
| ∂V {z } | ∂V {z }
heat flux Poynting flux
Like the viscous stress Π, the heat flux q must be calculated kinetically (in a plasma) or
tabulated (in an arbitrary complicated substance). In a gas, q = −κ∇T , but it is more
complicated in a magnetised plasma (see, e.g., Braginskii 1965; Helander & Sigmar 2005;
Parra 2017a).
Note that the magnetic energy and the work done by the Lorentz force are not included
in the first two terms on the right-hand side of Eq. (11.50) because all of that must already
be correctly acounted for by the Poynting flux. Indeed, since cE = −u × B + η∇ × B
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 93
[Eq. (11.13), with η renamed as in Eq. (11.24)], we have
B2
Z Z Z 2
c B BB
(E × B) · dS = u · dS + I− · u · dS
∂V 4π 8π 8π 4π
| ∂V {z } | ∂V {z }
magnetic energy work done by Maxwell stress
flow
(∇ × B) × B
Z
+ η · dS . (11.51)
4π
| ∂V {z }
resistive slippage accounting
for field not being precisely
frozen into flow
ρu2 B2 ρu2
∂ p γ
+ + =−∇· u+ pu + Π · u + q
∂t 2 γ−1 8π 2 γ−1
B 2 I − BB
(∇ × B) × B
+ ·u+η . (11.52)
4π 4π
It remains to separate the evolution equation for p by using the fact that we know the
equations for ρ, u and B and so can deduce the rates of change of the kinetic and
magnetic energies.
∂ ρu2 u2 ∂ρ ∂u
= + ρu ·
∂t 2 2 ∂t ∂t
u2 u2 B2
BB
= − ∇ · (ρu) − ρu · ∇ − u · ∇ · p + I− +Π
2 2 8π 4π
" #
ρu2 B2
BB
= −∇ · u + pu
+ I−
· u + ·
Π u
2 8π 4π
2
B BB
+ p∇ · u + I− : ∇u + Π : ∇u . (11.53)
8π 4π
| {z } | {z } | {z }
compressional energy exchange viscous
work with magnetic field dissipation
The flux terms (energy flows and work by stresses on boundaries) that have been crossed
out cancel with corresponding terms in Eq. (11.52) once Eq. (11.53) is subtracted from
it.
94 A. A. Schekochihin
11.11.2. Magnetic Energy
Using the induction equation (11.48),
∂ B2 B
· −u · ∇B + B · ∇u − B∇ · u + η∇2 B
=
∂t 8π 4π " #
B2 B2 |∇ × B|2
(∇ ×
B) ×
B
BB
= −∇ · u + η − I− : ∇u − η .
8π 4π 8π 4π 4π
| {z } | {z }
energy exchange Ohmic
with velocity field dissipation
(11.54)
Again, the crossed out flux terms will cancel with corresponding terms in Eq. (11.52).
The metamorphosis of the resistive term into a flux term and an Ohmic dissipation
term is a piece of vector algebra best checked by expanding the divergence of the flux
term. Finally, the u-to-B energy exchange term (penultimate on the rght-hand side)
corresponds precisely to the B-to-u exchange term in Eq. (11.53) and cancels with it if
we add Eqs. (11.53) and (11.54).
d p γ |∇ × B|2
= −∇ · q − p∇ · u − Π : ∇u + η . (11.55)
dt γ − 1 γ−1 4π
| {z } | {z } | {z } | {z } | {z }
advection of heat flux compressional viscous Ohmic
internal heating heating heating
energy
A further rearrangement and the use of the continuity equation (11.46) to express ∇·u =
−d ln ρ/dt turn Eq. (11.55) into
|∇ × B|2
d p γ−1
ln γ = −∇ · q − Π : ∇u + η . (11.56)
dt ρ p 4π
This form of the thermal-energy equation has very clear physical content: the left-hand
side represents advection of the entropy of the MHD fluid by the flow—each fluid ele-
ment behaves adiabatically, except for the sundry non-adiabatic effects on the right-hand
side. The latter are the heat flux in/out of the fluid element and the dissipative (viscous
and resistive) heating, leading to entropy production. Note that the form of the viscous
stress Π ensures that the viscous heating is always positive [see, e.g., Eq. (11.6)]. In these
Lectures, we will, for the most part, focus on ideal MHD and so use the adiabatic version
of Eq. (11.56), with the right-hand side set to zero.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 95
Let us reiterate the equations of ideal MHD, now complete:
∂ρ
+ ∇ · (ρu) = 0, (11.57)
∂t
∂u (∇ × B) × B
ρ + u · ∇u = −∇p + , (11.58)
∂t 4π
∂B
= ∇ × (u × B) , (11.59)
∂t
∂ p
+u·∇ = 0. (11.60)
∂t ργ
In what follows, we shall study various solutions and symptotic regimes of these rather
nice equations.
We will turn to more nontrivial equilibria in due course, but first we shall study this one
carefully—because it is very generic in the sense that many other, more complicated,
equilibria locally look just like this.
∂ξ
ρ = ρ0 + δρ, p = p0 + δp, u= , B = B0 ẑ + δB, (12.2)
∂t
where we have introduced the fluid displacement field ξ†. To start with, we consider all
† Thinking in terms of displacements makes sense in MHD but not in (homogeneous) hydro-
dynamics because in the latter case, just displacing a fluid element produces no back reaction,
whereas in MHD, because magnetic fields are frozen into the fluid and are elastic, displacing
fluid elements causes magnetic restoring forces to switch on.
96 A. A. Schekochihin
where k and ⊥ denote projections onto the direction (z) of B 0 and onto the plane (x, y)
perpendicular to it, respectively. Eqs. (12.5) tell us that parallel displacements produce no
perturbation of the magnetic field—obviously not, because the magnetic field is carried
with the fluid flow and nothing will happen if you displace a straight uniform field parallel
to itself.
The physics of magnetic-field perturbations become clearer if we observe that
δB δ(Bb) δB
= = δb + ẑ . (12.6)
B0 B0 B0
δB ⊥ δBk δB
= δb, = . (12.7)
B0 B0 B0
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 97
∂2ξ
= c2s ∇∇ · ξ + vA
2
∇⊥ ∇⊥ · ξ ⊥ + ∇2k ξ ⊥ ,
2
(12.9)
∂t
where two special velocities have emerged:
r
γp0 B0
cs = , vA = √ , (12.10)
ρ0 4πρ0
the sound speed and the Alfvén speed, respectively. The former is familiar from fluid
dynamics, while the latter is another speed, arising in MHD, at which perturbations can
travel. We shall see momentarily how this happens.
Let us seek wave-like solutions of Eq. (12.9), ξ ∝ exp(−iωt + ik · r). For such pertur-
bations,
ω 2 ξ = c2s kk · ξ + vA
2
k⊥ k⊥ · ξ ⊥ + kk2 ξ ⊥ .
(12.11)
Without loss of generality, let k = (k⊥ , 0, kk ) (i.e., by definition, x is the direction of k⊥ ;
see Fig. 40). Then Eq. (12.11) becomes
ω 2 ξx = c2s k⊥ (k⊥ ξx + kk ξk ) + vA
2 2
k ξx , (12.12)
2 2 2
ω ξy = vA kk ξy , (12.13)
2
ω ξk = c2s kk (k⊥ ξx + kk ξk ). (12.14)
98 A. A. Schekochihin
(a) (b) k⊥ 6= 0
Figure 41. Alfvén waves.
ω = ±kk vA , (12.19)
The two “+” solutions are the “fast magnetosonic waves” and the two “−” ones are the
“slow magnetosonic waves”.
Since both sound and Alfvén speeds are involved, it is obvious that the key parameter
demarcating different physical regimes will be their ratio, or, conventionally, the ratio of
the thermal to magnetic energies in the MHD medium, known as the plasma beta:
p0 2 c2s
β= = 2 . (12.24)
B02 /8π γ vA
The magnetosonic waves can be conveniently summarised by the so-called Friedricks
diagram, a graph of Eq. (12.23) in polar coordinates where the phase speed ω/k is radius
and the angle is θ, the direction of propagation with respect to B 0 (Fig. 42).
Clearly, magnetosonic waves contain perturbations of both the magnetic field and of
the “hydrodynamic” quantities ρ, p, u, but working them all out for the case of general
oblique propagation (θ ∼ 1) is a bit messy. The physics of what is going on is best
understood via a few particular cases.
• another Alfvén wave, this time with perturbation in the x direction (which, however,
is not physically different from the y direction when k⊥ = 0):
ω 2 ξx = kk2 vA
2
ξx ⇒ ω = ±kk vA , (12.25)
ξ = ξx x̂, δρ = 0, δp = 0, δB = 0, δb = ikk ξx x̂ (12.26)
(at high β, this is the slow wave, at low β, this is the fast wave);
100 A. A. Schekochihin
This is just the magnetically enhanced sound wave that was considered in §12.1.4. The
small corrections to it due to kk /k are not particularly interesting.
The lower sign in Eq. (12.31) gives the slow wave
cs vA
ω = ±kk p 2
, (12.33)
c2s + vA
which is more interesting. Let us find the corresponding eigenvector: from Eq. (12.14),
ξx kk c2s
=− 2 1 . (12.35)
ξk k⊥ c2s + vA
Using this equation together with Eqs. (12.15–12.18), we find that the perturbations in
102 A. A. Schekochihin
calculations in §12.1.5 are also valid, with the following simplifications arising from vA
being negligible compared to cs .
The upper sign in Eq. (12.31) again gives us the fast wave, which, this time, is a pure
sound wave:
ω = ±kcs . (12.42)
ω = ±kk vA . (12.43)
Because the slow wave’s dispersion relation in this limit looks exactly like that of an
Alfvén wave [Eq. (12.19)], it is called the pseudo-Alfvén wave. The similarity is deceptive
as the nature of the perturbation (the eigenvector) is completely different. Substituting
ω 2 = kk2 vA
2
into Eq. (12.14), we find
2
vA
k⊥ ξx + kk ξk = kk ξk kk ξk . (12.44)
c2s
This just says that, to lowest order in 1/β, ∇ · ξ = 0, i.e., the perturbations are incom-
pressible. In contrast to the anisotropic case [Eq. (12.35)], the perpendicular and parallel
displacements are now comparable (at least in general, assuming kk ∼ k⊥ ):
ξx kk
=− . (12.45)
ξk k⊥
Also in contrast to the anisotropic case, the density and pressure perturbations are now
104 A. A. Schekochihin
vanishingly small, but the field can be bent as well as compressed:
δρ v2
= −i(k⊥ ξx + kk ξk ) = −i A kk ξk → 0, (12.46)
ρ0 c2s
δp δρ
=γ → 0, (12.47)
p0 ρ0
kk
δb = ikk ξx x̂ = −i kk ξk , (12.48)
k⊥
δB
= −ik⊥ ξx = ikk ξk . (12.49)
B0
The δB and δb perturbations are in counter-phase, as are ξk and ξx (Fig. 45). It is easy
to check that pressure balance (12.40) is again maintained by these perturbations.
In the more general case of oblique propagation (kk ∼ k⊥ ) and finite beta (β ∼ 1),
the fast and slow magnetosonic waves generally have comparable frequencies and contain
perturbations of all relevant fields, with the fast waves tending to have the perturbations
of the thermal and magnetic pressure in phase and slow waves in counter-phase (Fig. 46).
kk uy uy p
|δb| ∼ kk ξy ∼ ∼ ∼ Ma 1 + β. (12.53)
ω vA
ω2 v2 kk uk
∇ · u ∼ ω(k⊥ ξx + kk ξk ) ∼ 2
ωξk ∼ 2 A 2 kk uk ∼ . (12.54)
kk cs cs + vA 1+β
Thus, the divergence of the flow velocity is small (the dynamics are incompressible) in
all three of our (potentially) small parameters:
∇·u kk 1
p ∼ Ma. (12.55)
2 2
k cs + vA k 1+β
106 A. A. Schekochihin
From this, we can immediately obtain an ordering for the density and pressure pertur-
bations: using Eqs. (12.3), (12.4), (12.33) and Eq. (12.54) [cf. Eqs. (12.36) and (12.46)],
δρ δp ∇·u Ma
∼ ∼∇·ξ ∼ ∼√ . (12.56)
ρ0 p0 ω β
The magnetic-field-strength (magnetic-pressure) perturbation is, using Eqs. (12.39) and
(12.33) [cf. Eq. (12.49)],
δB c2 kk uk p
∼ k⊥ ξx ∼ 2 s 2 ∼ β Ma, (12.57)
B0 cs + vA ω
or, perhaps more straightforwardly, from pressure balance (12.40) and using Eq. (12.56),
δB β δp p
=− ∼ β Ma. (12.58)
B0 2 p0
Finally, in a similar fashion, using Eqs. (12.17) and (12.57) [cf. Eqs. (12.38) and (12.48)],
we find
kk p
|δb| ∼ kk ξx ∼ β Ma (12.59)
k⊥
for slow-wave-like perturbations. Note that in all interesting limits this is overridden by
the Alfvénic perturbations (12.53).
u |δb| 1 δB p δρ p δp kk 1
Ma ≡ p ∼√ ∼√ ∼ β ∼ β ∼ √ 1 (12.62)
c2s + 2
vA 1+β β B ρ0 p0 k 1+β
† In the context of MHD turbulence theory, this principle, applied at each scale, is known as
the critical balance (see §12.4).
‡ Note that it is not absolutely necessary to work out detailed linear theory of a set of
equations in order to be able to construct such orderings: it is often enough to know roughly
where you are going and simply balance terms representing the physics that you wish to keep
(or expect to have to keep). An example of this approach is given in §12.2.9.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 107
and ω ∼ ku. The ordering can be achieved either in the limit of kk /k 1 or 1/β 1,
or both. Note that if one of these parameters is small, the other can be order unity or
even large (as long as it is not larger than the inverse of the small one).
The case of anisotropic perturbations and arbitrary β applies in a broad range of
plasmas, from magnetically confined fusion ones (tokamaks, stellarators) to space (e.g.,
the solar corona or the solar wind). We shall consider the implications of this ordering
in §12.3.
The case of high β applies, e.g., to high-energy galactic and extragalactic plasmas.
It is the direct generalisation to MHD of incompressible Navier–Stokes hydrodynamics,
i.e., in this case, all one needs to do is solve MHD equations assuming ρ = const and
∇ · u = 0. We shall consider this case now.
u ω 1 δB p δρ δp Ma
∼ ∼ √ ∼ Ma, |δb| ∼ ∼ β Ma ∼ 1, ∼ ∼ √ ∼ Ma2 .
cs kcs β B0 ρ0 p0 β
(12.63)
Thus, the density and pressure perturbations are minuscule, while magnetic perturbations
are order unity—magnetic fields are relatively easy to bend (i.e., subsonic motions can
tangle the field substantially in this regime). Because of this, it will not make sense to
split B into B 0 and δB explicitly, we will treat the magnetic field as a single field, with
no need for a strong mean component.
Let us examine the MHD equations (11.57–11.60) under the ordering (12.63).
Since ω ∼ ku, the convective derivative d/dt = ∂/∂t + u · ∇ survives intact in all
equations, allowing the advective nonlinearity to enter.
The continuity equation (11.57) simply reiterates our earlier statement that the velocity
field is divergenceless to lowest order:
1 dρ δρ
∇·u=− ∼ω ∼ Ma3 kcs → 0. (12.64)
ρ dt ρ0
The momentum equation (11.58) becomes
2
B2
δρ ∂u c δp B · ∇B
1 + + u · ∇u = −∇ s + + . (12.65)
ρ0 ∂t
|
γ p0
{z
8πρ0
}
4πρ0
≡ p̃
The density perturbation in the left-hand side is ∼ Ma2 and so negligible compared to
unity. The remaining terms in this equation are all the same order (∼ Ma2 kc2s ) and so
they must all be kept. The total “pressure” p̃ is determined by enforcing ∇ · u = 0
[Eq. (12.64)]. Namely, our equations are
∂u
+ u · ∇u = −∇p̃ + B · ∇B , (12.66)
∂t
where
∇2 p̃ = −∇∇ : (uu − BB) (12.67)
√
and we have rescaled the magnetic field to velocity units, B/ 4πρ0 → B.
In the induction equation, best written in the form (11.27), all terms are the same
order ∼ kuB ∼ Ma kcs B except the one containing ∇ · u, which is ∼ Ma3 kcs B and so
108 A. A. Schekochihin
must be neglected. We are left with
∂B
+ u · ∇B = B · ∇u . (12.68)
∂t
Finally, the internal-energy equation (11.60), which, keeping only the lowest-order
terms, becomes
∂ δp δρ
+u·∇ −γ = 0, (12.69)
∂t p0 ρ0
can be used to find δρ/ρ0 , once δp/p0 = γ(p̃ − B 2 /2)/c2s is calculated from the solution of
Eqs. (12.66–12.68). Note that δρ/ρ0 is merely a spectator quantity, not required to solve
Eqs. (12.66–12.68), which form a closed set.
Eqs. (12.66–12.68) are the equations of incompressible MHD (let us call it iMHD).
Note that while they have been obtained in the limit of β 1, all β dependence has
disappeared from them—basically, they describe subsonic dynamics on top of an infinite
heat bath.
Exercise 12.1. Show that iMHD conserves the sum of kinetic and magnetic energies,
2
B2
Z
d u
d3 r + = 0. (12.70)
dt 2 2
Exercise 12.2. Check that you can obtain the right waves, viz., Alfvén (§12.1.1) and pseudo-
Alfvén (§12.1.6), directly from iMHD.
∂Z ±
+ Z ∓ · ∇Z ± = −∇p̃ (12.72)
∂t
and, since ∇ · Z ± = 0,
∇2 p̃ = −∇∇ : Z + Z − . (12.73)
Thus, one can think of iMHD as representing the evolution of two incompressible “velocity
fields” advecting each other.
Let us restore the separation of the magnetic field into its mean and perturbed parts,
B = B 0 + δB = vA ẑ + δB (recall that B is in velocity units). Then
Z ± = ±vA ẑ + δZ ± (12.74)
and Eq. (12.72) becomes
∂δZ ±
∓ vA ∇k δZ ± + δZ ∓ · ∇δZ ± = −∇p̃ . (12.75)
∂t
12.2.8. Cross-Helicity
Eqs. (12.72) manifestly support two conservation laws:
|Z ± |2
Z
d
d3 r = 0, (12.78)
dt 2
i.e., the energy of each Elsasser field is individually conserved. This can be reformulated
as conservation of the total energy,
1 |Z + |2 |Z − |2
2
B2
Z Z
d d u
d3 r + = d3 r + = 0, (12.79)
dt 2 2 2 dt 2 2
and of a new quantity, known as the cross-helicity:
1 |Z + |2 |Z − |2
Z Z
d d
d3 r − = d3 r u · B = 0 . (12.80)
dt 2 2 2 dt
Exercise 12.3. To see why we neededR incompressibility to get this new conservation law, work
out the time evolution equation for d3 r u · B from the general (compressible) MHD equations
and hence the condition under which the cross-helicity is conserved.
kH 1 . (12.85)
After the equilibrium pressure balance is subtracted from Eq. (12.81), this equation becomes,
under any ordering in which δρ ρ0 ,
B2
du B · ∇B
ρ0 = −∇ δp + + − δρgẑ . (12.86)
dt 8π 4π
The last term is the buoyancy (Archimedes) force. In order for this new feature to give rise to
any nontrivial new physics, it must be ordered comparable to all the other terms in the equation:
using Eq. (12.83) to express g ∼ p0 /ρ0 H, we find
δρ δp δp
δρg ∼ kδp ⇒ ∼ kH , (12.87)
ρ0 p0 p0
2
kB δρ kH
δρg ∼ ⇒ ∼ 1 ⇒ β kH 1. (12.88)
4π ρ0 β
So we learn that the density perturbations must now be much larger than the pressure pertur-
bations, but, in order for the former to remain small and for the magnetic field to be in the
game, β must be high (it is in anticipation of this that we did not split B into B 0 and δB,
expecting them to be of the same order).
Let us now expand the internal-energy equation (11.60) in small density and pressure pertur-
bations. Denoting s = p/ργ = s0 (z) + δs (entropy density) and introducing the entropy scale
height
1 d ln s0 1 γ
≡ =− + (12.89)
Hs dz Hp Hρ
(assumed positive), we find†
d δs uz δs δp δρ δρ
=− , = −γ ≈ −γ . (12.90)
dt s0 Hs s0 p0 ρ0 ρ0
The last, approximate, expression follows from the smallness of pressure perturbations [Eq. (12.87)].
This then gives us
d δρ uz
= . (12.91)
dt ρ0 γHs
But, on the other hand, the continuity equation (11.57) is
d δρ uz 1 1 uz ∇·u 1
= −∇ · u + ⇒ ∇ · u = uz − = ⇒ ∼ 1.
dt ρ0 Hρ Hρ γHs γHp ku kH
(12.92)
Thus, the dynamics are incompressible again and the role of the continuity equation is to tell
us that we must find δp from the momentum equation (12.86) by enforcing ∇ · u = 0 to lowest
† We are able to take equilibrium quantities in and out of spatial derivatives because kH 1
and the perturbations are small.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 111
order. The difference with iMHD (§12.2.6) is that δρ/ρ now participates in the dynamics via the
buoyancy force and must be found self-consistently from Eq. (12.91).
Finally, we rewrite our newly found simplified system of equations for a stratified, high-β
atmosphere, in the following neat way:
∂u
+ u · ∇u = −∇p̃ + B · ∇B + aẑ, (12.93)
∂t
∂a
∇2 p̃ = −∇∇ : (uu − BB) + , (12.94)
∂z
∂a cs
+ u · ∇a = −N 2 uz , N= p , (12.95)
∂t γ Hs Hp
∂B
+ u · ∇B = B · ∇u, (12.96)
∂t
√
where we have rescaled B/ 4πρ0 → B and denoted the Archimedes acceleration
δρ δρ c2s
a=− g=− , (12.97)
ρ0 ρ0 γHp
a quantity also known as the buoyancy of the fluid. We shall call Eqs. (12.93–12.96) the equations
of stratified MHD (SMHD).
A new frequency N , known as the Brunt–Väisälä frequency, has appeared in our equations.†
In order for all the linear and nonlinear time scales that are present in our equations to coexist
legitimately within our ordering, we must demand that the Alfvén, Brunt–Väisälä and nonlinear
time scales all be comparable:
1 1
kvA ∼ N ∼ ku ⇒ √ ∼ ∼ Ma . (12.98)
β kH
This gives us a relative ordering between all the small parameters that have appeared so far,
including the new one, 1/kH. Using Eq. (12.91) and recalling Eq. (12.87), let us summarise the
ordering of the perturbations:
u δρ δp δB
∼ ∼ Ma, ∼ Ma2 , |δb| ∼ ∼ 1. (12.99)
cs ρ0 p0 B0
The difference with the iMHD high-β ordering (12.63) is that the density perturbations have
now been promoted to dynamical relevance, thankfully without jeopardising incompressibility
(i.e., still ordering out the sonic perturbations). The ordering (12.99) can be thought of as a
generalisation to MHD of the Boussinesq approximation in hydrodynamics.
u⊥ uk δB δρ δp ω kk
Ma ∼ ∼ ∼ |δb| ∼ ∼ ∼ ∼ ∼ 1 . (12.100)
cs cs B0 ρ0 p0 k⊥ cs k⊥
Starting again with the continuity equation (11.57), dividing through by ρ0 and order-
† N is real because we assumed Hs > 0 (a “stably stratified atmosphere”), otherwise the at-
mosphere becomes convectively unstable—this happens when the equilibrium entropy decreases
upwards (cf. §14.3, Q8 and Q5c).
112 A. A. Schekochihin
ing all terms, we get
∂ δρ
δρ
+ u⊥ · ∇⊥ + uk
· ∇k =− 1+ ∇⊥ · u⊥ + ∇k uk . (12.101)
∂t
|{z} | {z } ρ
| {z } |{z}0 ρ
0 | {z } | {z }
|{z}
Ma Ma Ma2
Ma Ma Ma2 Ma | {z }
| {z }
2 Ma
Ma
Thus, to lowest order, the perpendicular velocity field is 2D-incompressible:
O(Ma) : ∇⊥ · u⊥ = 0 . (12.102)
u⊥ = ẑ × ∇⊥ Φ . (12.104)
Similarly, for the magnetic field, we have
0 = ∇ · B = ∇⊥ · δB ⊥ +
∇k k ≈ ∇⊥ · δB ⊥ ,
δB (12.105)
| {z } | {z }
Ma Ma2
so δB ⊥ is also 2D-solenoidal and can be written in terms of a flux function:
δB
√ ⊥ = ẑ × ∇⊥ Ψ . (12.106)
4πρ0
√
Note that Ψ = −Ak / 4πρ0 , the parallel component of the vector potential.
Thus, Alfvénically polarised perturbations, u⊥ and δB ⊥ (see §12.1.1), can be described
by two scalar functions, Φ and Ψ. Let us work out the evolution equations for them.
∂Ψ ∂Φ
+ {Φ, Ψ} = vA . (12.112)
∂t ∂z
Turning now to the momentum equation (11.58), taking its perpendicular part and
dividing through by ρ ≈ ρ0 , we get
B2
2
du⊥ 1 B · ∇δB ⊥ cs δp 2 δB 2 δB ⊥
= −∇⊥ p + + = − ∇⊥ + vA + vA b·∇ .
dt ρ 0 8π 4π γ p 0 B 0 B
| {z } | {z } | {z 0 }
Ma2 Ma Ma2
(12.113)
To lowest order,
c2s δp v 2 δB
2 δB δp
O(Ma) : ∇⊥ + vA =0 ⇒ = −γ A . (12.114)
γ p0 B0 p0 c2s B0
This is a statement of pressure balance, which is physically what has been expected [see
Eq. (12.41)] and which will be useful in §12.3.2. In the next order, Eq. (12.113) contains
the perpendicular gradient of the second-order contribution to the total pressure. To
avoid having to calculate it, we take the curl of Eq. (12.113) and thus obtain
du⊥ δB ⊥
O(Ma2 ) : ∇⊥ × = vA2
∇⊥ × b · ∇ . (12.115)
dt B0
Finally, using again Eqs. (12.104), (12.106), (12.109) and (12.111) in Eq. (12.115), some
slightly tedious algebra leads us to
∂ 2 ∂
∇ Φ + {Φ, ∇2⊥ Φ} = vA ∇2⊥ Ψ + {Ψ, ∇2⊥ Ψ} . (12.116)
∂t ⊥ ∂z
Note that ∇2⊥ Φ is the vorticity of the flow u⊥ and so the above equation is the MHD
generalisation of the 2D Euler equation.
To summarise the equations (12.116) and (12.112) in their most compact form, we
have
d 2
∇ Φ = vA b · ∇∇2⊥ Ψ, (12.117)
dt ⊥
dΨ ∂Φ
= vA , (12.118)
dt ∂z
† Another easy route to this equation is to start from the induction equation in the form
(11.59), let B = ∇ × A, “uncurl” Eq. (11.59) and take the z component of the resulting
evolution equation for A.
114 A. A. Schekochihin
where the convective time derivative d/dt and the parallel spatial derivative b · ∇ are
given by Eqs. (12.109) and (12.111), respectively. Beautifully, these nonlinear equations
describing Alfvénic perturbations have decoupled completely from everything else: we do
not need to know δρ, δp, uk or δB in order to solve for u⊥ and δB ⊥ . Alfvénic dynamics
are self-contained.
Eqs. (12.117) and (12.118) are called the Equations of Reduced MHD (RMHD). They
were originally derived in the context of tokamak plasmas (Kadomtsev & Pogutse 1974;
Strauss 1976) and are extremely popular as a simple paradigm for MHD is a strong guide
field—not just in tokamaks, but also in space†.
where all terms are O(Ma2 ), δBk ≈ δB to leading order and we used Eq. (12.103) to
express ∇ · u. The derivatives d/dt and b · ∇ contain the nonlinearities involving Φ and
Ψ, which we already know from Eqs. (12.117) and (12.118).
To find an equation for uk , we take the z component of the momentum equation
(11.58):
B 2
B · ∇δBk
duk 1 ∂ duk 2 δB
= p
+
+ ⇒ = vA b·∇ . (12.120)
dt
|{z} ρ 0 ∂z
8π 4π dt B 0
| {z } | {z }
Ma2 Ma3 Ma2
The parallel pressure gradient is O(Ma3 ) because there is pressure balance [Eq. (12.114)]
to lowest order.
Finally, let us bring in the energy equation (11.60), as yet unused. To leading order,
it is
2
d δs d δp δρ d δρ vA δB
= −γ =0 ⇒ + 2 =0 , (12.121)
dt s0 dt p0 ρ0 dt ρ0 cs B0
where, to obtain the final version of the equation, we substituted Eq. (12.114) for δp/p0 .
Eqs. (12.119–12.121) are a complete set of equations for δB, uk and δρ, given Φ and
Ψ. These equations are linear in the Lagrangian frame associated with the Alfvénic
perturbations, provided the parallel distances are measured along perturbed field lines.
Physically, they tell us that slow waves propagate along perturbed field lines and are
passively (i.e., without acting back) advected by the perpendicular Alfvénic flows.
Exercise 12.4. Check that the linear relationships between various perturbations in a slow
wave derived in §12.1.5 are manifest in Eqs. (12.119–12.121).
† In the latter context, they are used most prominently as a description of Alfvénic tur-
bulence at small scales (see §12.4), for which the RMHD equations can be shown to be the
correct description even if the plasma is collisionless and in general requires kinetic treatment
(Schekochihin et al. 2009; Kunz et al. 2015).
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 115
In what follows, when we refer to RMHD, we will mean all five equations (12.117–
12.118) and (12.119–12.121).
∂δZk± c s vA ∂δZk±
∓p =
∂t 2
c2 + vA ∂z
"s ! ! #
1 1 1
− 1∓ p {ζ +
, δZk± } + 1± p {ζ −
, δZk± } .
2 2 /c2
1 + vA s
2 /c2
1 + vA s
(12.126)
Note the (expected) appearance of the slow-wave phase speed [cf. Eq. (12.33)] in the left-
hand side. Thus, slow waves interact only with Alfvénic perturbations—when vA cs ,
only with the counterpropagating ones, but at finite β, because the slow waves are slower,
a co-propagating Alfvénic perturbation can catch up with a slow one, have its way with
it in passing and speed on (it’s a tough world).
There is no energy exchange in these interactions: the “+” and “−” slow-wave energies
are individually conserved:
Z
d
d3 r |δZk± |2 = 0. (12.127)
dt
† At high β, vA cs , so we recover from Eqs. (12.125) and (12.122) the Elsasser fields as
defined for iMHD in Eq. (12.71).
116 A. A. Schekochihin
12.3.4. Entropy Mode
There are only two equations in Eq. (12.126), whereas we had three equations (12.119–
12.121) for our three compressive fields δB, uk and δρ. The third equation, Eq. (12.121),
was in fact for the entropy perturbation:
v 2 δB
dδs δs δρ
=0 , = −γ + A . (12.128)
dt s0 ρ0 c2s B0
ω=0 . (12.129)
This is the (famously often forgotten) 7th MHD mode, known as the entropy mode (there
are 7 equations in MHD, so there must be 7 linear modes: two fast waves, two Alfvén
waves, two slow waves and one entropy mode).
Exercise 12.5. Go back to §12.1 and find where we overlooked this mode.
Since the entropy mode is decoupled, its “energy” (variance) is individually conserved:
Z
d
d3 r |δs|2 = 0. (12.130)
dt
Thus, in RMHD, the (nonlinear) evolution
R 3 of ±all2 perturbations is constrained by 5 sepa-
rate conservation laws: d3 r |δZ ± 2
R R 3 2
⊥ | , d r |δZk | and d r |δs| are all invariants.
This simply says that the total pressure gradient is balanced by the tension force. A
general equilibrium for which this is satisfied is called a screw pinch.
118 A. A. Schekochihin
One simple particular case of this is the z pinch (Fig. 47a). This is achieved by letting
a current flow along the z axis, giving rise to an azimuthal field:
4π 1 r 0 0
Z
c 1 ∂
jθ = 0, jz = rBθ ⇒ Bθ = dr r jz (r0 ), Bz = 0. (13.9)
4π r ∂r c r 0
Eq. (13.8) becomes
∂p 1
= − jz Bθ . (13.10)
∂r c
The “pinch” comes from magnetic loops and is due to the curvature force: the loops
want to contract inwards, the pressure gradient opposes this and so plasma is confined
(Fig. 47b). This configuration will, however, prove to be very badly unstable (§14.4)—
which does not stop it from being a popular laboratory set up for short-term confinement
experiments (see, e.g., review by Haines 2011).
Another simple particular case is the θ pinch (Fig. 48a). This is achieved by imposing
a straight but radially non-uniform magnetic field in the z direction and, therefore,
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 119
azimuthal currents:
c ∂Bz
Bθ = 0, jz = 0, jθ = − . (13.11)
4π ∂r
Eq. (13.8) is then just a pressure balance, pure and simple:
B2
∂
p+ z =0 . (13.12)
∂r 8π
In this configuration, we can confine the plasma (Fig. 48c) or the magnetic flux (Fig. 48d).
The latter is what happens, for example, in flux tubes that link sunspots (Fig. 48b). The
θ pinch is a stable configuration (Q9).
The more general case of a screw pinch [Eq. (13.8)] is a superposition of z and θ
pinches, with both magnetic fields and currents wrapping themselves around cylindrical
flux surfaces.
The next step in complexity is to assume axial, but not cylindrical symmetry (∂/∂θ = 0,
∂/∂z 6= 0). This is explored in Q7.
−∇2 B = α∇ × B = α2 B ⇒ ∇2 + α 2 B = 0 ,
(13.15)
so the magnetic field satisfies a Helmholtz equation (to solve which, one must, of course,
specify some boundary conditions).
Thus, there is, potentially, a large zoo of MHD equilibria. Some of them are stable,
some are not, and, therefore, some are more interesting and/or more relevant than others.
120 A. A. Schekochihin
How does one tell? A good question to ask is as follows. Suppose we set up some initial
configuration of magnetic field (by, say, switching on some current-carrying coils, driving
currents inside plasma, etc.)—to what (stable) equilibrium will this system eventually
relax?
In general, any initially arranged magnetic configuration will exert forces on the plasma,
these will drive flows, which in turn will move the magnetic fields around; eventually, ev-
erything will settle into some static equilibrium. We expect that, normally, some amount
of the energy contained in the initial field will be lost in such a relaxation process because
the flows will be dissipating, the fields diffusing and/or reconnecting, etc.—the losses oc-
cur due to the resistive and viscous terms in the non-ideal MHD equations derived in
§11. Thus, one expects that the final relaxed static state will be a minimum-energy state
and so we must be able to find it by minimising magnetic energy:
B2
Z
d3 r → min . (13.16)
8π
Clearly, if the relaxation occured without any constraints, the solution would just be
B = 0. In fact, there are constraints. These constraints are topological: if you think
of magnetic field lines as a tangled mess, you will realise that, while you can change
this tangle by moving field lines around, you cannot easily undo linkages, knots, etc.—
anything that, to be undone, would require the field lines to have “ends”. This intuition
can be turned into a quantitative theory once we discover that the induction equation
(11.59) has an invariant that involves the magnetic field only.
13.2. Helicity
Magnetic helicity in a volume V is defined as
Z
H= d3 r A · B , (13.17)
V
The surface integral vanishes provided both u and B are parallel to the boundary (no
fields stick out and no flows cross). The resistive term in the surface integral can also
be ignored either by arranging V appropriately or simply by taking it large enough so
B → 0 on ∂V , or, indeed, by taking η → +0. Thus,
Z
dH
= −2η d3 r B · (∇ × B) , (13.24)
dt
Thus, as an initial magnetic configuration relaxes, while its energy can change quickly
(on dynamical times), its helicity changes only very slowly in the limit of small η. The
constancy of H (as η → +0) provides us with the constraint subject to which the energy
will need to be minimised.
Before we use this idea, let us discuss what the conservation of helicity means physically,
or, rather, topologically.
† The resistive term in the right-hand side of Eq. (13.24) is ∝ d3 r B · j, a quantity known
R
as the current helicity.
122 A. A. Schekochihin
By the same token, in general, in a system of many linked tubes, the helicity of tube i is
X
Hi = Φi Φthrough hole in tube i = Φi Φj Nij , (13.29)
j
where Nij is the number of times tube j passes through the hole in tube i. The total
helicity of the this entire assemblage of flux tubes is then
X
H= Φi Φj Nij . (13.30)
ij
Thus, H is the number of linkages of the flux tubes weighted by the field strength in
them. It is in this sense that it is a topological invariant.
∇ × B = αB ⇒ ∇2 B = −α2 B . (13.38)
Thus, our system will relax to a linear force-free state determined by Eq. (13.38) and
system-specific boundary conditions. Here α = α(H) depends on the (fixed by initial
conditions) value of H via the equation
Z Z
1
H(α) = d3 r A · B = d3 r B 2 , (13.39)
α
where B is the solution of Eq. (13.38) (since ∇ × B = αB = α∇ × A, we have
B = αA + ∇χ and the χ term vanishes under volume integration).
solution with the smallest energy must be the right one (if a system relaxed to a local
minimum, one can always imagine it being knocked out of it by some perturbation and
falling to a lower energy).
This gives us an interesting twisted field (Fig. 50), able to maintain itself in equilibrium
without help from pressure gradients.
Finally, we calculate its helicity according to Eq. (13.39): assuming that the length of
the cylinder is L, its radius R and so its volume V = πR2 L, we have
2πLB02 R
Z Z
1
d3 r B 2 = drr J02 (αr) + J12 (αr)
H=
α α 0
B02 V
2 2 2 2
= J 0 (αR) + 2J1 (αR) + J 2 (αR) − J 1 (αR)J 2 (αR) . (13.43)
α2 αR
If we solve this for α = α(H), our solution is complete.
Exercise 13.1. Work out what happens in the general case of ∂/∂θ 6= 0 and ∂/∂z 6= 0 and
whether the simple symmetric solution obtained above is the correct relaxed, minimum-energy
state (not always, it turns out). This is not a trivial exercise. The solution is in Taylor & Newton
(2015, §9), where you will also find much more on the subject of J. B. Taylor relaxation, relaxed
states and much besides—all from the original source.
There are other useful variational principles—other in the sense that the constraints that are
imposed are different from helicity conservation. The need for those arises when one consid-
ers magnetic equilibria in domains that do not completely enclose the field lines, i.e., at the
boundary, dS · B 6= 0. One example of such a variational principle, also yielding a force-free
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 125
field (although not necessarily a linear one), is given in Q1(e). A specific example of such a field
arises in Q7(f).
(Bernstein et al. 1958). Before we are in a position to prove this, we must do some
preparatory work.
14.1.1. Properties of the Force Operator F [ξ]
Since the right-hand side of Eq. (14.6) is symmetric with respect to swapping ξ ↔ η,
so must be the left-hand side:
Z Z
d3 r η · F [ξ] = d3 r ξ · F [η]. (14.9)
provided F has no complex coefficients (we shall confirm this explicitly in §14.2.1). Taking
the full scalar products (including integarting over space) of Eq. (14.11) with ξ ∗n and of
Eq. (14.12) with ξ n and subtracting one from the other, we get
Z Z Z
2 ∗
d r ρ0 |ξ n | = d r ξn · F [ξ n ] − d3 r ξn · F [ξ ∗n ] = 0
∗
2 3 2 3
− ωn − (ωn )
| {z }
>0
⇒ ωn2 = (ωn2 )∗ , q.e.d. (14.13)
This result implies that, if any MHD equilibrium is unstable, at least one of the eigenvalues
must be ωn2 < 0 and, since it is guaranteed to be real, any MHD instability will give rise
to purely growing modes (Fig. 51a), rather than growing oscillations (also known as
“overstabilities”; see Fig. 51b).
2) The eigenmodes {ξ n } are orthogonal.
Proof. Taking the full scalar products of Eq. (14.11) with ξ m (assuming m 6= n and
2
non-degeneracy of ωm,n ), and of the analogous equation
2
F [ξ m ] = −ρ0 ωm ξm (14.14)
with ξ n and subtracting them, we get†
Z Z Z
2 2
− (ωn − ωm ) d r ρ0 ξ n · ξ m = d r ξm · F [ξ n ] − d3 r ξn · F [ξ m ] = 0
3 3
| {z }
6= 0
Z Z
⇒ d3 r ρ0 ξ n · ξ m = δnm d3 r ρ0 |ξ n |2 , q.e.d.
(14.15)
Note
R that again δρ, δp and δB are all expressed as linear operators on ξ—and so δW =
δ d3 r B 2 /8π + p/(γ − 1) must also be some operator involving ξ and its gradients
but not ∂ξ/∂t (as we assumed in §14.1).
Finally, we deal with the momentum equation (to which we add gravity as this will
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 129
give some interesting instabilities):
∂u (∇ × B) × B
ρ + u · ∇u = −∇p + + ρg. (14.23)
∂t 4π
This gives us
∂2ξ (∇ × B 0 ) × δB (∇ × δB) × B 0
F [ξ] = ρ0 = −∇δp + + + + δρ g
∂t2 4π 4π
j × δB (∇ × δB) × B 0
= ∇ (ξ · ∇p0 + γp0 ∇ · ξ) − g∇ · (ρ0 ξ) + 0 + , (14.24)
c 4π
where j 0 = c(∇ × B 0 )/4π, we have used Eqs. (14.20) and (14.21) for δρ and δp, respec-
tively, and δB is given by Eq. (14.22).
(j 0 × δB) · ξ (∇ × δB) × B 0
− − ·ξ . (14.25)
| {zc } | 4π
{z }
j 0 · (ξ × δB) (∇ × δB) · (ξ × B 0 )
= =
c 4π
by parts
(δB × ∇) · (ξ × B 0 )
=
4π
δB · [∇ × (ξ × B 0 )]
=
4π
|δB|2
= by Eq. (14.22)
4π
Thus, we have arrived at a standard textbook (e.g., Kulsrud 2005) expression for the en-
ergy perturbation (this expression is non-unique because one can do various integrations
by parts):
Z
1
δW2 = d r (ξ · ∇p0 ) ∇ · ξ + γp0 (∇ · ξ)2 + (g · ξ) ∇ · (ρ0 ξ)
3
2
j 0 · (ξ × δB) |δB|2
+ + , (14.26)
c 4π
where δB = ∇ × (ξ × B 0 ). Note that two of the terms inside the integral (the second
and the fifth) are positive-definite and so always stabilising. The terms that are not sign-
definite and so potentially destabilising involve equilibrium gradients of pressure, density
and magnetic field (currents). It is perhaps not a surprise to learn that Nature might
dislike gradients—while it is of course not a rule that all such inhomogeneities render the
system unstable, we will see that they often do, usually when gradients exceed certain
critical thresholds.
All we need to do now is calculate δW2 according to Eq. (14.26) for any equilibrium
130 A. A. Schekochihin
that interests us and see if it can be negative for any class of perturbations (or show that
it is positive for all perturbations).
dp0
ρ0 = ρ0 (z) and p0 = p0 (z) satisfying = −ρ0 g (14.27)
dz
(gravity acts downward, against the z direction, g = −gẑ).
d ln s0
δW2 > 0 ⇔ >0 , (14.31)
dz
where s0 = p0 /ργ0 is the entropy function. Eq. (14.31) is the Schwarzschild criterion for
convective stability†. If this criterion is broken, there will be an instability, called the
interchange instability.
This calculation illustrates both the power and the weakness of the method:
—on the one hand, we have obtained a stability criterion quite quickly and without
having to solve the underlying equations,
—on the other hand, while we have established the condition for instability, we have
as yet absolutely no idea what is going on physically.
1/γ
ρ1 ρ01 p02 p02 p01
ρ1 < ρ02 ⇔ = <1 ⇔ < γ . (14.37)
ρ02 ρ02 p01 ργ02 ρ01
This is exactly the same as the Schwarzschild condition (14.31) for the interchange in-
stability (and this is why the instability is called that).
Note that, while this is of course a much simpler and more intuitive agrument than
the application of the Energy Principle, it only gives us a particular example of the
kind of perturbation that would be unstable under particular conditions, not any general
criterion of what equilibria might be guaranteed to be stable.
In Q8, we will explore how the above considerations can be generalised to an equilibrium that
also features a non-zero magnetic field.
As our second (also classic) example, we consider the stability of a z-pinch equilibrium
(§13.1.1, Fig. 47):
c 1 B0 (rB0 )0
B 0 = B0 (r)θ̂, j 0 = j0 (r)ẑ = (rB0 )0 ẑ, p00 (r) = − j0 B0 = − . (14.38)
4πr c 4πr
Since we are going to have to work in cylindrical coordinates, we must first write all
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 133
the terms in Eq. (14.26) in these coordinates and with the equilibrium (14.38):
!
0 1 ∂ @1 ∂ξθ ∂ξz
(ξ · ∇p0 )(∇ · ξ) = ξr p0 rξr + @ +
r ∂r r ∂θ
@ ∂z
2
ξ ∂ξr ∂ξz
= p00 r + p00 ξr + , (14.39)
r ∂r ∂z
2
1 ∂ 1 ∂ξθ ∂ξz
γp0 (∇ · ξ)2 = γp0 rξr + + , (14.40)
r ∂r r ∂θ ∂z
1 ∂ ∂ ∂ 1 ∂
δB = ∇ × (ξ × B 0 ) = r̂ ξr B0 + θ̂ − ξz B0 − ξr B0 + ẑ ξz B0 ,
r ∂θ ∂z ∂r r ∂θ
(14.41)
" #
B0
j 0 · (ξ × δB) j0 ∂ξz ∂ξr Z 1 ∂ξ
r
= (ξr δBθ − ξθ δBr ) = p00 ξr + + ξr 0 + ξθ Z ,
c c
|{z} ∂z ∂r B 0 r
∂θ
Z Z
p00
=−
B0
p0 B 0
∂ξz ∂ξr
= p00 ξr + + 0 0 ξr2 , (14.42)
∂z ∂r B0
" 2 2 # 2
2 2
B02 ∂ξz B00
|δB| B0 ∂ξr ∂ξz ∂ξr
= + + + + ξr . (14.43)
4π 4πr2 ∂θ ∂θ 4π ∂z ∂r B0
| {z }
2
B02 ∂ξz
∂ξr
= + +
4π ∂z ∂r
0
2B0 B0 ∂ξz ∂ξr B 02
ξr + + 0 ξr2
4π ∂z ∂r 4π
The terms that are crossed out have been dropped because they combine into a full
derivative with respect to θ and so, upon substitution into Eq. (14.26), vanish under
integration. Assembling all this together, we have
(
p00 rB00 rB002 ξr2 B0 B00
Z
1 3 0 0 ∂ξz ∂ξr
δW2 = d r p0 + + + 2 p0 + ξr +
2 B0 4π r 4π ∂z ∂r
| {z } | {z }
B02 B02
= 2p00 + =−
4πr 4πr
" 2 2 # 2 )
2 2
B ∂ξr ∂ξ z B ∂ξz ∂ξr
+ γp0 (∇ · ξ)2 + 0 2 + + 0 +
4πr ∂θ ∂θ 4π ∂z ∂r
( 2
ξ2 B 2 ∂ξz
Z
1 ∂ξr ξr
= d3 r 2p00 r + 0 + −
2 r 4π ∂z ∂r r
" 2 #)
2
2
B ∂ξr ∂ξ z
+ γp0 (∇ · ξ)2 + 0 2 + , (14.44)
4πr ∂θ ∂θ
where, in simplifying the first two terms in the integrand, we used the equilibrium equa-
tion (14.38):
We shall treat ξr and η as independent variables and minimise δW2 with respect to η:
B02
∂ integrand of ξr ξr 1 − γβ/2 ξr
=2 η− + 2γp0 η + =0 ⇒ η= ,
∂η Eq. (14.48) 4π r r 1 + γβ/2 r
(14.49)
where, as usual, β = 8πp0 /B02 . Putting this back into Eq. (14.48), we get
" 2 2 # 2
0
Z
rp 0 1 γβ γ 2 ξr
δW2 = d3 r p0 + +
p0 β 1 + γβ/2 2 1 + γβ/2 r2
Z 2
d ln p0 2γ ξr
= d3 r p0 r + . (14.50)
dr 1 + γβ/2 r2
There will be an instability (δW2 < 0) if (but not only if, because we are considering the
restricted set of axisymmetric displacements)
d ln p0 2γ
−r > , (14.51)
dr 1 + γβ/2
i.e., when the pressure gradient is too steep, the equilibrium is unstable.
What sort of instability is this? Recall that the perturbations that we have identified
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 135
as making δW2 < 0 are axisymmetric, have some radial and axial displacements and are
compressible: from Eq. (14.49),
ξr 2 ξr
∇·ξ =η+ = . (14.52)
r 1 + γβ/2 r
They are illustrated in Fig. 53. The mechanism of this aptly named sausage instability
is clear: squeezing the flux surfaces inwards increases the curvature of the azimuthal
field lines, this exerts stronger curvature force, leading to further squeezing; conversely,
expanding outwards weakens curvature and the plasma can expand further.
Exercise 14.1. Convince yourself that the displacements that have been identified cause mag-
netic perturbations that are consistent with the cartoon in Fig. 53.
∇ · ξ = 0, (14.53)
i.e., the most dangerous non-axisymmetric perturbations are incompressible (unlike for
the case of the axisymmetric sausage mode in §14.4.1: there we could not—and did not—
have such incompressible perturbations because we did not have ξθ at our disposal, to
be chosen in such a way as to enforce incompressibility).
To carry out further minimisation of δW2 , it is convenient to Fourier transform our
displacements in the θ and z directions—both are directions of symmetry (i.e., the equi-
librium profiles do not vary in these directions), so this can be done with impunity:
X
ξ= ξ mk (r) ei(mθ+kz) . (14.54)
m,k
Then Eq. (14.46) (with ∇ · ξ = 0) becomes, by Parseval’s theorem (the operator F [ξ]
being self-adjoint; see §14.1.1),
( " #)
∞ 2 2
|ξr | B2 m2
Z
1X ∂ ξr
2p00 + 0 + 2 |ξr |2 + |ξz |2
δW2 = 2πLz dr r r + ikξz .
2 0 r 4π ∂r r r
m,k
(14.55)
As ξz and ξz∗ only appear algebraically in Eq. (14.55) (no r derivatives), it is easy to
minimise δW2 with respect to them: setting the derivative of the integrand with respect
to either ξz or ξz∗ to zero, we get
m2 ikr3
∂ ξr ∂ ξr
−ik r + ikξz + ξz = 0 ⇒ ξz = . (14.56)
∂r r r2 m2 + k 2 r2 ∂r r
136 A. A. Schekochihin
Putting this back into Eq. (14.55) and assembling terms, we get
Z ∞ ( 0
m2 |ξr |2
X rp0
δW2 = πLz dr r 2p0 +
0 p0 β r2
m,k
" 2 # )
2
B02 k2 r2 m2 k 2 r2 ∂ ξr
+ 1− 2 + r . (14.57)
4π m + k2 r2 (m2 + k 2 r2 )2 ∂r r
| {z }
m2
=
m2 + k2 r2
The second term here is always stabilising. The most unstable modes will be ones with
k → ∞, for which the stabilising term is as small as possible. The remaining term will
allow δW2 < 0 and, therefore, an instability, if
d ln p0 m2
−r > . (14.58)
dr β
Again, the equilibrium is unstable if the pressure gradient is too steep. The most unstable
modes are ones with the smallest m, viz., m = 1.
What does this instability look like? The unstable perturbations are incompressible:
1 ∂ im
∇·ξ =0 ⇒ rξr + ξθ + ikξz = 0. (14.60)
r ∂r r
Setting m = 1 and using Eq. (14.56), we find
∂ k2 r4 ∂ ξr
iξθ = − rξr + 2 ≈ −2ξr and ξz ξr . (14.61)
∂r m + k 2 r2 ∂r r
| {z }
≈ r2
as k → ∞
The basic cartoon (Fig. 54) is as follows: the flux surfaces are bent, with a twist (to
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 137
remain uncompressed). The bending pushes the magnetic loops closer together and thus
increases magnetic pressure in concave parts and, conversely, decreases it in the convex
ones. Plasma is pushed from the areas of higher B to those with lower B, thermal pressure
in the latter (convex) areas becomes uncompensated, the field lines open up further, etc.
This is called the kink instability.
Similar methodology can be used to show that, unlike the z pinch, the θ pinch (§13.1.1, Fig. 48)
is always stable: see Q9.
(a) Argue that any magnetic field line can be described by the equations
α = const, β = const. (15.2)
This means that (α, β, `), where ` is the distance (arc length) along the field line, are a
good set of curvilinear coordinates, known as the Clebsch coordinates.
(b) Show that the magnetic flux through any area S in the (x, y) plane is
Z
Φ= dα dβ, (15.3)
S̃
where S̃ is the area S in new coordinates after transforming (x, y) → (α(x, y, 0), β(x, y, 0)).
(c) Show that if Eq. (15.1) holds at time t = 0 and α and β are evolved in time
according to
dα dβ
= 0, = 0, (15.4)
dt dt
where d/dt is the convective derivative, then Eq. (15.1) correctly describes the magnetic
field at all t > 0.
(d) Argue from the above that magnetic flux is frozen into the flow and magnetic field
lines move with the flow.
(e) Show that the field that minimises the magnetic energy within some domain subject
to the constraint that the values of α and β are fixed at the boundary of this domain
(i.e., that the “footpoints” of the field lines are fixed) is a force-free field.
A prototypical example of the kind of fields that arise from the variational principle in (e) is
the “arcade” fields describing magnetic loops sticking out of the Sun’s surface, with footpoints
anchored at the surface. One such field will be considered in Q7(f) and more can be found in
Sturrock (1994, §13).
(a) In the neighbourhood of the stagnation point (0, 0, 0), linearise the flow, assume
vertical magnetic field, B = (0, 0, B(t, x)) and derive an evolution equation for B(t, x),
including both advection by the flow and Ohmic diffusion. Suppose the field is initially
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 141
uniform, B(t = 0, x) = B0 = const. It should be clear to you from your equation that
magnetic field is being swept towards x = 0. What is the time scale of this sweeping?
Given the magnetic Reynolds number Rm = U L/η 1, show that flux conservation
holds on this time scale.
(b) Find the steady-state solution of your equation. Assume B(x) = B(−x) and use
flux conservation to determine the constants of integration (in terms of B0 and Rm).
What is the width of the region around x = 0 where the flux is concentrated? What is
the magnitude of the field there?
(c∗ ) Obtain the time-dependent solution of your equation for B and confirm that it
indeed converges to your steady-state solution. Find the time scale on which this happens.
√
Hint. The following changes of variables may prove useful: ξ = πRm x/L, τ = πU t/L,
X = ξeτ , s = (e2τ − 1)/2.
(d) Can you think of a quick heuristic argument based on the induction equation that
would tell you that all these answers were to be expected?
4. Zeldovich’s Antidynamo Theorem. Consider an arbitrary 2D velocity field: u =
(ux , uy , 0). Assume incompressibility. Show that, in a finite system (i.e., in a system that
can be enclosed within some volume outside which there are no fields or flows), this
velocity field cannot be a dynamo, i.e., any initial magnetic field will always eventually
decay.
Hint. Consider separately the evolution equations for Bz and for the magnetic field in
the (x, y)-plane. Show that Bz decays by working out the time evolution of the volume
integral of Bz2 . Then write Bx , By in terms of one scalar function (which must be possible
because ∂Bx /∂x + ∂By /∂y = 0) and show that it decays as well.
5. MHD Waves in a Stratified Atmosphere. The generalisation of iMHD to the case
of a stratified atmosphere is explained in §12.2.9. Convince yourself that you understand
how the SMHD equations and the SMHD ordering arise and then study them as follows.
(a) Work out all SMHD waves (both their frequencies and the corresponding eigenvec-
tors). It is convenient to choose the coordinate system in such a way that k = (kx , 0, kz ),
where z is the vertical direction (the direction of gravity). The mean magnetic field
B 0 = B0 b0 is assumed to be straight and uniform, at a general angle to z. We con-
tinue referring to the projection of the wave number onto the magnetic-field direction
as kk = k · b0 = kx b0x + kz b0z . Note that in the case of B0 = 0, you are dealing with
142 A. A. Schekochihin
stratified hydrodynamics, not MHD—the waves that you obtain in this case are the well
known gravity waves, or “g-modes”.
(b) Explain the physical nature of the perturbations (what makes the fluid oscillate)
in the special cases (i) kz = 0 and b0 = ẑ, (ii) kz = 0 and b0 = x̂, (iii) kx = 0, (iv)
kz 6= 0, kx 6= 0 and b0 = ẑ.
(c) Under what conditions are the perturbations you have found unstable? What is the
physical mechanism for the instability? What role does the magnetic field play (stabilising
or destabilising) and why? Cross-check your answers with §14.3 and Q8.
(d) Find the conserved energy (a quadratic quantity whose integral over space stays
constant) for the full nonlinear SMHD equations (12.93–12.96). Give a physical interpre-
tation of the quantity that you have obtained—why should it be conserved?
6. Electron MHD. In certain physical regimes (roughly realised, for example, in the
solar-wind and other kinds of astrophysical turbulence at scales smaller than the ion
Larmor radius; see Schekochihin et al. 2009, Boldyrev et al. 2013), plasma turbulence
can be described by an approximation in which the magnetic field is frozen into the
electron flow ue , while ions are considered motionless, ui = 0. In this approximation,
Ohm’s law becomes†
ue × B
E=− . (15.6)
c
Here ue can be expressed directly in terms of B because the current density in a plasma
consisting of motionless hydrogen ions (ni = ne ) and moving electrons is
j = ene (ui − ue ) = −ene ue , (15.7)
but, on the other hand, j is known via Ampère’s law. Here ne is the electron number
density and e the electron charge.
(a) Using this and Faraday’s law, show that the evolution equation for the magnetic
field in this approximation is
∂B
= −di ∇ × [(∇ × B) × B] , (15.8)
∂t
√
p B/ 4πmi ni → B,
where the magnetic field has been rescaled to Alfvénic velocity units,
and di = c/ωpi is the ion inertial scale (“ion skin depth”), ωpi = 4πe2 ni /mi . Eq. (15.8)
is the equation of Electron MHD (EMHD), completely self-consistent for B.
(b) Show that magnetic energy is conserved by Eq. (15.8). Is magnetic helicity con-
served? Does J. B. Taylor relaxation work and what kind of field will be featured in the
relaxed state? Is it obvious that this field is a good steady-state solution of Eq. (15.8)?
(c) Consider infinitesimal perturbations of a straight-field equilibrium, B = B0 ẑ + δB,
and show that they are helical waves with the dispersion relation
ω = ±kk vA kdi . (15.9)
Figure 56. A simple equilibrium from Q7, superficially resembling the poloidal cross-section
of a tokamak.
(f) Seek solutions to Eq. (15.15) that are linear force-free fields. Show that in this
case, Eq. (15.15) reduces to the Bessel equation (a substitution ψ = rf (r, z) will prove
useful). Set Bz (0, 0) = B0 . Find solutions of two kinds: (i) ones in a semi-infinite domain
z > 0, with the field vanishing exponentially at z → ∞; (ii) ones periodic in z. If you
also impose the boundary condition Br = 0 at r = R, how can this be achieved? Can
either of these solutions be the result of J. B. Taylor relaxation of an MHD system? If
so, how would one decide whether it is more or less likely to be the correct relaxed state
than the solution derived in §13.4?
You will find the solution of the type (i) in Sturrock (1994, §13) (who also shows how to construct
many other force-free fields, useful in various physical and astrophysical contexts). Think of this
solution in the context of Q1(e). The solution of type (ii) is a particular case of the general
(∂/∂θ 6= 0) equilibrium solution derived and discussed in Taylor & Newton (2015, §9). However,
the axisymmteric solution is not very useful because, as they show, depending on the values of
helicity and of R, the true relaxed state is either the cylidrically and axially symmetric solution
derived in §13.4 or one which also has variation in the θ direction.
If field-line bending is allowed (∂ξ/∂x 6= 0), another instability emerges, the Parker (1966)
instability. Do investigate.
9. Stability of the θ Pinch. Consider the following cylindrically and axially symmetric
equilibrium:
B02
c 0 d
B 0 = B0 (r)ẑ, j 0 = j0 (r)θ̂ = − B (r)θ̂, p0 + =0 (15.19)
4π 0 dr 8π
(a θ pinch; see §13.1.1, Fig. 48). Consider general displacements of the form
ξ = ξ mk (r)eimθ+ikz . (15.20)
Show that the θ pinch is always stable. Specifically, you should be able to show that
Z ∞ ( "
2
#)
2 B02 2 2 2
ξr ∂ξr imξθ
δW2 = 2πLz dr r γp0 |∇ · ξ| + k |ξr | + |ξθ | + + + > 0,
0 4π r ∂r r
(15.21)
where Lz is the length of the cylinder.
146 A. A. Schekochihin
Acknowledgments
I am grateful to many at Oxford who have given advice, commented, pointed out er-
rors, asked challenging questions and generally helped me get a grip. Particular thanks
to Andre Lukas and Fabian Essler, who were comrades in arms in the creation of the
MMathPhys programme, to Paul Dellar and James Binney, with whom I have collabo-
rated on teaching Kinetic Theory and MHD, and to Felix Parra and Peter Norreys, who
have helped design a coherent sequence of courses on plasma physics. I am also grateful to
M. Kunz and to the students who took the course, amongst them R. Cooper, J. Graham,
C. Hamilton, M. Hawes and M. Hardman, for critical comments and spotting lapses in
these notes. More broadly, I am in debt to the students who have soldiered on through
these lectures and, sometimes verbally, sometimes not, projected back at me some sense
of whether I was succeeding in communicating what I wanted to communicate. Section 5
appeared as a result of G. Lesur’s invitation to lecture on plasma stability at the 2017
Les Houches School on Plasma Physics—forcing me to face the challenge of teaching
something fundamental, beautiful, but not readily available in every plasma course. My
exposition was greatly influenced by the Princeton lectures by Russell Kulsrud (on MHD;
now published: Kulsrud 2005) and the UCLA ones by Steve Cowley (Cowley 2004), as
well as by the excellent books of Kadomtsev (1965), Krall & Trivelpiece (1973), Alexan-
drov et al. (1984) and Kingsep (2004).
REFERENCES
Abel, I. G., Plunk, G. G., Wang, E., Barnes, M., Cowley, S. C., Dorland, W. &
Schekochihin, A. A. 2013 Multiscale gyrokinetics for rotating tokamak plasmas: fluctu-
ations, transport and energy flows. Rep. Prog. Phys. 76, 116201.
Adams, D. 1979 The Hitch-Hiker’s Guide to the Galaxy. London: Pan Books.
Adkins, T. & Schekochihin, A. A. 2017 A solvable model of Vlasov-kinetic turbulence in a
two-dimensional phase space. In preparation.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 147
Alexandrov, A. F., Bogdankevich, L. S. & Rukhadze, A. A. 1984 Principles of Plasma
Electrodynamics. Berlin: Springer.
Arnold, V. I. & Khesin, B. A. 1999 Topological Methods in Hydrodynamics. Berlin: Springer.
Balbus, S. A. 2015 Astrophysical Fluid Dynamics. Lecture Notes for the Oxford MMathPhys
course; URL: https://www2.physics.ox.ac.uk/contacts/people/balbus.
Balbus, S. A. & Hawley, J. F. 1991 A powerful local shear instability in weakly magnetized
disks. I. Linear analysis. Astrophys. J. 376, 214.
Balescu, R. 1963 Statistical Mechanics of Charged Particles. New York: Interscience Press.
Bedrossian, J. 2016 Nonlinear echoes and Landau damping with insufficient regularity.
arXiv:1605.06841.
Beresnyak, A. & Lazarian, A. 2015 MHD turbulence, turbulent dynamo and applications.
In Magnetic Fields in Diffuse Media (ed. A. Lazarian, E. M. de Gouveia Dal Pino &
C. Melioli), Astrophysics and Space Science Library, vol. 407, p. 163. Berlin: Springer,
e-print arXiv:1406.1185.
Bernstein, I. B. 1958 Waves in a plasma in a magnetic field. Phys. Rev. 109, 10.
Bernstein, I. B., Frieman, E. A., Kruskal, M. D. & Kulsrud, R. M. 1958 An energy
principle for hydromagnetic stability problems. Proc. R. Soc. London A 244, 17.
Besse, N., Elskens, Y., Escande, D. F. & Bertrand, P. 2011 Validity of quasilinear theory:
refutations and new numerical confirmation. Plasma Phys. Control. Fusion 53, 025012.
Binney, J. J. 2016 Kinetic Theory of Self-Gravitating Systems. Lecture Notes
for the Oxford MMathPhys course on Kinetic Theory; URL: http://www-
thphys.physics.ox.ac.uk/user/JamesBinney/kt.pdf.
Bohm, D. & Gross, E. P. 1949a Theory of plasma oscillations. A. Origin of medium-like
behavior. Phys. Rev. 75, 1851.
Bohm, D. & Gross, E. P. 1949b Theory of plasma oscillations. B. Excitation and damping of
oscillations. Phys. Rev. 75, 1864.
Boldyrev, S., Horaites, K., Xia, Q. & Perez, J. C. 2013 Toward a theory of astrophysical
plasma turbulence at subproton scales. Astrophys. J. 777, 41.
Braginskii, S. I. 1965 Transport processes in a plasma. Rev. Plasma Phys. 1, 205.
Buneman, O. 1958 Instability, turbulence, and conductivity in current-carrying plasma. Phys.
Rev. Lett. 1, 8.
Case, K. M. 1959 Plasma oscillations. Ann. Phys. (N.Y.) 7, 349.
Chandrasekhar, S. 2003 Hydrodynamic and Hydromagnetic Stability. Dover.
Chapman, S. & Cowling, T. G. 1991 The Mathematical Theory of Non-Uniform Gases.
Cambridge: Cambridge University Press, 3rd Edition.
Chew, G. F., Goldberger, M. L. & Low, F. E. 1956 The Boltzmann equation and the one-
fluid hydromagnetic equations in the absence of particle collisions. Proc. R. Soc. London A
236, 112.
Connor, J. W., Hastie, R. J. & Taylor, J. B. 1979 High mode number stability of an
axisymmetric toroidal plasma. Proc. R. Soc. London A 365, 1.
Cowley, S. C. 2004 Plasma Physics. Lecture Notes for the UCLA course; URL: http://www-
thphys.physics.ox.ac.uk/research/plasma/CowleyLectures/.
Cowley, S. C. & Artun, M. 1997 Explosive instabilities and detonation in magnetohydrody-
namics. Phys. Rep. 283, 185.
Cowley, S. C., Cowley, B., Henneberg, S. A. & Wilson, H. R. 2015 Explosive instability
and erupting flux tubes in a magnetized plasma. Proc. R. Soc. London A 471, 20140913.
Davidson, P. A. 2004 Turbulence: An Introduction for Scientists and Engineers. Oxford: Oxford
University Press.
Davidson, P. A. 2016 Introduction to Magnetohydrodynamics (2nd edition). Cambridge: Cam-
bridge University Press.
Davidson, R. C. 1983 Kinetic waves and instabilities in a uniform plasma. In Handbook of
Plasma Physics, Volume 1 (ed. A. A. Galeev & R. N. Sudan), p. 519. Amsterdam: North-
Holland.
Dellar, P. J. 2016 Kinetic Theory of Gases. Lecture Notes for the Oxford MMathPhys course
on Kinetic Theory; URL: http://people.maths.ox.ac.uk/dellar/MMPkinetic.html.
Dellar, P. J. 2017 Hydrodynamics of Complex Fluids. Lecture Notes for
148 A. A. Schekochihin
the Oxford MMathPhys course on Advanced Fluid Dynamics; URL:
http://people.maths.ox.ac.uk/dellar/MMP Non Newt fluids.html.
Drummond, W. E. & Pines, D. 1962 Nonlinear stability of plasma oscillations. Nucl. Fusion
Suppl., Part 3 p. 1049.
Elsasser, W. M. 1950 The hydromagnetic equations. Phys. Rev. 79, 183.
Eyink, G. & Frisch, U. 2010 Robert H. Kraichnan. arXiv:1011.2383.
Fowler, T. K. 1963 Lyapunov’s stability criteria for plasmas. J. Math. Phys. 4, 559.
Fowler, T. K. 1964 Bounds on plasma instability growth rates. Phys. Fluids 7, 249.
Fowler, T. K. 1968 Thermodynamics of unstable plasmas. In Advances in Plasma Physics
(ed. A. Simon & W. B. Thompson), , vol. 1, p. 201. New York: Interscience Press.
Fowler, T. K. 2016 Speculations about plasma free energy, 50 years later. J. Plasma Phys.
82, 925820501.
Freiderg, J. 2014 Ideal MHD. Cambridge: Cambridge University Press.
Fried, B. D. & Conte, S. D. 1961 The Plasma Dispersion Function. New York: Academic
Press.
Frisch, U. 1995 Turbulence: The Legacy of A. N. Kolmogorov . Cambridge: Cambridge Univer-
sity Press.
Furth, H. P., Killeen, J. & Rosenbluth, M. N. 1963 Finite-resistivity instabilities of a
sheet pinch. Phys. Fluids 6, 459.
Gardner, C. S. 1963 Bound on the energy available from a plasma. Phys. Fluids 6, 839.
Goedbloed, H. & Poedts, S. 2004 Principles of Magnetohydrodynamics. Cambridge: Cam-
bridge University Press.
Goldman, M. V. 1984 Strong turbulence of plasma waves. Rev. Mod. Phys. 56, 709.
Gould, R. W., O’Neil, T. M. & Malmberg, J. H. 1967 Plasma wave echo. Phys. Rev. Lett.
19, 219.
Haeff, A. 1949 The electron-wave tube—a novel method of generation and amplification of
microwave energy. Proc. IRE 37, 4.
Haines, M. G. 2011 A review of the dense Z-pinch. Plasma Phys. Control. Fusion 53, 093001.
Helander, P. 2017 Available energy and ground states of collisionless plas-
mas. J. Plasma Phys., submitted; preprint available at http://www-
thphys.physics.ox.ac.uk/research/plasma/JPP/frontiers17.html .
Helander, P. & Sigmar, D. J. 2005 Collisional Transport in Magnetized Plasmas. Cambridge:
Cambridge University Press.
Howes, G. G., Cowley, S. C., Dorland, W., Hammett, G. W., Quataert, E. &
Schekochihin, A. A. 2006 Astrophysical gyrokinetics: basic equations and linear theory.
Astrophys. J. 651, 590.
Jackson, J. D. 1960 Longitudinal plasma oscillations. J. Nuclear Energy C 1, 171.
Kadomtsev, B. B. 1965 Plasma Turbulence. Academic Press.
Kadomtsev, B. B. & Pogutse, O. P. 1974 Nonlinear helical perturbations of a plasma in the
tokamak. Sov. Phys.–JETP 38, 575.
Kanekar, A., Schekochihin, A. A., Dorland, W. & Loureiro, N. F. 2015 Fluctuation-
dissipation relations for a plasma-kinetic Langevin equation. J. Plasma Phys. 81,
305810104.
Kingsep, A. S. 2004 Introduction to Nonlinear Plasma Physics. Moscow: MZ Press (in Russian).
Kingsep, A. S., Rudakov, L. I. & Sudan, R. N. 1973 Spectra of strong Langmuir turbulence.
Phys. Rev. Lett. 31, 1482.
Klimontovich, Yu. L. 1967 The Statistical Theory of Non-Equilibrium Processes in a Plasma.
MIT Press.
Kolmogorov, A. N. 1941 The local structure of turbulence in incompressible viscous fluid at
very large Reynolds numbers. Dokl. Acad. Nauk SSSR 30, 299.
Krall, N. A. & Trivelpiece, A. W. 1973 Principles of Plasma Physics. New York: McGraw-
Hill.
Krommes, J. A. 2015 A tutorial introduction to the statistical theory of turbulent plasmas, a
half-century after Kadomtsev’s Plasma Turbulence and the resonance-broadening theory of
Dupree and Weinstock. J. Plasma Phys. 81, 205810601.
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 149
Kruskal, M. D. & Oberman, C. R. 1958 On the stability of plasma in static equilibrium.
Phys. Fluids 1, 275.
Kulsrud, R. M. 1983 MHD description of plasma. In Handbook of Plasma Physics, Volume 1
(ed. A. A. Galeev & R. N. Sudan), p. 1. Amsterdam: North-Holland.
Kulsrud, R. M. 2005 Plasma Physics for Astrophysics. Princeton: Princeton University Press.
Kunz, M. W., Schekochihin, A. A., Chen, C. H. K., Abel, I. G. & Cowley, S. C. 2015
Inertial-range kinetic turbulence in pressure-anisotropic astrophysical plasmas. J. Plasma
Phys. 81, 325810501.
Landau, L. 1936 Transport equation in the case of Coulomb interaction. Zh. Eksp. Teor. Fiz.
7, 203.
Landau, L. 1946 On the vibration of the electronic plasma. Zh. Eksp. Teor. Fiz. 16, 574, English
translation: J. Phys. USSR 10, 25 (1946).
Landau, L. D. & Lifshitz, E. M. 1987 Fluid Mechanics (L. D. Landau and E. M. Lifshitz’s
Course of Theoretical Physics, Volume 6). Oxford: Pergamon Press.
Lifshitz, E. M. & Pitaevskii, L. P. 1981 Physical Kinetics (L. D. Landau and E. M. Lifshitz’s
Course of Theoretical Physics, Volume 10). Elsevier.
Loureiro, N. F., Schekochihin, A. A. & Cowley, S. C. 2007 Instability of current sheets
and formation of plasmoid chains. Phys. Plasmas 14, 100703.
Mallet, A. & Schekochihin, A. A. 2017 A statistical model of three-dimensional anisotropy
and intermittency in strong Alfvénic turbulence. Mon. Not. R. Astron. Soc. 466, 3918.
Mallet, A., Schekochihin, A. A. & Chandran, B. D. G. 2016 Disruption of sheetlike
structures in Alfvénic turbulence by magnetic reconnection. arXiv:1612.07604.
Malmberg, J. H. & Wharton, C. B. 1964 Collisionless damping of electrostatic plasma waves.
Phys. Rev. Lett. 13, 184.
Malmberg, J. H., Wharton, C. B., Gould, R. W. & O’Neil, T. M. 1968 Plasma wave
echo experiment. Phys. Rev. Lett. 20, 95.
Mazitov, R. K. 1965 Damping of plasma waves. J. Appl. Mech. Tech. Phys. 6, 22.
Moffatt, H. K. 1978 Magnetic Field Generation in Electrically Conducting Fluids. Cambridge:
Cambridge University Press.
Musher, S. L., Rubenchik, A. M. & Zakharov, V. E. 1995 Weak Langmuir turbulence.
Phys. Rep. 252, 177.
Nazarenko, S. 2011 Wave Turbulence. Berlin: Springer.
Newcomb, W. A. 1962 Lagrangian and Hamiltonian methods in magnetohydrodynamics. Nu-
clear Fusion: Supplement, Part 2 p. 451.
Nyquist, H. 1932 Regeneration theory. Bell Syst. Tech. J. 11, 126.
Ogilvie, G. I. 2016 Astrophysical fluid dynamics. J. Plasma Phys. 82, 205820301.
Ogilvie, G. I. & Proctor, M. R. E. 2003 On the relation between viscoelastic and magneto-
hydrodynamic flows and their instabilities. J. Fluid Mech. 476, 389.
O’Neil, T. 1965 Collisionless damping of nonlinear plasma oscillations. Phys. Fluids 8, 2255.
Parker, E. N. 1966 The dynamical state of the interstellar gas and field. Astrophys. J. 145,
811.
Parra, F. I. 2017a Collisional Plasma Physics. Lecture Notes for an Oxford MMathPhys
course; URL: http://www-thphys.physics.ox.ac.uk/people/FelixParra/
CollisionalPlasmaPhysics/CollisionalPlasmaPhysics.html.
Parra, F. I. 2017b Collisionless Plasma Physics. Lecture Notes for an Oxford MMathPhys
course; URL: http://www-thphys.physics.ox.ac.uk/people/FelixParra/
CollisionlessPlasmaPhysics/CollisionlessPlasmaPhysics.html.
Pelletier, G. 1980 Langmuir turbulence as a critical phenomenon. Part 1. Destruction of the
statistical equilibrium of an interacting-modes ensemble. J. Plasma Phys. 24, 287.
Penrose, O. 1960 Electrostatic instabilities of a uniform non-Maxwellian plasma. Phys. Fluids
3, 258.
Pfirsch, D. & Sudan, R. N. 1993 Nonlinear ideal magnetohydrodynamics instabilities. Phys.
Fluids B 5, 2052.
Pierce, J. R. & Hebenstreit, W. B. 1949 A new type of high-frequency amplifier. Bell Syst.
Tech. J. 28, 33.
Robinson, P. A. 1997 Nonlinear wave collapse and strong turbulence. Rev. Mod. Phys. 69, 507.
150 A. A. Schekochihin
Rudakov, L. I. & Tsytovich, V. N. 1978 Strong Langmuir turbulence. Phys. Rep. 40, 1.
Ruyer, C., Gremillet, L., Debayle, A. & Bonnaud, G. 2015 Nonlinear dynamics of the
ion Weibel-filamentation instability: An analytical model for the evolution of the plasma
and spectral properties. Phys. Plasmas 22, 032102.
Sagdeev, R. Z. & Galeev, A. A. 1969 Nonlinear Plasma Theory. W. A. Benjamin.
Schekochihin, A. A. 2016 Kinetic Theory and Statistical Physics. Lec-
ture Notes for the Oxford Physics course, Paper A1; URL: http://www-
thphys.physics.ox.ac.uk/people/AlexanderSchekochihin/A1/2014/A1LectureNotes.pdf.
Schekochihin, A. A. & Cowley, S. C. 2007 Turbulence and magnetic fields in astrophysi-
cal plasmas. In Magnetohydrodynamics: Historical Evolution and Trends (ed. S. Molokov,
R. Moreau & H. K. Moffatt), p. 85. Berlin: Springer, reprint: URL: http://www-
thphys.physics.ox.ac.uk/people/AlexanderSchekochihin/mhdbook.pdf.
Schekochihin, A. A., Cowley, S. C., Dorland, W., Hammett, G. W., Howes, G. G.,
Quataert, E. & Tatsuno, T. 2009 Astrophysical gyrokinetics: kinetic and fluid turbulent
cascades in magnetized weakly collisional plasmas. Astrophys. J. Suppl. 182, 310.
Schekochihin, A. A., Cowley, S. C., Rincon, F. & Rosin, M. S. 2010 Magnetofluid dy-
namics of magnetized cosmic plasma: firehose and gyrothermal instabilities. Mon. Not. R.
Astron. Soc. 405, 291.
Schekochihin, A. A., Parker, J. T., Highcock, E. G., Dellar, P. J., Dorland, W. &
Hammett, G. W. 2016 Phase mixing versus nonlinear advection in drift-kinetic plasma
turbulence. J. Plasma Phys. 82, 905820212.
Squire, J., Schekochihin, A. A. & Quataert, E. 2017 Amplitude limits and nonlinear
damping of shear-Alfvén waves in high-beta low-collisionality plasmas. New J. Phys., in
press [arXiv:1701.03175].
Strauss, H. R. 1976 Nonlinear, three-dimensional magnetohydrodynamics of noncircular toka-
maks. Phys. Fluids 19, 134.
Sturrock, P. A. 1994 Plasma Physics. An Introduction to the Theory of Astrophysical, Geo-
physical and Laboratory Plasmas. Cambridge: Cambridge University Press.
Taylor, J. B. & Newton, S. L. 2015 Special topics in plasma confinement. J. Plasma Phys.
81, 205810501.
Thornhill, S. G. & ter Haar, D. 1978 Langmuir turbulence and modulational instability.
Phys. Rep. 43, 43.
Tobias, S. M., Cattaneo, F. & Boldyrev, S. 2012 MHD dynamos and turbulence. In Ten
Chapters in Turbulence (ed. P. A. Davidson, Y. Kaneda & K. R. Sreenivasan), p. 351.
Cambridge University Press, e-print arXiv:1103.3138.
Tonks, L. & Langmuir, I. 1929 Oscillations in ionized gases. Phys. Rev. 33, 195.
Tsytovich, V. N. 1995 Lectures on Non-linear Plasma Kinetics. Berlin: Springer.
Uzdensky, D. A., Loureiro, N. F. & Schekochihin, A. A. 2010 Fast magnetic reconnection
in the plasmoid-dominated regime. Phys. Rev. Lett. 105, 235002.
van Kampen, N. G. 1955 On the theory of stationary waves in plasmas. Physica 9, 949.
van Kampen, N. G. & Felderhof, B. U. 1967 Theoretical Methods in Plasma Physics. Am-
sterdam: North-Holland.
Vedenov, A. A. 1963 Quasilinear plasma theory (theory of a weakly turbulent plasma). Plasma
Phys. (J. Nucl. Energy C) 5, 169.
Vedenov, A. A., Velikhov, E. P. & Sagdeev, R. Z. 1962 Quasilinear theory of plasma
oscillations. Nucl. Fusion Suppl., Part 2 p. 465 (in Russian).
Villani, C. 2014 Particle systems and nonlinear Landau damping. Phys. Plasmas 21, 030901.
Weibel, E. S 1958 Spontaneously growing transverse waves in a plasma due to an anisotropic
velocity distribution. Phys. Rev. Lett. 2, 83.
Wicks, R. T., Roberts, D. A., Mallet, A., Schekochihin, A. A., Horbury, T. S. &
Chen, C. H. K. 2013 Correlations at large scales and the onset of turbulence in the fast
solar wind. Astrophys. J. 778, 177.
Zakharov, V. E. 1972 Collapse of Langmuir waves. Sov. Phys.–JETP 35, 908.
Zakharov, V. E., Lvov, V. S. & Falkovich, G. 1992 Kolmogorov Spectra of
Turbulence I: Wave Turbulence. Berlin: Springer (updated online version URL:
https://www.weizmann.ac.il/chemphys/lvov/research/springer-series-nonlinear-
dynamics).
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 151
Zeldovich, Ya. B. 1956 The magnetic field in the two-dimensional motion of a conducting
turbulent liquid. Zh. Eksp. Teor. Fiz. 31, 154, English translation: Sov. Phys. JETP 4, 460
(1957).
152 A. A. Schekochihin
CONTENTS
1. Kinetic Description of a Plasma 1
1.1. Quasineutrality 1
1.2. Weak Interactions 1
1.3. Debye Shielding 2
1.4. Micro- and Macroscopic Fields 3
1.5. Maxwell’s Equations 4
1.6. Vlasov–Landau Equation 5
1.7. Collision Operator 6
1.8. Klimontovich’s Version of BBGKY 7
1.9. So What’s New and What Now? 8
2. Equilibrium and Fluctuations 10
2.1. Plasma Frequency 10
2.2. Slow vs. Fast 11
2.3. Multiscale Dynamics 11
2.4. Hierarchy of Approximations 13
3. Linear Theory: Plasma Waves, Landau Damping and Kinetic Instabli-
ties 14
3.1. Initial-Value Problem 14
3.2. Calculating the Dielectric Function: the “Landau Prescription” 17
3.3. Solving the Dispersion Relation: Slow-Damping/Growth Limit 20
3.4. Langmuir Waves 20
3.5. Landau Damping and Kinetic Instabilities 22
3.6. Physical Picture of Landau Damping 23
3.7. Hot and Cold Beams 25
3.8. Ion-Acoustic Waves 27
3.9. Damping of Ion-Acoustic Waves and Ion-Acoustic Instability 29
3.10. Ion Langmuir Waves 30
3.11. Summary of Electrostatic (Longitudinal) Plasma Waves 30
3.12. Plasma Dispersion Function: Putting Linear Theory on Industrial Basis 31
4. Energy, Entropy, Free Energy, Heating, Irreversibility and Phase Mix-
ing 33
4.1. Energy Conservation and Heating 34
4.2. Entropy and Free Energy 35
4.3. Structure of Perturbed Distribution Function 38
4.4. Landau Damping is Phase Mixing 39
4.5. Role of Collisions 40
4.6. Further Analysis of δf : the Case–van Kampen Mode 42
4.7. Free-Energy Conservation for a Landau-Damped Solution 43
5. General Kinetic Stability Theory 44
5.1. Linear Stability: Nyquist’s Method 44
5.2. Nonlinear Stability: Thermodynamic Method 52
6. Nonlinear Theory: Two Pretty Nuggets 57
6.1. Nonlinear Landau Damping 58
6.2. Plasma Echo 58
7. Quasilinear Theory 58
7.1. General Scheme of QLT 58
7.2. Conservation Laws 60
7.3. Quasilinear Equations for the Bump-on-Tail Instability in 1D 61
Oxford MMathPhys Lecture Notes: Plasma Kinetics and MHD 153
7.4. Resonant Region: QL Plateau and Spectrum 62
7.5. Energy of Resonant Particles 64
7.6. Heating of Non-Resonant Particles 65
7.7. Momentum Conservation 66
8. Kinetics of Quasiparticles 67
8.1. QLT in the Language of Quasiparticles 68
8.2. Weak Turbulence in the Language of Quasiparticles 71
8.3. General Scheme for Calculating Probabilities in WT 71
9. Langmuir Turbulence 71
9.1. Zakharov’s Equations 71
9.2. Secondary Instability of a Langmuir Wave 71
9.3. Weak Langmuir Turbulence 71
9.4. Solitons 71
9.5. Langmuir Collapse 71
9.6. Kingsep–Rudakov–Sudan Turbulence 71
9.7. Pelletier’s Equilibrium Ensemble 71
9.8. Theories Galore 71
10.Stochastic Echo and Phase-Space Turbulence 71
11.MHD Equations 82
11.1. Conservation of Mass 82
11.2. Conservation of Momentum 83
11.3. Electromagnetic Fields and Forces 84
11.4. Maxwell Stress and Magnetic Forces 85
11.5. Evolution of Magnetic Field 86
11.6. Magnetic Reynolds Number 86
11.7. Lundquist Theorem 87
11.8. Flux Freezing 88
11.9. Amplification of Magnetic Field by Fluid Flow 90
11.10.MHD Dynamo 91
11.11.Conservation of Energy 92
12.MHD in a Straight Magnetic Field 95
12.1. MHD Waves 95
12.2. Incompressible MHD 104
12.3. Reduced MHD 111
12.4. MHD Turbulence 116
13.MHD Relaxation 116
13.1. Static MHD Equilibria 117
13.2. Helicity 120
13.3. J. B. Taylor Relaxation 122
13.4. Relaxed Force-Free State of a Cylindrical Pinch 124
13.5. Parker’s Problem and Topological MHD 125
14.MHD Stability and Instabilities 125
14.1. Energy Principle 125
14.2. Explicit Calculation of δW2 128
14.3. Interchange Instabilities 130
14.4. Instabilities of a Pinch 132
15.Further Reading 137
15.1. MHD Instabilities 137
15.2. Resistive MHD 137
15.3. Dynamo Theory and MHD Turbulence 138
154 A. A. Schekochihin
15.4. Lagrangian MHD and MHD Action Principle 138
15.5. Hall MHD, Electron MHD, Braginskii MHD 138
15.6. Double-Adiabatic MHD and Onwards to Kinetics 139
Appendix A 146
A.1. Dimensional Theory of the Kolmogorov Cascade 146
A.2. Exact Laws 146
A.3. Intermittency 146
A.4. Turbulent Mixing 146