Plasma Physics Primer for Tokamak Control¶
A ground-up introduction to fusion plasma physics, written for engineers and
scientists who will use scpn-control but have no prior plasma background.
Every equation referenced here has a direct counterpart in the codebase.
1. What Is Plasma?¶
The Fourth State of Matter¶
Heat a solid and it melts. Heat the liquid and it boils. Heat the gas further -- past roughly 10,000 K -- and electrons strip from atoms. The result is plasma: a gas of free ions and electrons that responds collectively to electromagnetic fields.
Plasma is not merely "hot gas." Two properties distinguish it:
- Quasineutrality. On length scales larger than the Debye length, the electron and ion charge densities nearly cancel:
where \(Z_i\) is the ion charge state. Microscopically there are charge fluctuations, but macroscopically the plasma is electrically neutral.
- Collective behaviour. Each particle interacts simultaneously with many distant particles through their shared electromagnetic fields, not just with its nearest neighbour via binary collisions.
Debye Shielding¶
Drop a test charge into a plasma. Free electrons rush toward a positive charge (or flee a negative one), forming a cloud that shields the test charge's field. The characteristic shielding distance is the Debye length:
For a typical fusion plasma (\(T_e \sim 10\) keV, \(n_e \sim 10^{20}\) m\(^{-3}\)), \(\lambda_D \approx 70\;\mu\text{m}\) -- far smaller than the plasma radius (metres). This separation of scales is what makes quasineutrality valid.
A plasma must satisfy \(\lambda_D \ll L\) (system size) and \(N_D = n_e \frac{4}{3}\pi\lambda_D^3 \gg 1\) (many particles in a Debye sphere). Fusion plasmas have \(N_D \sim 10^9\).
Plasma Frequency¶
If electrons are displaced from their equilibrium, the restoring electric field produces oscillations at the plasma frequency:
with \(n_e\) in m\(^{-3}\). At \(n_e = 10^{20}\) m\(^{-3}\), \(\omega_{pe}/(2\pi) \approx 90\) GHz. Electromagnetic waves below this frequency cannot propagate into the plasma -- they reflect. This matters for microwave heating (Section 7).
Why Plasma Is Hard to Contain¶
Charged particles in a magnetic field gyrate around field lines with Larmor radius:
For a 10 keV deuteron in a 5 T field, \(\rho_L \approx 4\) mm. The particle is tied to its field line -- but free to stream along it at the thermal speed (\(\sim 10^6\) m/s for ions). A straight magnetic bottle leaks out both ends. A toroidal (doughnut-shaped) field solves end losses but introduces new drift-driven losses (Section 2).
2. Magnetic Confinement and the Tokamak¶
Toroidal Geometry¶
A tokamak confines plasma inside a torus. Two lengths define its shape:
- Major radius \(R_0\): distance from the torus axis of symmetry to the centre of the plasma cross-section.
- Minor radius \(a\): radius of the plasma cross-section.
The inverse aspect ratio \(\varepsilon = a/R_0\) characterises how tightly wound the torus is. Typical values: \(\varepsilon \approx 0.3\).
┌─────────────────────────┐
│ TF coil │
│ ┌─────────────────┐ │
│ │ plasma │ │
│ │ (hot core) │ │
│ │ R0 ←──→ a │ │
│ └─────────────────┘ │
└─────────────────────────┘
↑ toroidal direction (the long way around)
→ poloidal direction (the short way around)
Toroidal and Poloidal Fields¶
A purely toroidal field (produced by external TF coils) cannot confine plasma alone. The curvature and gradient of \(B_\phi\) cause charge-dependent vertical drifts: ions drift up, electrons drift down (or vice versa). The resulting vertical electric field drives an outward \(\mathbf{E} \times \mathbf{B}\) drift that expels the plasma.
The fix: add a poloidal magnetic field \(B_\theta\), so field lines wrap helically around the torus. Particles following these helical lines spend equal time above and below the midplane, averaging out the vertical drift.
In a tokamak, \(B_\theta\) is generated by a toroidal plasma current \(I_p\) (typically several megaamps), driven inductively by a central solenoid acting as the primary of a transformer.
Safety Factor¶
The safety factor \(q\) measures how many toroidal transits a field line makes per poloidal transit:
where \(\rho = r/a\) is the normalised minor radius. At the edge, \(q_{95}\) (evaluated at the surface enclosing 95% of the poloidal flux) is the standard stability metric. Typical operating range: \(q_{95} \in [3, 5]\).
If \(q\) drops below specific rational values \(m/n\) (integers), resonant perturbations can grow into islands and destroy confinement (Section 6).
Magnetic Surfaces and Flux Coordinates¶
In axisymmetric equilibrium, the poloidal magnetic flux function \(\psi(R, Z)\) is constant on nested toroidal surfaces. Particles and heat flow primarily along these surfaces, with slow cross-field transport setting the confinement time.
The coordinate \(\rho = \sqrt{(\psi - \psi_0)/(\psi_a - \psi_0)}\) maps linearly from the magnetic axis (\(\rho = 0\)) to the last closed flux surface (\(\rho = 1\)). Profiles of temperature \(T(\rho)\), density \(n(\rho)\), and current \(j(\rho)\) are functions of \(\rho\).
Reference Machines¶
| Parameter | ITER | SPARC | DIII-D |
|---|---|---|---|
| \(R_0\) (m) | 6.2 | 1.85 | 1.67 |
| \(a\) (m) | 2.0 | 0.57 | 0.67 |
| \(I_p\) (MA) | 15.0 | 8.7 | 1.5 |
| \(B_T\) (T) | 5.3 | 12.2 | 2.1 |
| \(\kappa\) | 1.7 | 1.97 | 1.8 |
| Status | Construction | Under construction | Operating |
These machines are preconfigured in scpn_control.core.tokamak_config.
3. Plasma Equilibrium: The Grad-Shafranov Equation¶
Force Balance¶
A magnetically confined plasma in steady state satisfies:
The magnetic force balances the pressure gradient. This is a necessary condition for confinement: without it, the plasma expands freely.
The Equation¶
In axisymmetric toroidal geometry with cylindrical coordinates \((R, \phi, Z)\), force balance reduces to a single scalar PDE for the poloidal flux \(\psi(R, Z)\):
where:
- \(\Delta^* \psi = R \frac{\partial}{\partial R}\!\left(\frac{1}{R} \frac{\partial \psi}{\partial R}\right) + \frac{\partial^2 \psi}{\partial Z^2}\) is the toroidal elliptic operator (Grad-Shafranov operator)
- \(p(\psi)\) is the pressure profile (a free function)
- \(F(\psi) = R B_\phi\) is the poloidal current function (a free function)
- Primes denote \(d/d\psi\)
The two free functions \(p'(\psi)\) and \(FF'(\psi)\) encode the profiles and must be supplied (from transport modelling or experiment).
What \(\psi\) Represents¶
\(\psi(R, Z)\) is the poloidal magnetic flux per radian. Contours of constant \(\psi\) are magnetic surfaces. The magnetic axis (where \(|\nabla\psi| = 0\) and \(\psi\) is extremal) sits at the plasma centre. The separatrix (in diverted plasmas) or limiter contact point defines the plasma boundary.
Boundary Conditions¶
Two configurations exist:
- Limiter (fixed boundary): The plasma touches a material wall. \(\psi\) is prescribed on a known boundary contour.
- Diverted (free boundary): External poloidal field coils create an X-point where \(\nabla\psi = 0\) but \(\psi \neq \psi_\text{axis}\). The separatrix passes through the X-point, and exhaust flows into the divertor (Section 8). The boundary shape is part of the solution.
How scpn-control Solves It¶
scpn_control.core.fusion_kernel.FusionKernel solves the Grad-Shafranov
equation on a \((R, Z)\) grid via Picard iteration with under-relaxation:
- Guess an initial \(\psi\).
- Evaluate the RHS: \(-\mu_0 R^2 p'(\psi) - FF'(\psi)\).
- Solve the resulting linear elliptic PDE (five-point stencil, direct or multigrid).
- Under-relax: \(\psi_\text{new} = \omega\,\psi_\text{solved} + (1-\omega)\,\psi_\text{old}\).
- Repeat until \(\|\psi_\text{new} - \psi_\text{old}\| < \epsilon\).
Both fixed-boundary and free-boundary variants are implemented. The free-boundary solver couples external coil currents to the Green's function response of the vacuum region.
Profiles support L-mode (linear) and H-mode (modified-tanh pedestal) via
scpn_control.core.pedestal.PedestalProfile.
4. Transport: Why Plasma Loses Energy¶
Three Transport Regimes¶
Plasma loses energy by cross-field transport of heat and particles. Three mechanisms operate, differing by orders of magnitude:
Classical transport. Binary Coulomb collisions cause particles to step across field lines by one Larmor radius per collision. Predicts diffusivities \(D \sim \rho_L^2 / \tau_\text{col} \sim 10^{-3}\) m\(^2\)/s. Far too weak to explain observations.
Neoclassical transport. In toroidal geometry, particles with low parallel velocity get trapped in local magnetic mirrors on the outboard side (where \(B\) is weaker). Their orbits project as banana-shaped trajectories in the poloidal plane, with width \(\Delta_b \approx q \rho_L / \sqrt{\varepsilon}\). Collisions scatter particles across these wider orbits, enhancing diffusion by a factor \(\sim q^2/\varepsilon^{3/2}\) over classical. Neoclassical theory predicts \(\chi_\text{neo} \sim 0.1\) m\(^2\)/s -- still too small.
Anomalous transport. Micro-instabilities (small-scale turbulence) drive fluctuations in the electric and magnetic fields. The resulting \(\tilde{E} \times B\) drifts carry particles and heat across field lines at rates \(\chi_\text{anom} \sim 1\text{--}10\) m\(^2\)/s. This is the dominant loss channel in all tokamaks.
Micro-Instabilities¶
Three principal turbulent modes drive anomalous transport:
| Mode | Full name | Driven by | Scale |
|---|---|---|---|
| ITG | Ion Temperature Gradient | \(\nabla T_i\) | \(k_\perp \rho_i \sim 0.3\) |
| TEM | Trapped Electron Mode | \(\nabla n_e\), \(\nabla T_e\) | \(k_\perp \rho_i \sim 0.5\) |
| ETG | Electron Temperature Gradient | \(\nabla T_e\) | \(k_\perp \rho_e \sim 0.3\) |
Each mode becomes unstable when its driving gradient exceeds a critical threshold. This produces stiff transport: the gradient is clamped near the threshold because any excess gradient drives turbulence that flattens the profile.
Critical Gradient Model¶
A practical model for turbulent diffusivity:
where \(\kappa_c\) is the critical gradient, \(\chi_s\) is the stiffness coefficient, and \(\alpha \sim 1\). Above threshold, transport increases sharply, preventing the gradient from rising much further.
Energy Confinement Time¶
The global energy confinement time is:
where \(W = \frac{3}{2} \int (n_e T_e + n_i T_i)\,dV\) is the stored thermal energy and \(P_\text{loss} = P_\text{heating} - dW/dt\) is the power crossing the separatrix.
Empirically, \(\tau_E\) scales with machine parameters. The standard H-mode scaling law IPB98(y,2) (ITER Physics Basis, Nuclear Fusion 39, 2175, 1999) is:
with \(I_p\) in MA, \(B_T\) in T, \(\bar{n}_{e,19}\) in \(10^{19}\) m\(^{-3}\), \(P_\text{loss}\) in MW, \(R\) in m, and \(M\) in AMU.
H-Mode vs L-Mode¶
At sufficiently high heating power, the plasma spontaneously transitions from L-mode (Low confinement) to H-mode (High confinement). The transition is triggered when the edge \(\mathbf{E} \times \mathbf{B}\) shear suppresses turbulence, forming a steep-gradient region called the pedestal at the plasma edge (\(\rho \approx 0.9\text{--}1.0\)). H-mode roughly doubles \(\tau_E\).
The pedestal height (pressure at the top of the barrier) sets a boundary condition that propagates inward through profile stiffness, determining core performance.
How scpn-control Models Transport¶
scpn_control.core.integrated_transport_solver.TransportSolver (subclass of
FusionKernel) evolves 1D radial profiles of \(T_e\), \(T_i\), \(n_e\) using the
diffusion equation with source terms:
Transport coefficients come from
scpn_control.core.gyrokinetic_transport.GyrokineticTransportModel, which
implements a reduced quasilinear dispersion model for ITG, TEM, and ETG modes
(parameters in GyrokineticsParams, fluxes in TransportFluxes).
Scaling law predictions use scpn_control.core.scaling_laws.ipb98y2_tau_e and
ipb98y2_with_uncertainty.
5. Fusion Power and the Burning Plasma¶
The DT Reaction¶
The most accessible fusion reaction is:
Total energy release: 17.58 MeV per reaction. The alpha particle (\(^4\)He) is charged and confined by the magnetic field, depositing its 3.52 MeV in the plasma (alpha heating). The neutron escapes and deposits its energy in the surrounding blanket/wall structure.
The volumetric fusion power density is:
where \(\langle\sigma v\rangle_\text{DT}\) is the Maxwellian-averaged DT reactivity (a function of ion temperature \(T_i\)). The reactivity peaks near \(T_i \approx 65\) keV and has a broad maximum above \(\sim 10\) keV.
Lawson Criterion¶
For a plasma to produce net energy, the fusion power must exceed losses. The triple product quantifies this:
This is the Lawson criterion. Achieving it requires simultaneously high density, temperature, and confinement time.
Q Factor¶
The fusion gain:
- \(Q = 1\): breakeven (fusion power = input power)
- \(Q = 5\): alpha heating equals external heating (burning plasma threshold)
- \(Q = 10\): ITER's target
- \(Q = \infty\): ignition (no external heating needed)
Alpha Heating and Self-Sustaining Burn¶
At \(Q \gtrsim 5\), the 3.52 MeV alpha particles deposit enough energy in the plasma to sustain the temperature without relying on external heating. This is the burning plasma regime. ITER aims to demonstrate \(Q = 10\) sustained for 300--500 s.
Alpha heating creates a positive feedback loop: more fusion \(\rightarrow\) more alphas \(\rightarrow\) higher \(T\) \(\rightarrow\) more fusion. Controlling this without thermal runaway requires active feedback on heating power and density.
Radiation Losses¶
Two radiation channels compete with heating:
Bremsstrahlung. Free electrons decelerated by ion Coulomb fields radiate:
This loss grows with impurity content (\(Z_\text{eff}\)).
Synchrotron radiation. Electrons gyrating in the magnetic field emit cyclotron radiation. At reactor-relevant temperatures (\(T_e > 20\) keV) and fields (\(B > 5\) T), synchrotron losses become non-negligible. Partial reflection from the wall reduces the net loss.
How scpn-control Implements It¶
Fusion reactivity: scpn_control.core.integrated_transport_solver computes
the Bosch-Hale parametrisation of \(\langle\sigma v\rangle_\text{DT}\) via
bosch_hale_reactivity. UQ-aware scaling law predictions are in
scpn_control.core.scaling_laws.ipb98y2_with_uncertainty.
6. Instabilities and Disruptions¶
MHD Instabilities¶
Magnetohydrodynamic (MHD) instabilities arise when the plasma current or pressure profile exceeds stability limits. Three categories matter most:
Tearing Modes / Neoclassical Tearing Modes (NTMs). At rational surfaces where \(q = m/n\), resistive reconnection tears the field into magnetic islands. The island width \(w\) evolves according to the Modified Rutherford Equation:
where \(\Delta'\) is the classical tearing stability index, \(a_\text{bs}\) is the bootstrap drive (neoclassical effect), \(w_d\) is the threshold island width, and \(f_\text{ECCD}\) is the ECCD stabilisation factor. NTMs can grow to large islands that degrade confinement or trigger disruptions.
Kink Modes. Current-driven instabilities that distort the plasma column. The internal kink (\(m=1, n=1\)) produces sawtooth oscillations. External kinks limit the total plasma current (Kruskal-Shafranov: \(q_\text{edge} > 1\)).
Ballooning Modes. Pressure-driven instabilities localised on the outboard (bad-curvature) side. They set the maximum pressure gradient:
Exceeding \(\alpha_\text{crit}\) triggers ideal ballooning, leading to an Edge-Localised Mode (ELM) in H-mode.
Vertical Displacement Events (VDEs)¶
Elongated plasmas (\(\kappa > 1\)) are vertically unstable: a small vertical displacement grows exponentially on the wall time \(\tau_w \sim 5\)--50 ms. Active feedback coils must counteract this in real time. Loss of vertical control leads to a VDE, where the plasma drifts into the wall, causing a rapid current quench and large halo currents that exert forces on the vessel.
The vertical instability growth rate for an elongated plasma scales roughly as:
where \(n_\text{stab}\) is the stability index (must remain \(< 1\) for passive stability).
Resistive Wall Modes (RWMs)¶
Above the no-wall beta limit \(\beta_{N,\text{nowall}}\), the plasma is stabilised only by eddy currents in the conducting wall. As these currents decay (on timescale \(\tau_\text{wall}\)), the mode grows at rate:
Active feedback with external coils can stabilise RWMs up to the ideal-wall limit \(\beta_{N,\text{wall}}\).
Disruptions¶
A disruption is an abrupt loss of confinement and plasma current. The sequence:
- Instability trigger: a locked tearing mode, density limit violation, or VDE.
- Thermal quench (\(\sim 1\) ms): stored thermal energy dumps onto the wall. Heat loads can reach \(\sim\)GW/m\(^2\).
- Current quench (\(\sim 10\) ms): the plasma current decays rapidly, inducing eddy currents and electromagnetic forces on the vessel (\(\sim\)MN for ITER).
- Halo currents: current flows through the scrape-off layer into the wall, producing asymmetric forces and torques.
Greenwald Density Limit¶
Empirically, the line-averaged density is limited by:
with \(I_p\) in MA and \(a\) in m. Exceeding this limit triggers a radiative collapse and disruption.
How scpn-control Handles It¶
- NTM dynamics:
scpn_control.core.ntm_dynamics.NTMIslandDynamicsevolves the Modified Rutherford Equation.NTMControlleruses ECCD for stabilisation.scpn_control.core.stability_mhdprovides ballooning and NTM seeding analysis. - Vertical stability:
scpn_control.control.sliding_mode_vertical.VerticalStabilizer(super-twisting sliding mode controller) andscpn_control.control.rzip_model.VerticalStabilityAnalysis. - RWM feedback:
scpn_control.control.rwm_feedback.RWMPhysics(growth rate model),RWMFeedbackController,RWMStabilityAnalysis. - Disruption prediction:
scpn_control.control.disruption_predictorprovides physics-based synthetic disruption generation (simulate_tearing_mode) and aDisruptionTransformer(sequence model for real-time prediction). - Disruption mitigation:
scpn_control.control.spi_mitigation.DisruptionMitigationController(Shattered Pellet Injection) andscpn_control.control.halo_re_physics(halo current and runaway electron models).
7. Plasma Heating and Current Drive¶
Ohmic Heating¶
The plasma current \(I_p\) heats the plasma resistively:
As the plasma heats up, resistivity drops and ohmic heating becomes inefficient above \(\sim\)3 keV. All reactor-grade plasmas require auxiliary heating.
Neutral Beam Injection (NBI)¶
High-energy neutral atoms (typically deuterium at 50--1000 keV) are injected into the plasma. Upon entering, they ionise through charge exchange and impact ionisation, becoming fast ions that transfer energy to the thermal plasma via Coulomb collisions.
NBI also drives current: the injected fast ions carry momentum preferentially in one direction. The current drive efficiency depends on beam energy, plasma density, and tangency radius.
Electron Cyclotron Current Drive (ECCD)¶
Microwave beams at the electron cyclotron frequency \(\omega_{ce} = eB/(m_e)\) (or harmonics) resonantly heat electrons. In a tokamak, \(B\) varies as \(1/R\), so different radial locations have different resonant frequencies. This allows localised deposition -- critical for NTM stabilisation (Section 6), where ECCD must be aimed precisely at the island O-point on a rational surface.
The electron cyclotron frequency at \(B = 5\) T:
ECCD also drives current through the Fisch-Boozer mechanism: preferential heating of electrons moving in one direction creates an asymmetric resistivity.
Bootstrap Current¶
In a tokamak with a pressure gradient, trapped particles (banana orbits) carry a net current without any external drive. This bootstrap current is:
Bootstrap current is "free" -- it reduces the transformer flux needed to sustain \(I_p\). In advanced scenarios, bootstrap current provides 50--90% of the total current, enabling steady-state operation. The Sauter model provides the standard neoclassical calculation.
How scpn-control Models It¶
scpn_control.core.current_drive.ECCDSource: Gaussian power deposition and driven current density profiles.scpn_control.core.current_drive.NBISource: beam heating and current drive with tangency radius and beam energy.scpn_control.core.current_drive.CurrentDriveMix: combines multiple sources (ECCD, NBI, bootstrap) into a total current profile.- Bootstrap current:
TransportSolver.calculate_bootstrap_currentwith both a simplified estimate and the full Sauter neoclassical model (calculate_sauter_bootstrap_current_full).
8. Scrape-Off Layer and Divertor¶
Plasma Boundary¶
The last closed flux surface (LCFS) separates the confined plasma from the scrape-off layer (SOL). Outside the LCFS, field lines are open: they intersect material surfaces (the divertor targets). Particles and heat stream along these open field lines to the targets.
Confined Region SOL Wall
┌──────────────────┐ ┌───────┐ ┌──────┐
│ nested flux │ │ open │ │ │
│ surfaces │ │ field │ │ │
│ T ~ 10 keV │ │ lines │ │ │
│ n ~ 10²⁰/m³ │ │ │ │ │
└──────────────────┘ └───┬───┘ └──────┘
LCFS │
┌────▼────┐
│ Divertor │
│ Target │
└─────────┘
Two-Point Model¶
The simplest SOL model connects two locations along a field line:
- Upstream (midplane): density \(n_u\), temperature \(T_u\).
- Target (divertor plate): density \(n_t\), temperature \(T_t\).
Parallel heat conduction along the field line (Spitzer-Harm):
where \(L_\parallel = \pi q_{95} R_0\) is the connection length and \(\kappa_0 \approx 2000\) W m\(^{-1}\) eV\(^{-7/2}\) is the electron parallel conductivity.
Heat Flux Width¶
The heat flux arriving at the divertor decays exponentially across the SOL:
where \(\lambda_q\) is the SOL heat flux width. The Eich scaling (Eich et al., Nuclear Fusion 53, 093031, 2013):
For ITER, \(\lambda_q \approx 1\) mm. The entire exhaust power (\(\sim\)100 MW) passes through a strip about 1 mm wide mapped to the outboard midplane. This is the exhaust problem: peak heat fluxes at the divertor can exceed 10 MW/m\(^2\), beyond what any material can sustain in steady state.
Detachment¶
The solution is detachment: increase the SOL density until the plasma in front of the divertor target cools below \(\sim\)5 eV. At these temperatures, impurity radiation (seeded nitrogen, neon, or argon) and volume recombination convert directed kinetic energy into isotropic radiation, spreading the heat load over a much larger wall area.
Operationally, detachment is controlled by gas puffing (increasing upstream density) and impurity seeding. The transition occurs when \(T_\text{target}\) drops below \(\sim\)5 eV, sharply reducing the ion flux to the target.
How scpn-control Models It¶
scpn_control.core.sol_model.TwoPointSOL implements the two-point model with:
- Eich scaling for \(\lambda_q\) via
eich_heat_flux_width - Parallel heat flux mapping from upstream to target
- Spitzer conduction along the connection length
- Target conditions (\(T_t\), \(n_t\), \(q_\text{peak}\))
peak_target_heat_fluxfor divertor load estimation
Further Reading¶
-
Wesson, J. Tokamaks, 4th ed. Oxford University Press, 2011. The standard reference. Covers equilibrium, stability, transport, heating, and SOL physics. Start here.
-
Freidberg, J. P. Plasma Physics and Fusion Energy. Cambridge University Press, 2007. Excellent on MHD theory and confinement physics, with derivations accessible to a graduate student.
-
Chen, F. F. Introduction to Plasma Physics and Controlled Fusion, 3rd ed. Springer, 2016. The best introductory text. Volume 1 covers single-particle motion, fluid theory, and waves. Good for building intuition.
-
Stacey, W. M. Fusion Plasma Physics, 2nd ed. Wiley-VCH, 2012. Comprehensive treatment of fusion-specific plasma physics including transport, heating, and reactor design.
-
Hazeltine, R. D. and Meiss, J. D. Plasma Confinement. Dover,
-
Rigorous treatment of drift orbits, neoclassical transport, and MHD stability. More mathematical than Wesson.
-
ITER Physics Basis. Nuclear Fusion 39, 2175 (1999) and updates in Nuclear Fusion 47 (2007). The definitive compilation of scaling laws, stability limits, and design parameters used by the fusion community.