Plasma Derivations

Alex Zylstra

February 15, 2014


1 General Plasma
 1.1 Definition of a plasma
 1.2 Debye Length
 1.3 Plasma Frequency
2 Single Particle Motion
 2.1 Cyclotron Motion
 2.2 E × B Drift
 2.3 F × B Drift
 2.4 B Drift
 2.5 Curvature/Gravitational Drift
 2.6 Polarization Drift
 2.7 Adiabatic Invariants
  2.7.1 μ
 2.8 Magnetic Mirror Machine
 2.9 Runaway Electrons
 3.1 2-fluid derivation
  3.1.1 Momentum Equation
  3.1.2 Continuity Equation
  3.1.3 Energy Equation
  3.1.4 Maxwell’s Equations
 3.2 Single-fluid MHD
 3.3 Pinches (MHD Equilibria)
  3.3.1 𝜃 Pinch
  3.3.2 Z Pinch
  3.3.3 Screw Pinch
 3.4 Diamagnetic Drift
 3.5 Drift Waves
 3.6 Two-stream Instability
 3.7 Flux Freezing
 3.8 MHD Instabilities
  3.8.1 Classification
  3.8.2 Qualitative Descriptions
4 Waves
 4.1 Gas Dynamics Sound Waves
 4.2 Electromagnetic Waves
 4.3 Polarization
 4.4 Electrostatic Plasma Oscillations
 4.5 Electron Plasma Waves
 4.6 Ion Acoustic Waves
 4.7 Electromagnetic Waves in Plasmas
 4.8 Electromagnetic Ion Waves
  4.8.1 Alfvén Wave
  4.8.2 Magnetosonic Wave
 4.9 Brief Summary of More Waves
  4.9.1 Upper Hybrid Waves
  4.9.2 Ion Cyclotron Waves
  4.9.3 Lower Hybrid Waves
  4.9.4 O/X Waves
  4.9.5 R/L Waves
5 Transport
 5.1 Diffusion in weakly-ionized gases
 5.2 Ambipolar Diffusion
 5.3 Diffusion in a Slab
 5.4 Diffusion in Magnetic Fields
 5.5 Diffusion in Fully Ionized Plasmas
 5.6 Plasma Resistivity
 5.7 Coulomb Collisions
 5.8 Coulomb Collision Frequencies
 5.9 Distribution Function
 5.10 Overview of Kinetic Theories
 5.11 Landau Damping
6 Diagnostics
 6.1 Refractive measurements of density
 6.2 Refractive Measurements of Astrophysical Length
 6.3 Faraday Rotation
 6.4 Magnetic Field Probes
  6.4.1 ˙B Probes
  6.4.2 Hall Probes
  6.4.3 Rogowski Coils
 6.5 Langmuir Probes
7 Tokamaks
 7.1 Introduction
 7.2 Lawson Criterion
 7.3 MHD Equilibrium in Tokamaks
  7.3.1 Screw Pinch Radial Balance
  7.3.2 Toroidal Force Balance - Hoop
  7.3.3 Toroidal Force Balance - Tire Tube
  7.3.4 Toroidal Force Balance - 1/R
  7.3.5 Toroidal Force Balance - Wall Stabilization
  7.3.6 Toroidal Force Balance - Vertical B Stabilization
 7.4 Heating
  7.4.1 Ohmic Heating
  7.4.2 Neutral Beams
  7.4.3 RF Heating
  7.4.4 Lower Hybrid Current Drive
 7.5 Banana Orbits
 8.1 Transport
  8.1.1 Bremsstrahlung Radiation
  8.1.2 Electron Heat Conduction
 8.2 Laser-Plasma Interactions
  8.2.1 Critical Density
  8.2.2 Inverse Bremsstrahlung
  8.2.3 Resonance Absorption
  8.2.4 Parametric Instabilities
  8.2.5 Hot Electrons
 8.3 Misc
  8.3.1 EM Field Generation
  8.3.2 Rayleigh-Taylor Instability
9 Back of the Envelope
 9.1 Plasma Frequency
  9.1.1 Electron
  9.1.2 Ion
 9.2 Debye Length
 9.3 Cyclotron Motion
  9.3.1 Gyroradii
  9.3.2 Cyclotron Frequencies
 9.4 Collision Frequencies
 9.5 Thermal Velocities
 9.6 Mean Free Path
 9.7 Ion Sound Speed
 9.8 E × B drift speed
 9.9 Plasma Parameter
 9.10 Plasma Beta
 9.11 Fermi Energy
 9.12 Plasma Degeneracy
 9.13 Plasma Coupling
 9.14 Critical Density

Creative Commons License
Plasma Derivations by Alex Zylstra is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.

Chapter 1
General Plasma

1.1 Definition of a plasma

A plasma is, at a basic level, an ionized gas. When electrons are stripped from ions, the gas becomes a quasi-neutral system of two or more interpenetrating fluids.

In a gas, microscopic behavior is governed by collisions between atoms or molecules. The macroscopic behavior is determined by taking volumes with large numbers of particles within and taking statistical averages over known particle distributions (i.e. a Maxwellian). A volume element of fluid must be small relative to the system size and contain a large number of particles for this to be valid.

In a plasma, electrons and ions still undergo collisions, but additionally they interact via the long-range Coulomb force. This allows the plasma to screen both DC and AC fields, as we will see in subsequent sections. This introduces two more conditions:

  1. λDe L
    The Debye length, or DC electric field screening length, must be much smaller than the system size for standard plasma behavior.
  2. ND = 4π 3 nλDe3 1
    The number of particles within a Debye sphere, as defined above, must be large. This plus 1 is analogous to the condition placed on ideal gases.
  3. ωpeτ 1
    If ωpe is a typical frequency for plasma oscillations, as governed by the Coulomb force, then the ’plasma frequency’ must be much greater than a typical frequency for hydrodynamic behavior ω = 1τ which could be taken as a collision frequency with neutral atoms. Thus we arrive at this condition.

1.2 Debye Length

Here we seek to define the typical length scale for DC electric field screening using simple arguments.

In an infinite homogenous plasma, we can take the electron distribution function as

fe(u) = Ae(mu22)kBT (1.1)

where A is a constant of proportionality, u is the velocity of a particle, Te is the temperature, and kB is Boltzmann’s constant. In the presence of an electrostatic potential, ϕ, this becomes

fe(u) = Ae(mu22+qϕ)kBT (1.2)

where q is the electron charge q = e. This is somewhat intuitive and also can be derived rigorously by taking the partition function for Boltzmann statistics and using qϕ as the ‘chemical potential’ μ. This derivation omitted here for simplicity.

Now consider the effect of a small deviation in density between electrons and ions, using the Poisson equation:

2ϕ = e 𝜖0(ni ne) (1.3)

where we have taken Z = 1, a hydrogen plasma, for simplicity. Next we assume that only the electrons move in response to the applied potential, which is reasonable given the small electron mass memi 1. In this case the ion density is simply the equilibrium value n0 and we can rewrite Poisson’s equation as

2ϕ = n0e 𝜖0 (1 nen0). (1.4)

We can calculate the perturbed electron density from the distribution functions by taking the ratio of Eq 1.2 to Eq 1.1:

ne n0 = eeϕkBT . (1.5)

Substituting into the Poisson equation,

2ϕ = n0e 𝜖0 (1 eeϕkBT ). (1.6)

In general this cannot be solved analytically, but in the limit eϕkBT 1 we can Taylor expand:

eeϕkBT 1 + eϕ kBT + ... (1.7)

which we put into the Poisson equation to get:

2ϕ = n0e 𝜖0 eϕ kBT (1.8)

2ϕ = n0e2 𝜖0kBTϕ (1.9)

which has solutions

ϕ(r) = ϕ0e|rr0|λDe , (1.10)


λDe = 𝜖0 kB T n0e2 (1.11)

which is the electron Debye length. Rigorously, the ion motion must be treated too, in which case

1 λD = 1 λDe + 1 λDi. (1.12)

Physically, the Debye length represents the distance an electrostatic field can penetrate into a plasma before it is screened. Electrons and ions move to oppose the imposed field. At zero temperature the screening is perfect, but at finite temperature the Debye length is finite due to thermal ions and electrons having enough energy to sample the screened potential. The density dependence is logical - if there are more available charges, the screening length with be short.

1.3 Plasma Frequency

In the previous section we considered the plasma DC response, now we consider the AC response. A simple intuitive approach is taken here, which is rigorously confirmed later.

If we consider a high-frequency oscillation, we can consider the ion inertia to be infinite.

A simple intuitive system is as follows. Consider semi-infinite slabs of rigid electron and ion fluids, with finite length L in the x direction. Take the ion slab as stationary and the electron slab as displaced by a distance x. There will be a restoring Coulomb force between the two slabs.

First we need the field of a charge slab. Using Poisson’s equation,

2ϕ = ρ𝜖 0 (1.13)

the charge density is ρ = ne with ni = ne = n taken for simplicity (Z=1 hydrogen plasma). Using Stokes’ Theorem, and E = ϕ,

SE dA = ρenc𝜖0 (1.14)

for a Gaussian pillbox, the field outside a single slab is,

E = nq 𝜖0 x. (1.15)

Now we consider the restoring force on the electron slab,

F = m = eE = ne2 𝜖0 x. (1.16)

which simplifies to

= ne2 me𝜖0x (1.17)

which of course has well-known oscillatory solutions of the form

x(t) = x0 + Aeiωpet (1.18)

with characteristic frequency

ωpe = ne2 me𝜖0. (1.19)

which is the electron plasma frequency.

At lower frequencies, the ion plasma frequency is also important, which is obtained by taking me mi.

Chapter 2
Single Particle Motion

2.1 Cyclotron Motion

First, we consider the motion of a single charged particle in constant uniform magnetic field - the cyclotron motion. So we take E = 0 and B = Bz where we can choose the magnetic field along the z axis without loss of generality. Then the Lorentz Force,

F = q(E + v ×B) = qv ×B (2.1)

leads to equations of motion

m = qvyBz (2.2) mÿ = qvxBz (2.3) mz̈ = 0 (2.4)

The z equation, of course, is trivial: constant motion in the direction. The other two axes can be addressed by taking the derivative of the y equation of motion to substitute in from the x:

mvy ̈ = qvx ̇Bz = qBzqvyBz m , (2.5)

which simplifies to

vy ̈ = qBz m 2v y = ωc2v y, (2.6)

where we have defined the cyclotron frequency ωc = qBzm. This clearly allows oscillatory solutions; without loss of generality we take

vy(t) = vcosωct (2.7)

combined with the original equations of motion, and using the initial conditions E = 1 2mv2 we arrive at the x velocity

vx(t) = vsinωct. (2.8)

Clearly the total perpendicular kinetic energy is conserved, vx2 + vy2 = v2, as expected since the magnetic force does no work.

Next we consider the particle position. We simply need to integrate the velocity equations once more, which gives:

x(t) = x0 v ωc cosωct (2.9) y(t) = y0 + v ωc sinωct (2.10)

The average position (x0,y0) is generally referred to as the ‘guiding center’ position. The quantity

rL = v ωc = mv qBz (2.11)

is the Larmor radius, which gives the size of gyrations due to the Lorentz force. These equations describe the motion of any charged particle in a constant and uniform magnetic field. The motion is generally helical gyrations about the magnetic field lines.

2.2 E ×B Drift

Next we consider the motion in constant and uniform E and B fields.

First, consider the situation where E B. Clearly this will lead to continuous acceleration along the magnetic field lines where v ×B is zero. This is not generally an interesting phenomenon, except in the case of runaway electrons, which will be considered later.

So for now, we consider the case E B. Without loss of generality we take E = Exx̂ and B = Bz. In this case the general Lorentz force,

F = q(E + v ×B), (2.12)

reduces to equations of motion,

m = qEx + qvyBz (2.13) mÿ = qvxBz (2.14) mz̈ = 0. (2.15)

Once again, the equation of motion along the magnetic field lines is trivial. Focusing on the x y plane, we start by taking the derivative of the y equation,

mvy ̈ = qBz (2.16)

which allows us to substitute the first equation of motion, giving,

mvy ̈ = qBz(qEx + qvyBz)m. (2.17)

Rearranging terms, we get that,

vy ̈ = q2Bz2 m2 ExBz + vy (2.18)

This clearly suggests a coordinate transformation,

vy = v y + ExBz (2.19)

in which case the equation of motion becomes

mvÿ = ω c2v y (2.20)

which is simply cyclotron motion in the primed coordinate system. But now that we have introduced the primed system, we note that it is drifting relative to the original reference frame with a constant velocity:

vD = ExBz (2.21)

More generally, if we allow E to lie anywhere in the xy plane, we get that

vE×B = E ×B B2 (2.22)

where this is called the E ×B drift for obvious reasons. We note the important result that this drift does not depend on the particle mass or the particle charge, which means that in a plasma electrons and ions will drift with the same velocity and the same direction. The E ×B drift can therefore cause net motion of the plasma, and in certain situations can be problematic for plasma confinement in magnetic fusion devices.

2.3 F ×B Drift

We observe that in the derivation of the E ×B drift, there is nothing ‘special’ about the electrostatic force

F = qE (2.23)

and we notice that we can substitute E = Fq in the previous result to get that

vF = 1 q F ×B B2 (2.24)

the general force drift. We note that in the general case there is a dependence on the particle charge q, meaning that electrons and ions will drift in opposite directions due to the ‘general force’.

2.4 B Drift

Up until now, we have only considered uniform fields. That will now change. Take E = 0 and non-uniform but constant B = Bz(x). The magnetic field points in the z direction but has a gradient in the x direction (taken without loss of generality). As usual we start with the Lorentz force,

F = q(E + v ×B) (2.25)

which gives equations of motion

m = qvyBz(x) (2.26) mÿ = qvxBz(x) (2.27) mz̈ = 0 (2.28)

Again, the z equation of motion is trivially solved. Focusing on the xy plane, the equations of motion are now more complicated than previously due to the position dependence of the magnetic field strength.

If we take the limit rL |B|B, or that the particle Larmor radius is much smaller than the gradient length scale for changes in B, then we can Taylor expand the magnetic field around the particle guiding center position:

Bz(x) Bz(x0) 1 + 1 B0 Bz x (x x0) (2.29)

substituting this into the y equation of motion,

mÿ = qvxBz(x0) 1 + 1 B0 Bz x (x x0) (2.30)

In general this would be very difficult to solve, but if we treat it perturbatively and use the general cyclotron motion as a 0th order solution, we can rewrite this as

mÿ = qvcos(ωct) Bz(x0) + Bz x rL cosωct. (2.31)

Now consider averaging over a gyration. The first term averages to 0, and the second term becomes (using cos2 12 over a whole gyration)

Fy = qvrL 2 Bz x (2.32)

using our previous definitions of rL,

Fy = qv2 2ωc Bz x (2.33)

Combined with our result for the general force drift, we can simply write the drift due to the magnetic field gradient as

vB = v2 2ωc B ×B B2 (2.34)

where we have made the generalization to arbitrary directions of the magnetic field gradient by intuition. We note that the B drift depends on both the particle charge and mass, which means that electrons and ions will drift in opposite directions at potentially different rates due to magnetic field gradients. The sign of the above equation is the sign of the charge. If the species temperatures are equal (Te = Ti) then mv2 is equal and the electrons and ions drift in opposite directions but at the same rate.

2.5 Curvature/Gravitational Drift

We can use our generalized force drift equation to examine two other interesting situations.

First, consider the effect of a gravitational field on the plasma. In terrestrial experiments there is an unavoidable force

F = mg (2.35)

which leads to a drift velocity

vg = m q g ×B B2 (2.36)

So a gravitational field will induce a drift, which depends on both the particle mass and charge. In real experiment, however, this is usually neglected due to the smallness of g relative to other forces (i.e. the electromagnetic force).

Now consider a curved plasma. If the magnetic field lines are bent with some curvature radius, how will that affect the particle motion?

The simple answer is to consider the particles as primarily streaming along the field lines. If the field lines are bent into a circle, as in a tokamak, then the curvature is essentially equivalent to considering the motion in a rotating reference frame. This induces a centrifugal force opposite the direction of curvature and perpendicular to the magnetic field:

Fcf = mΩ × (r ×Ω) = mv2Rc Rc2 (2.37)

this is written by inspection and intuition but could be derived rigorously. The velocity parallel to the magnetic field lines is denoted by v. If we substitute into the general force drift equation, then we get

vc = mv2 qRc2 Rc ×B B2 . (2.38)

The curvature is perpendicular to the radius of curvature vector and the magnetic field. Its magnitude depends on the radius of curvature, the magnetic field, and particle info. In particular we note that it is proportional to the ratio mq which means that electrons and ions will drift in opposite directions and differing rates (unless Te = Ti).

2.6 Polarization Drift

Consider a time-varying electric field given by

E = E0eiωtx̂ (2.39)

Also allow a constant uniform magnetic field B = Bz, so that the equations of motion become

m = qvyBz + qE0eiωt (2.40) mÿ = qvxBz (2.41) mz̈ = 0 (2.42)

as usual the z equation is trivial. Taking a derivative of the first equation,

vx ̈ = ωc2 v x iω ωc Ex ̃ Bz (2.43)

where we have used Fourier analysis Ė = iωE and the over-tilde denotes time varying quantities. If we define two drift velocities

vp = iω ωc Ex ̃ Bz (2.44) vE = Ex ̃ Bz . (2.45)

The latter drift is simply the E ×B drift but with a time-varying electric field; the drift velocity will vary sinusoidally with the imposed field as we might expect. The first drift is the so-called polarization drift. Rewriting the equations of motion,

vx ̈ = ωc2 v x vp , (2.46)

the motion in the x̂ direction is simply a superposition of the cyclotron motion and the polarization drift. In the ŷ direction,

vy ̈ = ωc2(v y vE) (2.47)

the motion is a superposition of cyclotron motion and the E ×B drift.

In the case of an electric field which cannot be simply decomposed into Fourier components, we can rewrite a more general expression by inspection,

vp = 1 ωcB E t (2.48)

2.7 Adiabatic Invariants

