Part 2, Chapter 5

Quasilinear Theory

Self-consistent evolution of waves and particles through resonant interactions

5.1 Introduction

Quasilinear theory describes the self-consistent evolution of wave fields and particle distribution functions when they interact through resonant wave-particle interactions. It goes beyond linear theory (which assumes fixed f₀) to account for how waves modify the distribution function, which in turn affects wave growth/damping.

The theory is called "quasilinear" because it retains terms that are second-order in the wave amplitude (δf₁ · E₁) while still assuming the perturbations are small. This allows us to describe:

Wave-induced diffusion: Resonant particles diffuse in velocity space due to waves
Plateau formation: Distribution function flattens at resonant velocities
Instability saturation: How growing waves eventually stabilize by modifying f₀
Turbulent heating: Energy transfer from waves to particles via diffusion

💡 Key Insight: Quasilinear theory bridges linear kinetic theory (Vlasov, Landau damping) and fully nonlinear phenomena (particle trapping, turbulence). It's valid when wave amplitudes are small but interactions persist long enough to change f₀.

5.2 Ordering and Assumptions

Expansion Scheme

We expand the distribution function and fields in powers of a small parameter ε (wave amplitude):

$$f(\mathbf{r}, \mathbf{v}, t) = f_0(\mathbf{v}, t) + f_1(\mathbf{r}, \mathbf{v}, t) + f_2(\mathbf{r}, \mathbf{v}, t) + \cdots$$

$$\mathbf{E}(\mathbf{r}, t) = \mathbf{E}_1(\mathbf{r}, t) + \mathbf{E}_2(\mathbf{r}, t) + \cdots$$

f₀ ~ O(1), f₁ ~ O(ε), f₂ ~ O(ε²) | E₁ ~ O(ε), E₂ ~ O(ε²)

Time Scale Separation

The theory assumes a separation of time scales:

Fast Time Scale

t_fast ~ ω⁻¹, k⁻¹v_th

• Wave oscillations
• Particle gyration
• Linear response (f₁ to E₁)

Slow Time Scale

t_slow ~ ε⁻² ω⁻¹

• Evolution of f₀
• Diffusion in velocity space
• Plateau formation

Validity Condition:

$$\frac{eE_1}{m\omega v_{th}} \ll 1 \quad \text{(weak turbulence)}$$

Wave electric field is small compared to particle thermal energy

5.3 Quasilinear Diffusion Equation

Derivation from Vlasov Equation

Starting from the Vlasov equation and averaging over fast oscillations, we obtain the quasilinear diffusion equation for the slowly evolving background distribution:

$$\frac{\partial f_0}{\partial t} = \frac{\partial}{\partial \mathbf{v}} \cdot \left( \mathbf{D} \cdot \frac{\partial f_0}{\partial \mathbf{v}} \right)$$

Quasilinear diffusion equation

where D is the diffusion tensor in velocity space:

$$D_{ij} = \frac{\pi q^2}{m^2} \sum_{\mathbf{k}} \delta(\omega_{\mathbf{k}} - \mathbf{k} \cdot \mathbf{v}) \, k_i k_j \, |E_{\mathbf{k}}|^2$$

Sum over all wave vectors; δ-function enforces resonance ω = k·v

Physical Interpretation

• Resonance condition: Only particles with v satisfying ω = k·v contribute to diffusion
• Random walk: Particles undergo random kicks from wave electric fields, leading to velocity diffusion
• Second-order effect: D ~ |E₁|², arises from correlations ⟨E₁ · δf₁⟩
• Irreversibility: Unlike linear Landau damping, quasilinear diffusion is irreversible (entropy increases)

1D Electrostatic Case

For electrostatic waves propagating in one direction (longitudinal waves along x):

$$\frac{\partial f_0}{\partial t} = \frac{\partial}{\partial v} \left( D(v) \frac{\partial f_0}{\partial v} \right)$$

$$D(v) = \frac{\pi e^2}{m^2} \sum_{k} k^2 |E_k|^2 \delta(\omega_k - kv)$$

Diffusion coefficient peaks at resonant velocities v_res = ω_k/k

5.4 Wave Energy Evolution

Wave Kinetic Equation

The waves also evolve due to their interaction with particles. The spectral energy density W_k = ε₀|E_k|²/2 satisfies:

$$\frac{\partial W_{\mathbf{k}}}{\partial t} = 2\gamma_{\mathbf{k}} W_{\mathbf{k}}$$

