Part 3, Chapter 7

The Closure Problem

Truncating the moment hierarchy with physical approximations

7.1 The Infinite Moment Hierarchy

Taking moments of the kinetic equation generates an infinite chain. Each equation for moment of order k contains the moment of order k+1:

$$\text{0th moment (continuity):} \quad \frac{\partial n}{\partial t} + \nabla\cdot(n\mathbf{u}) = 0 \quad\leftarrow\text{ contains } \mathbf{u}$$

$$\text{1st moment (momentum):} \quad mn\frac{d\mathbf{u}}{dt} = \ldots - \nabla\cdot\mathbf{P} \quad\leftarrow\text{ contains } \mathbf{P}$$

$$\text{2nd moment (energy):} \quad \frac{3}{2}n\frac{dT}{dt} = \ldots - \nabla\cdot\mathbf{q} \quad\leftarrow\text{ contains } \mathbf{q}$$

$$\text{3rd moment (heat flux):} \quad \frac{\partial \mathbf{q}}{\partial t} = \ldots - \nabla\cdot\boldsymbol{\mathcal{R}} \quad\leftarrow\text{ contains 4th moment}$$

Why does this happen?

The general moment of order k is:

$$M_{i_1 i_2 \cdots i_k} = m\int w_{i_1} w_{i_2}\cdots w_{i_k}\,f\,d^3v$$

When we take the moment of the advection term v dot nabla f, we get:

$$\int w_{i_1}\cdots w_{i_k}\,\mathbf{v}\cdot\nabla f\,d^3v \quad\longrightarrow\quad \nabla\cdot\int w_{i_1}\cdots w_{i_k}\,\mathbf{w}\,f\,d^3v$$

This is a moment of order k+1. The advection term always couples each moment equation to the next higher moment. This is a fundamental consequence of the velocity-space streaming in kinetic theory.

The Closure Problem: We have N equations but N+1 unknowns at every truncation level. To obtain a closed (solvable) system, we must express the highest-order moment in terms of lower-order moments. This approximation is called a closure.

7.2 Adiabatic Closure

The simplest physically motivated closure: assume zero heat flux (q = 0), meaning no thermal energy is conducted. This is valid when the dynamical timescale is much shorter than the heat conduction timescale.

Step 1: Start from the energy equation with q = 0

$$\frac{3}{2}n\frac{dT}{dt} + p\nabla\cdot\mathbf{u} = 0$$

We also set the viscous stress pi and collisional heating Q to zero.

Step 2: Use the continuity equation

From dn/dt + n nabla dot u = 0, we have nabla dot u = -(1/n)(dn/dt). Substituting:

$$\frac{3}{2}n\frac{dT}{dt} - nk_BT\cdot\frac{1}{n}\frac{dn}{dt} = 0$$

$$\frac{3}{2}\frac{dT}{T} = \frac{dn}{n}$$

Step 3: Integrate

$$\frac{3}{2}\ln T = \ln n + \text{const} \quad\Rightarrow\quad T \propto n^{2/3}$$

Since p = nk_BT, we get:

$$p \propto n \cdot n^{2/3} = n^{5/3}$$

Result: Adiabatic Equation of State

$$\boxed{p = p_0\left(\frac{n}{n_0}\right)^\gamma \quad\text{where}\quad \gamma = \frac{N+2}{N}}$$

For 3D (N=3 degrees of freedom): gamma = 5/3. For 2D: gamma = 2. For 1D: gamma = 3.

Validity condition

$$\tau_{\text{dynamic}} \ll \tau_{\text{heat conduction}} \quad\Leftrightarrow\quad \frac{L/c_s}{L^2/\chi} = \frac{\chi}{Lc_s} \ll 1$$

The process must be fast compared to heat diffusion. Valid for fast magnetosonic waves, shocks.

7.3 Isothermal Closure

The opposite limit: heat conduction is so fast that temperature equilibrates instantly.

Step 1: Set T = constant

$$T(\mathbf{r}, t) = T_0 = \text{const}$$

This means dT/dt = 0, so the energy equation is automatically satisfied if q adjusts to carry away the compression heating.

Step 2: The equation of state

$$p = nk_BT_0 \propto n$$

This is equivalent to gamma = 1 in the polytropic relation:

$$p = p_0\left(\frac{n}{n_0}\right)^1$$

Result: Isothermal Equation of State

$$\boxed{p = nk_BT_0, \quad T = \text{const}, \quad c_s = \sqrt{\frac{k_BT_0}{m}}}$$

Validity condition

$$\tau_{\text{heat conduction}} \ll \tau_{\text{dynamic}} \quad\Leftrightarrow\quad \frac{\chi}{Lc_s} \gg 1$$

