Magnetic confinement, MHD equilibrium and tokamaks

A plasma at n = 10²⁰ m⁻³ and T = 10 keV has pressure p = 2nT = 2 × 10²⁰ × 10 × 1.602 × 10⁻¹⁶ J = 3.2 × 10⁵ Pa — three atmospheres at a hundred million kelvin, and no wall can touch it. Ideal MHD holds it with the field: at rest the momentum equation is J × B = ∇p, so B⋅∇p = 0 and J⋅∇p = 0 — pressure is constant on magnetic surfaces and the current flows within them. The field's own pressure is B²/2μ₀ = 11.2 MPa at 5.3 T, so the plasma uses only β = 2μ₀p/B² ≈ 2.9 % of it. A straight cylinder loses particles out of its ends, so bend it into a torus — but a purely toroidal field is not an equilibrium: its 1/R gradient and curvature drive ions and electrons vertically in opposite directions, the resulting E × B drift pushes the whole column outward, and nothing balances the hoop force. The cure is a poloidal component that winds each field line helically, so a particle's upward drift on one pass is cancelled on the next. A tokamak makes that poloidal field with the plasma's own toroidal current, which is why Iₚ is the number everything else scales with.

Write the steady-state ledger per unit volume for a 50:50 D–T plasma with n_D = n_T = n/2 and Tₑ = Tᵢ = T. Alpha heating is (n²/4)⟨σv⟩E_α with E_α = 3.5 MeV; transport removes the stored energy 3nT once per confinement time, 3nT/τ_E; and the electrons radiate bremsstrahlung C_B n²√T with C_B = 5.35 × 10⁻³⁷ W m³ per √keV. Ignition — alpha heating alone sustaining the plasma — needs (n²/4)⟨σv⟩E_α ≥ 3nT/τ_E + C_B n²√T. At T = 15 keV the Bosch–Hale value is ⟨σv⟩ = 2.74 × 10⁻²² m³/s, so the heating coefficient is ¼ × 2.74 × 10⁻²² × 5.61 × 10⁻¹³ J = 3.84 × 10⁻³⁵ W m³ against bremsstrahlung 2.07 × 10⁻³⁶ W m³, a 5.4 % tax. Then nτ_E ≥ 3 × 15 keV/[(3.84 − 0.21) × 10⁻³⁵ W m³] = 7.21 × 10⁻¹⁵ J/(3.63 × 10⁻³⁵ W m³) = 1.99 × 10²⁰ s m⁻³, and nTτ_E ≥ 3.0 × 10²¹ keV s m⁻³. Because ⟨σv⟩ ∝ T² between 10 and 20 keV, the T² in the numerator cancels and the requirement is nearly flat there; below about 4.3 keV bremsstrahlung outruns alpha heating at any density — the ideal ignition temperature. Every term here is bookkeeping of ⟨σv⟩ from the last two lessons; nothing in it is confinement physics yet.

Axisymmetry does the work. Write B = ∇ψ × ∇φ + F∇φ, where ψ(R, Z) is the poloidal flux per radian and F = RB_φ; ∇⋅B = 0 is then automatic. Because B⋅∇p = 0, p is a function of ψ alone, and J⋅∇p = 0 forces F = F(ψ) too. Substituting into J × B = ∇p and keeping the component along ∇ψ leaves one scalar equation, Δ*ψ = −μ₀R² p′(ψ) − FF′(ψ), where Δ* = R∂_R(R⁻¹∂_R) + ∂²_Z is elliptic and carries the 1/R that distinguishes a torus from a cylinder. The unknown is ψ on the (R, Z) plane; the inputs are two profiles, p(ψ) and F(ψ), supplied by measurement, not by MHD; the boundary condition is ψ = ψ_b on the last closed flux surface — a limiter surface, or a separatrix whose X-point (∇ψ = 0) feeds the divertor. The right-hand side depends on the solution, so solve numerically: discretise Δ* on a grid, evaluate the source with the current ψ, invert the linear operator, repeat — the Picard loop equilibrium codes run between shots. One output is the Shafranov shift: the magnetic axis sits outboard of the geometric centre by about (a²/2R)(βₚ + lᵢ/2), roughly 0.35 m for ITER.

