Generalized Grad-Shafranov Equation (GGS)

• GGS is a generalization of the Grad-Shafranov equation that models MHD equilibria by incorporating arbitrary flows and pressure anisotropy.
• It combines both analytic and numerical methods to address complex plasma confinement, including applications to tokamaks, spheromaks, and stellarators.
• Solution techniques such as Lie symmetry analysis, conformal mapping, and finite-element solvers enable its use in advanced confinement and astrophysical plasma studies.

The Generalized Grad-Shafranov equation (GGS) describes axisymmetric or even non-axisymmetric magnetohydrodynamic (MHD) equilibria in plasmas with flows of arbitrary direction, pressure anisotropy, and, in specialized extensions, more exotic geometric or gravitational effects. It encapsulates equilibrium in the presence of incompressible flow, poloidal and toroidal velocity components, and pressure tensors deviating from isotropy, thus generalizing the classical Grad-Shafranov equation, which assumes static, isotropic conditions. The GGS underpins the analysis and control of plasma confinement in various magnetic confinement devices (e.g., tokamaks, spherical tokamaks, spheromaks) and develops the theoretical infrastructure essential for interpreting auxiliary heating, current drive, and advanced confinement regimes.

1. Mathematical Structure and Derivation

The GGS equation generalizes the static Grad-Shafranov equation by incorporating terms for incompressible flows and pressure anisotropy. In axisymmetric cylindrical coordinates $(R,\phi,z)$, and assuming a Chew-Goldberger-Low pressure tensor $$P = p_\perp I + (p_\parallel - p_\perp) \mathbf{b}\mathbf{b}$$ with $\mathbf{b} = \mathbf{B}/|\mathbf{B}|$, the governing equations are

$$\begin{aligned} \nabla \cdot (\rho \mathbf{v}) &= 0, \ \rho (\mathbf{v} \cdot \nabla) \mathbf{v} + \nabla \cdot \mathbf{P} &= \mathbf{J} \times \mathbf{B}, \ \nabla \times \mathbf{B} &= \mu_0 \mathbf{J}, \ \nabla \cdot \mathbf{B} &= 0, \end{aligned}$$

and the incompressible-flow condition enforces $\rho = \rho(\psi)$, with $\psi(R,z)$ labelling flux surfaces. Under axisymmetry, the magnetic and velocity fields are expanded as

$$\mathbf{B} = \frac{I(\psi)}{R} \mathbf{e}_\phi + \frac{1}{R} \nabla \psi \times \mathbf{e}_\phi,$$

$$\rho \mathbf{v} = \Theta(\psi) \nabla \phi + \nabla \phi \times \nabla F(\psi).$$

Surface functions $X(\psi)$, $M_p(\psi)$ (poloidal Alfvén Mach number), $\bar{p}_s(\psi)$ (static part of the effective pressure), $\Phi(\psi)$ (electrostatic potential), $\rho(\psi)$ (density), and $\sigma(\psi)$ (pressure anisotropy parameter) are constant on each flux surface.

The GGS equation for arbitrary incompressible flow with anisotropy becomes (using primes for $d/d\psi$): $$(1-\sigma-M_p^2)\,\Delta^*\psi + \frac{1}{2}(\sigma-M_p^2)'|\nabla \psi|^2 + \frac{1}{2}\left(\frac{X^2}{1-\sigma-M_p^2}\right)' + \mu_0 R^2 \bar{p}_s' + \mu_0 \frac{R^4}{2}\left[\frac{ \rho (\Phi')^2 }{1-\sigma-M_p^2}\right]' = 0,$$ where $\Delta^* \equiv R^2 \nabla \cdot ( \nabla / R^2 )$. The Bernoulli relation for the effective pressure is