Heat conduction must be fast enough to maintain uniform temperature. Valid for electrons along field lines in many fusion plasmas (kappa_parallel is very large).

Implied heat flux

The heat flux is not zero in the isothermal model -- it is whatever is needed to keep T constant. From the energy equation with dT/dt = 0:

$$\nabla\cdot\mathbf{q} = -p\nabla\cdot\mathbf{u}$$

The heat flux must exactly compensate the compressive heating.

7.4 CGL (Chew-Goldberger-Low) Closure

In a magnetized plasma, the pressure is naturally anisotropic because particles gyrate freely perpendicular to B but stream along B. The CGL closure treats the parallel and perpendicular pressures separately.

Step 1: Anisotropic pressure tensor

$$\mathbf{P} = p_\perp(\mathbf{I} - \hat{b}\hat{b}) + p_\parallel\hat{b}\hat{b}$$

where b-hat = B/|B| is the unit vector along the magnetic field. p_perp and p_parallel are the perpendicular and parallel pressures.

Step 2: Physical picture for perpendicular pressure

Perpendicular energy is the kinetic energy of gyromotion. The magnetic moment mu = mv_perp²/(2B) is an adiabatic invariant. For a fluid element:

$$\mu = \frac{m\langle v_\perp^2\rangle}{2B} = \frac{p_\perp}{nB} = \text{const (following the flow)}$$

Step 3: Derive the perpendicular CGL relation

From d(p_perp/(nB))/dt = 0:

$$\frac{d}{dt}\left(\frac{p_\perp}{nB}\right) = 0 \qquad\Rightarrow\qquad \frac{p_\perp}{nB} = \text{const}$$

Step 4: Physical picture for parallel pressure

The parallel motion is 1D adiabatic (along the field line). For 1D adiabatic compression with gamma_1D = 3:

$$p_\parallel \propto n_\parallel^3$$

The "1D density" along a flux tube of cross-section A is n_parallel = n times A. Since flux is conserved (BA = const), A = const/B. So n_parallel = n/B. Therefore:

$$p_\parallel \propto \left(\frac{n}{B}\right)^3 \cdot B^3 \cdot\frac{1}{n^3} \quad\Rightarrow\quad \frac{p_\parallel B^2}{n^3} = \text{const}$$

Result: CGL Double-Adiabatic Equations

$$\boxed{\frac{d}{dt}\left(\frac{p_\perp}{nB}\right) = 0 \qquad \frac{d}{dt}\left(\frac{p_\parallel B^2}{n^3}\right) = 0}$$

These two relations close the MHD equations for anisotropic pressure without needing to specify the heat flux.

Validity and limitations

CGL is valid when: (1) collision frequency is low (collisionless), (2) the magnetic moment is conserved (gradual field variations), and (3) there is no heat flux along field lines. It fails near resonant surfaces and in regions with strong parallel temperature gradients.

7.5 Braginskii Transport Closure

For collisional plasmas where the mean free path is much less than the system size (lambda_mfp much less than L), the distribution function is close to Maxwellian. The Chapman-Enskog expansion provides systematic transport coefficients.

Step 1: Chapman-Enskog expansion

Expand the distribution function about a local Maxwellian:

$$f = f_M\left(1 + \frac{\epsilon}{v_{th}}\,g_1(\mathbf{w}) + \frac{\epsilon^2}{v_{th}^2}\,g_2(\mathbf{w}) + \cdots\right)$$

where epsilon = lambda_mfp/L is the small parameter. The zeroth order f_M gives n, u, T. The first-order correction g₁ gives the transport fluxes (heat flux, viscous stress).

Step 2: Determine g₁ from the kinetic equation

Substituting f = f_M(1 + g₁) into the Boltzmann equation and linearizing gives:

$$C[g_1] = -\frac{Df_M}{Dt}\cdot\frac{1}{f_M}$$

where C is the linearized collision operator and D/Dt contains the streaming and force terms. Solving this integral equation for g₁ yields the transport coefficients.

Heat Flux (Braginskii result)

$$\mathbf{q} = -\kappa_\parallel \nabla_\parallel T - \kappa_\perp \nabla_\perp T + \kappa_\wedge \hat{b} \times \nabla T$$

Parallel conductivity

$$\kappa_\parallel \approx 3.16\frac{nT\tau_e}{m_e}$$

Dominant; unimpeded along B. Scales as T^5/2 through the collision time tau_e.

Perpendicular conductivity