2.7.1 μ

This is the first adiabatic invariant. We want to start from the action integral, in general form,

pdq (2.49)

which will be an adiabatic invariant of the motion (in certain limits). First we consider the cyclotron motion of a particle, so let us take p = mvrL and q = 𝜃 and integrate over one gyration

pdq = mvrLd𝜃 = 2πmvrL (2.50)

Now we need to rewrite this quantity as

4π m |q|×mv2 2B (2.51)

The first part is constant if the charge-to-mass ratio is not changing over the motion (a good assumption). In that case, the second part of the above equation is conserved, leading to

μ mv2 2B (2.52)

This quantity μ is an invariant of the motion. However, we have assumed that pure cyclotron motion is a good approximation of the particle motion over short time scales, or equivalently, we have assumed slow motion ω ωc. In cases where this is not true, then μ is not an invariant of the motion.

2.8 Magnetic Mirror Machine

The adiabatic invariant μ can be directly applied to a machine of fusion (historical) interest - the magnetic mirror.

Consider a linear cylindrical machine where the magnetic field is weakest at the center and reaches a peak field value near the ends. This could be arranged, for example, by a two-coil configuration.

A particle near the center of the machine, i.e. a low-field region, is heading towards the high-field region. The problem is to define when the particle is confined. The initial particle kinetic energy can be decomposed into parallel and perpendicular components which will satisfy

v = v0 sin𝜃 (2.53)

v = v0 cos𝜃 (2.54)

where 𝜃 is the pitch angle relative to the magnetic field. Based on the previous derivation, in the slow motion limit, we will have the adiabatic invariant of the motion

μ = mv2 2B . (2.55)

We can see from the definition of μ that as B increases, the particle’s perpendicular velocity must increase to conserve μ. For a given particle’s kinetic energy, then, there is a maximum value of B that can be reached before the particle is reflected. Consider the extremes of the motion, and setting μi = μf we get

vi2 Bi = vf2 Bf (2.56)

or equivalently,

Bf Bi = vf vi (2.57)

using vf = v0 and the pitch angle definition,

Bf Bi = 1 sin2𝜃. (2.58)

Next we define the ‘mirror ratio’ of the experimental machine as the ratio of minimum to maximum field strength, R BmaxBmin > 1. In this case there is a critical pitch angle for reflection,

sin2𝜃 c = Bmin Bmax. (2.59)

Any particles with 𝜃 > 𝜃c are confined, while particles with 𝜃 < 𝜃c are not confined. This leads to a ‘loss cone’ in phase space.

The mirror machine was originally proposed as a scheme for fusion energy. Unfortunately, the loss cone particle loss is too extreme for an efficient machine. Coulomb scattering continually scatters particles into the loss cone, and they are lost out the ends of the machine. A rigorous derivation finds that the maximum theoretical Q value for the mirror machine is 1.1. While this is actually slightly greater than unity, a real machine will not achieve the theoretical result, and even if it did a Q = 1.1 machine will never be economical for energy generation.

2.9 Runaway Electrons

Consider a single electron moving in a tokamak. In today’s machines, large transformers are used to inductively drive a toroidal current in the machine. In the (realistic) event that the plasma has finite resistivity, this will create a non-zero toroidal EMF. The value of this field is not important for the problem; we simply consider the case of a toroidal electric field.

The electron equation of motion will be

mv̇ = qE νve (2.60)

where ν is the collision frequency for momentum transfer, generally ν = νei + νee. It is convenient to introduce the form we = vevTe where the denominator is the electron thermal velocity vTe = 2Te me. If we rewrite the electron equation of motion as,

dwe dt = Ê νr we2, (2.61)

then the parts of the right hand side are normalized to

Ê = e|E| 2me Te (2.62)


νR = 3 82π ne4 logΛ 𝜖02meTe32, (2.63)

which we will not derive here. The important part is that with the velocity normalization we have done, neither Ê or νR depend on we. We can then immediately see that for

Ê νRwe2 (2.64)


we2 ν RÊ (2.65)

the electron will continuously accelerate. This is the runaway condition. In particular, we can write down a ‘critical velocity’

wc = νR E ̂, (2.66)

where all electrons starting with a velocity greater than wc will continuously accelerate. This phenomenon arises because the electron collision frequency for momentum transfer is proportional to 1v3 so that fast electrons collide very rarely. Runaway electrons are potentially problematic for tokamak systems.

Chapter 3

3.1 2-fluid derivation

3.1.1 Momentum Equation

This is more a physical argument than a strict derivation. We start from the well-known Navier-Stokes equation, i.e. the momentum equation for ordinary hydrodynamics:

ρ u t + (u )u = p ρν2u, (3.1)

where ρ is the mass density, u is the ordinary fluid velocity, and p is the scalar pressure. ν in this equation is the fluid kinematic viscosity.

Right off the bat, we know to make the substitution ρ mjnj for plasma species j. The kinematic viscosity term ρν2u is absorbed with the scalar pressure term into a tensor pressure term, giving us:

mjnj uj t + (uj )uj = P + ... (3.2)

where the bold-face P denotes a tensor, and the ellipses are included to show that we are still missing terms in this equation.

The first obvious addition that must be included when going from ordinary hydrodynamics to plasma magnetohydrodynamics is the effect of the Lorentz force. While ordinary fluids are composed of charge neutral particles, plasmas of course consist of electrons and ions and the fluid species j generally has non-zero charge qj. In this case each particle feels a Lorentz force due to the local fields:

F = qj E + v ×B (3.3)

which is generalized to a force on a fluid element by multiplying by the particle density nj and allowing the particle velocity to go to the average fluid velocity. This term is added to the momentum equation to obtain

mjnj uj t + (uj )uj = njqj E + uj ×B P + ... (3.4)

There is one more effect we must include. In ordinary hydrodynamics the fluid motion is determined by collisions within the fluid, and by external forces represented by the pressure term. But in a plasma, multiple fluids can be interpenetrating, most obviously the electron and ion fluids which must be co-located to preserve charge quasi-neutrality. This allows for a transfer of momentum between fluids via collisions, which will generally obey a relation

mjnjνji(uj ui) (3.5)

for collisions between fluids i and j, with collision frequency for momentum transfer νji = νijmimj. This is added to the momentum equation to give us its final two-fluid form:

mjnj uj t + (uj )uj = njqj E + uj ×B P mjnjνji(uj ui). (3.6)

If there are more than two fluid species (i.e. in a multi ion species plasma) then we must sum the last term with i going over all other species in the plasma, though primarily we will be interested in electron+ion two-fluid plasmas.

3.1.2 Continuity Equation

In ordinary hydrodynamics the continuity equation is:

ρ t + (ρu) = 0. (3.7)

If we simply make the same substitutions, u uj and ρ njmj we get that

nj t + (njuj) = 0. (3.8)

3.1.3 Energy Equation

This is not as often quoted in fluid mechanics. Basically, we must think about how the internal energy of a fluid element changes. It can change via the first-order flow of the fluid, via compressive (PdV ) work, and via heat conduction. Additionally, there may be any number of ‘external’ sources or sinks of energy.

The average internal energy of a particle in an ideal gas, equilibrated to a Maxwellian, can be given by (32)kBT. The overall internal energy is the particle density multiplied by this quantity. Changes in density are included in the PdV term; here we consider the change in temperature due to the flow. We can simply write this via the convective derivative as:

3 2nj t + uj Tj. (3.9)

Next, we consider the compressive work done on (or by) a fluid element, which is simply the pressure multiplied by the divergence of the fluid velocity:

pjuj (3.10)

Finally, we can simply write down the heat flux in terms of the thermal conductivity and temperature gradient,

q (3.11)

with q = κ T.

We combine these various terms together, with the total change in internal energy from the three equal to the external sources or sinks of energy:

3 2nj t + uj Tj + pjuj + q = Sj (3.12)

where Sj represents several sources and sinks of energy:

3.1.4 Maxwell’s Equations

For completeness we have to include Maxwell’s equations for electromagnetism to the fluid equations above, since the electric and magnetic fields must be determined self-consistently in MHD. In plasma physics we prefer to use the vacuum equations with plasma serving as the source terms. Anyways, in SI units the four Maxwell equations are:

E = 1 𝜖0ρ (3.13) B = 0 (3.14) ×E = B t (3.15) ×B = 1 c2 E t + μ0j (3.16)

where ρ is the total charge density:

ρ iniqi (3.17)

and j is the total current density:

j iniqiui (3.18)

3.2 Single-fluid MHD

In many plasmas, the generalized two-fluid system described above can be simplified considerably to a one-fluid model. The goal of this section is to present a derivation of the one-fluid model with a discussion of when it is applicable. For simplicity we will use a single-species Z = 1 plasma with the usual assumptions.

First, we need to define the single-fluid variables which will be used throughout this section:

ρ = mini + mene (3.19) v = 1 ρ(miniui + meneue) (3.20) j = e(niui neue) (3.21)

We start off by writing the momentum equations for both electron and ion species. For simplicity we take the pressure tensor as isotropic, meaning that the pressure appears only as the usual scalar pressure. The viscosity terms are small if the ion Larmor radius is much smaller than scale lengths for typical variations, for instance this condition can be written as rLi ρρ using the density gradient length scale.

We also drop the convective terms (uj )uj. Chen describes this as ‘hard to justify’. I think that the best justification is that in motions where the fluid velocity is ‘small’, and implicitly the 0th-order velocity is zero, then this convective term is second-order in the small velocity and can be ignored. A valid question is ‘small relative to what?’ If we consider the ratio of the convective term to the partial time derivative:

(u )u ut iku2 iωu u ωk (3.22)

If the plasma response is, say, acoustic (the slowest case), then ωk cs and the convective term is negligible if the fluid motion is sufficiently sub-sonic.

In any event, we can now write the ion momentum equation as

miniui t = nie(E + ui ×B) pi + miniνie(ui ue) (3.23)

and similarly for the electrons:

meneue t = nee(E + ue ×B) pe meneνei(ue ui) (3.24)

where in the collisional momentum transfer term (the last one) we note that νie = νeimemi, which leads to the physically intuitive result that momentum lost by one species is gained by the other (terms are equal and opposite).

We start off by taking the sum of the two momentum equations, which will lead to the single-fluid equation of motion. We get that,

n miui t + meue t = en(ui ue) ×B (pi + pe) (3.25)

The left-hand side is simply the time derivative of the single-fluid velocity multiplied by the density. The electric field terms have canceled. The magnetic field part of the Lorentz force reduces to using the single-fluid total current. We take the single-fluid total pressure as the sum of the individual ion and electron pressures, which is sensible. By taking the sum, the collision term has also canceled (what is lost by one fluid is gained by the other). We are thus left with the single-fluid equation of motion:

ρv t = j ×B p. (3.26)

In some cases a gravitational force term + ρg is added to the right hand side, which is not included in this derivation but can be seen via physical intuition.

To get to the various Ohm’s Laws, we must take the less-obvious step of calculating the difference of the ion and electron momentum equations. But in particular, we multiply the ion equation by the electron mass and vice versa, so that the difference becomes

mimen t(ui ue) = en(mi + me)E + en(meui + miue) ×B mepi + mipe (mi + me)nνei(ue ui), (3.27)

which is kind of a mess but simplifies considerably. The left hand side can be written in terms of the current,

mimen t(ui ue) = mimen e t j n . (3.28)

The electric field term can easily be simplified by recognizing that it is just eρE. The magnetic field term is trickier. We have to write:

meui + miue = miui + meue + mi(ue ui) + me(ui ue) = ρv n (mi me)j ne (3.29)

Then the entire magnetic field term is:

eρv ×B (mi me)j ×B (3.30)

The collision term (momentum transfer) is re-written in terms of the resistivity η = meνeine2 so that it becomes

(mi + me)neηj (3.31)

Combining all of these together we get

mimen e t j n = eρE + eρv ×B (mi me)j ×B mipe + mepi (mi + me)neηj (3.32)

An immediate and obvious simplification is to take the limit me mi which means that the terms (mi ± me) mi and min ρ. We also observe that, since pj = njkBTj the electron and ion pressures are pi pe for Ti Te in which case mepi mipe so we can drop the ion pressure gradient term. We can also rearrange terms and divide by ρe to get that

E + v ×B ηj = 1 eρ mimen e t j n mipe + mij ×B (3.33)

We then generally take the limit of slow motions, i.e. where inertial (cyclotron motion) effects are unimportant. Explicitly this is the limit va ωci where a is the plasma scale size. We can get this result by considering the ratio of the current time derivative term to the j ×B term within square brackets:

(mimene)t(jn) mij ×B meωj ejB ω ωce (3.34)

So, if ω ωce we can neglect the time derivative term.

In any event, this limit reduces the equation to:

E + v ×B = ηj + 1 en j ×B pe . (3.35)

This is known as the Generalized Ohm’s Law. The j ×B term is the Hall current term. In many physical systems it turns out that the Hall current and pressure gradient terms are small, in which case this expression reduces to the Resistive MHD Ohm’s Law:

E + v ×B = ηj (3.36)

In some plasmas the resistivity is small enough to be neglected (i.e. high temperature and low density). In this case the resistive Ohm’s law reduces to the Ideal MHD Ohm’s Law:

E + v ×B = 0 (3.37)

In summary, the single-fluid MHD model is:

ρv t = j ×B p (3.38) ρ t = (ρv) (3.39) E + v ×B = ηj + 1 ne j ×B + pe (Generalized) (3.40) = ηj(Resistive) (3.41) = 0(Ideal) (3.42)

with the various options given for Ohm’s Law (more complex ones can be derived but are not included).

3.3 Pinches (MHD Equilibria)

In this section we consider a few cases of MHD equilibria analysis, using the single-fluid MHD theory just derived. First, we note that in equilibrium several simplifications can be made to the single-fluid equations, in particular by taking time derivatives to zero (no change in an equilibrium situation). For simplicity we take ideal MHD.

The continuity equation can be rewritten as:

dρ dt + ρv = 0 (3.43)

in equilibria this implies that v = 0, which means also that j = 0. There are no particle or current sources/sinks in equilibrium solutions in our simple model (which does not include, for instance, current drive or neutral-beam heating).

The momentum equation becomes simply:

p = j ×B (3.44)

The complete set of equations which define simple MHD equilibrium systems is combined with two of Maxwells’ equations; in total we have that:

p = j ×B (3.45) ×B = μ0j (3.46) B = 0 (3.47)

One important qualitative feature of solutions to this set of equations is that both the lines of force (B field lines) and the current j must lie in a plane perpendicular to the pressure gradient. The rest of this section is devoted to identifying a few useful solutions of the MHD equilibrium problem.

3.3.1 𝜃 Pinch

Consider a cylindrical plasma which has a current flow in the azimuthal (𝜃̂) direction, thus the name. This can be created via the diamagnetic effect of plasmas by an imposed axial field. In any event, the unknowns are p = p(r), Bz = Bz(r) and j𝜃 = j𝜃(r) with the total field and current given by B = Bz(r) and j = j𝜃(r)𝜃̂.

By symmetry the condition B = 0 is automatically satisfied. Ampère’s Law reduces to

dBz dr = μ0j𝜃. (3.48)

And the pressure balance equation becomes

dp dr = j𝜃Bz = 1 μ0BzdBz dr = 1 2μ0 d dr Bz2 (3.49)

this can be rewritten, combining the derivatives with respect to r, to get that

d dr p + Bz2 2μ0 = 0 (3.50)

which can be trivially integrated, yielding the pressure balance equation:

p(r) + Bz2(r) 2μ0 = B02 2μ0 (3.51)

where we have automatically solved for the constant of integration by noting that the field Bz must become the vacuum value for radii outside the plasma where p = 0, i.e. r > a.

3.3.2 Z Pinch

We now consider the complementary system, the Z pinch. In this case, the current flows along the axis of a cylindrical plasma j = jz(r) which creates an azimuthal magnetic field B = B𝜃(r)𝜃̂ + B0. There is generally assumed to be an imposed magnetic field along the axis as well, given by the B0 term. The current in this pinch could be driven by a set of electrodes at either end of the plasma.

Once again we have the condition B trivially satisfied by the geometry. And Ampère’s Law becomes:

μ0jz = 1 r d dr(rB𝜃). (3.52)

Combining this with the pressure balance equation,

dp dr = jzB𝜃 = B𝜃 μ0r d dr(rB𝜃). (3.53)

If we expand the rB𝜃 derivative we get that:

dp dr + 1 2μ0 d dr(B𝜃2) + B𝜃2 μ0r = 0. (3.54)

which can be simplified by combining the ddr terms:

d dr p + B𝜃2 2μ0 + B𝜃2(r) μ0r = 0 (3.55)

At this point one must typically assume a current profile to solve analytically for the pressure profile.

3.3.3 Screw Pinch

A screw pinch is a linear combination of 𝜃 and Z- pinch, named because the magnetic field lines wrap around the cylindrical plasma like the threads on a screw (or the stripes on a barber’s pole). The defining pressure balance equation can be written by combining Eqs. 3.50 and 3.55:

d dr p + Bz2 2μ0 + B𝜃2 2μ0 + B𝜃2(r) μ0r = 0 (3.56)

The equation has three unknowns - p(r),Bz(r),B𝜃(r). In general solutions are obtained by assuming profiles for two of the three; MHD will then define the third.

3.4 Diamagnetic Drift

In our previous discussion of particle drifts, we left out one important drift - the diamagnetic drift. This is because the diamagnetic drift arises as a fluid effect and not from single-particle motion. We start out with the momentum equation for an arbitrary fluid (omitting subscripts for clarity):

mn u t + (u )u = nq(E + u ×B) p. (3.57)

First, consider the effect of the time derivative term. Let t iω and take the ratio of the time derivative term to the u ×B term considering only the motion perpendicular to the magnetic field:

iωmnv qnvB = ω ωc (3.58)

For drifts slow compared to the cyclotron motion, ω ωc we can neglect the time derivative term.

What about the term (u )u ? Observe that if the drift velocity we will derive ends up being vd u then this term will be 0. Therefore we will assume off the bat that this is the case, and at the end of the derivation we can check to make sure that this assumption is valid.

