Gyrokinetic Modeling in Burning Plasmas

• Gyrokinetic modeling in burning plasmas is a first-principles approach that simulates turbulence, transport, and multi-species interactions by solving the gyrokinetic Vlasov–Maxwell equations.
• Numerical methods, including both global and flux-tube models, accurately capture electromagnetic effects and zonal flow regulation, critical for predicting fusion reactor performance.
• Hybrid frameworks integrating surrogate optimization and extended MHD–gyrokinetic coupling enable efficient, high-fidelity predictions of confinement and operational scenarios in ITER and beyond.

Gyrokinetic modeling in burning plasmas provides a first-principles framework for simulating turbulence, transport, and energetic particle-driven phenomena in magnetically confined fusion devices. Core to this approach is the high-fidelity numerical solution of the gyrokinetic Vlasov–Maxwell system for multiple species (thermal ions, electrons, alpha particles, impurities) under reactor-relevant conditions characterized by high temperatures, low collisionality, strong electromagnetic effects, and intricate interactions between microturbulence and macroscopic instabilities. Fully predictive simulations of burning plasmas are essential for optimizing fusion performance, predicting confinement properties, and assessing operational scenarios in ITER and next-generation pilot plants.

1. Governing Gyrokinetic Equations in Burning Plasmas

The gyrokinetic model evolves the non-adiabatic component of the distribution function for each species, $g_s(\mathbf{R}, v_\parallel, \mu, t)$, in a five-dimensional phase space. For reactor-relevant regimes, both electrostatic ($\phi$) and electromagnetic ($A_\parallel, B_\parallel$) fluctuations are retained. A generic local $\delta f$ gyrokinetic equation takes the form

$$\frac{\partial g_s}{\partial t} + \left(v_\parallel \mathbf{b}^* + \mathbf{v}_D + \frac{c}{B} \mathbf{b}^* \times \nabla \langle \phi \rangle\right) \cdot \nabla g_s - C[g_s] = - \frac{Z_s e}{T_s} \frac{\partial \langle \phi \rangle}{\partial t} f_{Ms}$$

where $\mathbf{b}^*$ is the modified magnetic field, $\mathbf{v}_D$ is the drift velocity, and $C[\cdot]$ is the (often Sugama or Landau–Boltzmann) collision operator, implemented to conserve density, momentum, and energy. These equations are closed by the field equations: gyrokinetic quasineutrality for $\phi$ and Ampère’s law for $A_\parallel$ (Howard et al., 25 Apr 2024, Siena et al., 5 Nov 2025, Rodriguez-Fernandez et al., 2023).

In addition to the core system, hybrid models couple reduced MHD evolution of the bulk (single-fluid or compressible models) with the gyrokinetic treatment of energetic particles, capturing both thermal ion compressibility, diamagnetic effects, and kinetic drive of shear Alfvén and related modes (Wang et al., 2010, Liu et al., 22 Feb 2024).

2. Numerical Methods and Geometry: Global and Local Approaches

Gyrokinetic simulations employ both global and local (flux-tube) methodologies:

• Global modeling: The entire plasma radius (or a substantial fraction, e.g., $\rho_{\rm tor} \in [0, 0.6]$), is simulated with field-aligned coordinates $(x, y, z)$ covering flux label, binormal, and parallel coordinates, incorporating realistic magnetic equilibria (shaping, $q$-profile) and non-uniform meshes to resolve core, edge, and geometry-driven phenomena. Codes such as GENE-Tango and GYSELA implement high-order finite difference, semi-Lagrangian, or hybrid spatial discretizations, often with non-uniform polar or mapped meshes to represent D-shaped and X-point boundary plasmas (Bouzat et al., 2017, Siena et al., 5 Nov 2025).
• Flux-tube modeling: Small radial domains centered on a flux surface with periodic boundary conditions in $(y, z)$ and fixed local equilibrium. Flux-tube runs are crucial for linear benchmarking, turbulence studies, and scanning parameter dependencies, although they can overestimate profile stiffness and miss global rational-surface effects (Siena et al., 5 Nov 2025, Rodriguez-Fernandez et al., 2023).