$$\kappa_\perp \approx \frac{4.66\,nT}{m_e\Omega_e^2\tau_e} = \frac{\kappa_\parallel}{(\Omega_e\tau_e)^2}\cdot\text{const}$$

Strongly suppressed by the magnetic field. Ratio kappa_perp/kappa_parallel approximately (Omega_etau_e)^-2.

Diamagnetic (cross) conductivity

$$\kappa_\wedge \approx \frac{5nT}{2m_e\Omega_e}$$

Perpendicular to both B and nabla T. Intermediate magnitude.

Viscous stress tensor (Braginskii)

The viscous stress has five components in a magnetized plasma, reflecting the anisotropy:

$$\boldsymbol{\pi} = -\eta_0 \mathbf{W}_0 - \eta_1 \mathbf{W}_1 - \eta_2 \mathbf{W}_2 - \eta_3 \mathbf{W}_3 - \eta_4 \mathbf{W}_4$$

The dominant coefficient is the parallel viscosity:

$$\eta_0 \approx 0.96\,n_i T_i \tau_i \qquad\text{(ion viscosity along B)}$$

The perpendicular viscosity coefficients eta₁, eta₂ are suppressed by factors of (Omega_itau_i)^-2 and the gyroviscosity eta₃, eta₄ by (Omega_itau_i)^-1.

Validity

Braginskii transport is valid when: lambda_mfp much less than L (collisional regime), and rho_i much less than L (magnetized). The ordering is epsilon = lambda_mfp/L much less than 1. The result is accurate to first order in epsilon.

7.6 Advanced Closures

Gyrokinetic Closure

For low-frequency turbulence with k_perprho_i ~ 1, gyrokinetics provides an ordering:

$$\frac{\omega}{\Omega_i} \sim \frac{k_\parallel}{k_\perp} \sim \frac{\delta B}{B} \sim \frac{e\phi}{T_e} \sim \epsilon \ll 1$$

This reduces the 6D kinetic equation to 5D (averaging over the fast gyromotion), retaining finite Larmor radius effects while eliminating the cyclotron timescale.

Landau Fluid Closure (Hammett-Perkins)

For collisionless plasmas, Landau damping must be captured. The Hammett-Perkins closure models the heat flux as:

$$q_\parallel = -n\chi_\parallel \frac{\partial T}{\partial z} \cdot \text{sgn}(k_\parallel) \quad\text{where}\quad \chi_\parallel = \sqrt{\frac{8}{\pi}}\frac{v_{th}}{|k_\parallel|}$$

This non-local closure (involving |k| in Fourier space) mimics the phase-mixing physics of Landau damping within a fluid framework. The key insight is that Landau damping is equivalent to a velocity-dependent diffusion in phase space.

7.7 Comparison of Closures

Closure	Equation of State	Valid When	Sound Speed
Cold	p = 0	omega much greater than kv_th	c_s = 0
Isothermal (gamma=1)	p = nkT₀	Fast heat conduction	c_s = (kT/m)^1/2
Adiabatic 3D (gamma=5/3)	p proportional to n^5/3	No heat flux	c_s = (5kT/3m)^1/2
CGL	p_perp/(nB) = const	Collisionless, magnetized	Depends on angle
Braginskii	q = -kappa nabla T	lambda_mfp much less than L	Modified by transport
Landau fluid	Non-local q	Collisionless, k_parallelv_th ~ omega	Complex (damped)

Interactive Simulation

Compare Sound Wave Speeds Under Different Closures

Python

script.py194 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Compare wave propagation under different closures
# Solve the linearized 1D Euler equations: dn/dt + n0*du/dx = 0
# and n0*m*du/dt + dp/dx = 0
# with different p(n) relations

# Parameters (normalized)
n0 = 1.0       # background density
T0 = 1.0       # background temperature (in energy units, kT)
m = 1.0        # mass
p0 = n0 * T0   # background pressure

# Grid
Nx = 800
Lx = 20.0
dx = Lx / Nx
x = np.linspace(0, Lx, Nx)

# Time
dt = 0.005
Nt = 2000

# Sound speeds for different closures
gamma_values = {
    'Cold (gamma->0)': 0.001,
    'Isothermal (gamma=1)': 1.0,
    'Adiabatic 1D (gamma=3)': 3.0,
    'Adiabatic 3D (gamma=5/3)': 5.0/3.0,
    'Adiabatic 2D (gamma=2)': 2.0,
}

colors = {
    'Cold (gamma->0)': 'gray',
    'Isothermal (gamma=1)': 'cyan',
    'Adiabatic 1D (gamma=3)': 'red',
    'Adiabatic 3D (gamma=5/3)': 'orange',
    'Adiabatic 2D (gamma=2)': 'lime',
}