So what we have is that for slow drifts with the drift velocity perpendicular to the gradient,

0 = nq(E + u ×B) p. (3.59)

If we take the equation’s cross product with B,

0 = nq(E ×B + (u ×B) ×B) p ×B. (3.60)

Using the vector identity A × (B ×C) = (A C)B (A B)C for the middle term:

0 = nq E ×B + (u B)B uB2 p ×B. (3.61)

By definition of course uB = 0. If we take this equation, rearrange and solve for u we get that:

u = E ×B B2 p ×B nqB2 . (3.62)

The first term is the familiar E ×B drift. The second term is the Diamagnetic Drift:

vD = p ×B nqB2 . (3.63)

We observe that the diamagnetic drift is indeed perpendicular to the pressure gradient, as we had assumed before. Due to the dependence on q, electrons and ions will drift in opposite directions and thus the diamagnetic current can induce azimuthal currents in cylindrical plasmas. For Z = 1 and p = nkT, we can write the diamagnetic current as

jD = ne(uDi uDe) = (kTi + kTe)B ×n B2 . (3.64)

3.5 Drift Waves

The diamagnetic drift of the previous section can lead to an unstable situation in cylindrical or toroidal plasmas, leading to so-called drift waves.


Figure 3.1: Drift wave geometry

The basic idea is as follows. We have a plasma magnetically confined by the 0th order B field, shown as out of the page (or in the direction). Due to the density gradient as shown, the diamagnetic drift can set up electric fields if there is a perturbation in the isobaric surface (shown in red). This electric field causes E ×B drifts, shown as v1. A detailed analysis follows.

The zeroth-order drifts are due to the diamagnetic drift:

ve0 = vDe = kBTe eB0 n0 n0 ŷ (3.65)

and also for the ions. Since the electrons can flow along the magnetic field, and we assume that this wave is slow relative to the electron motion, we can use the Boltzmann relation for the electron density:

n1 n0 = eϕ1 kBTe (3.66)

which is to first order, in the limit eϕ kBTe. Where the plasma is perturbed to higher pressure (higher density), as in the top of the figure, the potential is positive. The electric field thus points from peak to trough of the wave, as shown.

We can now write the first-order velocity due to the E ×B drift as:

v1x = Ey B0 = ikyϕ1 B0 (3.67)

which can be obtained from E = ϕ = ikϕ. We can also use a sort of ‘continuity equation’ for guiding centers to write:

n1 t = v1xn0 x (3.68)

Doing some Fourier analysis on this continuity equation,

iωn1 = ikv1xn0 = ikyϕ1 B0 n0 (3.69)

Where n0 represents the density perturbation due to the warped isobar. Using the Boltzmann relation from earlier we can also write (separately) that

iωn1 = iωn0 eϕ1 kBTe (3.70)

Equation the RHS of these two gives:

kyϕ1n0 B0 = ωn0eϕ1 kBTe (3.71)

leading to:

ω ky = kBTe eB0 n0 n0 = vDe (3.72)

The waves travel at the electron diamagnetic drift velocity. This is the behavior in the azimuthal direction. There is a small component kz ky which causes these perturbations to propagate. Along the way we implicitly assumed that electron currents cannot simply flow along B0 to neutralize the electric field; this means that the plasma must be resistive (and these are sometimes called ‘resistive drift waves’).

We have also neglected the instability analysis. One can see via qualitative arguments that this situation is very similar to the gravitational instability, and thus it can be argued that the perturbations in isobars tend to grow. The full analysis requires treating the polarization drift and non-uniform E drift as well.

3.6 Two-stream Instability

We wish to consider the following scenario: an infinite uniform plasma, where the ions are stationary in the lab frame but the electrons have a non-zero uniform ue0. The plasma is unmagnetized (B0 = 0) and cold (Ti = Te = 0) for simplicity. The fluid momentum equation is:

mn u t + (u )u = nq(E + u ×B) p mnν(u u0) (3.73)

In a 1-D flow the u ×B first-order terms are zero for the electrons (ions have no first-order contribution from this term). Further, the p terms are zero because we assumed a cold plasma. The collisional momentum transfer term also goes to 0 because we note that ν T32.

For the ions, (u )u contains no first-order terms since u0 = 0. For the electrons, a term (u0 )u1 remains. For simplicity we take Z = 1. We thus get an ion equation:

minui1 t = neE1, (3.74)

and an electron equation:

men ue1 t + (ue0 )ue1 = enE1. (3.75)

Assuming electrostatic waves of the form:

E1 = E1ei(kxωt)x̂, (3.76)

we can use Fourier analysis to simplify the ion equation to:

iωminui1 = neE1, (3.77)

ui1 = ie ωmiE1x̂. (3.78)

And the electron equation similarly:

iωmenue1 + ikmenue0ue1 = enE1 (3.79)

ue1 = ie me(ω ku0)E1x̂. (3.80)

Next we use the fluid equation of continuity, which is generally:

n t + (nu) = 0. (3.81)

So for the ions, employing our previous assumptions and Fourier analysis,

iωni1 + ikn0ui1 = 0 (3.82)

ni1 = kn0 ω ui1 = iekn0 miω2 E1 (3.83)

where we have used the previously-derived relation between ui1 and E1. Following through the same steps for the electrons:

ne1 t + n0ue1 + (ue0 )ne1 = 0 (3.84)

iωne1 + ikn0ue1 + ikne1ue0 = 0 (3.85)

ne1 = kn0 ω kue0ue1 = iekn0 me(ω kue0)2E1. (3.86)

We need one more equation to close these relations. We note that this is a high-frequency oscillations so we can use Poisson’s equation to eliminate the electric field:

E = e 𝜖0(ni ne). (3.87)

Using the first order density perturbations we have already derived on the right side, and Fourier analysis again on the left side, we get:

ikE1 = e 𝜖0 iekn0 1 E1 × 1 miω2 1 me(ω kue0)2 (3.88)

We can now eliminate the electric field, and recognizing the plasma frequency ωpe = ne2me𝜖0 on the right hand side,

1 = ωpe2 memi ω2 1 (ω kue0)2 . (3.89)

An important question is whether the plasma is stable or unstable to perturbations. If we define the dimensionless variables x = ωωpe and y = kue0ωpe we can define:

F(x,y) memi x2 + 1 (x y)2 = 1 (3.90)

There are singularities at x = 0 and x = y. But we note that solutions to F = 1 are solutions to the dispersion relation. If y is sufficiently large then there are 4 real roots for x and thus ω, which all correspond to stable though oscillatory solutions. On the other hand, if there are only two real solutions for ω then there are necessarily two complex, one of which will be a damped oscillation and the other will be unstable (growing). This occurs for sufficiently small values of y, or more physically, for small kue0 or large-wavelength perturbations.

3.7 Flux Freezing

This is an important question in the dynamics of magnetized ideal plasmas. We want to know how plasma elements move with respect to the magnetic field lines, or vice versa. Consider the magnetic flux through a surface of plasma:

ψ(t) =B n̂dS (3.91)

where n̂ is the unit vector normal to the surface S. We now want to explore how the magnetic flux through the plasma changes with time. The time derivative of the previous equation can be written as:

dψ dt =B t n̂dS + B ×u dl. (3.92)

The left hand side is a simple time derivative. The right hand side has been decomposed: flux can change due to a changing magnetic field (the first term) or due to motion of the plasma surface relative to the magnetic field (the second term).

Using Faraday’s Law,

×E = B t , (3.93)

and the ideal Ohm’s Law,

E = v ×B, (3.94)

dψ dt = [(vu) ×B] dl. (3.95)

The important result is that for v = u then the flux is unchanging with time. Physically this occurs when the magnetic field lines are moving with the plasma, and we call this ”frozen in flux”. The obvious analogy is a superconductor. It is a well-known result that magnetic flux cannot penetrate a superconducting volume. If it did, and the flux was changing in time (as it must) then it would induce an EMF in the superconductor, which would cause an infinite current due to the zero resistivity.

A real plasma will have finite resistivity but in many scenarios the resistivity is small enough that ideal MHD is a decent approximation. This is particularly true in certain high-temperature magnetically confined plasmas and in very low-density space plasmas.

3.8 MHD Instabilities

First we discuss the classification schemes generally used for MHD instabilities, then discuss qualitatively a few important examples.

3.8.1 Classification

First an instability is generally described as internal or external mode. An internal mode instability is one in which the plasma surface does not move. Generally these affect transport and impose operational limits but do not impact confinement. Conversely, an external mode instability is one in which the surface of the plasma moves. These can cause loss of confinement.

Next we classify the source of the instability. Generally a non-equilibrium current flows, which modifies the MHD. If a current flows perpendicular the magnetic field, these instabilities are often called pressure-driven since p = J×B. On the other hand, if a parallel current flows to drive the instability we call it a current-driven mode.

Finally, the wall characteristics are often important. We can see this simply as follows: If the plasma surface moves towards the wall in an ideal MHD system (frozen flux), the magnetic flux between the plasma surface and wall will be compressed. This will affect the plasma motion and feed back. If the wall is resistive then the compressed flux will ohmically dissipate, but if it is superconducting then this cannot happen. We can therefore discuss no wall, conducting wall, or superconducting wall configurations (though the last is operationally difficult).

3.8.2 Qualitative Descriptions

Z-pinch interchange / sausage


Figure 3.2: The Z-pinch interchange or ‘sausage’ instability.

We start with a qualitative description of the Z-pinch interchange or ‘sausage’ instability, which is depicted in Fig. 3.2. As a reminder, the 0th order current flows axially and thus left to right (or vice versa) in the figure. Consider a perturbation where the surface of the plasma is rippled such that the radius varies with axial position, in particular shown for r1 > r2 > r3. In this case we know that B1 < B2 < B3 which implies that the magnetic confinement is higher at 3 than at 1. This causes the plasma to expand at 1 and contract further at 3.

Z-pinch twisting instability


Figure 3.3: The Z-pinch twisting instability.

Continuing with the Z-pinch equilibrium confinement scheme, consider the scenario in Fig. 3.3. In this case the cylindrical symmetry is broken by allowing the axis (dotted line) to be twisted. Since the current is flowing axially in the plasma, this creates a magnetic field scheme as shown in blue. At the top of the indicated twist, the azimuthal magnetic field is stronger than the equilibrium value. Below the twist it is weaker. This creates a force on the plasma that tends to reinforce the kink because the confinement is coming from the MHD equilibrium term p = j ×B.

Screw pinch kink instability


Figure 3.4: The screw pinch twisting or ‘kink’ instability.

This configuration is very similar to the previous. We simply note that a screw pinch suffers from the same sensitivity to kinking or twisting as the Z-pinch. See Fig. 3.4. The mechanism is exactly the same as in the previous section, but now we have to remember that there is a potentially large axial/toroidal field Bϕ as well as the B𝜃 field.

Favorable vs unfavorable curvature


Figure 3.5: Favorable versus unfavorable curvature.

The previous sections raise a general question - when is curvature of the plasma surface favorable versus unfavorable? Consider Fig. 3.5. The equilibrium surface of the plasma is denoted in black. We follow the general coloring scheme (plasma is red, magnetic field is blue) used in the previous sections. Consider a perturbation from the equilibrium, as shown.

The general picture follows again from considering the p = j ×B confinement in MHD equilibria systems. When the plasma surface is perturbed as shown on the left of Fig. 3.5, it will tend to grow because the confinement is not good at the peak. More generally, it turns out that when the surface of the plasma is curved towards the plasma then the surface is unstable, so we call that ‘unfavorable’ curvature. Conversely, when the curvature is away from the plasma the system is well-confined and we thus call that ‘favorable’ curvature.

Gravitational instability

This instability is named the ‘gravitational’ instability but really can occur for any non-electromagnetic force applied to the plasma.


Figure 3.6: Equilibrium state for the gravitational instability.

Consider the plasma configuration shown in Fig. 3.6. The plasma is confined by a magnetic field, shown as out of the page, with a well-defined surface horizontal on the page. There is a vertical force away from the plasma surface denoted by g. We take g to be any generalized non-electromagnetic force. Gravity is one potential application.

As we saw in the single particle motion derivations, there is a generalized force drift which can be written:

vD = 1 q g ×B B2 . (3.96)

We can immediately see that the electrons and ions drift horizontally on the page due to the force g, and in opposite directions as denoted by vDi and vDe in the figure.


Figure 3.7: Perturbed state for the gravitational instability.

Now we have to consider what will happen if the surface of the plasma is perturbed. Consider Fig. 3.7, in which we have assumed an imposed sinusoidal horizontal perturbation in the plasma surface. The imposed force g is still vertically down and B0 is still out of the page. Because of the g ×B0 drift, the ions tend to ‘pile up’ on the left side of a plasma protrusion while the electrons pile up on the right hand side of a protrusion. This creates a first-order electric field E1 which is horizontal in this scheme.

On the right where the plasma is protruding from its equilibrium position, the first-order drift E1 ×B0 points down, as shown. As we know the E × B drift affects ions and electrons equally, and thus this creates a net force on the plasma that reinforces the perturbation. In cases where the plasma is depressed from its equilibrium location the first-order electric field is reversed in direction, and thus the E1 ×B0 drift also reinforces the perturbation. Therefore the plasma configuration explored in this section is unstable to perturbations in the surface as shown.

It is worth noting that this can be thought of as a magnetically-confined plasma version of the well-known hydrodynamic Rayleigh-Taylor instability. In this case the plasma is a ‘heavy’ fluid supported by the magnetic field, a ‘light’ fluid.

Drift instability / wave

Qualitatively one can arrive at the same sort of instability without an imposed external force g if instead there is a density gradient. Recall that the diamagnetic drift is defined as:

vD = p ×B nqB2 . (3.97)

If instead of the external force, the plasma in Fig. 3.6 had a vertical density gradient (high density at the top) then the zero order g ×B0 are instead replaced by the electron and ion diamagnetic drifts. The remaining analysis is exactly the same.

This is a very important instability because in cylindrical systems, there is necessarily a radial pressure gradient and thus the plasma is susceptible to drift instabilities. In a cylindrical plasma this is sometimes referred to as the flute instability instead of the drift instability since the resulting plasma surface looks like a fluted column.

In toroidal systems the same instability occurs. The result of detailed analysis is that there is a small component along the toroidal axis. In this case the perturbations curve slightly and wrap around the plasma like screw threads or stripes on a barber’s pole. This causes the perturbations to propagate and they are called true ‘drift waves’.

Weibel Instability


Figure 3.8: The Weibel instability.

The final instability we consider is the Weibel instability. This situation is illustrated in Fig. 3.8 and described here. Consider a plasma where the electron temperature is much hotter in one direction than the other two. This can arise in magnetically confined systems (where and directions have different behavior) or in laser plasmas, particularly in direct laser illumination between the ablation surface and the critical surface.

Anyways, in the figure we take Ty Tx,Tz. as shown. Consider a randomly-arising magnetic field in the x z plane as shown in blue. Due to the field, electrons will tend to curve as shown for a few different locations. A detailed analysis shows that the electrons tend to curve to form current sheets in the vertical (ŷ) direction. These current sheets end up reinforcing the imposed field. The plasma is therefore unstable to randomly-arising magnetic fields in the x z plane.

Chapter 4

4.1 Gas Dynamics Sound Waves

Before we get into plasma physics, we start with ordinary gas dynamics. Starting with the Navier-Stokes equation:

ρ u t + (u )u = p ρν2u, (4.1)

we can derive sound waves as follows. Neglecting kinematic viscosity (ν 0), assuming u0 = 0 and taking small waves (so keep only first order terms), we get that

ρu t = p. (4.2)

Using Fourier analysis,

iωρu = ikp. (4.3)

We know that the group velocity of waves is ωk so we can rewrite this as

ω = kpuρ. (4.4)

We will also need to use the continuity equation:

ρ t + ρu = 0, (4.5)

ikρu iωρ = 0, (4.6)

u = ωk. (4.7)

Using this with the momentum equation derived dispersion relation, we get that

w = kpρ, (4.8)

which directly implies

vg = cs = p ρ. (4.9)

We implicitly assumed γ = 1 earlier, a more general result is that

cs = γp ρ (4.10)

4.2 Electromagnetic Waves

Before we introduce plasma effects, it is useful to derive electromagnetic waves without plasma first. First we should write down Maxwell’s equations:

E = ρ𝜖0 (4.11) B = 0 (4.12) ×E = B t (4.13) ×B = μ0j + 1 c2 E t (4.14)

In a vacuum we can immediately drop the source terms

E = 0 (4.15) B = 0 (4.16) ×E = B t (4.17) ×B = 1 c2 E t (4.18)

For plane waves the first two equations will be automatically satisfied. Taking the curl of Faraday’s law,

× (×E) = t(×B), (4.19)

and combining with Ampère’s law,

× (×E) = 1 c2 2E t2 . (4.20)

If we assume plane wave Fourier solutions of the form

E = Eei(ωtkz)x̂ (4.21)

and use the vector identity

× (×A) = (A) 2A (4.22)

We get that

2E = 1 c2 2E t2 (4.23)

which reduces to, using Fourier analysis,

w = ck (4.24)

which leads to the group and phase velocities:

vg ω k = c (4.25)

vp ω k = c (4.26)

In media we generally define the index of refraction,

n = ckω (4.27)

and the group / phase velocities become

vg = cn (4.28)

vp = cn (4.29)

4.3 Polarization

In the derivation of waves we neglected the geometry of the system, which leads to a discussion of the wave polarization. Starting off with the Fourier component solution for the electric field,

E = Eei(ωtkz)x̂ (4.30)

we note that the solutions here are waves with E x̂ but the wave propagation direction is . We can also find the magnetic field for this configuration, it is obtained from

×B = 1 c2 E t (4.31)

and the solution is

B = Eω c2kei(ωtkz)ŷ = Bei(ωtkz)ŷ (4.32)

with B = Ec. So to summarize, in this nice plane polarized solution we have E x̂, B ŷ, and finally k . The choice of coordinate system does not cause loss of generality. It is important to note that we will always have E ×B k in a vacuum.

