Part 3, Chapter 2

MHD Equations

Ideal and resistive magnetohydrodynamics derived from two-fluid theory

2.1 The MHD Approximation

MHD treats plasma as a single conducting fluid. This is obtained from the two-fluid equations by summing over species and making several key approximations:

Quasi-neutrality: n_e = Zn_i = n (charge separation negligible on scales L much greater than the Debye length)
Low frequency: omega much less than Omega_ci (slower than ion cyclotron frequency)
Large scales: L much greater than rho_i, lambda_D, c/omega_pi
Non-relativistic: u much less than c (displacement current negligible)
Small electron mass: m_e/m_i much less than 1 (electron inertia negligible)

2.2 Mass Density and Center-of-Mass Velocity

We define the single-fluid (MHD) variables by summing over the two species:

Step 1: Mass density

$$\rho = m_i n_i + m_e n_e = m_i n + m_e n \approx m_i n$$

Since m_e/m_i is approximately 1/1836, the electron mass contribution is negligible.

Step 2: Center-of-mass velocity

$$\rho\mathbf{v} = m_i n_i \mathbf{u}_i + m_e n_e \mathbf{u}_e \approx m_i n \mathbf{u}_i$$

So the MHD velocity is essentially the ion velocity: v approximately equals u_i.

Step 3: Current density

$$\mathbf{J} = en(\mathbf{u}_i - \mathbf{u}_e)$$

This allows us to express the electron velocity as u_e = u_i - J/(en) = v - J/(en).

Step 4: Total pressure

$$p = p_i + p_e = n_i k_B T_i + n_e k_B T_e = n k_B(T_i + T_e)$$

In MHD we work with a single scalar pressure that is the sum of electron and ion pressures.

2.3 MHD Continuity Equation

Add the two species continuity equations, weighted by their masses:

Step 1: Ion continuity

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

Step 2: Electron continuity

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

Step 3: Sum and use definitions

Adding these and using rho = m_in_i + m_en_e and rho v = m_in_iu_i + m_en_eu_e:

Result: MHD Mass Continuity

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

2.4 MHD Momentum Equation

