Part 3, Chapter 1

Moment Equations

Deriving fluid equations from kinetic theory through velocity moments

1.1 From Kinetic to Fluid Description

The kinetic description via the Boltzmann equation contains complete information through the distribution function f(r, v, t). By taking velocity moments, we derive fluid equations for macroscopic quantities like density, velocity, and pressure.

The Boltzmann equation for species s reads:

$$\frac{\partial f_s}{\partial t} + \mathbf{v}\cdot\nabla f_s + \frac{q_s}{m_s}(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\nabla_v f_s = \left(\frac{\delta f_s}{\delta t}\right)_{\text{coll}}$$

The general procedure is to multiply the Boltzmann equation by a velocity-dependent quantity g(v) and integrate over all of velocity space:

$$\int g(\mathbf{v}) \left[\frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla f + \frac{q}{m}(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\nabla_v f\right] d^3v = \int g(\mathbf{v})\left(\frac{\delta f}{\delta t}\right)_{\text{coll}} d^3v$$

Choosing g = 1, v, and mv²/2 yields the continuity, momentum, and energy equations respectively. We will derive each one step by step.

1.2 Macroscopic Quantities

Before deriving the moment equations, we define the macroscopic quantities as velocity-space integrals of the distribution function. We also introduce the random (thermal) velocity w = v - u.

Number Density (Zeroth Moment)

$$n(\mathbf{r}, t) = \int f(\mathbf{r}, \mathbf{v}, t) \, d^3v$$

Mean (Fluid) Velocity (First Moment)

$$\mathbf{u}(\mathbf{r}, t) = \frac{1}{n} \int \mathbf{v} \, f \, d^3v = \langle \mathbf{v} \rangle$$

Pressure Tensor (Second Moment)

$$P_{ij} = m \int (v_i - u_i)(v_j - u_j) f \, d^3v = m \int w_i w_j f \, d^3v$$

The scalar pressure is p = (1/3)Tr(P) = nk_BT. The traceless part is the viscous stress tensor pi_ij.

Heat Flux Vector (Third Moment)

$$\mathbf{q} = \frac{m}{2} \int |\mathbf{w}|^2 \mathbf{w} \, f \, d^3v$$

Energy flux carried by random thermal motion relative to the mean flow.

1.3 Zeroth Moment: Continuity Equation

Set g(v) = 1 and integrate the Boltzmann equation over velocity space. We handle each term individually.

Step 1: Time derivative term

$$\int \frac{\partial f}{\partial t}\,d^3v = \frac{\partial}{\partial t}\int f\,d^3v = \frac{\partial n}{\partial t}$$

We can interchange the integration and time derivative since the velocity integration limits are fixed (all of velocity space).

Step 2: Advection term

$$\int \mathbf{v}\cdot\nabla f\,d^3v = \nabla\cdot\int \mathbf{v} f\,d^3v = \nabla\cdot(n\mathbf{u})$$

The spatial gradient passes through the velocity integral. We used the definition of the mean velocity nu = integral of vf d3v.

Step 3: Lorentz force term

$$\int \frac{q}{m}(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\nabla_v f\,d^3v$$

Integrate by parts: since f vanishes as |v| tends to infinity,

$$= -\frac{q}{m}\int f\,\nabla_v\cdot(\mathbf{E}+\mathbf{v}\times\mathbf{B})\,d^3v$$

Now E is independent of v, so its divergence in velocity space vanishes. For the magnetic term:

$$\nabla_v\cdot(\mathbf{v}\times\mathbf{B}) = \sum_i \frac{\partial}{\partial v_i}(\mathbf{v}\times\mathbf{B})_i = 0$$

This vanishes because (v x B)_i does not depend on v_i (the cross product is perpendicular). Therefore the entire Lorentz force term contributes zero.

Step 4: Collision term

$$\int \left(\frac{\delta f}{\delta t}\right)_{\text{coll}} d^3v = 0$$

Collisions conserve particle number (they can redistribute particles in velocity space but cannot create or destroy them), so this integral vanishes.

Result: Continuity Equation

Combining all four terms:

$$\boxed{\frac{\partial n}{\partial t} + \nabla\cdot(n\mathbf{u}) = 0}$$

This is the conservation of particle number. In terms of mass density rho = mn, it becomes the familiar mass continuity equation.

1.4 First Moment: Momentum Equation

Set g(v) = mv and integrate over velocity space. We decompose v = u + w where w is the random velocity.

Step 1: Time derivative term

$$\int m\mathbf{v}\frac{\partial f}{\partial t}\,d^3v = m\frac{\partial}{\partial t}\int \mathbf{v}f\,d^3v = m\frac{\partial(n\mathbf{u})}{\partial t}$$

Step 2: Advection term

We write v = u + w and expand:

$$\int m\mathbf{v}(\mathbf{v}\cdot\nabla f)\,d^3v = m\nabla\cdot\int \mathbf{v}\mathbf{v}\,f\,d^3v$$

The dyadic vv = (u+w)(u+w) = uu + uw + wu + ww. Taking the integral:

$$m\int \mathbf{v}\mathbf{v}\,f\,d^3v = mn\mathbf{u}\mathbf{u} + m\int \mathbf{w}\mathbf{w}\,f\,d^3v = mn\mathbf{u}\mathbf{u} + \mathbf{P}$$

The cross terms uw and wu vanish because the average of w is zero by definition. Therefore this term gives nabla dot (mnuu + P).

Step 3: Electric field term

$$\int m\mathbf{v}\cdot\frac{q}{m}\mathbf{E}\cdot\nabla_v f\,d^3v = q\mathbf{E}\int \nabla_v f\cdot d^3v \cdot\mathbf{v}$$

Integrating by parts (boundary terms vanish as f goes to zero at infinity):

$$= -q\int f\nabla_v\cdot(\mathbf{v}\mathbf{E}^T)\,d^3v \quad\Rightarrow\quad \text{applying } \frac{\partial v_j}{\partial v_i} = \delta_{ij}$$

$$= -q\int f\,\mathbf{E}\,d^3v = -qn\mathbf{E}$$

With the sign convention in the Boltzmann equation, this gives +qnE on the right-hand side.

Step 4: Magnetic field term

Similarly integrating by parts for the magnetic force:

$$\int m v_j \frac{q}{m}(\mathbf{v}\times\mathbf{B})_i \frac{\partial f}{\partial v_i}\,d^3v = -q\int f\frac{\partial}{\partial v_i}\left[v_j(\mathbf{v}\times\mathbf{B})_i\right]d^3v$$

Using the identity that the divergence of v x B in velocity space vanishes (shown in Step 3 of zeroth moment), plus the term from the derivative of v_j:

$$= -q\int f\,(\mathbf{v}\times\mathbf{B})_j\,d^3v = -qn(\mathbf{u}\times\mathbf{B})_j$$

This gives +qn(u x B) on the right-hand side.

Step 5: Collision term

$$\int m\mathbf{v}\left(\frac{\delta f}{\delta t}\right)_{\text{coll}} d^3v = \mathbf{R}_s$$

Collisions conserve total momentum but not individual species momentum. R_s is the friction force from collisions with other species. For electron-ion collisions: R_e = -R_i.

Step 6: Simplify with the continuity equation

Combining Steps 1 and 2 gives the left-hand side m d(nu)/dt + nabla dot P. Using the product rule on the time derivative:

$$m\frac{\partial(n\mathbf{u})}{\partial t} + m\nabla\cdot(n\mathbf{u}\mathbf{u}) = mn\frac{\partial \mathbf{u}}{\partial t} + m\mathbf{u}\frac{\partial n}{\partial t} + mn(\mathbf{u}\cdot\nabla)\mathbf{u} + m\mathbf{u}\nabla\cdot(n\mathbf{u})$$

The second and fourth terms combine to give u times the continuity equation, which is zero:

$$m\mathbf{u}\left[\frac{\partial n}{\partial t} + \nabla\cdot(n\mathbf{u})\right] = 0$$

Result: Momentum Equation

$$\boxed{mn\left(\frac{\partial\mathbf{u}}{\partial t} + \mathbf{u}\cdot\nabla\mathbf{u}\right) = nq(\mathbf{E}+\mathbf{u}\times\mathbf{B}) - \nabla\cdot\mathbf{P} + \mathbf{R}}$$

The left side is the convective derivative of the fluid momentum. The right side has the Lorentz force, the divergence of the pressure tensor (including both scalar pressure gradient and viscous stress), and the collisional friction force.

1.5 Second Moment: Energy Equation

Set g(v) = mv²/2 and integrate. The derivation is more involved; we sketch the key steps.

Step 1: Identify the kinetic energy density

The total kinetic energy density splits into bulk flow and thermal parts using v = u + w:

$$\frac{m}{2}\int v^2 f\,d^3v = \frac{1}{2}mnu^2 + \frac{m}{2}\int w^2 f\,d^3v = \frac{1}{2}mnu^2 + \frac{3}{2}nk_BT$$

We identify the thermal energy density as (3/2)nk_BT = (3/2)p where p = nk_BT is the scalar pressure.

Step 2: Time derivative of total energy

$$\frac{\partial}{\partial t}\left(\frac{1}{2}mnu^2 + \frac{3}{2}p\right)$$

Step 3: Advection of energy and the heat flux

The advection term yields:

$$\nabla\cdot\left[\left(\frac{1}{2}mnu^2 + \frac{5}{2}p\right)\mathbf{u} + \boldsymbol{\pi}\cdot\mathbf{u} + \mathbf{q}\right]$$

The 5/2 factor (rather than 3/2) arises because the pressure does pdV work on the flow. The heat flux q is the third moment in the random velocity.

Step 4: Electromagnetic work

The Lorentz force contributes through work done on the fluid:

$$\int \frac{mv^2}{2}\frac{q}{m}(\mathbf{E}+\mathbf{v}\times\mathbf{B})\cdot\nabla_v f\,d^3v = -nq\mathbf{u}\cdot\mathbf{E}$$

The magnetic force does no work (v x B is perpendicular to v). Only the electric field does work on the plasma.

Step 5: Subtract the bulk kinetic energy equation

Taking the dot product of u with the momentum equation gives the bulk kinetic energy equation. Subtracting this from the total energy equation yields the internal (thermal) energy equation:

$$\frac{3}{2}\frac{\partial p}{\partial t} + \frac{3}{2}\nabla\cdot(p\mathbf{u}) + p\nabla\cdot\mathbf{u} + \nabla\cdot\mathbf{q} + \boldsymbol{\pi}:\nabla\mathbf{u} = Q$$

Where Q represents collisional heating (energy exchange between species).

Result: Energy Equation (in terms of temperature)

Using p = nk_BT and the continuity equation to simplify:

$$\boxed{\frac{3}{2}n\frac{dT}{dt} + p\nabla\cdot\mathbf{u} + \nabla\cdot\mathbf{q} + \boldsymbol{\pi}:\nabla\mathbf{u} = Q}$$

Each term has a clear physical meaning: (1) rate of change of thermal energy following the fluid, (2) compressive heating/cooling (pdV work), (3) divergence of heat flux, (4) viscous heating, and (5) collisional energy transfer.

Closure Problem: The energy equation contains the heat flux q (a third-order moment) and the viscous stress tensor pi (part of the second moment). The evolution equation for q would introduce the fourth moment, and so on ad infinitum. This infinite hierarchy must be truncated with a closure approximation.

1.6 Summary: The Moment Hierarchy

The three moment equations form a hierarchy. We collect them here for reference:

Continuity (zeroth moment: multiply by 1)

$$\frac{\partial n}{\partial t} + \nabla \cdot (n\mathbf{u}) = 0$$

1 equation, introduces u (from first moment)

Momentum (first moment: multiply by mv)

$$mn\left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) = qn(\mathbf{E} + \mathbf{u} \times \mathbf{B}) - \nabla \cdot \mathbf{P} + \mathbf{R}$$

3 equations (vector), introduces P (from second moment)

Energy (second moment: multiply by mv²/2)

$$\frac{3}{2}n\frac{dT}{dt} + p\nabla\cdot\mathbf{u} + \nabla\cdot\mathbf{q} + \boldsymbol{\pi}:\nabla\mathbf{u} = Q$$

1 equation (scalar for isotropic case), introduces q (from third moment)

Counting unknowns: We have n, u (3 components), P (6 independent components of symmetric tensor), q (3 components) = 13 unknowns, but only 1 + 3 + 6 = 10 equations (if we take the full tensor energy equation). We always need a closure relation for the highest moment.

1.7 Common Closure Approximations

Cold Plasma (T = 0)

$$p = 0, \quad \mathbf{q} = 0$$

Truncates after continuity+momentum. Valid when thermal speed is much less than phase speed.

Isothermal

$$p = nk_BT, \quad T = \text{const}$$

Valid when thermal equilibration is much faster than dynamics.

Adiabatic

$$\mathbf{q} = 0 \quad\Rightarrow\quad p = p_0\left(\frac{n}{n_0}\right)^\gamma$$

No heat flux; valid for fast processes where thermal conduction is negligible.

Braginskii

$$\mathbf{q} = -\kappa_\parallel \nabla_\parallel T - \kappa_\perp\nabla_\perp T$$

Full transport theory for collisional magnetized plasmas.

Interactive Simulation

Verify Moment Equations: Integrate a Maxwellian Distribution

Python

script.py128 lines

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

# Define Maxwellian distribution function (1D for clarity)
# f(v) = n * (m / (2 pi kT))^(1/2) * exp(-m v^2 / (2 kT))
# We use shifted Maxwellian with drift velocity u0

n0 = 1e19       # density [m^-3]
T0 = 100.0      # temperature [eV]
kT = T0 * 1.6e-19  # in Joules
m = 1.67e-27    # proton mass [kg]
u0 = 1e5        # drift velocity [m/s]

vth = np.sqrt(kT / m)  # thermal speed
print(f"Thermal speed: vth = {vth:.2e} m/s")
print(f"Drift velocity: u0 = {u0:.2e} m/s")

# Velocity grid: wide enough to capture the tail
v = np.linspace(-8*vth + u0, 8*vth + u0, 10000)
dv = v[1] - v[0]

# 1D shifted Maxwellian
f = n0 * np.sqrt(m / (2 * np.pi * kT)) * np.exp(-m * (v - u0)**2 / (2 * kT))

# ---- Compute moments by numerical integration ----
# Zeroth moment: density
n_calc = np.trapezoid(f, v)
print(f"\nZeroth moment (density):")
print(f"  Analytical: n = {n0:.4e}")
print(f"  Numerical:  n = {n_calc:.4e}")
print(f"  Relative error: {abs(n_calc - n0)/n0:.2e}")

# First moment: mean velocity
u_calc = np.trapezoid(v * f, v) / n_calc
print(f"\nFirst moment (mean velocity):")
print(f"  Analytical: u = {u0:.4e}")
print(f"  Numerical:  u = {u_calc:.4e}")
print(f"  Relative error: {abs(u_calc - u0)/u0:.2e}")

# Second moment: pressure (using random velocity w = v - u)
w = v - u_calc
p_calc = m * np.trapezoid(w**2 * f, v)
p_analytical = n0 * kT
print(f"\nSecond moment (pressure):")
print(f"  Analytical: p = n*kT = {p_analytical:.4e} Pa")
print(f"  Numerical:  p = {p_calc:.4e} Pa")
print(f"  Relative error: {abs(p_calc - p_analytical)/p_analytical:.2e}")

# Temperature from pressure
T_calc = p_calc / (n_calc * 1.6e-19)  # convert back to eV
print(f"  Temperature: T = {T_calc:.2f} eV (expected {T0:.2f} eV)")

# Third moment: heat flux (should be zero for Maxwellian)
q_calc = 0.5 * m * np.trapezoid(w**3 * f, v)
print(f"\nThird moment (heat flux):")
print(f"  Analytical: q = 0 (Maxwellian is symmetric)")
print(f"  Numerical:  q = {q_calc:.4e} W/m^2")

# ---- Verify continuity equation numerically ----
# For a static Maxwellian, dn/dt = 0 and div(nu) = 0
# Let's verify that integral of v * df/dx = d(nu)/dx for a spatially varying density
print(f"\n--- Verify continuity: spatial gradient ---")
# n(x) = n0 * (1 + 0.1*sin(kx)), at x=0: dn/dx = 0.1*n0*k
kx = 2.0  # wavenumber
# df/dx at x=0 for n(x)*Maxwellian = (dn/dx)*Maxwellian_normalized
dn_dx = 0.1 * n0 * kx
flux_gradient = dn_dx * u0  # d(nu)/dx at x=0
advection_integral = np.trapezoid(v * f, v) * (0.1 * kx)  # scale by dn/dx / n0
print(f"  d(nu)/dx = {flux_gradient:.4e}")
print(f"  integral(v * df/dx) = {advection_integral:.4e}")
print(f"  These should match (continuity equation verified).")

# ---- Plot the distribution and moments ----
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Plot 1: Distribution function
axes[0].plot(v/vth, f * vth / n0, 'c-', linewidth=2)
axes[0].axvline(u0/vth, color='orange', linestyle='--', label=f'u = {u0:.1e} m/s')
axes[0].set_xlabel('v / v_th', fontsize=12)
axes[0].set_ylabel('f(v) * v_th / n', fontsize=12)
axes[0].set_title('Shifted Maxwellian', fontsize=14, color='white')
axes[0].legend(fontsize=10)
axes[0].set_facecolor('#0a0a0a')
axes[0].tick_params(colors='white')
axes[0].spines['bottom'].set_color('white')
axes[0].spines['left'].set_color('white')

# Plot 2: Integrands for moments
axes[1].plot(v/vth, f/n0*vth, 'c-', linewidth=2, label='f (density)')
axes[1].plot(v/vth, v*f/(n0*u0*vth) * vth, 'orange', linewidth=2, label='v*f (momentum)')
axes[1].plot(v/vth, (v-u_calc)**2*f/(n0*vth**2) * vth, 'lime', linewidth=2, label='w^2*f (pressure)')
axes[1].set_xlabel('v / v_th', fontsize=12)
axes[1].set_ylabel('Integrand (normalized)', fontsize=12)
axes[1].set_title('Moment Integrands', fontsize=14, color='white')
axes[1].legend(fontsize=9)
axes[1].set_facecolor('#0a0a0a')
axes[1].tick_params(colors='white')
axes[1].spines['bottom'].set_color('white')
axes[1].spines['left'].set_color('white')

# Plot 3: Moment values as bar chart
moments = ['n (density)', 'u (velocity)', 'p (pressure)', 'q (heat flux)']
analytical = [n0, u0, p_analytical, 0.0]
numerical = [n_calc, u_calc, p_calc, q_calc]
# Normalize for display
norm_ana = [1.0, 1.0, 1.0, 0.0]
norm_num = [n_calc/n0, u_calc/u0, p_calc/p_analytical, 0.0]
x_pos = np.arange(len(moments))
axes[2].bar(x_pos - 0.15, norm_ana, 0.3, color='cyan', alpha=0.7, label='Analytical')
axes[2].bar(x_pos + 0.15, norm_num, 0.3, color='orange', alpha=0.7, label='Numerical')
axes[2].set_xticks(x_pos)
axes[2].set_xticklabels(['n', 'u', 'p', 'q'], fontsize=12)
axes[2].set_ylabel('Value / Analytical', fontsize=12)
axes[2].set_title('Moment Verification', fontsize=14, color='white')
axes[2].legend(fontsize=10)
axes[2].set_ylim(0, 1.3)
axes[2].set_facecolor('#0a0a0a')
axes[2].tick_params(colors='white')
axes[2].spines['bottom'].set_color('white')
axes[2].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("\nPlot saved to output.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

-- Fluid equations are velocity moments of the Boltzmann/Vlasov equation
-- Zeroth moment (multiply by 1) gives the continuity equation for particle conservation
-- First moment (multiply by mv) gives the momentum equation with Lorentz force, pressure gradient, and friction
-- Second moment (multiply by mv^2/2) gives the energy equation with compressive heating, heat flux, and viscous dissipation
-- Each moment equation introduces the next higher moment, creating an infinite hierarchy
-- The system must be closed by approximating the highest moment (cold, isothermal, adiabatic, or Braginskii closure)

← Part III Overview Next: MHD Equations →

Moment Equations

1.1 From Kinetic to Fluid Description

1.2 Macroscopic Quantities

Number Density (Zeroth Moment)

Mean (Fluid) Velocity (First Moment)

Pressure Tensor (Second Moment)

Heat Flux Vector (Third Moment)

1.3 Zeroth Moment: Continuity Equation

Step 1: Time derivative term

Step 2: Advection term

Step 3: Lorentz force term

Step 4: Collision term

Result: Continuity Equation

1.4 First Moment: Momentum Equation

Step 1: Time derivative term

Step 2: Advection term

Step 3: Electric field term

Step 4: Magnetic field term

Step 5: Collision term

Step 6: Simplify with the continuity equation

Result: Momentum Equation

1.5 Second Moment: Energy Equation

Step 1: Identify the kinetic energy density

Step 2: Time derivative of total energy

Step 3: Advection of energy and the heat flux

Step 4: Electromagnetic work

Step 5: Subtract the bulk kinetic energy equation

Result: Energy Equation (in terms of temperature)

1.6 Summary: The Moment Hierarchy

Continuity (zeroth moment: multiply by 1)

Momentum (first moment: multiply by mv)

Energy (second moment: multiply by mv2/2)

1.7 Common Closure Approximations

Cold Plasma (T = 0)

Isothermal

Adiabatic

Braginskii

Interactive Simulation

Verify Moment Equations: Integrate a Maxwellian Distribution

Key Takeaways

Energy (second moment: multiply by mv²/2)