We could have chosen a wave solution in which the electric field vector rotated in the xy plane as it propagated along z. In this case the magnetic field vector will also have to rotate. These are circularly or elliptically polarized solutions. We can characterize these solutions by taking the ratio

iExEy (4.33)

when this = 1, the wave is right-hand circularly polarized. When it is = 1 the wave is left-hand circularly polarized instead.

4.4 Electrostatic Plasma Oscillations

We start off with the MHD picture of an electrostatic plasma oscillation. We will see at the end why it is an ’oscillation’ and not a wave. We need to start with the electron momentum equation:

mene ue t + (ue )ue = nee(E + ue ×B) pe mnν(ue ui). (4.34)

We need to make several simplifications and approximations in this section to make the problem tractable. First, we take Te = 0 meaning that pe = 0. Next, we neglect collisions, ν = 0. Next we take ue0 = E0 = B0 = 0 and keep only first-order quantities:

mn0ue t = n0eE (4.35)

for convenience we have dropped the e subscripts. If we assume sinusoidally oscillatory solutions, then we can use Fourier analysis to get that

iωmn0u1 = n0eE1 (4.36)

ωmu1 = ieE1 (4.37)

We have also assumed plane waves. This is the basis of our dispersion relation.

Next we can invoke the continuity equation:

n t + (nu) = 0. (4.38)

Using the same set of simplifications, we get:

iωn1 + ikn0u1 = 0. (4.39)

We can use this to rewrite the dispersion relation to eliminate u:

u1 = ωn1 kn0 (4.40)

So that the dispersion relation becomes:

ω2mn1 kn0 = ieE1. (4.41)

Next we need to eliminate both E1 and n1. We can do this by the use of Poisson’s equation. As we know we have to be careful in plasma physics as to when we can invoke Poisson’s equation, but in this case since we assumed high frequency waves and that the ion inertia is infinite, there must be a resulting electric field.

E = e 𝜖0(ni ne) (4.42)

Using ne = n0 + n1 and ni = n0 plus our usual set of assumptions, this simplifies to

ikE1 = en1𝜖0 (4.43)

E1 = en1 ik𝜖0 (4.44)

Using this result the dispersion relation becomes

ω2mn1 kn0 = ie2n1 ik𝜖0 (4.45)

which can be rewritten as

ω2 = n0e2 me𝜖0 = ωpe2 (4.46)

These disturbances do not propagate, k = 0, but instead are stationary oscillatory solutions with a given frequency, known as the plasma frequency:

ωpe = ne e2 me𝜖0 (4.47)

4.5 Electron Plasma Waves

The electron oscillations we derived in the previous result can become true waves (i.e. k0) when there is a finite electron temperature. They are called Electron Plasma Waves, or sometimes Bohm-Gross Waves or Langmuir Waves. Revisiting the electron momentum equation, with e subscripts omitted:

mn u t + (u )u = ne(E + u ×B) p mnν(u ui) (4.48)

Again we neglect collisions, and take plane waves with small perturbations (1st order terms only) but keep the pressure gradient, and take u0 = E0 = B0 = 0:

mnu1 t = neE1 p (4.49)

Generally we must make an assumption about the plasma equation of state, for example adiabatic:

p ργ p p = γn n (4.50)

if we then substitute the ideal gas law p = nkT into the above, and in this case, if we take 1-D isothermal compression and expansion of the plasma then γ = 3.

p = γkTn = 3kTn1 (4.51)

since there is no gradient in the equilibrium density n0. Plugging this into the momentum equation, taking kBT T for simplicity of notation,

mnu1 t = neE1 3Tn1, (4.52)

which we can use Fourier analysis to simplify further to:

iωmnu1 = neE1 3ikTn1. (4.53)

We need two relations between u1,E1,n1 to eliminate them. First we use the continuity equation:

n t + (nu) = 0 (4.54)

which, with first-order terms Fourier analyzed:

iωn1 + ikn0u1 = 0. (4.55)

Solving for u1:

u1 = ωn1 kn0 (4.56)

Next, we can use the Poisson equation (valid because these are fast oscillations, and ion momentum can be approximated as infinite), to get a relation between n and E. To start:

E = e 𝜖0(ni ne) (4.57)

with ne = n0 + n1 and ni = n0, and using ik,

ikE1 = en1𝜖0 (4.58)

or solving for E1:

E1 = en1 ik𝜖0 (4.59)

Using these expressions for u1 (continuity eq) and E1 (Poisson eq) the momentum equation becomes:

iωmn0ωn1 kn0 = ne en1 ik𝜖0 + 3ikTn1 (4.60)

ω2 = n0e2 m𝜖0 + 3n0T m k2 (4.61)

using the definition of the thermal velocity vth = 2Tm, we get the electron plasma wave dispersion relation:

ω2 = ω pe2 + 3 2vth2k2. (4.62)

We can see immediately that as T 0, then vth 0 and this dispersion relation reduces the electron oscillation found previously. Furthermore, we note that propagation requires ω2 > ωpe2 for the electron plasma wave.

4.6 Ion Acoustic Waves

Next we consider the ion acoustic wave, which as we will see is the analog of traditional hydrodynamic sound waves. Since the ions are involved, it is by necessity a slow wave. We therefore take n = ne = ni, since the electron response will be very fast compared to the wave. This requires that we do not use Poisson’s equation in the derivation!

Starting with the ion momentum equation:

min ui t + (ui )ui = ne(E + ui ×B) pi (4.63)

I note that we have omitted the collision term. Taking u0 = E0 = B0 = 0, keeping only first-order terms, and using E = ϕ, we simplify this to:

minu1 t = neϕ1 γTin1. (4.64)

Doing our traditional simplification with Fourier analysis t iω and ik, and taking plane waves, we get that:

iωmin0u1 = ikn0eϕ1 ikγTin1. (4.65)

To simplify this further requires that we obtain relations between u1, ϕ1, and n1 so that they can be eliminated.

First, we use the continuity equation as always:

n t + (nu) = 0, (4.66)

iωn1 + ikn0u1 = 0, (4.67)

u1 = ωn1 kn0 (4.68)

Next we note that in the presence of an electric potential, the electrons can be written as obeying the Boltzmann distribution:

ne = n0eeϕkT = n 0 1 + eϕ T + ... (4.69)

since the electron response is much faster than changes in ϕ. In the last step we used the assumption eϕ T to Taylor expand the exponential. Recognizing the contribution from the perturbation in density, n1, we write:

n1 = n0eϕ1 Te . (4.70)

Or equivalently, we can eliminate ϕ using:

ϕ1 = n1Te en0 (4.71)

Now we can use these two results to solve the momentum equation for the dispersion relation.

iωmin0ωn1 kn0 = ikn0en1Te en0 ikγTin1. (4.72)

We can eliminate n1 now and further simplify algebraically to get that:

ω2 = k2 Te + γTi mi = cs2k2 (4.73)

with the equivalent of the sound speed in a plasma being:

cs = Te + γTi mi . (4.74)

We note that the electrons essentially have γe = 1 while the ion contribution is the typical one and actually will reduce to the hydrodynamic result in that limit.

4.7 Electromagnetic Waves in Plasmas

Next we turn to the critical question of electromagnetic waves in plasmas. In vacuum, of course, this would be a typical light wave propagating at the speed of light, but in a plasma the E and B fields can induce plasma ‘sources’ of charge and current that feed back to the wave. The plasma response can be written in terms of the induced charge and current densities:

ρ = e 𝜖0(ni ne) (4.75) j = e(niui neue) (4.76)

and thus Maxwell’s equations are:

E = e 𝜖0(ni ne) (4.77) B = 0 (4.78) ×E = B t (4.79) ×B = μ0e(niui neue) + 1 c2 E t (4.80)

We start by taking the curl of Faraday’s Law:

× (×E) = t(×B), (4.81)

and substitute in Ampère’s Law to get the dispersion relation:

× (×E) = t μ0e(niui neue) + 1 c2 E t . (4.82)

Of course we know our vector identity:

× (×E) = (E) 2E. (4.83)

Unlike the vacuum scenario, E0. We use this result with our proto-dispersion-relation equation, collecting ‘vacuum’ terms on the left and ‘plasma’ terms on the right to get:

1 c2 2E t2 2E = (E) μ 0e t(niui neue). (4.84)

At this point it is useful to note that if we set ne = ni = 0 we recover the vacuum result ω2 = c2k2. If we assume transverse waves, then E = ik E = 0, and

1 c2 2E t2 2E = μ 0e t(niui neue). (4.85)

We know that light waves are very fast, so we can treat the ion inertia as infinite and consider the current as coming from the electron response only. In this case we have to solve the electron equation of motion, i.e. the momentum equation:

mene ue t + (ue )ue = nee(E + ue ×B) pe (4.86)

At this point we assume Te = 0 pe = 0, and take the case u0 = E0 = B0 = 0. Furthermore, assuming that the wave-induced perturbation is small and keeping only first-order terms, we get that:

men0u1 t = n0eE1. (4.87)

Using Fourier analysis:

iωmen0u1 = n0eE1, (4.88)

u1 = e iωmeE1. (4.89)

We can substitute this result into the dispersion relation to get:

1 c2 2E t2 2E = μ 0en0 e iωme E t . (4.90)

Now we have to use Fourier analysis on the whole thing, t iω and ik to get:

ω2 c2 E + k2E = μ0n0e2 me E (4.91)

We can now cancel the electric field, use 𝜖0μ0 = c2 and ωpe = ne2 me 𝜖0, and rearrange to get the dispersion relation for EM waves in plasmas:

ω2 = c2k2 + ω pe2. (4.92)

The most important thing to note from this dispersion relation is that if ω2 < ωpe2, then k is imaginary and the wave is evanescent. We call this a cutoff for the wave. It is also useful to calculate the group velocity:

vg ω k (4.93)

If we directly take the derivative k of the dispersion relation, we get that:

2ωvg = 2c2k, (4.94)


vg = c2kω(= nc) (4.95)

If we do some algebraic manipulation of the original dispersion relation,

c2k2 = ω2 1 ω pe2ω2 (4.96)


k ω = 1 c1 ωpe 2 ω2 (4.97)

using this we get a group velocity of:

vg = c1 ωpe 2 ω2 (4.98)

For propagating waves, we must have ω > ωpe which nicely gives us the result that vg < c. We also note that, since ωpe ne, the wave propagation depends on the plasma density in this case.

4.8 Electromagnetic Ion Waves

In this section we consider two electromagnetic ion waves: the Alfvén wave and the magnetosonic wave.

4.8.1 Alfvén Wave

This is a fundamental plasma wave. We wish to consider low-frequency ion oscillations in the presence of a magnetic field. We must consider a magnetized plasma with non-zero B0, which we take as B0 = B0. The geometry of this wave is k B0, E1,j1 B0, and v,B1 perpendicular to both B0 and j1. By convention we take E1 x̂ without loss of generality.

Starting with Maxwell:

× (×E) = k(k E1) + k2E 1 = ω2 c2 E1 + iω 𝜖0c2j1 (4.99)

Since k E1 = 0 by geometry, the only non-trivial equation is for the x direction:

𝜖0(ω2 c2k2)E 1 = iωn0e(uix uex) (4.100)

We have to calculate the ion response from the momentum equation.

mini ui t + (ui )ui = nie(E + ui ×B) p (4.101)

For simplicity we take the zero temperature limit, and use E = ϕ. We also neglect second-order terms. This gives us

miui1 t = eϕ1 + eui ×B0. (4.102)

Splitting into x and y components, and using Fourier analysis, we get

iωmiuix = ikeϕ1 + euiyB0 iωmiuiy = euixB0 (4.103)

Starting with the x equation, we can write that

ωmi ieB0 eB0 iωmi uix = keϕ1, (4.104)

uix = ke ωmiϕ1 1 ωci2 ω2 1. (4.105)

We can immediately use this result to get the form for y:

uiy = ieB0 ωmiuix = iωci ω uix, (4.106)

uiy = iωci ω ke ωmiϕ1 1 ωci2 ω2 1. (4.107)

In summary, then, the ion motion in the xy plane is given by

uix = ke ωmiϕ1 1 ωci2 ω2 1 (4.108) uiy = iωci ω ke ωmiϕ1 1 ωci2 ω2 1 (4.109)

We can immediately obtain the electron equations of motion by substituting mi me, ωci ωce, and taking the limit ωceω 1.

uex = ke ωmiϕ1 ω2 ωce2 0 (4.110) uey = iωce ω ke ωmiϕ1 ω2 ωce2 E1 B0 (4.111)

There is no electron motion in the x direction because the cyclotron motion overpowers it; however, the electrons have the usual E × B drift which is ŷ in this geometry. Since the current of interest is in the x̂ direction, it must be from the ion motion only.

Going back to the dispersion relation,

𝜖0(ω2 c2k2)E 1 = iωn0euix, (4.112)

𝜖0(ω2 c2k2)E 1 = iωn0e ke ωmiϕ1 1 ωci2 ω2 1. (4.113)

Since E1 = ikϕ1, we can eliminate them from the dispersion relation and simplify to:

(ω2 c2k2)E 1 = n0e2 𝜖0mi 1 ωci2 ω2 1. (4.114)

We recognize the plasma frequency, and simplify further to:

(ω2 c2k2) = ω pi2 1 ωci2 ω2 1. (4.115)

In the limit ω ωci, then:

ω2 c2k2 = ωpi2ω2 ωci2 = ω2 ρ 𝜖0B02 (4.116)

Rearranging, after some algebra one obtains the dispersion relation for Alfvén waves:

ω2 k2 = c2 1 + (ρμ0B02)c2 (4.117)

We recognize the Alfvén velocity:

vA = B0 2 ρμ0 . (4.118)

In the limit vA c then:

ω = vAk (4.119)

Basically, this wave is one where the lines of force (magnetic field lines) and the plasma move together in the plane perpendicular to the initial magnetic field.

4.8.2 Magnetosonic Wave

In the last section we considered waves along the initial magnetic field; we now consider low-frequency electromagnetic waves that propagate across B0. Again we take B0 = B0, E1 = E1x̂, but now k = kŷ to satisfy the above. At the beginning of this derivation, we note that E1 ×B0 k, implying that the E × B drifts will compress and expand the plasma as the wave propagates. This means that we have no choice but to keep the pressure gradient term in the ion momentum equation:

mini ui t + (ui )ui = nie(E + ui ×B) p (4.120)

However we let E0 = u0 = 0, keep only first-order terms, and use p = γTini,

miniui1 t = nieE1 + nieui1 ×B0 γTini (4.121)

employing Fourier analysis gets us to:

iωmin0ui1 = n0eE1 + n0eui1 ×B0 ikγTin1 (4.122)

ui1 = ie ωmi(E1 + ui1 ×B0) + kγTi miω n1 n0ŷ (4.123)

Next we need to split this into x and y components, which will make it tractable:

uix = ie ωmi(E1 + uiyB0) (4.124) uiy = ie ωmiuixB0 + kγTi miω n1 n0 (4.125)

We now need to use the continuity equation:

n t + (nu) = 0 (4.126)

iωn1 + ikn0uiy = 0 (4.127)

n1 n0 = k ωuiy (4.128)

We use this with the y equation of motion to get

uiy = ieB0 ωmi 1 k2γTi ω2mi 1u ix (4.129)

We can now plug this into the equation of motion for x:

uix = ie ωmi E1 ieB02 ωmi 1 k2γTi ω2mi 1u ix (4.130)

1 ωci2 ω2 1 k2γTi ω2mi 1 u ix = ie ωmiE1 (4.131)

uix = ie ωmiE1 1 ωci2 ω2 1 k2γiTi ω2mi 1 1 (4.132)

For the final current, we also need to know the electron velocity, which can be obtained by mi me and ωci ωce. We also take the limit ω ωce in which case the term in square brackets can be simplified by the binomial approximation, and we get:

uex = ie ωmeE1 ω2 ωce2 1 k2γeTe ω2me ie ωme ω2 ωce2 k2γeTe ω2me E1 (4.133)

Similarly to the last section, we know that the dispersion relation here will be:

𝜖0(ω2 c2k2)E x = iωn0e(uix uex) (4.134)

using the fluid velocities we just derived,

𝜖0(ω2 c2k2) = n0e2 me mi me 1 ωci2 ω2 1 k2γiTi ω2mi 1 1 k2 ωce2 γeTe me (4.135)

the plasma frequency,

ω2 c2k2 = ω pe2 mi me 1 ωci2 ω2 1 k2γiTi ω2mi 1 1 k2 ωce2 γeTe me (4.136)

taking the limit ω ωci,

ω2 c2k2 = ω pe2 mi me ω2 ωci2 1 k2γiTi ω2mi k2 ωce2 γeTe me (4.137)

After a bit of algebra, we can rearrange this to be

ω2 c2k2 1 + γeTe mivA2 + ωpi2 ωci2 ω2 k2γiTi mi = 0 (4.138)

where we have used the definition of the Alfvén speed, vA B0μ0 ρ = cωciωpi. From the derivation of the ion acoustic wave, recall the sound speed is:

cs = γe Te + γi Ti mi (4.139)

so the dispersion relation becomes, after some algebraic hammering,

ω2 1 + c2 vA2 = c2k2 1 + cs2 vA2 (4.140)

or after rearranging, the magnetosonic wave dispersion relation:

ω2 k2 = c2vA2 + cs2 vA2 + c2 . (4.141)

This wave is essentially an acoustic (sound) wave but the compression/expansion is created by E × B drifts. In the limit B0 0, this wave becomes the ion acoustic wave. In the limit T 0 the pressure gradient forces drop out, cs 0, and this wave becomes a modified Alfvén wave.

Because the phase velocity of this wave is almost always higher than vA, it is sometimes called the ‘fast’ hydromagnetic wave.

4.9 Brief Summary of More Waves

We give a brief summary of the rest of the menagerie of plasma waves here, with dispersion relations given from general plasma physics references (e.g. Chen).

4.9.1 Upper Hybrid Waves