# Initial perturbation: Gaussian density pulse
sigma = 0.3
n_pert = 0.1 * np.exp(-((x - Lx/2)**2) / (2*sigma**2))

results = {}

for label, gamma in gamma_values.items():
    cs = np.sqrt(gamma * T0 / m)

# Initialize
    n1 = n_pert.copy()   # density perturbation
    u1 = np.zeros(Nx)    # velocity perturbation

# Time evolution using Lax-Wendroff scheme
    for step in range(Nt):
        # Pressure perturbation depends on closure
        if gamma < 0.01:
            dp = np.zeros(Nx)  # cold plasma
        else:
            dp = gamma * T0 * n1  # linearized: dp = gamma * (p0/n0) * dn = gamma*T*dn

# Spatial derivatives (central differences)
        dn_new = np.zeros(Nx)
        du_new = np.zeros(Nx)
        for i in range(1, Nx-1):
            # dn/dt = -n0 * du/dx
            du_dx = (u1[i+1] - u1[i-1]) / (2*dx)
            dn_new[i] = n1[i] - dt * n0 * du_dx
            # n0*m*du/dt = -dp/dx
            dp_dx = (dp[i+1] - dp[i-1]) / (2*dx)
            du_new[i] = u1[i] - dt * dp_dx / (n0 * m)

# Add numerical diffusion for stability
        nu_diff = 0.001
        for i in range(1, Nx-1):
            dn_new[i] += nu_diff * (n1[i+1] - 2*n1[i] + n1[i-1])
            du_new[i] += nu_diff * (u1[i+1] - 2*u1[i] + u1[i-1])

n1 = dn_new
        u1 = du_new

results[label] = (n1.copy(), cs)

# ---- Plotting ----
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Final density perturbation for all closures
ax = axes[0, 0]
for label in gamma_values:
    n1_final, cs = results[label]
    ax.plot(x, n1_final / 0.1, color=colors[label], linewidth=2, label=f'{label}, cs={cs:.2f}')
ax.plot(x, n_pert / 0.1, 'w--', linewidth=1, alpha=0.3, label='Initial')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('delta_n / delta_n_0', fontsize=12)
ax.set_title('Density Perturbation at t=10', fontsize=14, color='white')
ax.legend(fontsize=8, loc='upper right')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# Plot 2: Sound speed vs gamma
ax = axes[0, 1]
gamma_range = np.linspace(0.01, 4.0, 200)
cs_range = np.sqrt(gamma_range * T0 / m)
ax.plot(gamma_range, cs_range, 'c-', linewidth=2)
for label, gamma in gamma_values.items():
    cs = np.sqrt(gamma * T0 / m)
    ax.plot(gamma, cs, 'o', color=colors[label], markersize=10, label=label)
ax.set_xlabel('gamma (adiabatic index)', fontsize=12)
ax.set_ylabel('Sound speed cs', fontsize=12)
ax.set_title('Sound Speed vs Closure', fontsize=14, color='white')
ax.legend(fontsize=8)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# Plot 3: CGL anisotropic pressure under compression
ax = axes[1, 0]
B_ratio = np.linspace(0.5, 3.0, 200)  # B/B0
n_ratio = np.linspace(0.5, 3.0, 200)  # n/n0

# For a 1D compression along B: n increases, B stays constant
# p_perp/p0 = (n/n0)*(B/B0) and p_par/p0 = (n/n0)^3 / (B/B0)^2
# If B ~ n (perpendicular compression): p_perp ~ n^2, p_par ~ n

# Case 1: perpendicular compression (B increases with n)
n_r = np.linspace(0.5, 3.0, 200)
p_perp_perp_comp = n_r**2  # p_perp/p0 when B ~ n
p_par_perp_comp = n_r      # p_par/p0 when B ~ n

# Case 2: parallel compression (B constant)
p_perp_par_comp = n_r       # p_perp/p0 when B = B0
p_par_par_comp = n_r**3     # p_par/p0 when B = B0

# Isotropic adiabatic for comparison
p_iso = n_r**(5.0/3.0)

ax.plot(n_r, p_perp_perp_comp, 'c-', linewidth=2, label='CGL p_perp (perp compress)')
ax.plot(n_r, p_par_perp_comp, 'c--', linewidth=2, label='CGL p_par (perp compress)')
ax.plot(n_r, p_perp_par_comp, 'orange', linewidth=2, label='CGL p_perp (par compress)')
ax.plot(n_r, p_par_par_comp, 'orange', linewidth=2, linestyle='--', label='CGL p_par (par compress)')
ax.plot(n_r, p_iso, 'w--', linewidth=1.5, label='Isotropic gamma=5/3')
ax.set_xlabel('n / n_0', fontsize=12)
ax.set_ylabel('p / p_0', fontsize=12)
ax.set_title('CGL vs Isotropic Pressure', fontsize=14, color='white')
ax.legend(fontsize=8)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_ylim(0, 6)

