🔬 physics

Gradient-Based Construction of Collisionless Steady-State Guiding-Center Distributions in Tokamaks and Stellarators

This paper presents a matrix-free, JAX-accelerated residual-minimization method for constructing collisionless steady-state guiding-center distribution functions on finite grids in tokamaks and stellarators, validated through axisymmetric and three-dimensional test cases while explicitly noting that the results provide finite-grid evidence rather than establishing continuous equilibria or transport predictions.

Original authors: Jingyi Yu, Chang Liu

Published 2026-09-18

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Jingyi Yu, Chang Liu

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Inside the heart of a fusion reactor, a magnetic cage holds a seething cloud of superheated gas. To keep this reaction going, scientists must understand how energetic particles move within that cage. These particles do not simply drift in straight lines; they follow complex, winding paths dictated by the shape of the magnetic field. If the magnetic field is not perfectly symmetrical, these paths can twist and turn in ways that make it difficult to predict where the particles will be a moment later. This creates a problem for computer simulations: if scientists start a simulation with a guess about where the particles are, that guess might be wrong from the very first second. The particles would immediately begin to shift and settle, creating an artificial "relaxation" that hides the true behavior researchers are trying to study. For decades, finding a starting point that is perfectly stable in these complex magnetic cages has been a significant hurdle, especially in devices that lack perfect symmetry.

A team of researchers at Peking University has developed a new way to solve this starting problem. They created a method to adjust a guessed distribution of particles until it fits the magnetic cage perfectly, ensuring the particles stay put without artificial shifting. Their approach treats the problem as a search for the closest possible match between a starting guess and a state of perfect balance. They use a powerful mathematical tool to measure how far off the guess is, then iteratively nudge the particle positions to reduce that error. Crucially, they do this without building a massive, unwieldy map of every possible interaction, which would be too large for even the fastest computers to handle. Instead, they calculate the necessary adjustments on the fly, processing the data in small, independent slices that can be handled simultaneously by modern graphics processors.

The researchers tested their method on several different magnetic configurations, including designs that are perfectly symmetrical and others that are intentionally twisted to improve stability. In the symmetrical cases, they could compare their results against known mathematical truths to verify that their method worked correctly. They found that their approach could take a rough guess and refine it until the particles were essentially stationary, reducing the initial error by nearly 99 percent in some cases. When they applied the same technique to more complex, three-dimensional magnetic shapes, the results were similarly successful. The particles, once adjusted, showed significantly less variation as they moved along their paths, indicating that the starting state was much closer to a true equilibrium than before.

However, the researchers are careful to define the limits of what they have achieved. Their work proves that they can find a better starting point on a specific computer grid, but it does not prove that a perfect, continuous state of balance exists for every possible magnetic shape. The method works within the boundaries of the grid used for the calculation, and it does not guarantee that the particle density remains strictly positive everywhere, though the negative values found were extremely small. The study demonstrates that the new method reduces the artificial drift of particles over a finite period of time, but it stops short of claiming to solve the broader problem of predicting long-term transport or proving that a continuous equilibrium exists for all scenarios.

The power of this new technique lies in its efficiency and its ability to handle complex data without getting bogged down in massive calculations. By running the process on a single high-performance graphics card, the researchers completed a benchmark test in just over one and a half seconds, a task that took more than two minutes on a standard high-end processor. This speed allows them to test many different scenarios quickly. The method also includes a way to check its own work by tracking individual particles along their paths to see if the distribution remains stable. In their tests, the adjusted distributions showed a dramatic reduction in variation compared to the original guesses, confirming that the particles were indeed settling into a more natural state.

This work provides a practical tool for scientists designing the next generation of fusion reactors. By offering a reliable way to generate stable starting conditions, it removes a layer of uncertainty from simulations, allowing researchers to focus on the physics that truly matters. The team has made their code and data available, ensuring that others can verify their findings and apply the method to their own designs. While the quest for a perfect fusion reactor continues, this new method offers a clearer, more stable foundation from which to begin the journey, turning a difficult guess into a precise, calculated starting point.

Technical Summary: Gradient-Based Construction of Collisionless Steady-State Guiding-Center Distributions

Problem Statement
In collisionless guiding-center simulations for tokamaks and stellarators, a prescribed distribution function $F$ is stationary only if it satisfies $\dot{Z} \cdot \nabla_Z F = 0$ along unperturbed characteristics. Standard local Maxwellian distributions often fail this condition when finite-width orbits cross flux surfaces, leading to artificial initial relaxation that obscures physical responses in global energetic-particle or full- $f$ simulations. While axisymmetric systems allow the construction of stationary states using invariants ( $E, P_\phi, \mu$ ), general stellarator geometries lack a third exact toroidal invariant, making direct prescription of a stationary distribution impossible. Existing methods like direct orbit averaging suffer from interpolation dependencies and long periods near trapped-passing boundaries. This work addresses the numerical problem of finding a nearby state that minimizes a specified discrete guiding-center residual for a given input distribution and fixed discretization, without relying on a third invariant.