Advanced interpolation and gyroaveraging schemes enable accurate treatment of the $r=0$ singularity and efficient evaluation of field couplings on non-uniform grids, delivering high-order convergence and substantial memory savings (Bouzat et al., 2017).

3. Physics of Turbulence and Transport in Burning Plasmas

The dominant microinstabilities in burning plasma regimes include:

• Ion-temperature-gradient (ITG) and trapped-electron-mode (TEM) turbulence: Primary drivers of ion and electron thermal transport.
• Electromagnetic modes: Microtearing modes (MTMs), kinetic ballooning modes (KBMs), and Alfvénic ion temperature gradient (AITG) modes, all of which contribute notably to fluxes at low $k_y\rho_s$ and must be captured with high resolution and full electromagnetic modeling (Siena et al., 5 Nov 2025).
• Electron temperature gradient (ETG) modes: ETG-driven transport is typically suppressed in ITER-relevant regimes by strongly growing zonal flows which, under extremely low collisionality, quench electron heat flux on long timescales.

Burning plasma gyrokinetics consistently finds "stiff" profile response: modest changes (e.g., ±10%) in normalized temperature or density gradients produce order-unity shifts in heat flux, reflecting turbulence pinned to instability thresholds and strong nonlinear regulation by zonal flows (Siena et al., 5 Nov 2025, Rodriguez-Fernandez et al., 2023). This stiffness sets challenges and opportunities for integrated predictive modeling.

Alpha particle and fast ion stabilization of turbulence is essential in burning plasmas. Gyrokinetic studies demonstrate that non-Maxwellian (slowing-down or bi-Maxwellian) distributions must be used for accurate modeling, as Maxwellian approximations overpredict turbulence suppression and can lead to spurious confinement predictions (Siena et al., 2018).

4. Surrogate Modeling and Computational Strategies

First-principles nonlinear profile prediction is computationally intensive—global, high-dimensional simulations require $\mathcal{O}(10^5)$ GPU-hours for full-fidelity transport iteration. The PORTALS surrogate-based optimization framework accelerates this process by constructing Gaussian-process surrogates for turbulent and neoclassical fluxes at key flux surfaces ($r/a\sim0.35, 0.55, 0.75, 0.875, 0.9$), reducing the number of expensive nonlinear runs by factors of $3$–$10$. This enables efficient multi-channel (electron temperature, ion temperature, density, impurities, rotation) flux-matched profile prediction (Rodriguez-Fernandez et al., 2023, Howard et al., 25 Apr 2024).

Key aspects of PORTALS methodology include:

• Training GPs on initial sample points generated by direct CGYRO+NEO runs,
• Minimizing composite residuals to match target power and particle sources,
• Automatic uncertainty quantification by mapping GP variances to fusion power bands,
• Reusing trained surrogates for scenario scans (e.g., edge density, auxiliary power),
• Integration with highly parallel GPU-accelerated CGYRO runs ($\mathcal{O}(10^4)$ GPUs, $400$M cell-updates/sec) (Rodriguez-Fernandez et al., 2023).

5. Energetic Particle-Driven MHD and Hybrid Modeling

Energetic-particle physics, especially fast alpha-driven Alfvénic activity, is handled in extended hybrid MHD–gyrokinetic codes (e.g., GMEC, XHMGC):

• Reduced MHD equations are solved for the bulk, incorporating perturbed vorticity, magnetic potential, and pressure. The model absorbs fluid compressibility and includes kinetic and diamagnetic effects due to thermal ions and fast ions (Liu et al., 22 Feb 2024, Wang et al., 2010).
• Gyrokinetic equations for energetic particles evolve $\delta f$ in a particle-in-cell (PIC) scheme, providing pressure-coupling terms to the MHD solver and capturing resonance effects (drift resonance, Landau damping, and drive).
• Field-aligned coordinate systems and mesh twist–shift periodicity are critical for efficient resolution of high-$n$ Toroidal Alfvén Eigenmodes (TAEs) and Reversed-Shear Alfvén Eigenmodes (RSAEs).
• The generalized fishbone-like dispersion relation (GFLDR) encapsulates the interplay between field-line bending, fluid potential energy, and kinetic resonant drive, capturing the discrete gap mode structure (e.g., kinetic BAE, TAE) and the critical effect of thermal-ion compressibility in opening gaps and reducing continuum damping (Wang et al., 2010).

