Part 3, Chapter 5

MHD Instabilities

Stability analysis, energy principles, and common instability modes

5.1 Linearized MHD and the Force Operator

To analyze stability, we linearize the ideal MHD equations about an equilibrium. Let every quantity Q = Q₀ + Q₁ where Q₀ is the equilibrium and Q₁ is a small perturbation.

Step 1: Introduce the displacement vector

The fluid displacement xi relates to the perturbed velocity by:

$$\mathbf{v}_1 = \frac{\partial\boldsymbol{\xi}}{\partial t}$$

All perturbed quantities can be expressed in terms of xi.

Step 2: Express perturbed quantities

$$\rho_1 = -\nabla\cdot(\rho_0\boldsymbol{\xi})$$

$$p_1 = -\boldsymbol{\xi}\cdot\nabla p_0 - \gamma p_0\nabla\cdot\boldsymbol{\xi}$$

$$\mathbf{B}_1 = \mathbf{Q} \equiv \nabla\times(\boldsymbol{\xi}\times\mathbf{B}_0)$$

Step 3: The linearized momentum equation

$$\rho_0\frac{\partial^2\boldsymbol{\xi}}{\partial t^2} = \mathbf{F}(\boldsymbol{\xi})$$

where F(xi) is the MHD force operator:

$$\mathbf{F}(\boldsymbol{\xi}) = \nabla(\boldsymbol{\xi}\cdot\nabla p_0 + \gamma p_0\nabla\cdot\boldsymbol{\xi}) + \frac{1}{\mu_0}(\nabla\times\mathbf{Q})\times\mathbf{B}_0 + \frac{1}{\mu_0}(\nabla\times\mathbf{B}_0)\times\mathbf{Q}$$

Step 4: Self-adjointness

The force operator F is self-adjoint (Hermitian), meaning:

$$\int \boldsymbol{\eta}^*\cdot\mathbf{F}(\boldsymbol{\xi})\,d^3r = \int \boldsymbol{\xi}\cdot\mathbf{F}(\boldsymbol{\eta})^*\,d^3r$$

This guarantees that eigenvalues omega² are real: modes are either purely oscillatory (omega² greater than 0, stable) or purely growing/decaying (omega² less than 0, unstable). There are no overstable modes in ideal MHD.

5.2 The Energy Principle

Since F is self-adjoint, stability can be determined without solving the full eigenvalue problem. The equilibrium is stable if and only if the potential energy change is positive for all trial displacements.

Step 1: Define the energy functional

$$\delta W = -\frac{1}{2}\int \boldsymbol{\xi}^*\cdot\mathbf{F}(\boldsymbol{\xi})\,d^3r$$

The negative sign ensures that delta W greater than 0 corresponds to stability (restoring force opposes displacement).

Step 2: Expand into physical terms

After integration by parts and using the equilibrium relations, the energy functional becomes:

$$\delta W = \frac{1}{2}\int \left[\frac{|\mathbf{Q}|^2}{\mu_0} + \gamma p_0|\nabla\cdot\boldsymbol{\xi}|^2 + (\boldsymbol{\xi}\cdot\nabla p_0)(\nabla\cdot\boldsymbol{\xi}^*) + \mathbf{J}_0\cdot(\boldsymbol{\xi}^*\times\mathbf{Q})\right]d^3r$$

Step 3: Identify each term

Stabilizing Terms (always positive)

|Q|²/mu₀: Energy to bend field lines. Proportional to (k dot B)² -- vanishes when perturbation is perpendicular to B.
gamma p|nabla dot xi|²: Energy to compress the plasma. Always positive.

Potentially Destabilizing Terms

(xi dot nabla p)(nabla dot xi*): Pressure-driven instabilities. Unfavorable when displacement is along the pressure gradient with compression.
J₀ dot (xi* x Q): Current-driven instabilities (kink modes).

Stability Criterion

$$\boxed{\delta W[\boldsymbol{\xi}] > 0 \quad\text{for all admissible } \boldsymbol{\xi} \quad\Leftrightarrow\quad \text{Stable}}$$

If there exists any trial displacement that makes delta W negative, the equilibrium is unstable.

5.3 Rayleigh-Taylor Instability in MHD

A heavy fluid supported against gravity by a light fluid (or magnetic field) is unstable. In MHD, magnetic tension provides a stabilizing effect.