Next we consider electron plasma waves. It turns out that the electron plasma wave along B0 is unaffected by the magnetic field. The perpendicular case, however, is the upper hybrid wave. The easiest situation to consider is the zero temperature limit. If we consider longitudinal waves, E1 k and k B0 we can immediately see that there will be an E × B drift for the electron motion in the wave. This affects the electron equation of motion. At the end of the derivation, we would arrive at the upper hybrid frequency:

ωuh2 = ω pe2 + ω ce2 (4.142)

Once again, we immediately see that the B0 0 limit corresponds to the usual electron plasma oscillation.

4.9.2 Ion Cyclotron Waves

In the derivation of the ion acoustic wave we assumed unmagnetized plasma. Instead we consider the ion acoustic wave in a magnetized plasma. It is tempting to set k B0 = 0, but the problem for this situation is that the electrons are unable to move between wave fronts because of their small gyroradii. Instead, if k and B0 are almost perpendicular but not quite, then this wave can propagate. It turns out that the critical angle is χ me mi which is indeed small.

When one goes through the derivation, analogously to the ion acoustic wave, one arrives at the electrostatic ion cyclotron wave dispersion relation:

ω2 = ω ci2 + k2c s2 (4.143)

we can see that this reduces to the ion acoustic wave in the limit B0 0 (equivalently ωci 0).

4.9.3 Lower Hybrid Waves

So what happens to the ion acoustic wave when the propagation angle is exactly π2, i.e. the wave is exactly perpendicular to the magnetic field? In this case, keeping finite electron mass, it turns out that the compression/expansion of the normal ion acoustic wave is unimportant because the electron motion is constrained. Starting with the electron and ion equations of motion, it turns out that there is an oscillation at the lower hybrid frequency:

ωlh = ωci ωce (4.144)

4.9.4 O/X Waves

We now move on to considering the electromagnetic wave in magnetized plasmas, first for the waves with k B0. First, if E1 B0 then there is no change to the plasma response, and we have the Ordinary (‘O’) Wave:

ω2 = ω pe2 + c2k2 (4.145)

which is the same as the unmagnetized result.

On the other hand, if E1 B0 then there is a drift which alters the plasma response. Generally the wave’s electric field is taken as elliptically polarized in the plane perpendicular to B0. Working through the math is somewhat involved, but at the end of the day the dispersion relation is

c2k2 ω2 = 1 ωpe2 ω2 ω2 ωpe2 ω2 ωuh2 (4.146)

which includes the previously-encountered upper-hybrid frequency.

4.9.5 R/L Waves

The last set of waves are the electromagnetic waves with k B0. In general the wave’s electric field can lie anywhere in the plane perpendicular to the initial magnetic field. This results in two solutions, which are elliptically polarized waves with right-hand or left-hand polarization, and which are respectively the R wave

n2 = c2k2 ω = 1 ωpe2ω2 1 ωceω (4.147)

and the L wave

n2 = c2k2 ω = 1 ωpe2ω2 1 + ωceω (4.148)

Chapter 5

5.1 Diffusion in weakly-ionized gases

If we consider a weakly-ionized gas, then the diffusion of plasma species is dominated by collisions with the neutral atoms. This is the simplest diffusion case, so we start with it. We begin with the momentum equation (for an arbitrary plasma species):

mn u t + (u )u = nq(E + u ×B) p mnν(u un) (5.1)

if we consider a steady-state equilibrium system, then the term ut is zero. We take u0 = B0 = 0. If u is small or ν is large then we can likewise neglect the convective derivative term as a small term, in which case,

v = u un = q mνE T mν n n (5.2)

where we have also used p = nT. We can simplify this expression by defining the mobility:

μ |q| mν (5.3)

and diffusion coefficients:

D T mν (5.4)

so that the total flux of a plasma species j is a generalized Fick’s law:

Γj = ±μjnjE Djnj (5.5)

of course if q = 0 or E = 0 then this expression reduces to the typical Fick’s law from classical physics.

5.2 Ambipolar Diffusion

To start off with, we rewrite the continuity equation in terms of the particle flux:

n t + Γ = 0 (5.6)

It is clear that the particle flux can cause time-varying concentrations of plasma species. If we consider a hydrogen ion - electron plasma, with just two species, the immediate question is if the diffusion process violates quasi-neutrality. We note that both the mobility and diffusion coefficients depend inversely on the mass, and so naïvely one might ask what stops the electrons from all diffusing away and leaving an ion-only plasma. The answer is that an imposed positive potential on the high-density regions of plasma creates an electric field which enhances ion diffusion and suppresses electron diffusion. The plasma will quickly charge to essentially force Γe = Γi. This process is known as ambipolar diffusion.

We can derive the spontaneously-generated ambipolar electric field by setting the diffusion rates for electrons and ions equal to each other:

μinEa Din = μenEa Den (5.7)

where we have taken ne = ni = n. We can solve this equation for the ambipolar electric field:

Ea = Di De μi + μe n n (5.8)

with definitions for the (electron and ion) mobility and diffusion coefficients given in the previous section.

We can also derive the common resulting diffusion as:

Γ = μinEa Din (5.9)

Γ = Di De μi + μe μiDi n (5.10)

Γ = μiDe + μeDi μi + μe n (5.11)

which is simply Fick’s Law with a new ambipolar diffusion coefficient:

Da = μiDe + μeDi μi + μe (5.12)

5.3 Diffusion in a Slab

We consider the simple application to diffusion in a slab geometry. If the ambipolar diffusion coefficient is constant, then the continuity equation becomes:

n t = Da2n. (5.13)

We can solve this via a separation of variables technique. Take solutions of the form:

n(r,t) = T(t)X(x) (5.14)

plugging this into the continuity equation, we get that:

XdT dt = TDad2X dx2 (5.15)

which we can rearrange to get:

1 T dT dt = Da X d2X dx2 . (5.16)

We note that the left-hand side depends only on t and the right-hand side depends only on x, so by the separation of variables technique they must both be equal to a common constant. We can therefore solve them individually. Starting with the time,

1 T dT dt = const = 1 τ (5.17)

where we have set the constant equal to 1τ for reasons that will be apparent later. This equation can be simply solved:

dT dt = T τ (5.18)

T(t) = T0etτ (5.19)

The plasma decays exponentially in time, with time constant τ.

Next, we consider the solution in the spatial (x) direction. The equation to solve is:

Da X d2X dx2 = 1 τ (5.20)

d2X dx2 = 1 τDaX (5.21)

which has sinusoidal solutions:

X(x) = Asin xτDa + Bcos xτDa (5.22)

We require that, for a plasma extending to x = ±L, the solution must satisfy X(±L) = 0 so we can immediately set A = 0 and B = X0. We also now have the requirement that:

L τDa = π 2 (5.23)


τ = 2L π 2 1 Da (5.24)

is the timescale for the plasma to decay by ambipolar diffusion. Combining these results, the solution for the plasma density as a function of time and space is:

n(x,t) = n0etτ cos πx 2L (5.25)

5.4 Diffusion in Magnetic Fields

We revisit the problem of diffusion in weakly-ionized plasmas but now introduce a magnetic field, which will have the effect of reducing the diffusion. Of course, this is the main goal of magnetic confinement fusion!

In the direction parallel to B, there is no effect on the diffusion rates. If we take B , then:

Γz = ±μnEz Dn z (5.26)

For the perpendicular direction, we have to return to the fluid momentum equation:

mndv dt = ±en(E + v×B) p mnνv (5.27)

Assuming that the collisions are fast enough that the velocity derivative term (left hand side) can be neglected, and also assuming an isothermal plasma, the two components are:

mnνvx = ±enEx ± envyB0 T dn dx (5.28) mnνvy = ±enEy envxB0 T dn dy (5.29)

which can be simplified immediately to:

vx = ±μEx ±ωc ν vy D n dn dx (5.30) vy = ±μEy ωc ν vx D n dn dy (5.31)

We can solve these for vy via substitution of vx. This gives us:

(1 + ωc2τ2)v y = ±μEy D n dn dy ωc2τ2Ex B0 ± ωc2τ2 T eBn dn dx (5.32)

and similarly for vx:

(1 + ωc2τ2)v x = ±μEx D n dn dx + ωc2τ2Ey B0 ωc2τ2 T eBn dn dy (5.33)

We notice that the last two terms are the E × B drift and the diamagnetic drift, respectively, which are important generally but not for this problem. The first two terms are the standard mobility and diffusion terms, but because of the factor (1 + ωc2τ2) on the left-hand side the coefficients are reduced from the no-field case:

μ = μ 1 + ωc2τ2 (5.34)

D = D 1 + ωc2τ2 (5.35)

5.5 Diffusion in Fully Ionized Plasmas

We give a brief qualitative discussion of the difference between diffusion in partially- and fully-ionized plasmas. The previous sections dealt with diffusion in partially-ionized plasmas, in which both electrons and ions diffused via collisions with the background neutral particles, and we did not care about the concentrations of neutrals.

In the case of a fully-ionized plasma, the collisions which matter are ion-ion, electron-electron, or ion-electron / electron-ion. These can be lumped into like-particle and unlike-particle collisions.

In the case of a Coulomb collision between two identical (like) particles, it turns out that the individual particle velocities can be changed quite significantly in the collision - i.e. they are reversed in a 180 collision, and a 90 collision changes the velocity vector direction by 90. However, in a magnetized plasma, the important quantity is actually the particle’s guiding center. It turns out that the ‘center of mass of the guiding centers’ does not change in like-particle collisions. This means that like-particle collisions cannot lead to diffusion.

All diffusion in fully-ionized plasmas therefore comes from unlike-particle collisions. In an electron colliding with an ion, for instance, there can be net momentum exchange between the species which leads to diffusion. The ion itself is much less affected because of its large mass, of course, but at the end of the day momentum conservation requires that both species will diffuse based on these unlike-particle collisions.

5.6 Plasma Resistivity

We present a simple and intuitive derivation of the plasma resistivity, which is of course an important quantity. From the momentum equations, we know that the rate of momentum transfer between electrons and ions is:

Pie = Pei = mnνei(ui ue) (5.36)

We can also write this another way as follows. The plasma resistivity for two interpenetrating fluids with relative flow will be proportional to the relative velocity, (ui ue), the scattering species density (ni), the electron density ne, and the scattering force which is proportional to the two charges (e2). There is also going to be a constant of proportionality, so in total:

Pie = ηe2n2(u i ue) (5.37)

setting these two expressions for the momentum exchange equal to each other,

mnνei(ui ue) = η2e2n2(u i ue), (5.38)

which reduces to:

η = m ne2νei. (5.39)

The constant of proportionality η is the plasma resistivity. We now need to come up with an expression for the electron-ion collision rate.

The simple derivation is as follows. Consider a 90 collision between electron and ion. We can define this as the case where the initial kinetic energy is equal to half the potential energy at the initial impact parameter:

1 2mev2 = 1 2 e2 4π𝜖0b90 (5.40)

There is some degree of arbitrariness in this definition by about a factor of 2. More rigorous derivations are done but are much lengthier. Anyways, we can write the 90 impact parameter as:

b90 = e2 4πme𝜖0v2. (5.41)

A simple expression for the scattering cross section is thus:

σ = πb902 = e4 16π2me2𝜖02v4 (5.42)

and the electron-ion collision rate can be written:

νei = nσv = ne4 16π2me2𝜖02v3. (5.43)

It turns out that in many real plasmas, small-angle collisions actually dominate Coulomb interactions, or are at the very least important and cannot be neglected. We will see this more rigorously later. The relative importance of small- and large-angle scattering is characterized by the Coulomb logarithm, and the electron-ion collision rate is enhanced by lnΛ:

νei = nσv = ne4 lnΛ 16π2me2𝜖02v3. (5.44)

We can now write down an expression for the plasma resistivity:

η = m ne2νei = e2 lnΛ 16π2me𝜖02v3 (5.45)

Further, if the electrons have a Maxwellian distribution then v = Te me and this becomes:

η = mee2 lnΛ 16π2𝜖02Te32 lnΛ Te32. (5.46)

Since the Coulomb logarithm is very insensitive to plasma conditions, we can see that the plasma resistivity is mostly dependent on the electron temperature. We also note that the resistivity decreases (conductivity increases) at high temperature, which places limits on processes like Ohmic heating.

5.7 Coulomb Collisions

We must treat the collision of two charged particles; the classical Coulomb collision, in which the two particles interact only through the electrostatic potential:

E = ϕ, (5.47)

ϕ = q1q2 4π𝜖0r. (5.48)

In the lab frame we can define the kinematics by the positions and velocities of the two individual particles:

m1dv1 dt = q1q2 4π𝜖0 r2 r1 |r2 r1|3 (5.49) m2dv2 dt = q1q2 4π𝜖0 r1 r2 |r1 r2|3 (5.50) dr1 dt = v1 (5.51) dr2 dt = v2 (5.52)

Of course each of the previous four equations is in three dimensions, generally, and we must simplify from 12 equations before this is a tractable problem.

The first is to notice that conservation of momentum:

d dt m1v1 + m2v2 = 0, (5.53)

suggests a new coordinate system moving with the center of mass velocity:

V = m1v1 + m2v2 m1 + m2 (5.54)

which is a constant of the motion, plus the relative velocity:

v = v1 v2. (5.55)

Of course we can transform back to the original lab frame coordinates via:

v1 = V + m2 m1 + m2v, (5.56) v2 = V + m1 m1 + m2v. (5.57)

The motion in the center-of-mass frame is simply the relative motion:

dv dt = q1q2 4π𝜖0mr r r3 (5.58)

where we have introduced the reduced mass:

mr = m1m2 m1 + m2. (5.59)

The problem is now equivalent to a single particle of mass mr moving around a central potential at the center of mass. At this point we will need to introduce additional constraints via conserved quantities, in particular energy and angular momentum. One can derive these rigorously, but for simplicity:

1 2mrv2 + q1q2 4π𝜖0 1 r = E0 = const (5.60)

mrr ×v = L0 = const (5.61)

One nice result of conservation of angular momentum is that we can see that the motion is entirely in a plane, therefore we can simplify from three to two dimensions for the remainder of the derivation.

We can write the conserved quantities in terms of the original velocity v0 and impact parameter b (relative to the center of mass).

E0 = 1 2mrv02 (5.62)

L0 = mrbv0 (5.63)

Finally, we use cylindrical coordinates for this problem:

x = rcos𝜃,y = rsin𝜃 (5.64)

and we remember that:

v = r̂ + r𝜃̇𝜃̂ (5.65)

so that the conservation of energy and momentum, respectively, give us that:

2 + r2𝜃̇2 + q1q2 2π𝜖0mrr = v02 (5.66)


r2𝜃̇ = bv 0 (5.67)

the second equation can be used to eliminate 𝜃 in the first to get:

1 v022 = 1 b2 r2 2b90 r (5.68)

where we have defined the 90 impact parameter as:

b90 = q1q2 4π𝜖0mrv02 (5.69)

The next task is to solve for the scattering angle χ as a function of the impact parameter and relative velocity. If we define the point of closest approach by an angle 𝜃min and distance rmin we can see via geometric arguments that

χ = 2𝜃min π (5.70)

We will solve for 𝜃min using our previously-derived equation of motion, starting with:

𝜃min π 𝜃minπd𝜃 = π rmind𝜃 drdr (5.71)

with a change of variables from 𝜃 to r in the last step. Next we need to use:

d𝜃 dr = 𝜃̇ = bv0 r2 (5.72)

we can plug our previous relation for into this equation and then into the expression for 𝜃min to get that:

𝜃min = π brmin dr rr2 b2 2b90 r (5.73)

Actually evaluating this integral requires a non-obvious trick substitution (see, e.g. Friedberg):

b r = b90 b + 1 + b90 2 b2 sinα (5.74)

𝜃min = π α1α2 dα = π (α2 α1) (5.75)

with α2 = π2 and

sinα1 = b90 b2 + b90 2 (5.76)

What we end up with, omitting some algebra,

cot χ 2 = b b90 = 4π𝜖0mr q1q2 v02b (5.77)

5.8 Coulomb Collision Frequencies

Now that we have the general result in the previous section for Coulomb collisions, the immediate question is how often does a particle in a plasma, a ‘test particle’, collide with other particles. In this derivation, which follows along with Friedberg Section 9.3, we explicitly consider test electrons but the results can be applied generally.

First, we can write simply the change in the test electron’s momentum as:

d dt(meve) = (Δmeve)niσve = νei(meve) (5.78)

where in the last step we have essentially defined the electron-ion collision frequency νei as a rate of momentum loss. In general the cross section is not a constant and we in fact must rewrite this as:

νei = 1 meve(Δmeve)fi(vi)|ve vi|bdbdαdvi (5.79)

where the momentum loss must be integrated over all possible impact parameters b, scattering angles α in the plane perpendicular to the scattering, and also over the scattering body distribution, in this case the ions.

The first task is to evaluate Δmeve and then we must tackle the integral. If the initial velocity is in the x direction in the center-of-mass frame, then we can write the change in momentum as:

Δmeve = mr(vi vf) x̂ (5.80)

In terms of the scattering angle χ used previously we can see by intuition that the change in x-directed momentum will be

Δmeve = mr(ve vx)(1 cosχ) (5.81)

and based on what we previously derived for Coulomb scattering:

Δmeve = mr(ve vx) × b902 b2 + b902 (5.82)

Now we can tackle the integrations over impact parameter and α:

0αdα0bmax b902 b2 + b902bdb = πb902 ln 1 + bmax2 b902 (5.83)

There is an important point to be made here. In a naïve sense one might expect to take the integral over the impact parameter b to infinity since the Coulomb interaction has indefinite range, and we have previously argued that long-range interactions are fundamental aspects of a plasma. However, in this integration taking the b integration to infinity would lead to a logarithmic divergence, which cannot be allowed.

We can make a simple physical argument that as the impact parameter is taken to the large limit, it cannot exceed the Debye length because at that point the plasma will screen the Coulomb interaction. We therefore typically take bmax λDe. We typically use the definition Λ = bmaxb90 such that

πb902 ln 1 + bmax2 b902 2πb902 lnΛ (5.84)