Benchmarking against analytic results and comprehensive eigenmode databases shows that these hybrid models quantitatively reproduce Alfvénic mode frequencies and growth rates under realistic burning plasma parameters.

6. Scenario Predictions, Sensitivities, and Operational Implications

Surrogate-accelerated nonlinear gyrokinetic profile predictions for ITER's 15 MA baseline scenario show:

• Fusion gain $Q=12.2$ with self-consistent $\alpha$-heating, external power, and transport fluxes, in line with ITER's design objectives (Siena et al., 5 Nov 2025, Howard et al., 25 Apr 2024).
• Ion to electron heat flux ratio $Q_i/Q_e > 1$ everywhere, indicating ITG-dominated turbulence.
• Density peaking (e.g., $n_e(0.2)/\langle n_e\rangle\approx1.3$) and impurity transport consistent with H-mode database values.
• Weak or negligible turbulent isotope effect: $\tau_E$ for D-T, D, and H plasmas essentially identical under matched heating, contrary to naive scaling expectations (Howard et al., 25 Apr 2024).
• Electromagnetic stabilization (especially inclusion of MTMs and KBMs) is paramount for density peaking and overall confinement.
• Sensitivity studies reveal strong dependence of transport on the safety-factor profile: flattened $q(r)$ profiles with low magnetic shear destabilize KBMs, degrading confinement by orders of magnitude (Siena et al., 5 Nov 2025).
• External toroidal rotation is insufficient, at expected ITER values, to notably reduce turbulent fluxes.

7. Numerical Innovations, Convergence, and Accuracy Benchmarks

Realistic geometry handling, including D-shaped cross-section mapping, non-uniform polar grids, and robust high-order interpolation and gyroaverage operators, is essential for representing edge effects, peaking, and global turbulence features with moderate computational resources (Bouzat et al., 2017). Recent convergence studies demonstrate that, even with 10–30% reduced grid points, global accuracy is retained, enabling high-fidelity predictions with tractable numerical cost.

Verification and validation exercises—including comparison of global and local (flux-tube) runs, direct comparison to analytic theory (e.g., kinetic Alfvén wave dispersion), and cross-code benchmarking (e.g., MEGA, M3D-C1-K, QRH)—confirm the predictive capabilities of state-of-the-art gyrokinetic and hybrid codes across the dominant dynamic range of burning plasma operation (Siena et al., 5 Nov 2025, Wang et al., 2010, Liu et al., 22 Feb 2024, Rosen et al., 2021).

---

• (Wang et al., 2010) Wang et al., "An extended hybrid magnetohydrodynamics gyrokinetic model for numerical simulation of shear Alfvén waves in burning plasmas"
• (Bouzat et al., 2017) Grandgirard et al., "Targeting realistic geometry in Tokamak code Gysela"
• (Siena et al., 2018) Di Siena et al., "Non-Maxwellian fast particle effects in gyrokinetic GENE simulations"
• (Rosen et al., 2021) Rosen et al., "An E & B Gyrokinetic Simulation Model for Kinetic Alfvén Waves in Tokamak Plasmas"
• (Rodriguez-Fernandez et al., 2023) Rodriguez-Fernandez et al., "Enhancing predictive capabilities in fusion burning plasmas through surrogate-based optimization in core transport solvers"
• (Liu et al., 22 Feb 2024) White et al., "Development of a gyrokinetic-MHD energetic particle simulation code Part II: Linear simulations of Alfvén eigenmodes driven by energetic particles"
• (Howard et al., 25 Apr 2024) Howard et al., "Prediction of Performance and Turbulence in ITER Burning Plasmas via Nonlinear Gyrokinetic Profile Prediction"
• (Siena et al., 5 Nov 2025) Citrin et al., "First global gyrokinetic profile predictions of ITER burning plasma"