# Plot 4: Braginskii anisotropy kappa_par/kappa_perp vs Omega*tau
ax = axes[1, 1]
Omega_tau = np.logspace(-1, 3, 200)
kappa_ratio = Omega_tau**2  # kappa_par/kappa_perp ~ (Omega*tau)^2
kappa_cross_ratio = Omega_tau  # kappa_wedge/kappa_perp ~ Omega*tau

ax.loglog(Omega_tau, kappa_ratio, 'cyan', linewidth=2, label='kappa_par / kappa_perp')
ax.loglog(Omega_tau, kappa_cross_ratio, 'orange', linewidth=2, label='kappa_wedge / kappa_perp')
ax.loglog(Omega_tau, np.ones_like(Omega_tau), 'w--', linewidth=1, alpha=0.5)
ax.axvline(1.0, color='red', linestyle=':', alpha=0.5, label='Omega*tau = 1')
ax.fill_between(Omega_tau, 0.01, 1e8, where=Omega_tau > 1, alpha=0.05, color='cyan')
ax.set_xlabel('Omega_c * tau (magnetization parameter)', fontsize=12)
ax.set_ylabel('Conductivity ratio', fontsize=12)
ax.set_title('Braginskii Transport Anisotropy', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_ylim(0.01, 1e7)

for ax in axes.flat:
    ax.spines['bottom'].set_color('white')
    ax.spines['left'].set_color('white')

fig.patch.set_facecolor('#111111')
plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#111111')

# Print summary
print("=== Closure Comparison: Sound Speeds ===")
for label, gamma in sorted(gamma_values.items(), key=lambda x: x[1]):
    cs = np.sqrt(gamma * T0 / m)
    print(f"  {label:30s}: cs = {cs:.4f} (gamma = {gamma:.3f})")

print(f"\n=== CGL Double-Adiabatic Relations ===")
print(f"  Perpendicular: p_perp / (n*B) = const")
print(f"  Parallel:      p_par * B^2 / n^3 = const")

print(f"\n=== Braginskii Anisotropy ===")
for Ot in [0.1, 1.0, 10.0, 100.0]:
    print(f"  Omega*tau = {Ot:6.1f}: kappa_par/kappa_perp ~ {Ot**2:.0f}")

print("\nPlot saved to output.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

-- The moment hierarchy is infinite: each equation introduces the next higher moment
-- Adiabatic closure (q=0) gives p proportional to n^gamma with gamma = (N+2)/N; valid for fast processes
-- Isothermal closure (T=const) gives p proportional to n; valid when heat conduction is dominant
-- CGL double-adiabatic closure handles anisotropic pressure in collisionless magnetized plasmas
-- Braginskii closure provides full transport coefficients for collisional plasmas via Chapman-Enskog expansion
-- The choice of closure profoundly affects the wave speeds and stability properties of the plasma model

← Two-Fluid Theory Next: Magnetohydrodynamics →

The Closure Problem

7.1 The Infinite Moment Hierarchy

Why does this happen?

7.2 Adiabatic Closure

Step 1: Start from the energy equation with q = 0

Step 2: Use the continuity equation

Step 3: Integrate

Result: Adiabatic Equation of State

Validity condition

7.3 Isothermal Closure

Step 1: Set T = constant

Step 2: The equation of state

Result: Isothermal Equation of State

Validity condition

Implied heat flux

7.4 CGL (Chew-Goldberger-Low) Closure

Step 1: Anisotropic pressure tensor

Step 2: Physical picture for perpendicular pressure

Step 3: Derive the perpendicular CGL relation

Step 4: Physical picture for parallel pressure

Result: CGL Double-Adiabatic Equations

Validity and limitations

7.5 Braginskii Transport Closure

Step 1: Chapman-Enskog expansion

Step 2: Determine g1 from the kinetic equation

Heat Flux (Braginskii result)

Parallel conductivity

Perpendicular conductivity

Diamagnetic (cross) conductivity

Viscous stress tensor (Braginskii)

Validity

7.6 Advanced Closures

Gyrokinetic Closure

Landau Fluid Closure (Hammett-Perkins)

7.7 Comparison of Closures

Interactive Simulation

Compare Sound Wave Speeds Under Different Closures

Key Takeaways

Step 2: Determine g₁ from the kinetic equation