Free-Boundary Grad-Shafranov Evolutive (FGE)

• FGE is a numerical framework that self-consistently computes tokamak plasma equilibrium by solving the Grad–Shafranov equation alongside external circuit dynamics.
• It uses advanced techniques like finite element discretization, Newton-based solvers, and deflation methods to track multiple equilibrium branches and ensure convergence.
• Real-time reconstruction and integrated modeling in FGE frameworks enable robust scenario planning, control applications, and dynamic analysis in fusion research.

The Free-Boundary Grad–Shafranov Evolutive (FGE) framework encompasses a class of numerical methods and software packages for solving the nonlinear axisymmetric equilibrium of a tokamak plasma, where the plasma boundary is not imposed but self-consistently determined in concert with external conductors. The FGE paradigm unifies the solution of the Grad–Shafranov equation, passive/active circuit dynamics, and, where necessary, plasma transport and profile evolution. Codes including FGE (Heiß et al., 7 Dec 2025), FreeGSNKE (Pentland et al., 7 Mar 2025), NSFsim (Clark et al., 2024), Equinox (0909.4474), TokaMaker (Hansen et al., 2023), and FEQIS (Fable et al., 2024) all implement distinct FGE approaches, ranging from real-time reconstruction to high-fidelity scenario modeling and control-oriented simulation.

1. Mathematical Structure of the Free-Boundary Grad–Shafranov Problem

FGE codes are centered around the Grad–Shafranov (GS) equation, a second-order nonlinear elliptic PDE for the poloidal flux $\psi(R, Z)$: $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$ where $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$ is the linear GS operator, $\mu_0$ is the magnetic permeability, and $J_\phi$ is the toroidal current density, defined within the unknown plasma region $\Omega_p$ via

$$J_\phi = R\frac{dp}{d\psi} + \frac{1}{\mu_0 R} F \frac{dF}{d\psi}$$

with $p(\psi)$ the pressure and $F(\psi)=RB_\phi$ the toroidal field function. For free-boundary equilibria, $\psi$ on the computational boundary $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$0 is not prescribed but satisfies an integral Dirichlet condition reflecting the vacuum field from plasma and conductor sources: $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$1 where $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$2 is the axisymmetric Green's function of $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$3.

Passive and active external circuits introduce coupled ODEs for conductor currents, with mutual and self-inductance matrices linking the plasma and circuit dynamics. For example (Fable et al., 2024, Heiß et al., 7 Dec 2025): $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$4 with the plasma-induced flux $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$5 entering as a "back-emf" term.

2. Numerical Methods and Iterative Schemes