Step 1: Setup -- density interface with gravity

Consider a plasma with density rho(z) in a gravitational field g = -g z-hat, supported by a horizontal magnetic field B = B₀ x-hat. A perturbation proportional to exp(ik_xx + ik_yy + gamma t) is applied.

Step 2: Linearize the momentum equation

The linearized vertical momentum equation for incompressible perturbations gives:

$$\rho_0\gamma^2\xi_z = -g\frac{d\rho_0}{dz}\xi_z + \frac{(\mathbf{k}\cdot\mathbf{B}_0)^2}{\mu_0}\xi_z$$

Step 3: For a sharp interface (Atwood number)

At a sharp density jump from rho₁ (heavy, on top) to rho₂ (light, below), with Atwood number A = (rho₁ - rho₂)/(rho₁ + rho₂):

Result: MHD Rayleigh-Taylor growth rate

$$\boxed{\gamma^2 = gk\mathcal{A} - \frac{(k_x B_0)^2}{\mu_0(\rho_1 + \rho_2)/2} = gk\mathcal{A} - \frac{(\mathbf{k}\cdot\mathbf{B})^2}{\mu_0\bar{\rho}}}$$

The first term drives instability (gravity + density gradient). The second term stabilizes through magnetic tension. Modes with k perpendicular to B (so k dot B = 0) are most unstable -- identical to the hydrodynamic case. Modes with k parallel to B can be fully stabilized.

Critical wavenumber for stabilization

Setting gamma² = 0 for k along B:

$$k_{\text{crit}} = \frac{g\mathcal{A}\mu_0\bar{\rho}}{B_0^2}$$

Modes with k greater than k_crit (short wavelengths) are stabilized by magnetic tension. Only long wavelengths remain unstable if k is along B.

5.4 Kelvin-Helmholtz Instability

Velocity shear across a magnetic surface can drive the Kelvin-Helmholtz instability, even in the absence of gravity.

Setup

Two layers with flow velocities V₁ and V₂, densities rho₁ and rho₂, and magnetic fields B₁ and B₂ along the interface. The dispersion relation for the sharp interface is:

Dispersion relation

$$\gamma^2 = \frac{\rho_1\rho_2(\mathbf{k}\cdot\Delta\mathbf{V})^2}{(\rho_1+\rho_2)^2} - \frac{(\mathbf{k}\cdot\mathbf{B}_1)^2 + (\mathbf{k}\cdot\mathbf{B}_2)^2}{\mu_0(\rho_1+\rho_2)}$$

where Delta V = V₁ - V₂ is the velocity shear. The first term drives instability; the second stabilizes through magnetic tension on both sides of the interface.

Stability condition

$$\boxed{(\mathbf{k}\cdot\Delta\mathbf{V})^2 < \frac{(\rho_1+\rho_2)}{\mu_0\rho_1\rho_2}\left[(\mathbf{k}\cdot\mathbf{B}_1)^2 + (\mathbf{k}\cdot\mathbf{B}_2)^2\right]}$$

A strong magnetic field component along the flow stabilizes against K-H. This is why the magnetopause is often stable despite large velocity shears.

5.5 Sausage and Kink Instabilities in Z-Pinch

Consider a cylindrical plasma column (Z-pinch) with axial current I generating azimuthal field B_theta, and an applied axial field B_z. Perturbations proportional to exp(imtheta + ikz) are characterized by azimuthal mode number m.

Sausage Instability (m = 0)

Axisymmetric pinching: the column alternately thins and bulges. The growth rate for a Z-pinch of radius a is:

$$\gamma^2 = \frac{k^2 v_{A\theta}^2}{1 + k^2 a^2}\left(\frac{B_z^2}{B_\theta^2(a)} - 1 - k^2 a^2\right)$$

where v_{A theta} = B_theta(a)/sqrt(mu₀ rho).

Stability condition: B_z² greater than B_theta²(a)(1 + k²a²). For long wavelengths (ka much less than 1), this requires B_z greater than B_theta.

Kink Instability (m = 1)

Helical displacement of the entire column. The key stability criterion is the Kruskal-Shafranov condition:

$$q(a) = \frac{2\pi a B_z}{\mu_0 I} = \frac{a B_z}{R B_\theta(a)} > 1$$

