Grad-Shafranov Equation: MHD Equilibria

• Grad-Shafranov Equation is a fundamental elliptic PDE in MHD equilibria that determines the poloidal magnetic flux balance in plasma systems.
• It is formulated in cylindrical coordinates and captures nonlinearity from pressure and current profiles, enabling both analytic and numerical solutions.
• Applications span tokamak, stellarator, astrophysical magnetospheres, and machine learning-driven equilibrium reconstructions for real-world plasma devices.

The Grad-Shafranov Equation (GSE) is the central elliptic partial differential equation governing magnetohydrodynamic (MHD) equilibria in magnetically confined plasmas. It determines the poloidal magnetic flux, encapsulating the balance between pressure gradients and Lorentz forces in axisymmetric and certain classes of non-axisymmetric configurations, and forms the backbone of equilibrium modeling in fusion devices, astrophysical disks, and generalized MHD systems.

1. Mathematical Formulation and Physical Setting

The GSE emerges from the static MHD equilibrium equations: $$(\nabla\times\mathbf{B})\times\mathbf{B} = \nabla p,\quad\quad \nabla\cdot\mathbf{B}=0$$ where $\mathbf{B}$ is the magnetic field, and $p$ is the plasma pressure. In axisymmetric toroidal geometry, introducing the poloidal flux $\psi(R,Z)$ (in cylindrical coordinates $(R,\phi,Z)$), the GSE reads

$$\Delta^* \psi(R,Z) \equiv R\,\frac{\partial}{\partial R}\left(\frac{1}{R}\frac{\partial\psi}{\partial R}\right) + \frac{\partial^2\psi}{\partial Z^2} = -\mu_0 R^2 \frac{dp}{d\psi} - F(\psi)\frac{dF}{d\psi}$$

with $F(\psi)=RB_\phi$ the poloidal current function and $B_\phi$ the toroidal field component. The left-hand side, $\Delta^*\psi$, is the so-called toroidal Laplacian, encoding the curvature and poloidal flux diffusion properties specific to toroidal geometry (Ding et al., 24 Nov 2025, Ciro et al., 2014).

Boundary conditions are typically of Dirichlet type, reflecting the imposed value of $\psi$ on the last closed flux surface (LCFS), with shape parameterized via D-shaped, X-point, or more general contours. The equation is nonlinear unless $p(\psi)$ and $F^2(\psi)$ are linear or quadratic functions, leading to the classical Solov'ev solution and related families (Ciro et al., 2014, Kuiroukidis et al., 2018).

2. Generalizations Beyond Axisymmetry and Physical Effects

The structure of the equilibrium can be generalized along several dimensions:

• Non-axisymmetric (3D) Equilibria: The generalized Grad-Shafranov equation (GGS) incorporates hidden volume-preserving symmetries, such that any smooth 3D equilibrium off the magnetic axis admits a one-parameter family of volume-preserving diffeomorphisms generated by a divergence-free vector field $\mathbf{u}$ (Burby et al., 2020). A stream function $\psi$ is introduced via $\mathbf{u}\times\mathbf{B} = \nabla\psi$, and the GGS equation for $\psi$ becomes:

$$-\bar{\nabla}\cdot(R^{-2}\rho_{\bar{g}}\bar{\nabla}\psi) + C(\psi)\frac{\mathbf{u}}{R^2}\overline{\cdot}\bar{\nabla}\bar{\times}\left(\frac{\mathbf{u}}{R^2}\right) = \frac{p'(\psi) + R^{-2}C(\psi)C'(\psi)}{\rho_{\bar{g}}}$$

where $\bar{g}$ is an $S^1$-averaged metric, and all operations are appropriately averaged (Burby et al., 2020). In axisymmetry, the additional term vanishes, regaining the standard GSE.

• Pressure Anisotropy and Flow: With anisotropy and incompressible plasma flow, the GSE is further generalized to

$$(1-\sigma_d-M_p^2)\,\Delta^*\psi + \cdots + \mu_0 R^2\overline{p}_s' + \frac{\mu_0 R^4}{2}\left[\frac{(1-\sigma_d)\rho(\Phi')^2}{1-\sigma_d-M_p^2}\right]' = 0$$

where $M_p$ is the poloidal Mach number, $\sigma_d$ encodes the pressure anisotropy, and $\Phi$ is the electric potential function (Evangelias et al., 2016).