Add the ion and electron momentum equations. The collision terms cancel (R_ei = -R_ie by Newton's third law).

Step 1: Write out both momentum equations

$$m_i n\left(\frac{\partial \mathbf{u}_i}{\partial t} + \mathbf{u}_i\cdot\nabla\mathbf{u}_i\right) = en(\mathbf{E}+\mathbf{u}_i\times\mathbf{B}) - \nabla p_i + \mathbf{R}_{ie}$$

$$m_e n\left(\frac{\partial \mathbf{u}_e}{\partial t} + \mathbf{u}_e\cdot\nabla\mathbf{u}_e\right) = -en(\mathbf{E}+\mathbf{u}_e\times\mathbf{B}) - \nabla p_e + \mathbf{R}_{ei}$$

Step 2: Add the equations

On the left side, we neglect electron inertia (m_e much less than m_i), giving rho dv/dt. On the right side:

$$en(\mathbf{u}_i - \mathbf{u}_e)\times\mathbf{B} = \mathbf{J}\times\mathbf{B}$$

The electric field terms en(E) - en(E) cancel due to quasi-neutrality. The collision terms cancel: R_ie + R_ei = 0.

Step 3: Combine pressures

$$-\nabla p_i - \nabla p_e = -\nabla(p_i + p_e) = -\nabla p$$

Result: MHD Momentum Equation

$$\boxed{\rho\frac{d\mathbf{v}}{dt} = \rho\left(\frac{\partial\mathbf{v}}{\partial t} + \mathbf{v}\cdot\nabla\mathbf{v}\right) = \mathbf{J}\times\mathbf{B} - \nabla p}$$

The MHD momentum equation has only the J x B force and pressure gradient -- the electric field has been eliminated! This is a profound simplification.

2.5 Ohm's Law from Electron Momentum

The MHD Ohm's law is derived from the electron momentum equation alone (not the sum). Since m_e is small, the electron inertia is negligible.

Step 1: Electron momentum with negligible inertia

$$0 \approx -en(\mathbf{E}+\mathbf{u}_e\times\mathbf{B}) - \nabla p_e + \mathbf{R}_{ei}$$

We set the left side to zero because m_en du_e/dt is negligible compared to the electromagnetic forces.

Step 2: Express u_e in terms of MHD variables

Using u_e = v - J/(en):

$$\mathbf{u}_e\times\mathbf{B} = \mathbf{v}\times\mathbf{B} - \frac{\mathbf{J}\times\mathbf{B}}{en}$$

Step 3: Model the friction force

The electron-ion friction is proportional to the relative drift:

$$\mathbf{R}_{ei} = -m_e n \nu_{ei}(\mathbf{u}_e - \mathbf{u}_i) = m_e n\nu_{ei}\frac{\mathbf{J}}{en} = \frac{m_e\nu_{ei}}{ne^2}\mathbf{J} \equiv \eta_{res}\mathbf{J}$$

Here eta is the Spitzer resistivity, arising from electron-ion collisions.

Step 4: Solve for E

Rearranging the electron momentum equation:

$$\mathbf{E} + \mathbf{v}\times\mathbf{B} = \eta\mathbf{J} + \frac{1}{en}\left(\mathbf{J}\times\mathbf{B} - \nabla p_e\right)$$

This is the generalized Ohm's law. The J x B/(en) term is the Hall term, and the pressure gradient term drives diamagnetic effects.

Ideal MHD Ohm's Law (drop Hall, pressure, resistive terms)

$$\boxed{\mathbf{E} + \mathbf{v}\times\mathbf{B} = 0}$$

Resistive MHD Ohm's Law (keep only resistivity)

$$\boxed{\mathbf{E} + \mathbf{v}\times\mathbf{B} = \eta\mathbf{J}}$$

2.6 The Induction Equation

We derive the evolution equation for B by combining Faraday's law with Ohm's law.

Step 1: Start from Faraday's law

$$\frac{\partial \mathbf{B}}{\partial t} = -\nabla\times\mathbf{E}$$

Step 2: Substitute resistive Ohm's law

From E = -v x B + eta J:

$$\frac{\partial \mathbf{B}}{\partial t} = -\nabla\times(-\mathbf{v}\times\mathbf{B} + \eta\mathbf{J})$$

$$= \nabla\times(\mathbf{v}\times\mathbf{B}) - \nabla\times(\eta\mathbf{J})$$

Step 3: Use Ampere's law to eliminate J

In the MHD limit (no displacement current): J = (nabla x B)/mu₀. Assuming uniform resistivity:

$$\nabla\times(\eta\mathbf{J}) = \frac{\eta}{\mu_0}\nabla\times(\nabla\times\mathbf{B}) = \frac{\eta}{\mu_0}\left[\nabla(\nabla\cdot\mathbf{B}) - \nabla^2\mathbf{B}\right]$$

Since div B = 0 (always), the first term vanishes:

$$= -\frac{\eta}{\mu_0}\nabla^2\mathbf{B}$$

Result: Resistive Induction Equation

$$\boxed{\frac{\partial\mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B}) + \frac{\eta}{\mu_0}\nabla^2\mathbf{B}}$$

The first term is advection (field frozen to fluid) and the second is diffusion (field slips through fluid).

Ideal limit (eta = 0):

$$\frac{\partial\mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B})$$

This implies the frozen-in flux theorem: magnetic field lines are "frozen" into the fluid.

2.7 Magnetic Reynolds Number and Flux Freezing

The ratio of the advection term to the diffusion term defines the magnetic Reynolds number:

Dimensional analysis

$$|\nabla\times(\mathbf{v}\times\mathbf{B})| \sim \frac{vB}{L}, \qquad \left|\frac{\eta}{\mu_0}\nabla^2\mathbf{B}\right| \sim \frac{\eta B}{\mu_0 L^2}$$

$$\boxed{R_m = \frac{\mu_0 v L}{\eta} = \frac{\text{advection rate}}{\text{diffusion rate}}}$$

R_m much greater than 1

Ideal MHD valid. Field frozen to plasma. Typical of astrophysical plasmas (R_m ~ 10¹⁰ in the solar corona).

R_m ~ 1 or less

Resistive effects important. Field diffuses through plasma. Laboratory plasmas, reconnection layers.

