Plasma Physics from One Particle to Kinetic Theory#

A plasma can be described at three levels: as individual charged particles gyrating in fields, as a conducting fluid (magnetohydrodynamics), or as a distribution function in phase space (kinetic theory). This tutorial works through all three with physicskit.plasma. It checks each numerical result against the formula it should reproduce, and ends with Landau damping, a collisionless effect that only the kinetic description captures.

Note

Units. The single-particle and fluid functions take SI inputs and return SI outputs. QE, MP and ME are the elementary charge and the proton and electron masses in SI. The particle-in-cell (PIC) code instead uses the standard normalized units of plasma simulation: time in \(1/\omega_{pe}\), length in Debye lengths \(\lambda_D\), and velocity in thermal speeds \(v_{th} = \omega_{pe}\lambda_D\). Keep the two conventions apart.

Gyration: the cyclotron frequency and Larmor radius#

A charge in a uniform magnetic field circles it at the cyclotron frequency \(\omega_c = qB/m\), on a circle of Larmor radius \(r_L = m v_\perp/(|q|B)\). The Boris pusher (boris_integrate(); Boris, 1970), the standard integrator in PIC codes, splits each step into half electric kicks and a pure rotation. The magnetic force does no work, and a pure rotation can’t change the particle’s speed either, so the kinetic energy is conserved to rounding error:

import numpy as np
from physicskit import plasma

B0, v_perp = 1.0, 1.0e5  # 1 T, 100 km/s proton
omega_c = plasma.cyclotron_frequency(plasma.QE, plasma.MP, B0)
r_L = plasma.larmor_radius(v_perp, plasma.QE, plasma.MP, B0)
print(f"omega_c = {omega_c:.4e} rad/s, period = {2 * np.pi / omega_c * 1e9:.2f} ns, r_L = {r_L * 1e3:.4f} mm")
# omega_c = 9.5788e+07 rad/s, period = 65.59 ns, r_L = 1.0440 mm

pos, vel = plasma.boris_integrate(
    np.zeros(3), np.array([v_perp, 0.0, 0.0]), plasma.QE, plasma.MP,
    E=np.zeros(3), B=np.array([0.0, 0.0, B0]), dt=1e-10, steps=2000,
)
print(f"orbit radius = {np.ptp(pos[:, 1]) / 2 * 1e3:.4f} mm, "
      f"max |speed error| = {np.abs(np.linalg.norm(vel, axis=1) - v_perp).max():.1e} m/s")
# orbit radius = 1.0440 mm, max |speed error| = 1.5e-10 m/s

The simulated orbit has exactly the Larmor radius. Over three gyrations the speed changes by \(10^{-10}\) m/s out of \(10^5\), which is floating-point rounding.

Drifts: E×B is the same for every species#

Add a perpendicular electric field and the gyration centre drifts with

\[\mathbf v_E = \frac{\mathbf E\times\mathbf B}{B^2},\]

which contains neither the charge nor the mass. Electrons and ions drift together, so an \(\mathbf E\times\mathbf B\) drift carries the whole plasma without driving a current. To measure the drift cleanly, sample the position after a whole number of gyro-periods, where the gyration has returned to its starting phase:

E, B = np.array([0.0, 1.0e3, 0.0]), np.array([0.0, 0.0, B0])
print(plasma.exb_drift(E, B))
# [1000.    0.    0.]

for name, q, m in [("proton", plasma.QE, plasma.MP), ("electron", -plasma.QE, plasma.ME)]:
    period = 2 * np.pi / abs(plasma.cyclotron_frequency(q, m, B0))
    dt = period / 1000
    pos, _ = plasma.boris_integrate(np.zeros(3), np.array([v_perp, 0.0, 0.0]), q, m, E, B, dt=dt, steps=20 * 1000)
    print(f"{name:<8} drift over 20 gyro-periods: {pos[-1, 0] / (20 * period):.1f} m/s")

# proton   drift over 20 gyro-periods: 999.7 m/s
# electron drift over 20 gyro-periods: 999.7 m/s

The proton gyrates 1836 times more slowly than the electron and in the opposite direction, yet both guiding centres move at \(E/B = 1000\) m/s in \(+x\). A drift that does depend on charge, grad_b_drift(), separates the species and drives the ring current in Earth’s magnetosphere.

Collective scales: plasma frequency and Alfvén speed#

Collective behaviour starts once there are many particles. Displace the electrons and they oscillate about the ions at the plasma frequency \(\omega_{pe} = \sqrt{n e^2/(\varepsilon_0 m_e)}\), the natural frequency of a plasma. In the fluid (MHD) picture, magnetic field lines behave like strings under tension \(B^2/\mu_0\) loaded with the plasma’s mass density, and waves travel along them at the Alfvén speed \(v_A = B/\sqrt{\mu_0\rho}\) (Alfvén, Nature 150, 405, 1942):