• Relativistic and General Relativistic Regimes: For force-free Kerr magnetospheres and non-Riemannian spacetimes (e.g., with torsion):

$$\nabla\cdot\left[\frac{\alpha}{w^2}\Lambda(\psi)\nabla\psi\right] + \frac{w^2}{\alpha^2}(\Omega_F-\omega)\frac{d\Omega_F}{d\psi}|\nabla\psi|^2 + \frac{32\pi^2}{\alpha^2w^2c^2}I\frac{dI}{d\psi} = 0$$

where $\alpha$ is the lapse, $w$ measures cylindrical radius, $\Omega_F$ is the field line angular velocity, and the inclusion of torsion fields adds pseudoscalar corrections (Cirilo-Lombardo, 2018, Mahlmann et al., 2018).

• Quasisymmetric Stellarators: Asymptotic expansions around vacuum fields yield

$$\Delta^*\Psi + \Delta^*(I-a) = -\frac{\mu_0}{B_v^2}\frac{dp}{d\Psi} - H(\Psi)$$

to $O(\epsilon)$, with coordinates and constraints arising from quasisymmetry and vacuum structure (Nikulsin et al., 20 Jan 2025).

3. Analytic, Semi-Analytic, and Similarity Solutions

Several analytic reductions of the GSE are foundational:

• Solov'ev Solutions: For linear $p(\psi)$ and quadratic $F^2(\psi)$, the general solution decomposes into particular polynomial solutions and a homogeneous part involving Bessel, Legendre, or Whittaker functions, with the eigenvalue parameter connected to safety factor, $\beta$, and the Shafranov shift (Ciro et al., 2014, Evangelias et al., 2016).
• Similarity Reductions and Symmetry Methods: Classical and nonclassical Lie symmetry analysis provides systematic reductions, yielding invariant solutions, scaling (self-similar) solutions, and explicit multipole expansions (Nadjafikhah et al., 2011, Kuiroukidis et al., 17 Jan 2024, Ferreira, 2018).
• Generalized Similarity Reductions: Direct ("nonclassical") similarity ansätze lead to exact ordinary differential equations representing D-shaped, X-point, diverted, and planetary-magnetospheric equilibria via suitable variable transformations and superposition of different scaling modes (Kuiroukidis et al., 17 Jan 2024).
• Quasi-Analytic and Conformal Mappings: For specialized profiles and boundary geometries, separation of variables after conformal transformation produces spectral expansions and rapidly converging series for complex configurations (Kuiroukidis et al., 2018, Pataki et al., 2012).

4. Numerical Solution Methodologies

A variety of high-accuracy solvers have been developed:

Method	Geometry/Boundary	Nonlinearity	Key Strengths
Mimetic Spectral Element	Axisymmetric, fixed	Nonlinear	Arbitrary order, machine-precision, exact current continuity, robust on curved/X-point meshes (Palha et al., 2015)
Hybridizable Discontinuous	Axisymmetric, fixed	Semi-linear	High-order, handles piecewise-smooth/X-point boundaries, Anderson acceleration (Sánchez-Vizuet et al., 2017)
Adaptive Newton-FEM	Axisymmetric, free	Fully nonlinear	Adaptive mesh, shape calculus for free boundary, robust for Taylor states (Serino et al., 3 Jul 2024)
Cartesian Cut-Cell	Axisymmetric, free	Nonlinear	Parallelized, embedded boundary, shape control via coil current optimization, Aitken acceleration (Liu et al., 2020)
Physics-Informed Neural OP	Fixed, parametric	Nonlinear	Millisecond inference, end-to-end operator learning, semi-supervised physical loss, Transformer–KAN architecture (Ding et al., 24 Nov 2025, Rizqan et al., 29 Apr 2025, Jang et al., 2023)
Real-Time Least-Squares/FE	Axisymmetric, free	Fully nonlinear	Reduced basis, fixed-point, real-time for equilibrium reconstruction (0909.4474)
Data-free PINNs	Fixed/flexible, param	Nonlinear	Parameterized for various shapes/physics, hard/soft BCs, error $\ll1\%$, fast inference (Jang et al., 2023, Rizqan et al., 29 Apr 2025)

Spectral, finite element, and operator-learning techniques are now robustly available for both fixed-boundary and free-boundary settings, including adaptive refinement around separatrix, enforcement of global current constraints, and acceleration techniques such as Anderson or Aitken relaxation. Quasi-analytic (conformal mapping plus separation) and multipole constructions are used for analytic benchmarking and as the core of some high-order numerical solvers (Palha et al., 2015, Pataki et al., 2012).