Flux Freezing Theorem (Alfven's theorem)

When R_m is much greater than 1, the magnetic flux through any closed loop moving with the fluid is conserved:

$$\Phi = \int_S \mathbf{B}\cdot d\mathbf{A} = \text{const for ideal MHD}$$

Proof sketch: The rate of change of flux through a surface S(t) moving with velocity v is:

$$\frac{d\Phi}{dt} = \int_S \frac{\partial\mathbf{B}}{\partial t}\cdot d\mathbf{A} + \oint_C (\mathbf{v}\times\mathbf{B})\cdot d\boldsymbol{\ell}$$

Using Stokes' theorem on the line integral and substituting the ideal induction equation:

$$= \int_S \left[\frac{\partial\mathbf{B}}{\partial t} - \nabla\times(\mathbf{v}\times\mathbf{B})\right]\cdot d\mathbf{A} = 0$$

Magnetic diffusion time

When v = 0 (static plasma), the induction equation becomes a diffusion equation:

$$\frac{\partial\mathbf{B}}{\partial t} = \frac{\eta}{\mu_0}\nabla^2\mathbf{B} \qquad\Rightarrow\qquad \tau_d = \frac{\mu_0 L^2}{\eta}$$

For a tokamak with L ~ 1 m and eta ~ 10^-7 Omega m: tau_d ~ 10 seconds. For the Sun with L ~ 10⁹ m: tau_d ~ 10¹⁰ years.

2.8 MHD Force Balance: Magnetic Pressure and Tension

The Lorentz force J x B can be decomposed into magnetic pressure and tension:

Step 1: Use Ampere's law

$$\mathbf{J}\times\mathbf{B} = \frac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B}$$

Step 2: Apply the vector identity

$$(\nabla\times\mathbf{B})\times\mathbf{B} = (\mathbf{B}\cdot\nabla)\mathbf{B} - \nabla\left(\frac{B^2}{2}\right)$$

Result

$$\boxed{\mathbf{J}\times\mathbf{B} = -\nabla\left(\frac{B^2}{2\mu_0}\right) + \frac{(\mathbf{B}\cdot\nabla)\mathbf{B}}{\mu_0}}$$

Magnetic Pressure

$$p_B = \frac{B^2}{2\mu_0}$$

Acts isotropically, pushes perpendicular to field lines. Analogous to gas pressure.

Magnetic Tension

$$\frac{(\mathbf{B}\cdot\nabla)\mathbf{B}}{\mu_0} = \frac{B^2}{\mu_0 R_c}\hat{n}$$

Acts along curved field lines toward the center of curvature R_c. Like tension in a rubber band.

2.9 The Complete Ideal MHD System

Collecting all the equations together (with adiabatic closure), the full ideal MHD system is:

$$\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf{v}) = 0 \qquad\text{(mass)}$$

$$\rho\frac{d\mathbf{v}}{dt} = \mathbf{J}\times\mathbf{B} - \nabla p = \frac{1}{\mu_0}(\nabla\times\mathbf{B})\times\mathbf{B} - \nabla p \qquad\text{(momentum)}$$

$$\frac{d}{dt}\left(\frac{p}{\rho^\gamma}\right) = 0 \qquad\text{(energy/adiabatic)}$$

$$\frac{\partial\mathbf{B}}{\partial t} = \nabla\times(\mathbf{v}\times\mathbf{B}) \qquad\text{(induction)}$$

$$\nabla\cdot\mathbf{B} = 0 \qquad\text{(solenoidal constraint)}$$

Counting: 8 equations (1 scalar + 3 vector + 1 scalar + 3 vector) for 8 unknowns (rho, v_x, v_y, v_z, p, B_x, B_y, B_z). The system is closed. The div B = 0 constraint is maintained automatically if satisfied initially.

Interactive Simulation

Flux Freezing: 1D Magnetic Field Advection

Python

script.py147 lines

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

# Demonstrate flux freezing: B is advected with the flow
# Solve dB/dt = d/dx(v*B) + (eta/mu0)*d2B/dx2 in 1D
# Compare ideal (eta=0) vs resistive cases

# Grid parameters
Nx = 400
Lx = 10.0
dx = Lx / Nx
x = np.linspace(0, Lx, Nx)

# Time parameters
dt = 0.002
Nt = 1500
t_total = Nt * dt

# Physical parameters
mu0 = 1.0  # normalized
v0 = 1.0   # advection velocity (constant, uniform flow)

# Initial magnetic field: Gaussian pulse
B0_ideal = np.exp(-((x - 3.0)**2) / 0.5)
B0_resist_low = B0_ideal.copy()
B0_resist_high = B0_ideal.copy()

# Resistivities for three cases
eta_cases = [0.0, 0.005, 0.05]
labels = ['Ideal (eta=0)', 'Low resistivity (eta=0.005)', 'High resistivity (eta=0.05)']
colors = ['cyan', 'orange', 'red']