where q is the safety factor and R is the major radius (for a torus). Derivation: the kink is unstable when the field line pitch allows the perturbation to be resonant:

$$\mathbf{k}\cdot\mathbf{B} = \frac{m B_\theta}{r} + k B_z = 0 \quad\Rightarrow\quad k = -\frac{m B_\theta}{r B_z}$$

For m=1 and the perturbation wavelength fitting in the torus (k = 1/R), the resonance condition gives q = 1 at the edge.

Higher m modes

For m greater than or equal to 2 (fluting modes), the stability condition involves the Suydam criterion for cylindrical plasmas, requiring a minimum shear in the safety factor profile q(r) to stabilize against interchange-like modes.

5.6 Ballooning Modes in Tokamaks

Pressure-driven instabilities localized to the bad curvature (low-field) side of a torus.

Physical picture

In a tokamak, the outer midplane has unfavorable (bad) curvature: kappa points outward while nabla p points inward. This drives interchange-like instabilities. However, the perturbation must be consistent along the field line, which passes through both good and bad curvature regions. Ballooning modes maximize amplitude in the bad curvature region while maintaining periodicity.

The ballooning parameter

$$\alpha = -\frac{2\mu_0 R q^2}{B^2}\frac{dp}{dr} = -R q^2 \frac{d\beta}{dr}$$

This is the normalized pressure gradient. Instability occurs when alpha exceeds a critical value alpha_crit that depends on the magnetic shear s = (r/q)(dq/dr).

The s-alpha stability diagram

In the (s, alpha) plane, the first stability region is bounded approximately by:

$$\alpha_{\text{crit}} \approx 0.6\,s \quad\text{(for large aspect ratio)}$$

At low shear, a second stability region can open up at high alpha, relevant for advanced tokamak scenarios.

5.7 Tearing Mode and Magnetic Reconnection

The tearing mode is a resistive instability that breaks and reconnects magnetic field lines at resonant surfaces where k dot B = 0.

Step 1: Setup -- current sheet

Consider a sheared magnetic field B = B₀ tanh(x/a) y-hat (Harris sheet). At x = 0, the field reverses direction. The equilibrium current is:

$$J_z = \frac{B_0}{\mu_0 a}\text{sech}^2(x/a)$$

Step 2: Outer region (ideal MHD)

Away from x = 0, ideal MHD applies. The perturbed flux function psi₁ (where B₁ = nabla psi₁ x z-hat) satisfies:

$$\frac{d^2\psi_1}{dx^2} - k^2\psi_1 - \frac{B_0''}{B_0}\psi_1 = 0$$

This has a logarithmic singularity at the resonant surface. The jump in d(psi₁)/dx across x = 0 defines the tearing mode stability parameter:

$$\Delta' = \left[\frac{1}{\psi_1}\frac{d\psi_1}{dx}\right]_{0^-}^{0^+}$$

Step 3: Inner region (resistive layer)

Near x = 0, resistivity is essential. The resistive layer has width delta:

$$\delta \sim \left(\frac{\eta\gamma}{\mu_0 (B_0'/\sqrt{\mu_0\rho})^2}\right)^{1/4} a$$

Matching the inner and outer solutions determines the growth rate.

Result: Tearing mode growth rate

$$\boxed{\gamma \sim \left(\frac{\eta}{\mu_0}\right)^{3/5}\left(\frac{k B_0'}{\sqrt{\mu_0\rho}}\right)^{2/5} \sim \tau_A^{-2/5}\tau_R^{-3/5}}$$

where tau_A = a/v_A is the Alfven time and tau_R = mu₀a²/eta is the resistive time. The growth rate is a geometric mean -- much faster than pure diffusion but slower than the Alfven time.

Magnetic island width

The nonlinear saturation creates magnetic islands with width:

$$w = 4\sqrt{\frac{\psi_1}{B_0'}} \propto \psi_1^{1/2}$$

Stability criterion: the tearing mode is unstable when Delta prime greater than 0 (the outer solution has a positive jump in the derivative of psi across the resonant surface).

Interactive Simulation

Z-Pinch Instability Growth Rates: Sausage (m=0) and Kink (m=1)

Python

script.py162 lines

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

# Z-pinch stability analysis
# Plasma column of radius a, axial field Bz, azimuthal field Btheta(a)
# Current I -> Btheta(a) = mu0*I/(2*pi*a)

# Normalized units: a = 1, mu0 = 1, rho = 1
a = 1.0
mu0 = 1.0
rho = 1.0

# Vary the ratio Bz/Btheta to explore stability
Btheta = 1.0  # B_theta at r = a

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

# ===== Plot 1: Sausage (m=0) growth rate vs ka for different Bz/Btheta =====
ax = axes[0, 0]
ka = np.linspace(0.01, 5.0, 300)

Bz_ratios = [0.5, 0.8, 1.0, 1.2, 1.5, 2.0]
colors_bz = plt.cm.plasma(np.linspace(0.1, 0.9, len(Bz_ratios)))

for i, ratio in enumerate(Bz_ratios):
    Bz = ratio * Btheta
    vA_theta = Btheta / np.sqrt(mu0 * rho)

# gamma^2 = k^2 * vA_theta^2 / (1 + k^2*a^2) * (Bz^2/Btheta^2 - 1 - k^2*a^2)
    gamma_sq = (ka**2 * vA_theta**2) / (1 + ka**2) * (Bz**2/Btheta**2 - 1 - ka**2)

# Only plot where gamma^2 > 0 (unstable)
    gamma = np.where(gamma_sq > 0, np.sqrt(gamma_sq), 0)
    ax.plot(ka, gamma, color=colors_bz[i], linewidth=2, label=f'Bz/Bth={ratio:.1f}')

ax.set_xlabel('k*a (normalized wavenumber)', fontsize=12)
ax.set_ylabel('Growth rate gamma / (vA_theta/a)', fontsize=12)
ax.set_title('Sausage Instability (m=0)', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_ylim(0, 1.5)

# ===== Plot 2: Stability boundary in (ka, Bz/Btheta) space =====
ax = axes[0, 1]
ka_grid = np.linspace(0.01, 5.0, 300)
Bz_Bt_grid = np.linspace(0.01, 3.0, 300)
KA, BZ_BT = np.meshgrid(ka_grid, Bz_Bt_grid)

# Sausage: unstable when Bz^2/Bt^2 < 1 + (ka)^2
sausage_unstable = BZ_BT**2 < 1 + KA**2

# Kink (m=1): unstable when q < 1, i.e., Bz/Bt < 1 (simplified)
kink_unstable = BZ_BT < 1.0

ax.contourf(KA, BZ_BT, sausage_unstable.astype(float), levels=[0.5, 1.5],
            colors=['red'], alpha=0.3)
ax.contourf(KA, BZ_BT, kink_unstable.astype(float), levels=[0.5, 1.5],
            colors=['blue'], alpha=0.3)

# Stability boundary
ka_bound = np.linspace(0.01, 5.0, 300)
Bz_Bt_crit = np.sqrt(1 + ka_bound**2)
ax.plot(ka_bound, Bz_Bt_crit, 'r-', linewidth=2, label='Sausage boundary')
ax.axhline(1.0, color='blue', linewidth=2, linestyle='--', label='Kink boundary (q=1)')
ax.set_xlabel('k*a', fontsize=12)
ax.set_ylabel('Bz / B_theta', fontsize=12)
ax.set_title('Stability Diagram', fontsize=14, color='white')
ax.legend(fontsize=10)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax.set_ylim(0, 3)
ax.text(2.5, 0.5, 'UNSTABLE\n(both)', color='white', fontsize=11, ha='center')
ax.text(2.5, 2.5, 'STABLE', color='lime', fontsize=11, ha='center')

# ===== Plot 3: Rayleigh-Taylor growth rate with B stabilization =====
ax = axes[1, 0]
k = np.linspace(0.01, 10.0, 300)
g = 1.0  # gravity
A_atwood = 0.5  # Atwood number

B0_values = [0.0, 0.3, 0.5, 0.8, 1.0]
colors_b = plt.cm.cool(np.linspace(0.1, 0.9, len(B0_values)))

for i, B0 in enumerate(B0_values):
    # gamma^2 = g*k*A - (k*B0)^2 / (mu0 * rho_avg)
    rho_avg = 1.0
    gamma_sq = g * k * A_atwood - (k * B0)**2 / (mu0 * rho_avg)
    gamma = np.where(gamma_sq > 0, np.sqrt(gamma_sq), 0)
    ax.plot(k, gamma, color=colors_b[i], linewidth=2,
            label=f'B={B0:.1f} (k_crit={g*A_atwood*mu0*rho_avg/B0**2:.1f})' if B0 > 0 else f'B=0 (hydro)')

ax.set_xlabel('Wavenumber k', fontsize=12)
ax.set_ylabel('Growth rate gamma', fontsize=12)
ax.set_title('Rayleigh-Taylor with Magnetic Field', fontsize=14, color='white')
ax.legend(fontsize=9)
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')

# ===== Plot 4: Tearing mode scaling =====
ax = axes[1, 1]
S = np.logspace(2, 10, 200)  # Lundquist number S = tau_R / tau_A

# Growth rate scaling: gamma * tau_A ~ S^{-3/5}
gamma_tearing = S**(-3/5)
# Resistive layer width: delta/a ~ S^{-2/5}
delta_over_a = S**(-2/5)
# Island growth time: tau_island/tau_A ~ S^{3/5}
tau_island = S**(3/5)

ax2 = ax.twinx()
ax.loglog(S, gamma_tearing, 'cyan', linewidth=2, label='gamma * tau_A ~ S^{-3/5}')
ax.loglog(S, delta_over_a, 'orange', linewidth=2, label='delta/a ~ S^{-2/5}')

# Mark typical values
typical_S = {'Lab plasma': 1e4, 'Tokamak': 1e7, 'Solar corona': 1e10}
for name, S_val in typical_S.items():
    ax.axvline(S_val, color='gray', linestyle=':', alpha=0.5)
    ax.text(S_val, 0.5, name, color='gray', fontsize=8, rotation=90, va='bottom')

ax.set_xlabel('Lundquist number S = tau_R / tau_A', fontsize=12)
ax.set_ylabel('Normalized growth rate / layer width', fontsize=12)
ax.set_title('Tearing Mode Scaling', fontsize=14, color='white')
ax.legend(fontsize=10, loc='upper right')
ax.set_facecolor('#0a0a0a')
ax.tick_params(colors='white')
ax2.tick_params(colors='white')
ax2.set_ylabel('')

for ax_item in axes.flat:
    ax_item.spines['bottom'].set_color('white')
    ax_item.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("=== Z-Pinch Stability Summary ===")
print(f"Sausage (m=0): stable when Bz/Btheta > sqrt(1 + k^2*a^2)")
print(f"  -> For long wavelengths (ka << 1): Bz > Btheta")
print(f"  -> For short wavelengths: need stronger Bz")
print(f"Kink (m=1): stable when q(a) = a*Bz/(R*Btheta) > 1")
print(f"  (Kruskal-Shafranov criterion)")

print(f"\n=== Rayleigh-Taylor with B ===")
print(f"gamma^2 = g*k*A - (k.B)^2/(mu0*rho)")
print(f"Modes with k perp to B: no stabilization")
print(f"Modes with k || B: stabilized for k > k_crit")

print(f"\n=== Tearing Mode ===")
print(f"Growth rate: gamma ~ tau_A^{{-2/5}} * tau_R^{{-3/5}}")
print(f"Layer width: delta ~ a * S^{{-2/5}}")
print(f"Island width: w ~ 4*sqrt(psi/B')")
for name, S_val in typical_S.items():
    gamma_val = S_val**(-3/5)
    print(f"  {name} (S={S_val:.0e}): gamma*tau_A = {gamma_val:.2e}")

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 MHD force operator F(xi) is self-adjoint, so eigenvalues omega² are real (no overstability in ideal MHD)
-- Energy principle: delta W greater than 0 for all displacements implies stability, without solving the full eigenvalue problem
-- Rayleigh-Taylor: magnetic tension stabilizes short-wavelength modes with k along B, but modes with k perpendicular to B remain unstable
-- Kelvin-Helmholtz: magnetic field along the flow provides stabilization against velocity shear
-- Sausage (m=0): requires B_z greater than B_theta for long wavelengths; kink (m=1): requires q greater than 1
-- Ballooning modes set the beta-limit in tokamaks via the s-alpha stability boundary
-- Tearing mode: resistive reconnection with gamma scaling as S^-3/5; creates magnetic islands with width proportional to psi^1/2

← MHD Waves Next: Two-Fluid Theory →