In magnetically-confined plasmas we typically have lnΛ 15 20, and this quantity is called the Coulomb logarithm. It is essentially a metric of the relative importance of small- and large-angle collisions. In the regime lnΛ > 10 the behavior is dominated by small-angle collisions, the classical plasma regime. If lnΛ < 10 then we are in a moderately- or strongly-coupled plasma, which are more applicable to ICF plasmas.

We now use this in the expression for the collision frequency:

νei = 4πmr meve lnΛ(ve vx)b 902f i(vi)|ve vi|dvi (5.85)

we now have to tackle the integration over target velocities, which is best done in spherical coordinates:

vx = v e + vcos𝜃 (5.86) vy = vsin𝜃sinϕ (5.87) vz = vsin𝜃cosϕ (5.88)

in which case we can write that:

dvi = v2 sin𝜃dvd𝜃dϕ (5.89)

|ve vi| = v (5.90)

b902 = e2 4π𝜖0mr 2 1 v4 (5.91)

and returning to the expression for the collision frequency, some algebra reveals that it reduces to

νei = 2e4 8π32𝜖02 ni Ti32 mi32 memr lnΛI(we) (5.92)

with velocities normalized to the ion thermal value: we = vevTi and wi = vvTi. We have also introduced a dimensionless integral:

I(we) = ewe2 we 0dw0πd𝜃sin𝜃cos𝜃ew22wew cos 𝜃 (5.93)

This can be evaluated in terms of the error function, and then taking appropriate limits. We skip the details and go straight to the result:

I(we) = 1 we2 π 2weerf(we) ewe2 π 2 1 we3 + 3π4 (5.94)

Using this to complete the expression for the collision frequency, we get:

νei = e4 4π𝜖02 ni Ti32 1 memr lnΛ 1 ve3 + 1.3vTi3 (5.95)

For electron-ion collisions in particular we can take the limits mr me and ve vTe vTi, in which case:

νei = e4 4π𝜖02me2 ni Ti32 lnΛ 1 ve3 (5.96)

One can also examine the expression for the collision frequency and come up with a hierarchy of various species colliding with each other, in terms of the mass ratio μ = memi:

νee νei (5.97) νii μν ei (5.98) νie μνei (5.99)

5.9 Distribution Function

In the previous sections we have often explicitly or implicitly assumed that the particles obey a Maxwellian distribution, but in many plasmas of interest this is not the case and leads to full-fledged kinetic theory. First we discuss some properties of the distribution function.

In general the particle distribution can be a function of full phase space: position r, velocity v, and time t:

f = f(r,v,t) (5.100)

One can obtain functions such as the density via integration over part of phase space:

n(r,t) =f(r,v,t)dv (5.101)

In this particular case we sometimes use a normalized distribution function:

f̂ = fn (5.102)

Of course the most famous and commonly-used distribution function is the Maxwellian:

fm ̂ = m 2πkBT32ev2vth2 (5.103)

with vth = 2kB Tm and v = |v|.

5.10 Overview of Kinetic Theories

The basic problem of kinetic theory is to define how the distribution function of a plasma changes with time. One can simply derive this by taking the chain rule total time derivative of f:

df dt = f t + f x dx dt + f y dy dt + f z dz dt + f vx dvx dt + f vy dvy dt + f vz dvz dz (5.104)

The first term is simply the explicit time dependence of f. The second through fourth terms are observed to be:

v f (5.105)

and the fifth through seventh terms are, using F = mv̇, simply

F m f v (5.106)

Overall the total change in time of the distribution function is due to external interactions, which we call ‘collisions’ and now we can write the Boltzmann Equation:

f t + v f + F m f v = f t c. (5.107)

Of course in plasmas one often has the situation that the forces present are simply due to the particle’s Coulomb interaction, and furthermore, collisions can often be neglected in hot plasmas. In this limit we obtain the Vlasov Equation:

f t + v f + q m(E + v ×B) f v = 0. (5.108)

If there are collisions with neutral atoms, then one can use the Krook collision term:

f t c = fn f τ (5.109)

And another familiar limit is with Coulomb collisions, we get the Fokker-Planck Equation:

df dt = v (fΔv)1 2 2 vv : (fΔvΔv) (5.110)

which is too complex to delve into at this level.

5.11 Landau Damping

One of the most important and famous results of kinetic theory is Landau damping, or more generally wave-plasma collisionless coupling. In this case we consider electron plasma oscillations in an initially uniform plasma with E0 = B0 = 0. We also want to treat the wave perturbatively, i.e. let:

f(r,v,t) = f0(v) + f1(r,v,t) (5.111)

the first-order Vlasov equation for electrons is:

f1 t + v f1 e mE1 f0 v = 0 (5.112)

we will assume infinite inertia ions and use plane wave perturbations in the x direction, so that:

f1 ei(kxωt) (5.113)

so the Vlasov equation becomes, to first order,

iωf1 + ikvxf1 = e mEx f vx (5.114)

f1 = ieEx m f0vx ω kvx (5.115)

Now we have to use Poisson’s Equation:

𝜖0E1 = ik𝜖0Ex = in1 = ef1d3v (5.116)

Using these two results we can eliminate Ex to get:

1 = e2 km𝜖0f0vx ω kvx d3v (5.117)

If f is a Maxwellian, then we can easily separate the three coordinates and reduce this to a one-dimensional integral:

1 = wpe2 k2 f̂0(vx)vx vx (ωk) dvx (5.118)

We notice immediately the pole at vx = ωk. While the pole is not necessarily on the contour of integration if ωk is complex, it matters nevertheless. In the simplest limit of large phase velocity and weak damping, the pole lies near the real axis and we can take the contour of integration as along the axis except for a small semi-circle around the pole. In this case, by the residue theorem,

1 = ωpe2 k2 Pf̂0v v (ωk)dv + iπ f̂0 v v=ωk (5.119)

The first term in square brackets ends up being the electron plasma wave dispersion relation:

ω2 = ω pe2 + 3T 2mk2 (5.120)

see Chen section 7.4 for a derivation of this. The interesting part is the second term, which serves as a small and imaginary correction to the dispersion relation in the limit of small k:

ω2 1 iπωpe2 k2 f̂0 v v=vϕ = ωpe2 (5.121)

Or rearranging, and using again the smallness of the second term in parentheses,

ω = ωpe 1 + iπ 2 ωpe2 k2 f̂0 v v=vϕ (5.122)

Essentially, we can see that this derivation has introduced damping to the electron plasma wave. The damping depends on the wavenumber, and most critically on the shape of the distribution function where the particle velocity matches the wave’s phase velocity. At this point there is coupling predicted between the wave and the plasma particles. If (f0v) < 0 then we can see that the wave is damped and transfers energy to the particles. On the other hand, if (f0v) > 0 the the wave is actually unstable and grows in amplitude.

Chapter 6

6.1 Refractive measurements of density

One can use the refractive index of a plasma, or dispersion, as a density measurement via interferometry. Consider the following scheme, shown in Fig. 6.1, where we stick an unmagnetized plasma in one leg of an interferometer (Mach-Zender scheme shown). In this case, the detectors will measure the phase shift between the two legs, so we must calculate the phase shift due to plasma.

Consider the source S as a monochromatic source of electromagnetic radiation that obeys:

ω ωpe,f fpe (6.1)

Remember the back-of-the-envelope rule fpe = 9n Hz, with [n] = m3.


Figure 6.1: Density measurement.

The dispersion relation for electromagnetic waves in a plasma is

ω2 = ω pe2 + c2k2, (6.2)

ck = ω1 ωpe 2 ω2 (6.3)

which leads to a refractive index

n ck ω = (1 + ωpe2ω2)12 (6.4)

We know that we can write the phase shift for propagating waves as

ϕ =kd =nω c d (6.5)

So the difference between plasma and vacuum (with n = 1) is:

Δϕ =0L(n 1)ω c d, (6.6)

Δϕ = 0L (1 ω pe2ω2)12 1 ω c d (6.7)

Based on our initial assumption that ω ωpe we can use the binomial expansion:

Δϕ =0L1 2 ωpe2 ω2 ω c d = 1 2ωc0Lω pe2d (6.8)

Using the definition of the plasma frequency:

ωpe2 = nee2 𝜖0me (6.9)

we can write the phase shift as:

Δϕ = 1 2ωc e2 𝜖0me0Ln ed (6.10)

The phase shift is directly proportional to the line-integrated electron number density.

6.2 Refractive Measurements of Astrophysical Length

Now consider a somewhat different system, where we have a broadband source of radiation S that traverses a plasma of scale length L and known density ne to a detector. Given some information, for example the difference in arrival time between radiation of frequencies f1 and f2, what is the length L?


Figure 6.2: Length measurement.

We start off by stating:

L = vgt (6.11)

so for two different frequencies:

Δt = L(vg11 v g21) (6.12)

The electromagnetic radiation propagates at the group velocity:

vg ω k (6.13)

so we need the dispersion relation. We know that for electromagnetic radiation in an unmagnetized plasma:

ω2 = c2k2 + ω pe2 (6.14)

taking the derivative k gives us:

2ωvg = 2c2k, (6.15)

vg = c2k ω . (6.16)

Using the original dispersion relation to substitute for ω,

vg = c2k ck 1 + ωpe2 c2k2 12 (6.17)

We typically take the limit ωpe ω, which can be checked via the back-of-the-envelope relation fpe = 9n with f in Hz and [n] = m3. In that limit we can binomial approximate the last term,

vg = c 1 ωpe2 2c2k2 c 1 1 2 fpe2 f2 (6.18)

returning to the original expression between L and Δt, we can rewrite as:

L = Δt(vg11 v g21)1 (6.19)

since the terms fpef are small, when we substitute in for the group velocities,

L = cΔt fpe2 2f22 fpe2 2f12 1 (6.20)

at which point one must take some values for f1, f2, Δt, and ne and plug them into this expression.

6.3 Faraday Rotation

Next we consider using electromagnetic radiation as a probe of magnetic field conditions. This topic is similar to the first considered in Diagnostics, in that we will calculate a phase shift due to the plasma and consider that it is being compared to the vacuum value in an interferometer.

We consider the electromagnetic source of radiation as producing monochromatic plane-polarized waves, which propagate in the plasma such that k B0 . In this case the two fundamental solutions are the R and L waves, or right- and left-hand circularly polarized waves. The initial wave can be written as a superposition of these two, following Hutchinson’s notation:

E(z = 0) = E0 2 (1,i) + (1,i) (6.21)

where we are denoting the polarizations by ± i from the basic:

Ex Ey = ±i (6.22)

for circular polarization. Also note that we have taken E(z = 0) x̂ without loss of generality. In general, some distance later, the electric field will be given by:

E = E0 2 eik+z + eikz (6.23)

where the two circular polarizations can have different wave numbers. Using k = nωc,

E = E0 2 ein+ωzc + einωzc (6.24)

which we can rewrite using our (x,y) notation as:

E(z) = E0 exp[i(n+ + n)ωz2c] cos Δϕ 2 ,sin Δϕ 2 (6.25)


Δϕ = (n+ n)ωz c (6.26)

there is a phase difference between the RHCP and LHCP components, which leads to a rotation in the plane polarization as the wave propagates along z.

In general, for a wave propagating in an plasma with small 𝜃 between k and B0, the dispersion relation can be written as:

n2 = 1 X ± XY cos𝜃 (6.27)


X = ωpe2 ω2 ,Y = ωce ω (6.28)

the trick here is to write

n±2 1 X = 1 ±XY cos𝜃 1 X (6.29)

taking the square root, and approximating the second term as small (X is small):

n± 1 X = 1 ±XY cos𝜃 2(1 X) (6.30)

taking the difference between the two polarizations:

(n+ n) = XY cos𝜃 1 X . (6.31)

We now plug in for X and Y :

Δϕ = ωz c ωpe2ωce cos𝜃 ω31 ωpe 2 ω2 (6.32)

in the limit of small X, i.e. ω ωpe:

Δϕ = ωz c ωpe2ωce cos𝜃 ω3 (6.33)

we note that the polarization plane rotation is z. The Faraday rotation angle is typically defined as

α Δϕ 2 = ωz 2c ωpe2ωce ω3 cos𝜃 (6.34)

The Faraday rotation angle is linearly proportional to the propagation distance z, the plasma density (through plasma frequency), and the initial magnetic field strength through the cyclotron frequency.

In the common event that the plasma is not perfectly uniform, a WKB-style analysis can be adopted:

α =0Lϕ z ω c d (6.35)

some algebraic manipulation and our previous definition of Δϕ leads to:

α = e 2mc0Lne ncB d (6.36)

the Faraday rotation depends on both density and magnetic field profiles. To measure one, the other must be known ab initio.

6.4 Magnetic Field Probes

6.4.1 Probes


Figure 6.3: probe.

One of the most common and simplest probes for magnetic field is just a coil. The schematic is shown in Fig. 6.3. The magnetic field is sampled via flux through a coil, which induces an EMF and voltage. We know from Faraday’s law that

E d =SḂ dA (6.37)


V = NAḂ (6.38)

where V is the induced voltage, N is the number of coils, and A is the area of one coil. Typically these probes are used with an circuit integrator, in which case the post-integrator voltage is:

V = NAB RC (6.39)

6.4.2 Hall Probes


Figure 6.4: Hall probe.

The Hall probe is shown schematically in Fig. 6.4. The concept is that a current, j, is passed through the probe. There is a j ×B force on the current carriers. Since this is typically electrons, there is an induced charge separation and electric field as shown in the figure. A Hall probe utilizes this effect by measuring the potential induced across the probe, which can be related to the magnetic field and current. In most plasma experiments, however, there is too much pickup for this to be a useful technique.

6.4.3 Rogowski Coils


Figure 6.5: Rogowski coil.

The last simple diagnostic considered in this section is actually a current diagnostic, but the concept is very similar to the probe. Consider the schematic Fig. 6.5. In this case the coil surfaces are perpendicular to the magnetic field induced by the current I. By

B d = μ0Ienc, (6.40)

the induced flux on the probe is

Φ = NAμ0I (6.41)

where N is the number of coils and A is the area of each coil. We can thus write the voltage induced as

V = Φ̇ = NAμ0İ. (6.42)

Once again we typically would pair a Rogowski coil with an integrating circuit, in which case the final signal is:

V = NAμ0I RC (6.43)

6.5 Langmuir Probes

The last diagnostic considered is the important and basic Langmuir probe. Consider an electrically conducting probe placed into the plasma at arbitrary potential. Fig. 6.6 shows a schematic of the probe potential. A sheath develops around the probe; over the sheath the probe’s potential is Debye screened. In this example the probe potential V is less than the plasma potential.


Figure 6.6: Langmuir probe sheath potential.

The first important case is when the probe potential is much less than the plasma. This could arise, for example, if the probe is grounded and the plasma is positively charged due to ambipolar diffusion. If the difference in potentials is greater than a few times kTee, then the thermal electrons do not have enough energy to overcome the potential and the probe collects only ions. This is called the ion saturation current. We can estimate this as:

ISI = nsevA (6.44)

where ns is the sheath density, v is the ion drift velocity, and A is the probe area. If we take the drift velocity as

v kB Tmi (6.45)

and the sheath density as given by Debye screening:

ns = n0eeϕkBTe n0e12 = 0.6n 0 (6.46)

where we have approximated the sheath boundary as eϕ kBTe2. In this case the ion saturation current is given by:

ISI = 0.6n0eAkB Te mi (6.47)

By measuring the ion saturation current, one can determine the plasma density if the electron temperature is known.


Figure 6.7: Langmuir probe I V curve.

Of course one must also have a way to measure the electron temperature. This can actually be done with the same probe! Consider qualitatively what happens as the probe potential is increased. When it reaches the plasma potential there is no potential difference to drive a current, so the current must go to zero. As the potential is further increased, the probe potential exceeds the plasma potential. Now ions are repulsed by the probe potential and electrons are preferentially collected instead. The rate at which this transition occurs depends on the temperature, since the probe is essentially sampling an increasing fraction of the electron distribution as the voltage is increased. The derivation is omitted here, but the qualitative behavior of the I V curve is shown in Fig. 6.7.

Chapter 7

7.1 Introduction

In this section we attempt to give a brief motivation of the basic design of a tokamak and its general properties. In previous sections we have considered only linear machines. If one wishes to build a fusion reactor, then end losses in a linear pinch (Z or 𝜃) are too great to overcome. The magnetic mirror was proposed as a way around this limit, yet it turns out that the mirror scheme is not enough to overcome the end losses and the absolute theoretical maximum Q for a mirror machine is 1.1. For magnetic confinement, we must design a scheme with no end losses.

In the 1950s Soviet physicists Sakharov and Tamm proposed a scheme in which a linear pinch plasma is bent into a torus, thus eliminating the end losses. The basic geometry is shown in Fig. 7.1; here we call R the major radius of the torus and a the minor radius. If one imagines using coils around the torus to create a toroidal magnetic field, the confinement scheme will be very similar to the 𝜃 pinch MHD equilibrium.


Figure 7.1: Tokamak geometry.

The basic field geometry thus far is shown in Fig. 7.2. In the top half of the figure we show the toroidal current-carrying coils which create a toroidal field B. This field configuration is equivalent to a single infinite current-carrying wire on the origin for this simple analysis. We know that the magnetic field due to an infinite wire can be derived from Ampère’s Law:

B d = μ0Ienc (7.1)

|B| = μ0I 2πr (7.2)

The important point to take away from this is that |B| 1r.


Figure 7.2: Tokamak with toroidal fields only; similarity to current carrying wire.

So there is a magnetic field gradient. We know from Section 2.3 that a gradient in the magnetic field causes particle drifts:

vB = v2 2ωc B ×B B2 (7.3)

In our simple torus geometry, the result of the B drift is shown in Fig. 7.3. Because the drift velocity is 1q the electrons and ions drift in opposite directions. Working through the cross product with the right hand rule, the reader can verify that ions drift up and electrons drift down in this geometry.


Figure 7.3: B drift.

After some time has passed, the electron and ion drifts will necessarily induce charge separation. Positive charge will gather at the top of the machine and negative charge at the bottom. This is shown in Fig. 7.4. Naïvely one might expect that this will simply continue until the induced electric field cancels out the B drift and thus quiescence can be achieved.