# Store snapshots
snapshots = {0: [], 1: [], 2: []}
snap_times = [0, Nt//4, Nt//2, 3*Nt//4, Nt-1]

# Time evolution using finite differences (Lax-Wendroff + diffusion)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

for case_idx, eta in enumerate(eta_cases):
    B = B0_ideal.copy()

for n in range(Nt):
        # Advection: Lax-Wendroff scheme
        Bnew = np.zeros_like(B)
        for i in range(1, Nx-1):
            # Advection (upwind for stability)
            advection = -v0 * (B[i] - B[i-1]) / dx
            # Diffusion
            diffusion = (eta / mu0) * (B[i+1] - 2*B[i] + B[i-1]) / dx**2
            Bnew[i] = B[i] + dt * (advection + diffusion)
        # Boundary conditions (periodic-ish)
        Bnew[0] = Bnew[1]
        Bnew[-1] = Bnew[-2]
        B = Bnew

if n in snap_times:
            snapshots[case_idx].append((n*dt, B.copy()))

# Plot 1: All three cases at final time
ax = axes[0, 0]
ax.plot(x, B0_ideal, 'w--', linewidth=1, alpha=0.5, label='Initial B(x)')
for case_idx in range(3):
    t_final, B_final = snapshots[case_idx][-1]
    ax.plot(x, B_final, color=colors[case_idx], linewidth=2, label=labels[case_idx])
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('B(x)', fontsize=12)
ax.set_title('Final State: Ideal vs Resistive', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# Plot 2: Time evolution of ideal case
ax = axes[0, 1]
for t_val, B_snap in snapshots[0]:
    alpha = 0.3 + 0.7 * (t_val / t_total)
    ax.plot(x, B_snap, 'c-', linewidth=1.5, alpha=alpha, label=f't={t_val:.1f}')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('B(x)', fontsize=12)
ax.set_title('Ideal MHD: Frozen-in Flux', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# Plot 3: Time evolution of resistive case
ax = axes[1, 0]
for t_val, B_snap in snapshots[2]:
    alpha = 0.3 + 0.7 * (t_val / t_total)
    ax.plot(x, B_snap, 'r-', linewidth=1.5, alpha=alpha, label=f't={t_val:.1f}')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('B(x)', fontsize=12)
ax.set_title('Resistive MHD: Diffusion + Advection', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# Plot 4: Peak B vs time (conservation check)
ax = axes[1, 1]
for case_idx, eta in enumerate(eta_cases):
    B = B0_ideal.copy()
    peak_B = []
    flux = []
    times = []
    for n in range(Nt):
        peak_B.append(np.max(B))
        flux.append(np.trapezoid(B, x))
        times.append(n * dt)
        Bnew = np.zeros_like(B)
        for i in range(1, Nx-1):
            advection = -v0 * (B[i] - B[i-1]) / dx
            diffusion = (eta / mu0) * (B[i+1] - 2*B[i] + B[i-1]) / dx**2
            Bnew[i] = B[i] + dt * (advection + diffusion)
        Bnew[0] = Bnew[1]
        Bnew[-1] = Bnew[-2]
        B = Bnew
    ax.plot(times, np.array(peak_B)/peak_B[0], color=colors[case_idx], linewidth=2, label=labels[case_idx])

ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Peak B / B_max(0)', fontsize=12)
ax.set_title('Peak Field Decay', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_ylim(0, 1.1)

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 Rm values
print("Magnetic Reynolds Numbers:")
print(f"  Ideal:  Rm = infinity (perfect flux freezing)")
L = 0.5  # characteristic length of the pulse
for eta in [0.005, 0.05]:
    Rm = mu0 * v0 * L / eta
    tau_d = mu0 * L**2 / eta
    print(f"  eta={eta}: Rm = mu0*v*L/eta = {Rm:.1f}, tau_diff = {tau_d:.1f}")
print(f"\nAt Rm >> 1: field advected with flow (frozen-in)")
print(f"At Rm ~ 1: significant diffusion broadens and weakens the pulse")
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

-- MHD is derived from two-fluid equations by summing over species and exploiting m_e much less than m_i
-- The MHD momentum equation contains only J x B and pressure gradient (no explicit E field)
-- Ohm's law comes from the electron momentum equation with negligible electron inertia
-- The induction equation governs B evolution: advection (frozen flux) vs diffusion
-- Magnetic Reynolds number R_m = mu₀vL/eta determines the regime
-- J x B decomposes into magnetic pressure B²/(2mu₀) and magnetic tension (B dot nabla)B/mu₀

← Moment Equations Next: MHD Equilibrium →