The actual choice of method depends crucially on the treatment of nonlinearity, geometry complexity, need for real-time performance, and requirements for current conservation and boundary accuracy.

5. Machine Learning and Physics-Informed Approaches

Recent work has established neural operator and PINN-based models for the GSE:

• Physics-Informed Neural Operators (PINO): Learn the mapping from plasma boundary shape parameters to the solution operator $\mathcal{G}[\alpha]$, incorporating PDE residuals, boundary conditions, and total current constraints into the loss function, with architectures such as Transformer-KAN NO (TKNO) achieving high accuracy $\sim0.25\%$ and sub-millisecond inference (Ding et al., 24 Nov 2025).
• PINNs and FNOs: PINNs and Fourier Neural Operators (FNOs) enforce physics constraints (residuals, BCs) either fully data-free or in semi-supervised regimes, generalizing across boundary geometries and parameter space. PINNs have demonstrated robust accuracy (error $\lesssim1\%$), rapid evaluation (ms), and successful integration with verification tools (Marabou) for formal property validation (Rizqan et al., 29 Apr 2025, Jang et al., 2023).
• Parametric/Inverse Problem Solving: Parameterized PINNs facilitate rapid gradient-based inverse design or equilibrium inference for arbitrary shapes, pressure, and current profiles, enabling, e.g., shape optimization, uncertainty quantification, and robust extrapolation (Jang et al., 2023, Rizqan et al., 29 Apr 2025).

Hybrid methods, combining traditional high-fidelity solvers with neural surrogates, have shown promise for both acceleration and high-accuracy equilibrium reconstruction in fusion control systems.

6. Applications in Plasma Physics, Astrophysics, and Beyond

The GSE is structurally foundational to equilibrium calculations in:

• Tokamaks and Spherical Tokamaks: Used for real-time equilibrium reconstruction, stability, and transport modeling, including free-boundary and X-point (divertor) configurations (Ciro et al., 2014, Palha et al., 2015, Rizqan et al., 29 Apr 2025).
• Stellarators and Quasisymmetric Devices: Quasi-symmetric versions of the GSE enable rapid exploration and optimization of hybrid axisymmetric/non-axisymmetric equilibria, providing a coarse-grained approach to 3D equilibrium design (Nikulsin et al., 20 Jan 2025).
• Astrophysical Magnetospheres: The relativistic and torsion-generalized GSE systems govern magnetar magnetospheres, black hole magnetospheres in Kerr spacetime, and astrophysical jet launching mechanisms (Cirilo-Lombardo, 2018, Mahlmann et al., 2018).
• Space Physics and Planetary Magnetospheres: Generalized similarity reductions and analytic solutions are applied to planetary field structures and flux rope modeling (Kuiroukidis et al., 17 Jan 2024).

The equation's scope is thus considerably broader than tokamak modeling and extends to any system characterized by MHD equilibrium, including complex flows, anisotropy, and relativistic effects.

7. Open Problems and Research Directions

Key frontiers include:

• Full 3D Equilibria: The search for exact 3D MHD equilibria is now recast as finding optimal hidden volume-preserving symmetries, subsuming axisymmetric, helical, and fully non-axisymmetric equilibria within a single variational or operator framework (Burby et al., 2020).
• Efficient Real-Time Control: Operator-learning and accelerated PINO/PINN surrogates are critical enablers for next-generation control, response, and optimization in experimental devices and fusion reactors (Ding et al., 24 Nov 2025, Rizqan et al., 29 Apr 2025).
• Analytic–Numerical Hybridization: Advances in quasi-analytic/spectral solvers, together with similarity and group-invariant solutions, provide templates for benchmarking, code-verification, and as basis sets for more complex equilibrium models (Kuiroukidis et al., 2018, Kuiroukidis et al., 17 Jan 2024).
• Mathematical Structure: The relationships between symmetry (classical and nonclassical), similarity reduction, and the overdetermination arising in quasisymmetric and general 3D cases are active areas of research, with implications for uniqueness, existence, and optimization (Nadjafikhah et al., 2011, Nikulsin et al., 20 Jan 2025).

The Grad-Shafranov equation thus remains a central object of plasma theory, computational physics, and applied mathematics, continuing to evolve in scope and utility with advances in both analysis and computation.