Methodology
The authors formulate the problem as a matrix-free residual-minimization task in a physical Jacobian metric ( $J$ -metric).

Discrete Operator: The conservative stationary residual is discretized on a tensor grid. The objective function minimizes a weighted $L_2$ norm of the augmented residual, which includes a core term and a penalty for $v_\parallel$ edge nodes. Rate weighting ( $\omega_i$ ) is applied to rescale nodes by their local characteristic rates, ensuring the objective is approximately time-dimensionless without altering the zero-residual set.
Optimization Algorithm: The solver employs Polak–Ribière+ (PR+) conjugate gradient iterations with exact quadratic line searches.
- Unconstrained Branch: Solves for the nearest state in the $J$ -metric, projecting the input onto the kernel of the discrete operator.
- Constrained Branch: An optional extension enforcing non-negativity ( $F \ge 0$ ) via projection after each line step. This changes the feasible set but leaves the residual operator unchanged.
Implementation (JAX/XLA): The method is implemented in Python using JAX/XLA. It utilizes a matrix-free approach where the operator and its adjoint are evaluated via reverse-mode automatic differentiation without assembling a four-dimensional matrix. Independent magnetic-moment ( $\mu$ ) slices are batched and compiled into a single execution graph, enabling efficient parallelism on GPUs (tested on NVIDIA RTX 5090) and CPUs (AMD EPYC 9655).
Verification Hierarchy: The workflow includes an axisymmetric verification branch (using $R-Z$ coordinates and NUBEAM/EPCoM inputs) to validate invariant references and remapping chains before applying the method to three-dimensional stellarator geometries (precise-QA, HSX, W7-X).

Key Contributions

Explicit J-Metric Formulation: The paper provides a conditional "nearest-state" interpretation of the solution, explicitly defining the metric in which the distance is minimized.
Matrix-Free PR+ Implementation: A scalable solver using JAX/XLA that avoids normal matrix assembly, handling independent $\mu$ slices via vmap and compiled loops (lax.scan).
Axisymmetric Verification: A hierarchy of tests comparing NUBEAM (local coordinates) and EPCoM (constants-of-motion) representations, demonstrating that input fidelity and remapping accuracy significantly impact the final endpoint.
Matched 3D Calculations: The first application of this specific residual-minimization framework to matched precise-QA, HSX, and W7-X configurations, providing corrected finite-time orbit diagnostics.

Results

Axisymmetric Tests: In the DIII-D axisymmetric case, the constrained (non-negative) branch successfully reduced the residual. The EPCoM representation (constructed in constants-of-motion variables) required significantly less density correction than the native NUBEAM representation, highlighting the importance of input representation fidelity.
3D Stellarator Cases: For precise-QA, HSX, and W7-X, the unconstrained PR+ solver reduced the augmented objective function by factors ranging from $10^{-5}$ to $10^{-4}$ (e.g., W7-X reduced from $5.24 \times 10^3$ to $7.54 \times 10^{-2}$ ).
Orbit Diagnostics: Independent orbit tracking (using RK4 integrators) on the resulting endpoints showed reduced variation of the distribution along characteristics compared to the input.
- Precise-QA: 79.1% reduction in variation.
- HSX: 67.5% reduction.
- W7-X: 36.3% reduction.
Limitations Observed: The unconstrained solutions exhibited small negative mass ratios (e.g., $4.3 \times 10^{-5}$ for W7-X), indicating the method is not strictly positivity-preserving without the constrained branch. Grid convergence studies were inconclusive; fine-grid runs failed to meet the target residual within the fixed step budget, leaving grid convergence and product-rule refinement unresolved.

Claims and Significance
The authors explicitly frame their contribution as a reproducible numerical construction of improved discrete initial states rather than a proof of continuous equilibrium.

Scope: The results provide finite-grid, finite-time evidence that the constructed endpoints reduce the discrete residual and the variation of the distribution along independently integrated characteristics.
Modesty of Claims: The paper explicitly states it does not establish:
- A continuous guiding-center equilibrium.
- Grid convergence or product-rule convergence.
- Strict positivity of the distribution.
- A transport prediction.
Significance: The work demonstrates that for a fixed augmented discrete operator, a "nearest state" can be audibly and reproducibly computed. It highlights that a small discrete residual does not automatically guarantee stationarity along characteristics, and that input representation fidelity is a critical factor in the success of the remapping and minimization process. The method serves as a robust tool for preparing initial conditions in global simulations where a third invariant is unavailable.

Technical Summary: Gradient-Based Construction of Collisionless Steady-State Guiding-Center Distributions

More like this