γ_k is the growth/damping rate from linear theory

The growth rate depends on the current distribution function f₀(v,t):

$$\gamma_{\mathbf{k}} = -\frac{\pi \omega_{pe}^2}{2k} \left. \frac{\partial f_0}{\partial v} \right|_{v=\omega/k}$$

Growth rate proportional to slope of f₀ at resonance (Landau formula)

Energy Conservation

Total energy (particles + waves) is conserved in quasilinear theory:

$$\frac{d}{dt}\left( \int \frac{1}{2}mv^2 f_0 \, d^3v + \sum_{\mathbf{k}} W_{\mathbf{k}} \right) = 0$$

Energy lost by waves → gained by particles (or vice versa)

Physical Picture: Growing waves extract energy from gradients in f₀. The diffusion flattens f₀ at resonant velocities, reducing ∂f₀/∂v → reducing growth rate → saturation. The wave energy is eventually absorbed by particles through the diffusive heating process.

5.5 Plateau Formation

Relaxation to Marginal Stability

Consider an unstable initial distribution f₀(v,0) with ∂f₀/∂v > 0 in some velocity range (e.g., bump-on-tail). Waves grow and cause diffusion at resonant velocities. The system evolves toward a state where:

$$\frac{\partial f_0}{\partial v} \approx 0 \quad \text{at resonant velocities}$$

Plateau formation (marginal stability)

This is called a plateau in velocity space. Once formed, the plateau ismarginally stable: any small perturbation that creates ∂f₀/∂v ≠ 0 immediately drives diffusion to restore the plateau.

Time Scale for Plateau Formation

$$t_{plateau} \sim \frac{1}{\gamma_{QL}} \sim \frac{m^2 v_{th}^2}{\pi e^2 k^2 |E_k|^2}$$

γ_QL ~ D/v_th² is the quasilinear relaxation rate

Before Plateau

• Positive slope ∂f₀/∂v > 0
• Waves grow exponentially
• Wave energy W_k increases
• Strong diffusion at resonance

After Plateau

• Flat distribution ∂f₀/∂v ≈ 0
• Waves saturated (γ → 0)
• Wave energy constant
• Marginal stability maintained

5.6 Applications and Limitations

Applications

• Bump-on-tail instability: Saturation via plateau formation
• Current-driven instabilities: Quasilinear diffusion limits runaway electrons
• Solar wind heating: Wave-induced diffusion heats minor ions
• Radiation belts: Pitch-angle diffusion by whistler waves
• Anomalous transport: Enhanced collision rate in fusion plasmas

When Quasilinear Theory Fails

• Large amplitudes: eE/(mωv_th) ~ 1 → particle trapping
• Wave-wave coupling: Nonlinear three-wave interactions
• Coherent structures: BGK modes, phase-space holes

Interactive Simulations

Quasilinear Velocity Diffusion: Plateau Formation

Python

script.py57 lines

import numpy as np

# Quasilinear diffusion: df/dt = d/dv [D_ql * df/dv]
# D_ql = (e/m)^2 * pi * |E_k|^2 * delta(omega - k*v)
# Result: distribution flattens to form a plateau at v ~ omega/k

m_e = 9.109e-31
e = 1.602e-19

# Initial Maxwellian + beam
n0 = 1e18
T_eV = 10.0
v_th = np.sqrt(T_eV * e / m_e)

# Beam parameters
n_beam = 0.01 * n0
v_beam = 5 * v_th
dv_beam = 0.5 * v_th

v = np.linspace(-2*v_th, 8*v_th, 50)

# Initial distribution: Maxwellian + beam
f0 = n0 / (np.sqrt(2*np.pi) * v_th) * np.exp(-v**2 / (2*v_th**2))
f_beam = n_beam / (np.sqrt(2*np.pi) * dv_beam) * np.exp(-(v-v_beam)**2 / (2*dv_beam**2))
f_initial = f0 + f_beam

print("Quasilinear Diffusion: Plateau Formation")
print("=" * 55)
print(f"Background: n_0 = {n0:.0e} m^-3, T_e = {T_eV:.0f} eV")
print(f"Beam: n_b/n_0 = {n_beam/n0:.2f}, v_beam = {v_beam/v_th:.1f} v_th")
print()

# Simulate plateau formation (simplified)
# Plateau forms between v_min and v_beam where df/dv > 0
# Find the positive slope region
df_dv = np.gradient(f_initial, v)
unstable = df_dv > 0