The equilibrium fixes how tightly each field line winds. On a surface of minor radius r a line makes q = rB_φ/(RB_θ) toroidal turns per poloidal turn in the large-aspect-ratio circular limit — q is the pitch, and it is the coordinate stability speaks in. At the edge, Ampère's law gives B_θ(a) = μ₀Iₚ/2πa, so qₐ = 2πa²B_φ/(μ₀RIₚ): more current unwinds the pitch, more toroidal field tightens it. For ITER, 15 MA in a = 2.0 m gives B_θ(a) = 1.5 T and q_cyl = 2.0 × 5.3/(6.2 × 1.5) = 1.14; the D-shaped section with κ = 1.7 multiplies this by about (1 + κ²)/2 = 1.95, and toroidal corrections take q₉₅ to 3. Inside, q follows the current profile: for j ∝ (1 − r²/a²)^ν the enclosed current gives q(r) = qₐ(r²/a²)/[1 − (1 − r²/a²)^(ν+1)], so q₀ = qₐ/(ν + 1) — a parabolic profile puts the axis at half the edge value. Surfaces where q = m/n are rational: a helical perturbation with m poloidal and n toroidal periods lines up with the field there and bends no field lines, which is why those surfaces are where instabilities live. The figure below draws this profile.

Displace the equilibrium by ξ(r, t) and linearise. With p₁ = −ξ⋅∇p − γp∇⋅ξ and B₁ = ∇×(ξ×B) ≡ Q, the dynamics collapse to ρ∂²ξ/∂t² = F(ξ) with F(ξ) = J×Q + μ₀⁻¹(∇×Q)×B + ∇(ξ⋅∇p + γp∇⋅ξ). With n⋅ξ = 0 on a conducting wall, F is self-adjoint: ∫η⋅F(ξ) dV = ∫ξ⋅F(η) dV. For a normal mode ξe^(−iωt), −ρω²ξ = F(ξ), so ω² = −∫ξ*⋅F(ξ) dV/∫ρ|ξ|² dV is real, exactly as ⟨ψ|H|ψ⟩/⟨ψ|ψ⟩ is real for Hermitian H. Two cases remain: ω² > 0 oscillates, ω² < 0 grows as e^(γt) with γ = √(−ω²). Define δW = −½∫ξ*⋅F(ξ) dV; the plasma is stable if and only if δW ≥ 0 for every admissible ξ, and, by Rayleigh–Ritz, one trial ξ with δW < 0 proves instability outright. Written out, δW has three positive terms — field-line bending |Q_⊥|²/μ₀, field compression, plasma compression γp|∇⋅ξ|² — against two that can be negative, the pressure gradient in bad curvature and the parallel current: interchange and ballooning modes tap the first, kinks the second. Numerically, discretising −F gives a symmetric matrix A and ρ a positive-definite M; Aξ = ω²Mξ has real eigenvalues, and the sign of the lowest is the test.

The theorem's hypothesis is the ideal Ohm's law, which freezes field lines into the fluid and forbids any change of magnetic topology. Add resistivity η and the constraint fails in a thin layer at each rational surface, where field lines reconnect into islands. That is the tearing mode: it grows on the hybrid timescale τ_R³⁄⁵τ_A²⁄⁵, between the Alfvén time τ_A = a/v_A and the resistive time τ_R = μ₀a²/η, and it is judged by Δ′, a matching condition across the layer, not by δW, because a topology-changing ξ was never in ideal MHD's trial space. Neoclassical tearing modes, driven by the bootstrap current lost inside an island, set the β actually reachable below the Troyon value. Two more limits are empirical: the Greenwald density limit n_G = Iₚ/πa² (10²⁰ m⁻³ for Iₚ in MA, a in m), and τ_E itself, which comes from the IPB98(y,2) scaling of many machines, not from any equation here. The triple product is thus a target the equilibrium can hold and the ideal test can protect; whether τ_E reaches the 2 s the balance asks for at 10²⁰ m⁻³ is answered by turbulence and by experiment.