Spatial discretization uses either finite elements (TokaMaker (Hansen et al., 2023), Equinox (0909.4474)), finite differences (FreeGSNKE (Pentland et al., 7 Mar 2025), FEQIS (Fable et al., 2024)), or hybrid approaches, often relying on unstructured meshes to faithfully capture device geometry. Boundary conditions are imposed via explicit integral coupling (Green's function methods, e.g., FEQIS (Fable et al., 2024)) or inductance-matrix techniques.

Most FGE codes employ Newton-based nonlinear solvers. The global unknown—which may be $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$6 or $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$7—is updated at each step via Krylov or direct methods, with Jacobians computed analytically or using Jacobian-free Newton–Krylov (JFNK) approaches (Pentland et al., 7 Mar 2025, Heiß et al., 7 Dec 2025).

Time evolution, when included, is handled by implicit (backward Euler) schemes for stiff resistive or MHD timescales. The solution at each time slice involves nested loops, alternately updating plasma profiles, solving the GS equation with updated currents, and evolving circuit ODEs (Clark et al., 2024, Fable et al., 2024, Heiß et al., 7 Dec 2025).

Typical FGE iterative steps include:

• Calculate $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$8 given current profiles and coil currents.
• Update plasma boundary by identifying the last closed flux surface (LCFS) from $\Delta^* \psi = -\mu_0 R J_\phi[\psi; R, Z]$9.
• Restrict $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$0 to the updated plasma region and repeat until convergence.
• Couple with external circuits, evolving $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$1 using the latest $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$2 and mutual inductance matrices.
• For dynamic evolution, update transport profiles and repeat.

3. Handling Multiple Solutions and Branch Tracking

Recent findings conclusively demonstrate that the free-boundary GS problem admits multiple physically distinct equilibrium solutions under realistic geometry and integral boundary coupling (Pentland et al., 7 Mar 2025). Detection and tracking of these secondary (disconnected) solution branches utilize deflated continuation methods, where a deflation operator modulates the nonlinear residual to "repel" roots already found: $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$3 and the deflated residual becomes $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$4. This approach enables robust exploration of bifurcations and termination points as control parameters (e.g., total current, profile coefficients, coil setpoints) are scanned.

Empirically, typical scenarios display two dominant branches:

• Deeply-confined (diverted) equilibria with the LCFS well inside the wall and strong core current.
• Shallowly-confined (limited) equilibria where the edge of the LCFS contacts the wall, with smaller total flux.

Tracking diagnostics such as magnetic-axis flux $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$5 and boundary flux $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$6 yields bifurcation diagrams that characterize branch merging, stability, and existence limits imposed by the global integral boundary constraint.

4. Integration of Current Diffusion and Transport Physics

A distinguishing feature of recent FGE solvers is the self-consistent integration of current diffusion and transport equations, allowing simultaneous evolution of poloidal flux, plasma pressure, and toroidal-field profiles (Heiß et al., 7 Dec 2025, Clark et al., 2024). The resistive evolution of the current profile is governed by variants of the axisymmetric Ohm's law, formulated as either bulk (0D) or flux-surface-resolved (1D) current diffusion equations (CDEs). Weak-form and test-function-based implementations support profile-dependent resistivity, non-inductive currents, and feedback with the equilibrium solver.

The interplay between MHD equilibrium and transport evolution is enforced via operator splitting or tightly coupled iterative schemes. For scenario modeling and controller studies, this makes it possible to capture the effects of $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$7 evolution, bootstrap current redistribution, and edge current drive in a unified free-boundary framework (Clark et al., 2024, Heiß et al., 7 Dec 2025).

5. Physical and Computational Implications

The integral free-boundary condition introduces a global coupling between internal current sources and the vacuum region, restricting solution multiplicity compared with local Dirichlet or artificially-imposed boundaries (Pentland et al., 7 Mar 2025, Fable et al., 2024). Nevertheless, non-uniqueness remains: multiple branches are physically possible and, unless properly detected, secondary solutions may not be revealed by standard equilibrium codes, which typically converge only to a single root per parameter set (e.g., EFIT++, Fiesta, EFIT-TENSOR).

A plausible implication is that scenario modeling, transport, wave and MHD-stability simulations may be compromised if equilibrium branch switching—such as at a fold point in parameter space—is not properly captured. Skipping between branches during time-marching through parameter space can induce unphysical discontinuities or convergence failures in downstream applications (Pentland et al., 7 Mar 2025).

To mitigate these risks, embedding deflation or other multiroot-finding strategies within FGE solvers is advocated, enabling comprehensive exploration and characterization of all equilibrium branches for experimental or control-discriminated selection.

6. Performance, Validation, and Code Intercomparison

Recent FGE implementations have demonstrated:

• Millisecond-per-step performance at high spatial resolution (e.g., FGE: <10 ms/step for grid 28×65; FreeGSNKE: 65×65 grid) (Heiß et al., 7 Dec 2025, Pentland et al., 7 Mar 2025).
• Predictive fidelity and consistency versus EFIT, GSevolve, and SPIDER for reconstructed axis, LCFS shape, and magnetic diagnostics (Clark et al., 2024, Fable et al., 2024, Hansen et al., 2023).
• Robust real-time operation with sub-100 ms equilibrium reconstructions in Equinox, enabling nonlinearity source identification and high-throughput scenario modeling (0909.4474).
• Full state-space linearization for control-theoretic analysis, reducing the dynamical state to $$\Delta^* \equiv R \partial_R (1/R \partial_R) + \partial_Z^2$$8 (external conductors) independent of grid size, suitable for controller synthesis and model-based control studies (Heiß et al., 7 Dec 2025).

Empirically, FGE codes such as FGE, FreeGSNKE, and NSFsim reproduce both dynamic and stationary transport benchmarks (RAPTOR, KINX) to better than experimental uncertainty for profile evolution and vertical growth rates (Heiß et al., 7 Dec 2025, Clark et al., 2024).

7. Outlook, Limitations, and Future Directions

Current open problems in FGE methodology include the rigorous identification of all possible equilibria allowed by the global integral constraint, systematic cross-code benchmarking (especially in multi-axis/doublet topologies), and the extension of self-consistent transport coupling for arbitrary profiles (including turbulent and pedestal models).

Technical extensions in progress include improved handling of singularities in Green's function computation at high order, expanded parallelism, rigorous adjoint and optimization interfaces for design/control studies, and scalable integration with diagnostic data streams for experimental real-time operation (Hansen et al., 2023, Fable et al., 2024, Heiß et al., 7 Dec 2025).

The evolution of FGE is expected to underpin both integrated modeling workflows and real-time digital twins of plasma equilibrium, enabling robust scenario planning, active control, and systematic stability analysis across device generations and advanced configurations.