# After quasilinear relaxation: plateau in unstable region
f_final = f_initial.copy()
# Find plateau region
v_min_plateau = 3.0 * v_th  # approximate
v_max_plateau = v_beam + dv_beam
mask = (v >= v_min_plateau) & (v <= v_max_plateau)
if np.any(mask):
    plateau_val = np.mean(f_initial[mask])
    f_final[mask] = plateau_val

print(f"{'v/v_th':>8} {'f_initial':>12} {'f_final':>12} {'df/dv':>12} {'Status':>10}")
print("-" * 56)
for vi, fi, ff, dfi in zip(v[::3], f_initial[::3], f_final[::3], df_dv[::3]):
    status = "unstable" if dfi > 0 and vi > 2*v_th else "stable"
    print(f"{vi/v_th:8.2f} {fi:12.3e} {ff:12.3e} {dfi:12.3e} {status:>10}")

print(f"\nPlateau forms in range v/v_th = [{v_min_plateau/v_th:.1f}, {v_max_plateau/v_th:.1f}]")
print(f"Energy released goes into wave spectrum")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Quasilinear Wave Energy Evolution

Fortran

program.f9067 lines

program ql_wave_energy
  implicit none
  ! Quasilinear wave energy evolution:
  ! dW_k/dt = 2 * gamma_k * W_k
  ! where gamma_k = Landau growth/damping rate
  ! gamma_k ~ -(pi/2) * omega_pe^2/k^2 * (df/dv)|_{v=omega/k}
  real(8) :: e, m_e, eps0, k_B, pi
  real(8) :: n0, T_eV, v_th, omega_pe, lambda_D
  real(8) :: k, v_phase, gamma_k, W_k
  real(8) :: dt, t, f_slope
  real(8) :: n_beam, v_beam, dv_beam
  integer :: i, j, nsteps

e = 1.602d-19
  m_e = 9.109d-31
  eps0 = 8.854d-12
  pi = 4.0d0 * atan(1.0d0)

n0 = 1.0d18
  T_eV = 10.0d0
  v_th = sqrt(T_eV * e / m_e)
  omega_pe = sqrt(n0 * e**2 / (eps0 * m_e))
  lambda_D = v_th / omega_pe

! Beam parameters
  n_beam = 0.01d0 * n0
  v_beam = 5.0d0 * v_th
  dv_beam = 0.5d0 * v_th

write(*,'(A)') "Quasilinear Wave Energy Spectrum Evolution"
  write(*,'(A)') "==========================================="
  write(*,'(A,ES10.3,A)') "omega_pe = ", omega_pe, " rad/s"
  write(*,'(A,F6.1,A)') "v_beam/v_th = ", v_beam/v_th, ""
  write(*,*)

write(*,'(A12,A12,A14,A14,A14)') &
    "k*lam_D", "v_ph/v_th", "df/dv", "gamma/w_pe", "tau_grow(s)"
  write(*,'(A)') repeat("-", 68)

do i = 1, 15
    k = dble(i) * 0.03d0 / lambda_D
    v_phase = omega_pe / k  ! approximate for Langmuir wave

! Slope of f at resonant velocity
    ! Maxwellian part: df/dv = -v/(v_th^2) * f
    ! Beam part: df/dv = -(v-v_beam)/(dv_beam^2) * f_beam
    f_slope = -v_phase/v_th**2 * n0/(sqrt(2*pi)*v_th) &
              * exp(-v_phase**2/(2*v_th**2))
    f_slope = f_slope + (-(v_phase-v_beam)/dv_beam**2) &
              * n_beam/(sqrt(2*pi)*dv_beam) &
              * exp(-(v_phase-v_beam)**2/(2*dv_beam**2))

! Growth rate
    gamma_k = -pi/2.0d0 * omega_pe**2/k**2 * f_slope

if (abs(gamma_k) > 0.0d0) then
      write(*,'(F12.4,F12.2,ES14.3,ES14.3,ES14.3)') &
        k*lambda_D, v_phase/v_th, f_slope, &
        gamma_k/omega_pe, 1.0d0/abs(gamma_k)
    end if
  end do

write(*,*)
  write(*,'(A)') "Positive gamma: wave grows (beam-driven)"
  write(*,'(A)') "Negative gamma: Landau damping"
end program ql_wave_energy

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Previous: Fokker-Planck Next: Particle Trapping →