omega_pe = plasma.plasma_frequency(n=1.0e12)       # ionospheric F-layer, m^-3
print(f"f_pe = {omega_pe / (2 * np.pi) / 1e6:.3f} MHz")
# f_pe = 8.979 MHz

v_A = plasma.alfven_speed(B=5.0e-9, rho=5.0e6 * plasma.MP)   # solar wind at 1 AU
print(f"v_A = {v_A / 1e3:.1f} km/s")
# v_A = 48.8 km/s

The first number is the familiar rule of thumb \(f_{pe} \approx 8.98\sqrt{n}\) Hz. It explains why the ionosphere reflects shortwave radio below about 9 MHz and is transparent above. The second is the typical Alfvén speed in the solar wind, around ten times slower than the wind itself, so the solar wind is super-Alfvénic. cold_plasma_dispersion() and stix_parameters() extend these two scales to the full set of cold-plasma waves (Stix, Waves in Plasmas, 1992).

Landau damping: collisionless decay in a particle simulation#

The fluid picture predicts that a Langmuir wave (an electron density oscillation) with wavenumber \(k\) oscillates forever at the Bohm-Gross frequency \(\omega^2 = \omega_{pe}^2 + 3k^2v_{th}^2\). Landau (1946) showed that in a collisionless plasma the wave nevertheless decays. Electrons moving at nearly the wave’s phase velocity exchange energy with it, and a Maxwellian has slightly more slower electrons, which gain energy, than faster ones, which lose it. In the weak-damping limit the rate is

\[\gamma = -\omega_{pe}\sqrt{\frac{\pi}{8}}\frac{1}{(k\lambda_D)^3} \exp\left(-\frac{1}{2(k\lambda_D)^2} - \frac32\right).\]

A 1D electrostatic PIC simulation (pic_simulate()) contains no damping term at all. It only pushes particles through the self-consistent field, yet it reproduces the decay. Seed a small density perturbation with landau_damping_ic() (a “quiet start” that avoids random sampling noise in the seed) and follow the field energy, which oscillates at \(2\omega\) and decays at \(2\gamma\):

from scipy.signal import find_peaks

k, v_th = 0.5, 1.0                      # k * lambda_D = 0.5, normalized units
L = 2 * np.pi / k                       # one wavelength, periodic
print(f"gamma (Landau)  = {plasma.landau_damping_rate(k, v_th):.4f}")
# gamma (Landau)  = -0.1514

x0, v0 = plasma.landau_damping_ic(100000, L=L, k_mode=k, alpha=0.05, v_th=v_th, seed=1)
run = plasma.pic_simulate(x0, v0, L=L, ng=64, dt=0.05, steps=500)
t, W = run["t"], run["field_energy"]

peaks, _ = find_peaks(W)
t_peaks, W_peaks = t[peaks][:6], W[peaks][:6]
print("peak times", t_peaks.round(2))
# peak times [ 2.4   4.65  6.9   9.   11.3  12.6 ]

gamma_fit = 0.5 * np.polyfit(t_peaks, np.log(W_peaks), 1)[0]
omega_fit = np.pi / np.mean(np.diff(t_peaks[:5]))
print(f"gamma (PIC fit) = {gamma_fit:.4f}, omega (PIC fit) = {omega_fit:.3f}, "
      f"Bohm-Gross = {np.sqrt(1 + 3 * k**2):.3f}")
# gamma (PIC fit) = -0.1377, omega (PIC fit) = 1.412, Bohm-Gross = 1.323

Two things come out of the simulation:

  • The damping rate is within about 10% of Landau’s formula, from 100 000 particles and nothing but Newton’s and Gauss’s laws. The formula is itself the weak-damping approximation; the exact root of the kinetic dispersion relation at \(k\lambda_D = 0.5\) is \(\gamma = -0.1533\). The remaining spread comes from particle noise, which scales as \(1/\sqrt{N}\). Different seed values give fits between about \(-0.14\) and \(-0.17\), and after roughly six peaks the field energy reaches the noise floor, which is why the fit stops there.

  • The frequency is 1.41, not the Bohm-Gross fluid value 1.32. It agrees with the real part of the exact kinetic root, \(\omega = 1.4156\,\omega_{pe}\). At \(k\lambda_D = 0.5\) the fluid approximation is no longer accurate, and the particle simulation, which makes no fluid assumption, shows the kinetic answer.

Where to go next#