$$\bar{p}(R, \psi) = \bar{p}_s(\psi) - \rho(\psi) \left[ \frac{v^2}{2} - \frac{R^2 (\Phi')^2}{1-M_p^2} \right].$$

Specializing to parallel flows ($\Phi' = 0$), a reparametrization

$$u(\psi) = \int_0^\psi \sqrt{1-\sigma(f)-M^2(f)}\,df$$

eliminates the $|\nabla\psi|^2$ term. In $u$-variables, the GGS reduces to a form isomorphic to the static Grad-Shafranov equation (Poulipoulis et al., 2020, Evangelias et al., 2016).

2. Physical Interpretation and Parameter Space

The GGS equation contains five or six arbitrary surface functions, depending on the treatment of pressure anisotropy:

• $M_p(\psi)$: poloidal Mach number, reflects the amplitude and shear of plasma flow.
• $X(\psi)$: sets the toroidal field/current, often chosen to align with particular shaping.
• $\bar{p}_s(\psi)$: surface function for pressure.
• $\rho(\psi)$: mass density (for incompressible flows, a flux function).
• $\Phi(\psi)$: electrostatic potential, entering for non-parallel flows.
• $\sigma(\psi)$ or $\sigma_d(\psi)$: anisotropy parameter ($\sigma_d = \mu_0 (p_\parallel - p_\perp)/B^2$).

Free function choices (e.g., linear, polynomial, power-law) control equilibrium features such as pressure peaking, current-density distribution, rotation shear, and the location and structure of magnetic X-points. Pressure anisotropy enters linearly, flow quadratically, granting significant flexibility for shaping (Poulipoulis et al., 2020, Kuiroukidis et al., 2024, Evangelias et al., 2016).

The role of these parameters is summarized:

Parameter	Physical Effect	Scaling in GGS Equation
$M_p(\psi)$	Flow amplitude/shear	Quadratic in $M_p$
$X(\psi)$	Toroidal field/current	Linear/quadratic in $X$
$\bar{p}_s(\psi)$	Static pressure	Linear in $\bar{p}_s$
$\rho(\psi)$	Density (incompressibility)	Linear
$\Phi(\psi)$	Non-parallel electric field	Quadratic in $\Phi'$
$\sigma(\psi)$	Anisotropy	Linear

Anisotropy acts more strongly than flow on most equilibrium quantities, especially on current-density structure and the shaping of $p_\parallel$, $p_\perp$; Mach-number profiles control barrier-like steepening in the pressure and current profiles (Poulipoulis et al., 2020, Kuiroukidis et al., 2024, Evangelias et al., 2016).

3. Analytic, Quasi-Analytic, and Numerical Solutions

Analytic

Linear ansätze for the free functions (generalized Solov’ev, Hernegger-Maschke types) linearize the GGS, enabling explicit polynomial or Frobenius-type solutions for the poloidal flux in $r$ and $z$ (Kaltsas et al., 2019, Kuiroukidis et al., 2024, Kuiroukidis et al., 2024). The generalized Solov’ev seed solution forms the basis for constructing D-shaped, negative-triangularity, and diverted equilibria (Kuiroukidis et al., 2024, Kuiroukidis et al., 2024). Also, similarity reductions exploit scaling and conditional symmetries to reduce the GGS to ODEs for special families of group-invariant solutions (Kuiroukidis et al., 2024, Kuiroukidis et al., 2015, Cicogna et al., 2015, Nadjafikhah et al., 2011).

Quasi-Analytic

Conformal mapping transforms map the GGS onto separable coordinate systems (e.g., $(u, v)$), allowing construction of up-down asymmetric D-shaped and diverted equilibria through Sturm–Liouville ODEs in auxiliary variables. This approach enables high-fidelity quasi-analytic solutions under certain linear ansätze for the free functions (Kuiroukidis et al., 2018).

Numerical

For strongly nonlinear or general-form free functions, solutions are computed using finite-element codes (e.g., the extended HELENA solver for parallel rotation and anisotropy (Poulipoulis et al., 2015, Poulipoulis et al., 2020)) or via mesh-free deep neural networks that minimize the GGS residual (physics-informed neural networks) (Kaltsas et al., 2021). The extension of HELENA via the $u$-$\psi$ flux remapping directly leverages the structurally identical nature of the parallel-flow GGS to the standard Grad–Shafranov equation, admitting direct numerical computation of advanced-confinement equilibria.

Nonlinear translation-invariant solutions for the slab GGSE (with quartic free-function expansion) yield equilibria sensitive to localized sheared flow and equilibrium nonlinearity, capturing the onset and stabilization mechanisms pertinent to the L–H transition (Kuiroukidis et al., 2012).

4. Physical Effects of Flow and Anisotropy

Both pressure anisotropy and flows shape plasma equilibria, with the following main consequences (Poulipoulis et al., 2020, Evangelias et al., 2016):

• Current Density Modulation: Anisotropy alters $J_{\parallel}$ and $J_\phi$, with off-axis or peaked profiles introducing strong current shears or profile hollowing. For $\sigma>0$, $p_\parallel$ rises and $p_\perp$ falls; profile shear in $\sigma(\psi)$ exerts a larger effect than its amplitude.
• Pressure and $\beta$: Anisotropy and parallel flow act to modulate the effective pressure profile; $\bar{p}$ changes modestly, but $p_\parallel$/$p_\perp$ ratios can increase significantly (e.g., $p_\parallel/p_\perp \sim 1.2$ on axis for relevant tokamak cases). Toroidal $\beta_t$ is sensitive to the sign and profile of $\sigma$.
• Shaping and X-points: Additional flow and anisotropy terms provide the degrees of freedom needed for forming D-shaped, negative-triangularity, or multi-X-point equilibria. Poloidal flow steepens the pressure profile and can increase elongation and triangularity; electric-field driven (non-parallel) flow creates sheared flow layers and additional X-point structures (Kuiroukidis et al., 2024, Kuiroukidis et al., 2024).
• Stability Landscape: The classic sufficient stability criterion ($A\ge0$ everywhere throughout the plasma), which combines current-driven, shear, and flow-shear terms, is not generally satisfied for typical rotating equilibria—especially due to the dominant destabilizing role of $A_1$ (current-driven) relative to stabilizing magnetic shear $A_2$ (Poulipoulis et al., 2015, Kuiroukidis et al., 2012). However, increased sheared flow and nonlinear equilibrium shaping are synergistic in stabilizing edge-dominated H-mode-like states (Kuiroukidis et al., 2012).

5. Non-Axisymmetric and Extended Geometries

The GGS formalism generalizes beyond axisymmetry to cover non-axisymmetric equilibria (stellarators, optimized 3D tori) through an average over hidden volume-preserving symmetries. Here, the GGS equation is formulated in terms of a flux function $\psi$ associated with an $S^1$ action, and recovers the axisymmetric Grad–Shafranov equation as a special case. Force balance is satisfied only in an averaged sense; the construction thus serves as a framework for forming candidate 3D equilibria to be refined by further optimization (Burby et al., 2020).

Non-Riemannian extensions to the GGS—embedding dynamical spacetime torsion—introduce complex stream functions capturing axion–helicity couplings in magnetospheres (e.g., magnetars), with direct implications for charge and energy storage as well as gravitational-wave emission in extreme astrophysical settings (Cirilo-Lombardo, 2018).

6. Symmetry Analysis and Solution-Generating Techniques

Continuous Lie symmetries (scaling, exceptional “weak” symmetries) underpin the construction of large analytic families of GGS solutions. The symmetry structure determines possible similarity variables and invariant subspaces, leading to reductions of the GGS to ODEs and to the systematic expansion of the known classes of analytical equilibria (D-shaped, up–down asymmetric, hollow-shell, multi-X-point) (Cicogna et al., 2015, Kuiroukidis et al., 2015, Nadjafikhah et al., 2011, Kuiroukidis et al., 2024). Conditional and weak symmetries allow for symmetry reductions even when classical invariance is lost, thus providing quasi-group-invariant solution spaces with rich geometric content.

7. Applications and Future Directions

The GGS equation provides the basis for equilibrium reconstruction, controller design, and stability assessment in magnetically confined plasmas, including ITER- and NSTX-class tokamaks, spherical tokamaks, and spheromaks. Its analytic, semi-analytic, and numerical solution methodologies—spanning Lie symmetry techniques, conformal mapping, finite-element solvers, and neural networks—enable direct modeling of pressure, current, flow, and magnetic-shear shaping in modern devices.

Extension of GGS approaches to:

• free-boundary equilibria,
• highly nonlinear flow/anisotropy profiles,
• 3D stellarator optimization,
• general relativistic and non-Riemannian effects (for astrophysical plasmas), remains an active field, with open questions regarding the existence, smoothness, and stability of such generalized equilibria (Burby et al., 2020, Cirilo-Lombardo, 2018).

The GGS formalism, through both its mathematical generality and physical flexibility, underwrites the ability to explore and optimize advanced plasma confinement scenarios in theory and computational practice (Poulipoulis et al., 2020, Kuiroukidis et al., 2024, Evangelias et al., 2016, Kaltsas et al., 2021).