Figure 7.4: Induced electric field and thus E × B drift due to the B drift charge separation.

There is, however, a major problem for this design. As shown in Fig. 7.4, this induced electric field is perpendicular to the toroidal field. There is therefore a E × B drift, which is outwards as shown in the figure. As we know from Section 2.2, both electrons and ions will drift in the same direction at the same velocity due to this drift. The entire plasma is therefore shoved radially outwards. Numerical estimates show that this process will happen very fast (sub-ms) relative to necessary energy confinement times in tokamaks (s).


Figure 7.5: Tokamak field configuration with toroidal current and thus poloidal field.

One must therefore prevent the B drift from causing charge separation. The solution is shown in Fig. 7.5. If we add a toroidal current within the plasma itself, it will create a poloidal field and the overall field lines will wrap around the torus like stripes on a barber’s pole (or screw threads - indeed this is the toroidal analog of the screw pinch). In magnetic confinement machines we are typically limited to β 1, for instance β 0.1, in which case we know that the transport properties along the field lines are much faster than perpendicular to them. In this case then, we know that the electron conductivity along the lines of force is very high. With the field lines wrapped around the plasma, the top and bottom of the machine are essentially electrically shorted and no potential difference can be generated between them. The B drift can therefore generate some currents within the plasma but it cannot cause a charge separation, and we therefore avoid the disastrous E × B drift.

7.2 Lawson Criterion

The immediate question for any confinement scheme is how well it can perform relative to the demands set by fusion ignition and energy gain. We therefore follow a simple derivation for tokamaks by J. Lawson [Proc. Phys. Soc. B70, 6 (1957)].

We obviously consider a DT plasma confined in a tokamak. The fusion power which can be used for self-heating is generated by the alpha particles from the DT fusion reaction:

D + T n(14.1MeV) + α(3.5MeV) (7.4)

The total fusion power is:

Sα = 1 4Eαni2σv (7.5)

where the factor of 14 comes from nD = nT = ni2, Eα = 3.5 MeV, and σv is the fusion reactivity. Using that the total density n = (ne + ni) 2ni and that the total pressure is p = nT,

Sα = 1 16Eα p2 T2σv (7.6)

The fusion power must overpower loss mechanisms. In this simple analysis, we consider only Bremsstrahlung and heat conduction losses. The bremsstrahlung power derivation won’t be treated until later, so we simply give the result:

SB = CBZeffne2T 1 4CBZeff p2 T32 (7.7)

where CB is a constant of proportionality.

The last term to derive is a heat conduction loss term. In general, this would be

1 V q dA = 2κ r T r r=a (7.8)

Unfortunately the thermal conductivity κ is not generally well-known, and neither is the temperature gradient at the edge of the plasma. It is therefore convenient to introduce / define a general ‘energy confinement time’ such that this loss term is

3 2 p τE (7.9)

where τE characterizes the loss timescale for thermal energy due to thermal conduction, or other loss mechanisms not considered.

So, in general the power balance equation is:

Sα + Sh SB + Sκ (7.10)

where we also include an arbitrary heating term Sh, which is generally taken to be 0 in a steady-state ignited plasma. Substituting in the terms we have derived above:

Eασv 16 p2 T2 CB 4 p2 T32 + 3 2 p τE (7.11)

We can further simplify the constant coefficients:

Kασvp2 T2 KB p2 T32 + Kκ p τE (7.12)

this can be simplified algebraically to obtain a requirement on the pressure-confinement product:

pτE KκT2 Kασv KBT (7.13)

Since the fusion reactivity is just a function of the temperature, the right-hand side is a function of temperature only. This is the Lawson Criterion. For a given temperature that we can achieve in a fusion tokamak, the Lawson Criterion tells us what pressure-confinement product pτE we must achieve. For example, a detailed calculation reveals that if one can operate a tokamak at 15 keV then one must have pτE 8.3 atm-s. Typical energy confinement times are of order one second, which means that typical pressures in reactors will have to be of order ten atmospheres.

A plot of the minimum pτE is shown in Fig. 7.6


Figure 7.6: Plot of the Lawson Criterion (minimum pτ) versus temperature. The shaded region corresponds to positive fusion energy gain.

7.3 MHD Equilibrium in Tokamaks

7.3.1 Screw Pinch Radial Balance

Recall the work done previously in section 3.3. Now that we have added the toroidal current, and thus a poloidal field, we can see that the field configuration in the tokamak is at a most basic level the toroidal analog of the screw pinch:

d dr p + Bz2 2μ0 + B𝜃2 2μ0 + B𝜃2(r) μ0r = 0 (7.14)

In general, with two of the three unknowns specified (e.g. Bz and B𝜃) then MHD determines the third (e.g. p). However, because of the toroidal geometry there are several complications that arise, and are qualitatively discussed in the following sections.

7.3.2 Toroidal Force Balance - Hoop

First we consider the hoop force on a tokamak plasma. The situation is shown schematically in Fig. 7.7. The field of interest in this case is poloidal, and this is therefore applicable to Z-pinch or Screw-pinch plasmas. The plasma is shown in red in the figure, the current flows in or out of the page, and the poloidal field thus wraps around the plasma as shown.


Figure 7.7: Hoop force on a tokamak plasma.

Because of the toroidal geometry, the magnetic field induced has a 1r dependence and is therefore stronger on the inside of the plasma relative to the outside. On the other hand, the surface area on the inner half, S1, is smaller than the outer surface area since they are going like r. But the total pressure on one half surface of the plasma due to the magnetic pressure will be B2S and the quadratic dependence on the magnetic field wins. The force on the inside of the plasma is thus greater than on the outside, and this effect creates a net ‘Hoop’ force directed outwards. The name comes from the fact that this is analogous to the tension in a circular current-carrying wire loop.

7.3.3 Toroidal Force Balance - Tire Tube

The next situation to consider is analogous to the pressure exerted on the walls of a tire. Consider that the plasma surface is an isobar. In that case the force exerted on the surface of the plasma can be thought of as pS.Since S2 > S1 and p1 = p2 by assumption this creates a net force, which is also directed outwards. This is shown schematically in Fig. 7.8. This force does not depend on what the magnetic field is doing, and thus will occur in toroidal plasmas of all types (Z-pinch, 𝜃-pinch, Screw-pinch).


Figure 7.8: The tire tube force on a toroidal plasma.

7.3.4 Toroidal Force Balance - 1/R

The next toroidal force balance problem comes from the 1R dependence of the toroidal field. Because this requires a toroidal field, it occurs in the 𝜃-pinch or screw-pinch problems only. This situation is shown in Fig. 7.9. The coils carry a current Ic and the plasma has an induced opposite current Ip since it is diamagnetic. We assume that the plasma current is carried in an infinitesimal sheath on the outer surface.


Figure 7.9: Schematic of a situation in which the 1/R force arises.

We remember Ampère’s Law in integral form, which we can apply over a toroidal loop:

B d = 2πrBϕ = μ0Ienc. (7.15)

Just outside the plasma, the toroidal field is:

Bϕa = μ0Ic 2πr (7.16)

and just inside the inner plasma surface:

Bϕi = μ0(Ic Ip) 2πr (7.17)

The net force on a plasma surface is given by:

F = (Bϕa Bϕi)S2 2μ0 (7.18)

depending on the jump in field and the surface area.

We now need to consider the force on the inside of the plasma relative to the outside. Inside, 1R is smaller and the surface area S is also smaller, but because of the quadratic dependence on B and the fact that B 1R the total force on the inside of the plasma is actually greater than on the outside. This is analogous to the hoop force. The net result is that the 1R effect described here results in an outward directed net force on the plasma.

7.3.5 Toroidal Force Balance - Wall Stabilization

No matter which type of pinch we choose, the above effects all act to push the plasma radially outwards in a toroidal machine or tokamak. Obviously this needs to be avoided in a confinement or fusion machine, since the plasma must remain isolated from the first wall. The problem of stabilization thus arises.


Figure 7.10: Schematic of a perfectly-conducting wall stabilization.

One way to do this is to consider the effect of the wall. This is shown schematically in Fig. 7.10. In particular, consider a perfectly conducting (superconductor) wall. In this case we know that no magnetic flux can penetrate the wall. The initial situation is shown in the top half of the figure. The poloidal field is greater on the inside edge of the plasma, and due to the effects described in the previous sections we know that there is a net force on the plasma directed radially outwards. After some time, the plasma will be moved closer to the wall (bottom half of the figure). Since magnetic flux cannot penetrate a superconductor, the flux piles up on the outside of the plasma, increasing the strength of the poloidal magnetic field until equilibrium is achieved since this counteracts the outwards force.

This seems like an excellent stabilization property. Of course, the problem is that a superconducting first wall is technologically infeasible. In a real machine the first wall will be resistive. In this case the magnetic flux can penetrate the wall and resistively dissipate. The wall effect might slow down the outwards drift of the plasma, but it will not indefinitely stop it and provide stable confinement.

7.3.6 Toroidal Force Balance - Vertical B Stabilization

We note that many of the above effects depend on the fact that the poloidal field is stronger on the inside edge of the plasma. This is shown schematically in Fig. 7.11. A simple way to counteract this is to add a constant vertical magnetic field, shown below the schematic. This adds to the poloidal field on the plane and tends to increase the field strength on the outer part of the plasma while decreasing it on the inner edge. This counteracts the toroidal force effects and can negate the radial force. If we need to counteract effects like the tire tube force, which do not depend on a poloidal field, the vertical field could be further increased until the field on the outer edge is actually greater than the field on the inner edge, thus providing confinement.


Figure 7.11: The effect of adding a vertical field.

This technique is used in many real machines and is part of the design for ITER.

7.4 Heating

This section gives a brief overview of potential heating mechanisms in tokamak systems. We recall that the Lawson criterion requires we heat the plasma to make keV before fusion α self-heating can take over.

7.4.1 Ohmic Heating

We have already seen that for several reasons we will need to have a toroidal current in the plasma. Since the plasma has finite resistivity, the toroidal current will provide some Ohmic heating:

P = IV = I2R (7.19)

the power is proportional to the square of the toroidal current and the plasma resistivity. This looks great but we recall that the plasma resistivity is T32. This means that Ohmic heating becomes less efficient as the plasma is heated. A detailed analysis reveals that Ohmic heating is great for initial heating of the plasma but cannot get a plasma to fusion-relevant temperatures on its own.

For initial start-up, most machines use a transformer scheme (with the plasma as the secondary winding) to drive current inductively for Ohmic heating.

7.4.2 Neutral Beams

Another method for heating a plasma is to use neutral beams of energetic ions. For example, one could use an accelerator to generate beams of energetic D which impinge on the plasma. As the atoms hit the plasma they are ionized, and then slow quickly via Coulomb collisions, transferring their energy to both plasma species. The required ion energy is set by the temperature E kBT and also by the areal density of the tokamak plasma - ideally the ions stop near the center of the plasma.

Neutral beams are required so that the particles do not undergo cyclotron motion until they are collisionally ionized within the plasma itself. The scheme to create a neutral beam is a typical electrostatic ion accelerator plus a beam neutralizer.

In current machines, beam energies are typically of order 100 keV. A reactor-scale machine would require energies of order 1 MeV. Neutral beam heating is expected to play a part in any magnetic fusion machine.

7.4.3 RF Heating

Another very common scheme is to use high-frequency electromagnetic waves to heat the plasma. Generally the wave frequency is chosen to match either the electron or ion cyclotron resonance. We know that the electron cyclotron frequency rule-of-thumb is fce 28 GHz/T. So for a toroidal field in the machine of 10 T, we must produce radiation at 280 GHz which is technically challenging to produce that radiation. Ion cyclotron heating, on the other hand, is at much lower frequency fci 14A MHz/T. So in a 10 T field with A = 2 (D resonance), the ICH heating system must produce 70 MHz radiation. The absorption is collisionless so these processes are still efficient at high temperature. This is simple, but it turns out that the launching antenna must be very close to the plasma surface which is technically challenging.

ECH and/or ICH is expected to play a role in any fusion machine.

7.4.4 Lower Hybrid Current Drive

A related question is how we drive current in a tokamak. In the Ohmic heating section we discussed the use of a pulsed transformer to induce large toroidal currents which can Ohmically heat the plasma early in the discharge. However, we must have a scheme to drive steady-state toroidal currents in a fusion machine. Neutral beams and ECH/ICH can be used to do this but it turns out the efficiency is low. The problem of ‘current drive’ is an area of investigation. One way to do this is to launch waves at the lower-hybrid resonance such that the k of the radiation is toroidally directed. It turns out that these waves tend to collisionlessly Landau damp on the electrons, which drives a toroidal current.

Detailed analysis shows that this is a relatively efficient way to induce a steady-state toroidal current but it cannot provide quite enough. A real fusion machine will depend on a large amount of ‘bootstrap’ current for steady-state operation.

7.5 Banana Orbits

We finish this section with a brief qualitative description of an effect of toroidal geometry on single-particle motion.


Figure 7.12: Single-particle motion in a tokamak: Banana orbits.

We recall that the field strength exhibits a 1R dependence, which is shown schematically at the top of Fig. 7.12. Consider a particle which is confined by cyclotron motion to one line of force. Because a tokamak is a screw-pinch style confinement scheme, as we have argued previously, the lines of force wrap around the plasma. If we project a single line of force into the poloidal plane, though, it will be circular. This is shown by the dashed line in Fig. 7.12. Naïvely, one might expect that the particles will simply gyrate around the line of force all the way around. However, as the particle moves from the outside of the plasma towards the inside of the machine the magnetic field strength is increasing. If the particle has non-zero transverse momentum, then it will exhibit a mirror machine style confinement and reflection. Recall Section 2.8, where we derived that:

sin2𝜃 c = Bmin Bmax (7.20)

any particles at the location of Bmin with pitch angle greater than 𝜃c will be confined to a banana-style orbit as shown in green.


Figure 7.13: Schematic illustration of particles with various degrees of mirror confinement in a banana orbit.

Of course particles have a distribution (typically Maxwellian) and thus there is a wide range of banana orbits. This is shown very schematically in Fig. 7.13. The green orbit is a highly confined particle, i.e. with high pitch angle. The blue orbit is moderately confined, and the red orbit is barely confined. Of course there are also unconfined orbits, which transit the entire machine.

In a real machine, the banana orbits will affect transport properties in the plasma. We also note that, depending on the initial direction of the velocity, the banana orbits are qualitatively different. See Fig. 7.14. This can also affect transport properties.


Figure 7.14: Depending on the direction of the initial velocity, particle banana orbits are qualitatively different.

Chapter 8

In this chapter we derive a few scattered plasma physics topics related to inertial confinement fusion.

8.1 Transport

8.1.1 Bremsstrahlung Radiation

For simplicity we begin with the cross section for emission of a photon with energy hν by an electron with velocity ve interacting with an ion via the Coulomb force. This is the Kramers cross section:

dσ dν = 32π2 33 Z2e6 me2c3ve2 1 hν. (8.1)

The next task is to calculate the spectral emissivity per unit mass. To do this we have to integrate the Kramers cross section over the electron distribution:

ην = hν 4πAmpvmind3vf(v)vdσ dν, (8.2)

where the lower limit of integration is the minimum electron velocity needed to create a photon of energy hν, given by hν = mevmin22. Using the Kramers cross section together with the Maxwellian distribution:

f(v) = 2 π m kT32emev22kBT , (8.3)

we can write the spectral emissivity as

ην = hν 4πAmp 32π2 33 Z2e6 me2c3 1 hν 2 π m kT32vmindvvexp[m ev22k BT] (8.4)

after doing the integral, only the lower limit remains and we substitute hν = mevmin22. Also with some algebraic cleanup, the emissivity expression becomes:

ην = 16π 36 e6 me2c3 Z2ne AmpkB Te me exp hν kBTe (8.5)

The total power radiated into 4π is obtained by integrating η over all frequencies:

Pbr = 4π0η νdν, (8.6)

Pbr = 32π 36 e6 mec2kB Te mec2 Z3ρ (Amp)2 (8.7)

We note that ignoring the constants gives a proportionality:

Pbr TeZ3ρA2 (8.8)

8.1.2 Electron Heat Conduction

In many ICF-relevant situations the electron heat conduction is critically important for the implosion dynamics. The classical heat conductivity is derived from the Fokker-Planck equation:

f t + v f eE m vf = Av vf v v(v vf) v3 + Cee(f) (8.9)

where we have introduced A = (2πnZe4m2)lnΛe, with lnΛe the Coulomb log for electron-ion collisions. Cee represents a formal treatment of electron-electron collisions, which are neglected in this derivation.

Here will will assuming a temperature gradient in the x̂ direction, and consider a perturbative approach in which the distribution function becomes:

f(v,μ) = f0(v) + μf1(v) (8.10)

where f0 is the isotropic equilibrium Maxwellian. μ = cos𝜃 where 𝜃 is the pitch angle between the electron velocity and the temperature gradient. We are also implicitly assuming that the electron mean free path is much less than the gradient length scale:

L dlnTdx1 λ e. (8.11)

The 0-th order distribution function is written:

f0(v) = n (2π)32vth3 exp(v22v th2) (8.12)

If we insert the perturbed distribution function f(v,μ) into the Fokker-Planck equation and keep only terms that are first order in μ, we get:

f1 t + vf0 x eE m f0 v = 2A v3 f1 (8.13)

rearranging gives

f1 = v4 2A f0 x eE mv f0 v (8.14)

We need another equation at this point. We use the plasma property of quasi-neutrality to require that the current due to f1 vanishes in steady state, or:

jx = 0 = ex̂ dvμvμf1 (8.15)

using the previous expression for f0 and f1 , glossing over some algebra, yields:

eE = 5 2kBdT dx (8.16)

It is an interesting aside that we have shown here that a gradient in temperature generates an electric field. Anyways, this can be replaced in the equation for f1 which depends only on derivatives of f0 and can be shown to be:

f1(v) = f0(v)vth4 4A dlnT dx v vth 4 v2 vth2 8 (8.17)

It follows from this that the heat flux is:

q =d3v(mv22)μvμf 1(v) = χdT dx (8.18)

where we have defined the thermal conductivity:

χS = 8 π32G(Z)(kBT)52kB Ze4mlnΛ (8.19)

where the function G(Z) takes into account effects of electron collisions, which are more important in low-Z plasmas. We note:

χS G(Z) Z T52 lnΛ (8.20)

The conductivity is primarily a function of temperature.

A detailed look at the velocity integral reveals that electrons with ve 3.7vth end up doing the brunt of the work in conducting thermal energy. It is therefore important to revisit our assumption that f1 was a small perturbation to f0. If we write

|f1| f0 10q vthnkBT 1 (8.21)

the real heat flux must satisfy this inequality. We can define a free-streaming flux:

qF vthnkBT (8.22)

physically, this represents the entire thermal energy of the plasma nkBT moving with the thermal velocity, and represents an upper limit on flux. The real conductivity thus must satisfy:

q qF 0.1 (8.23)

The heat flux cannot exceed 10% of the free-streaming limit. Generally in hydrodynamics simulations one takes:

q = min(qS,fLqF ) (8.24)

which is the minimum of the Spitzer conductivity and the free-streaming flux limit multiplied by a flux limiter which is generally taken in the range 0.03 0.1.

Detailed Fokker-Planck simulations reveal that the Spitzer conductivity is accurate when Lλe 102, i.e. in the long gradient length limit. As the gradient length scale is decreased, flux-limited conduction must be used. This is particularly important in calculations of direct-drive ICF targets, where the thermal conduction between the critical surface and ablation front is fundamental to the entire dynamics of the problem, and the temperature gradient is very steep between the hot plasma at the critical surface and cold plasma at the ablation front.

8.2 Laser-Plasma Interactions

8.2.1 Critical Density

In Section 4.2 we derived the dispersion relation for electromagnetic waves in unmagnetized plasmas, which we reproduce here:

ω2 = ω pe2 + c2k2 (8.25)

We can immediately see that there is a cutoff (k = 0, no propagation) when the wave frequency is equal to the electron plasma frequency. Since the latter depends on the plasma density, this presents a limit on density:

ω2 = ω pe2 = nce2 𝜖0me (8.26)


nc = ω2𝜖0me e2 (8.27)

where nc is the critical density. Electromagnetic waves can only propagate for densities n < nc.

8.2.2 Inverse Bremsstrahlung

If we consider a laser beam propagating from low-density plasma towards a higher density region (n > nc) at some point the laser energy must be either scattered or absorbed in the plasma. The best situation for ICF is when the laser wave is absorbed collisionally, or more colloquially via the inverse bremsstrahlung mechanism.

If we treat the laser electric field as:

E(r,t) = E(r)eiωt (8.28)

and write the electron momentum equation as:

mene ue t + ue (u) = neqe(E + ue ×B) meneνeiue (8.29)

In previous derivations of EM waves in plasmas we neglected collisions, i.e. we had assumed νei = 0. In this problem we keep finite collision frequency, since that is the effect we are looking for. If we assume that u0 = E0 = B0 = 0 and keep only first-order terms:

meneue t = neeE meneνeiue, (8.30)

ue t = e meE νeiue. (8.31)

Where we have implicitly assumed that the ion inertia is infinite over the timescales of this process. Since the field is harmonic, we can use Fourier analysis on this equation:

iωue = e meE νeiue (8.32)

we can rearrange this equation to acquire the plasma response:

ue = ieE me(ω + iνei) (8.33)

which leads to the current:

j = neeue = i𝜖0ωpe2 ω + iνeiE (8.34)

We now go back to Maxwell’s equations for a moment:

×E = B t (8.35) ×B = μ0j + 1 c2 E t (8.36)

Combining these and using Fourier analysis, we get that:

× (×E) = μ0 tj + 1 c2 2E t2 (8.37)

Using the vector identity for the left hand side:

× (×E) = 2E (E), (8.38)

and using Fourier analysis we get that

k2E = ω2 c2 iωμ0j (8.39)

Now we need to substitute into this equation for the current to get:

ω2 = c2k2 + ω pe2(1 + iν eiω)1 (8.40)

If the damping is weak, we can use the binomial approximation:

ω2 = c2k2 + ω pe2 1 iνei ω (8.41)

Clearly we now have wave damping via the imaginary part of the dispersion relation. If we take k as imaginary with k = kr + iki2 with weak damping (i.e. small ki), we can isolate the imaginary part of the above:

0 = c2k rki ωpe2ν eiω (8.42)

From the undamped dispersion relation we can approximate:

kr = 1 cω2 ωpe 2 (8.43)

we follow Kruer and write the damping coefficient as:

ki = ωpe2 ω2 νei vg (8.44)

the energy damping length for inverse bremsstrahlung is simply ki1.

8.2.3 Resonance Absorption


Figure 8.1: The resonance absorption process for p-polarized light.

We start off a series of qualitative descriptions of non-collisional absorption processes with resonance absorption. For p-polarized light, the schematic is shown in Fig. 8.1. At the top we show a schematic density profile where the plasma has a gradient from 0 past the critical density, with n . Next we show the p-polarized wave propagating in the plasma in green. The cross-hatches represent the electric field’s direction. The wave is propagating at an angle α relative to the density gradient. Due to refraction, the wave bends, and reaches a maximum in density at z = Lsinα where L is the gradient length scale. However, since the electric field has a component in the direction there will be an evanescent component in the z direction. This reaches the surface of critical density, and can resonantly couple to the plasma there. The bottom part of Fig. 8.1. The resonant field coupling to the plasma at the critical density can excite electron plasma waves, or ‘plasmons’. For more detail see Kruer.

8.2.4 Parametric Instabilities

Now we discuss a few parametric instabilities. The essential process is three-wave coupling:

ωL = ω1 + ω2 (8.45)

kL = k1 + k2 (8.46)

An incident laser photon decays into two daughter waves, 1 and 2. These can be several combinations of scattered photons, plasmons (electron plasma waves), and phonons (ion acoustic waves).

Two-Plasmon Decay


Figure 8.2: The two-plasmon decay.

First up is the two-plasmon decay, which is shown schematically in Fig. 8.2. The dispersion relation for light waves is shown in green and the dispersion relation for electron plasma waves, ‘plasmons’, is shown in blue. A single incident photon can decay into two plasmons via the frequency and k matching conditions. As shown, one plasmon propagates in the direction of the original photon and the other goes backwards. Because the electron plasma wave dispersion relation is so much flatter than the light wave, by a ratio vthc, we can approximate ωpe ωL2. Since ωpe n, this implies that the two-plasmon decay can occur occur very close to a plasma density of nc4,the ‘quarter-critical’.

Stimulated Raman Scattering


Figure 8.3: The Stimulated Raman Scattering process.

The next possibility we consider is called ‘Stimulated Raman Scattering’ (SRS). The schematic for this decay is shown in Fig. 8.3. Once again the green curve shows the light wave dispersion relation, and the blue curve shows the plasmon dispersion relation. Once incident photon, labelled ‘inc’, decays into a scattered photon and a plasmon. The situation shown is a backscattered photon, but this is not necessarily the case. SRS-scattered photons can go forwards, sideways, etc. The maximum density at which SRS can occur is quarter-critical. However, it can also occur at much lower densities, and is thus distinct from TPD.

Stimulated Brillouin Scattering


Figure 8.4: The Stimulated Brillouin Scattering.

The last parametric decay instability is Stimulated Brillouin Scattering (SBS). Once again, this is shown schematically in Fig. 8.4. Yet again green denotes light waves. SBS, consists of an incident laser photon decaying into a backscattered photon and an Ion Acoustic Wave (IAW), or ‘phonon’. These are shown in red. Because of ω k for IAWs, and they are low frequency, the wave matching conditions require that the scattered photon is backscattered and the phonon goes nearly in the same direction as the original wave to conserve momentum ki 2kL. The SBS process can occur anywhere in under-dense plasma.

8.2.5 Hot Electrons

It is important to have a brief discussion of hot electrons. In ICF implosions, it is important to keep the fuel adiabat low so that it is easily compressible to very high densities during the implosion. This means that any unintended heating of the fuel early on in the target illumination is undesirable, since it will raise the adiabat. The relevance here is that many of these processes can generate very high-energy electrons. Any time a laser-plasma instability generates plasmons, or electron plasma waves, those waves can Landau damp on the plasma electrons. When the wave is very strong, there is a non-linear wave-breaking process which efficiently couples wave energy to particle energy and generates a two-temperature electron distribution. The hot electrons can have a temperature of many tens of keV. These electrons are energetic enough to penetrate the fuel layer in an implosion target and preheat it. It is therefore important to understand and control LPI, in addition to controlling LPI to increase laser coupling to the target.

8.3 Misc

We now discuss a few totally random but important topics of relevance to ICF.

8.3.1 EM Field Generation

Field generation in plasmas is an important topic, particularly for spontaneous generation of fields in ICF. If we recall Ohm’s law from the discussion of magnetohydrodynamics, there is a pressure gradient term. This work will follow that of Haines. To first order we can write from Ohm’s Law:

E = pe nee (8.47)

Physically, if there is a pressure gradient in the plasma the electrons quickly leave until there is an electrostatic field set up to mitigate this. See also the previous discussion of ambipolar diffusion. Other terms in the generalized Ohm’s law are omitted for clarity. This is also ignoring any fields from incident lasers in an ICF scheme. Using Faraday’s law we can get the magnetic field:

B t = ×E (8.48)

and using p = nkBT,

B t = ×pe nee = kBTe ×ne nee (8.49)

This is the most basic mechanism for spontaneous magnetic field generation. If the temperature and density gradients are not parallel in a plasma, then a magnetic field is generated. This is also known as the Biermann Battery.

8.3.2 Rayleigh-Taylor Instability

In a situation where a heavy fluid is supported against a gravitational potential by a lighter fluid, then it is clearly energetically and thus dynamically favorable to exchange material of the heavier fluid for the lighter one. Doing this from first principles is non-trivial. However, if we consider two homogenous fluids with an interface at z = 0 with densities ρ1 and ρ2. Stealing a result from Drake, the equation of motion for perturbations of this interface is

z ρnw z = k2 ρnw + wg n ρ z (8.50)

We take the perturbation in each fluid as:

w1 = w0ekz (8.51) w2 = w0ekz (8.52)

plugging into the first equation, we get that:

w0 g n(ρ2 ρ1) = n k2(ρ2 + ρ1)kw0 (8.53)

and thus, the growth rate, which we (Drake) have called n0, is given by:

n0 = ρ2 ρ1 ρ2 + ρ1kg = An kg (8.54)

where g is the gravitational acceleration, k is the wave number, and An is the Atwood Number:

An ρ2 ρ1 ρ2 + ρ1 (8.55)

Obvious An is the stable situation, and larger An correspond to larger differences in density across the interface and thus higher growth rates.

Chapter 9
Back of the Envelope

In this chapter we simplify a few of the previously-derived relations to ‘back of the envelope’ or Formulary-style equations which can be used to quickly calculate important plasma properties. The NRL Plasma Formulary (available online) is also an excellent resource for this sort of calculation.

9.1 Plasma Frequency

We start off with the plasma frequency.

9.1.1 Electron

Way back in Section 1.3 we derived the general result for the electron plasma frequency:

ωpe = ne2 𝜖0me(mks) = 4πne2 me (cgs) (9.1)

where the first expression is in the MKS system of units (generally used in this work) and the second is the CGS system (commonly used in older physics materials, and also in ICF). We also note that often f is a more desired quantity than ω, e.g. if one wishes to compare to a source of electromagnetic radiation. We can therefore generate Table 9.1.

Table 9.1: Electron Plasma Frequencies




9.1.2 Ion

The ion plasma frequency is not as widely used, however it is useful to note that

ωp Z m (9.2)

and therefore we can write the ion plasma frequency in terms of the electron values:

ωpi = 1 42 Z Aωpe (9.3)

where the fraction out front comes from the square root of the ion to electron mass ratio.

9.2 Debye Length

Recall the result of Section 1.2:

λDe = 𝜖0 kB T ne2 (MKS) = kB T 4πne2(CGS) (9.4)

In addition to the constants, we note that the Debye length depends on both temperature and density. Once again, we can give the result in both systems of units for convenience, though in both equations we absorb Boltzmann’s constant into the temperature and take [T] =eV:

λDe = 7430Tn(MKS) (9.5)

λDe = 743Tn(CGS) (9.6)

9.3 Cyclotron Motion

In this section, refer back to Section 2.1 for derivations. In this section, for convenience, we always take [B] =T.

9.3.1 Gyroradii

We recall the result, in MKS units:

rL = v ωc = mv qBz (9.7)

which leads to the following:

rLe = 2.4TeBcm (9.8)

rLi = 100ATe(ZB)cm (9.9)

where the results are given in cm, but the conversion to meters is trivial since we have used ‘convenient’ units for both Te and B.

9.3.2 Cyclotron Frequencies

From the previous results we can immediately write that:

ωc = qB m (9.10)

The cyclotron frequency only depends on the particle type and the magnetic field. By far the easiest values to remember are:

fce = 28GHzT (9.11)


fci = 14(ZA)MHzT (9.12)

one simply has to multiply by the magnetic field strength in Tesla to get the gyro frequencies.

9.4 Collision Frequencies

Here we will use the results of 5.8.

νei = e4 4π𝜖02me2 ni Ti32 lnΛ 1 ve3 (9.13)

For a thermal distribution, this expression simplifies to:

νei = (3 × 106)ne lnΛ Te32 s1(CGS) (9.14)

νei = (3 × 1012)ne lnΛ Te32 s1(MKS) (9.15)

where of course [ν] = s1 and the only different between systems of units is [n].

Also recall the scaling result:

νee νei (9.16) νii μν ei (9.17) νie μνei (9.18)

where μ = memi, so based on the species mass ratio, and there is an assumption buried in here that we are working in a Z = 1 plasma. For non-hydrogenic ions, the ion collision rate is enhanced by a factor of Z4.

9.5 Thermal Velocities

For a Maxwellian distribution function, we know from basic physics that the mean thermal velocity is:

vth = 2kB Tm (9.19)

This leads to the simple expressions:

vTe = (4.2 × 107)T ecm/s (9.20)

vTi = (9.8 × 105)T iAcm/s (9.21)

which are obviously in CGS units. Again we have used [T] =eV. The extension to MKS is trivial.

9.6 Mean Free Path

We just derived the collision frequencies and thermal velocities. We can simply write that:

λ = vth ν (9.22)

which leads to the following for electrons:

λe = (4.2 × 107)T e ×(3 × 106)ne lnΛ Te32 1 (9.23)

λe = (1.2 × 1013) Te2 ne lnΛcm (9.24)

and for ions:

λi = (9.8 × 105)T iA ×(3 × 106)μZ4ni lnΛ Ti32 1 (9.25)

λi = (1.4 × 1013) Ti2 AZ4ni lnΛcm (9.26)

all of these expressions are in CGS units with our usual convention of [T] =eV.

9.7 Ion Sound Speed

We use the derivation presented in Sec. 4.6. In particular:

cs = Te + γTi mi . (9.27)

for the ion sound speed. We can simplify this somewhat by writing:

cs = γTe mi . (9.28)

where we have to be somewhat careful about what exactly is meant by γ. Anyways, plugging in some numbers gets us that

cs (9.8 × 105)ZAγT ecm/s (9.29)

which is obviously in CGS units. Once again, [Te] =eV.

9.8 E × B drift speed

The drift velocity is:

v = E ×B B2 = E B (9.30)

where we assume the electric and magnetic fields are perpendicular. A handy formula, with [E] = V/m and [B] =T is

v = 104E Bm/s (9.31)

9.9 Plasma Parameter

The plasma parameter is defined as the number of particles in a Debye sphere:

ND = 4π 3 nλDe3 (9.32)

this is useful to know, since for classical plasma behavior we typically require ND 1. For simple calculations,

ND (1.7 × 109)T32 n (9.33)

in CGS units.

9.10 Plasma Beta

The plasma beta is defined as the ratio of thermal to magnetic pressure in a plasma;

β nkBT B28π, (9.34)

and is an important parameter for magnetically-confined plasmas. In the CGS system of units:

β (4 × 1011)nT B2 (9.35)

where as usual [T] =eV and [B] =T.

9.11 Fermi Energy

From Statistical Mechanics we know that the Fermi Energy is:

EF = 2(3π2)(23) 2me n23 (9.36)

in CGS units, this reduces to:

EF (3.6 × 1015)n23eV (9.37)

we can also calculate the Fermi pressure:

pF 2 5nEf (2.3 × 1033)n53bar (9.38)

9.12 Plasma Degeneracy

We just discussed the Fermi energy. In some ICF-relevant plasmas, the thermal energy can be comparable to or less than the Fermi energy, in which case we call it a degenerate plasma. It is useful to characterize these plasmas by:

𝜃 Thermal   Energy Fermi   Energy = kBT EF (9.39)

so, for quick back of the envelope calculations in CGS units,

𝜃 (2.7 × 1014) Te n23 (9.40)

the result, 𝜃, is a dimensionless quantity.

9.13 Plasma Coupling

In a similar vein, it is useful to think of the plasma coupling. We want to compare the thermal energy of plasma ions to the strength of the coupling between them. We define:

Γ = Potential Energy Thermal Energy = Z2e2a kBT (9.41)

so, once again for CGS units,

Γ = (1.4 × 107)n13 T (9.42)

9.14 Critical Density

We previously discussed the fact that the cutoff or critical density is very important for laser-plasma interactions, and thus for ICF. Recall, from section 8.2.1,

nc = ω2𝜖0me e2 (MKS) = ω2me 4πe2 (CGS) (9.43)

It is most useful to write this in terms of the laser wavelength λ = 2πcω. Then, when getting rid of the constants, we get:

nc = 1.1 × 1021 λμm2 cm3 (9.44)

the result is in CGS, and we have written the wavelength in units of μm for convenience. Thus, for a Nd:glass laser at 1ω, the critical density is 1021. If we frequency triple it, then the critical density is 1