Excitation Pullbacks in Plasma Simulations

• Excitation pullbacks are a numerical technique that restructures the parallel magnetic potential to eliminate large-amplitude cancellations in gyrokinetic PIC simulations.
• The method splits the parallel potential into symplectic and Hamiltonian components and applies a pullback remapping to correct the perturbed distribution function.
• This scheme enhances simulation stability, allows for significantly larger time steps, and reduces computational cost in complex plasma turbulence studies.

Excitation pullbacks are a numerical algorithmic technique developed to mitigate the large-amplitude cancellation problem encountered in electromagnetic gyrokinetic particle-in-cell (PIC) simulations, especially within the context of global plasma turbulence codes such as ORB5. The method systematically reorganizes the representation of the parallel magnetic potential in the mixed-variable Hamiltonian/symplectic formalism and corrects the perturbed distribution function, facilitating robust and efficient gyrokinetic PIC computations for electromagnetic phenomena at realistic plasma parameters (Mishchenko et al., 2018).

1. Hamiltonian/Symplectic Splitting and the Pullback Principle

In the ORB5 code, the perturbed parallel magnetic potential $A_\|(\mathbf{x},t)$ is decomposed into a “symplectic” part $A_\|^{(s)}$ and a “Hamiltonian” part $A_\|^{(h)}$, such that

$A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$

This splitting separates $A_\|^{(s)}$, which modifies the symplectic structure, from $A_\|^{(h)}$, which enters the Hamiltonian. The single-particle gyrocenter Lagrangian then reads:

$$L = \left(q\,\mathbf{A}^* + p_\|\mathbf{b}\right)\cdot\dot{\mathbf{R}} - \left(\tfrac12 m v_\|^2 + \mu B + q\,\langle\phi\rangle - q v_\| \langle A_\|^{(h)}\rangle\right),$$

where $\mathbf{A}^*=\mathbf{A}+(m/q)v_\|\mathbf{b}$, $\langle\cdot\rangle$ is the gyro-average, and $(\mathbf{R},v_\|,\mu)$ are mixed-variable coordinates.

The corresponding perturbed equations of motion, to first order in the fields, are:

$A_\|^{(s)}$0

with $A_\|^{(s)}$1 and $A_\|^{(s)}$2.

2. The Cancellation Problem and Pullback Mitigation

The electromagnetic $A_\|^{(s)}$3 PIC approach leads to an Ampère equation of the form:

$A_\|^{(s)}$4

Here, a large-amplitude numerical cancellation occurs between the field operator and the particle current term, causing severe stability and accuracy issues, particularly when high time-steps or low collisionality are present.

The pullback scheme addresses this challenge through a time-step-wise reorganization (“pullback”) of the parallel potential:

1. Remap the splitting at each time $A_\|^{(s)}$5:

$A_\|^{(s)}$6

$A_\|^{(s)}$7

2. Correct the mixed-variable perturbed distribution function:

$A_\|^{(s)}$8

This procedure ensures that nearly all of $A_\|^{(s)}$9 resides in $A_\|^{(h)}$0, eliminating problematic cancellations from $A_\|^{(h)}$1. There is also a fully nonlinear variant, which involves shifting particle velocities at constant weight (Mishchenko et al., 2018).

3. Discretization, Marker Representation, and Algorithmic Workflow

The perturbed distribution for each species is represented via a marker-in-cell approach:

$A_\|^{(h)}$2

All fields, including $A_\|^{(h)}$3, $A_\|^{(h)}$4, and $A_\|^{(h)}$5, are expanded on three-dimensional tensor-product B-splines:

$A_\|^{(h)}$6

with similar expansions for $A_\|^{(h)}$7 and $A_\|^{(h)}$8. In toroidal geometry, a toroidal Fourier transform allows independent treatment of each $A_\|^{(h)}$9 mode.

The simulation time step advances as follows:

1. Push markers with current fields via equations of motion.
2. Scatter density $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$0 and parallel current $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$1 markers to the field grid.
3. Solve quasi-neutrality, Ohm, and Ampère equations.
4. Execute the pullback:
   • Remap $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$2, $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$3.
   • Update marker weights: $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$4.
5. Run diagnostics or energy checks as needed.

This cycle is repeated until the final simulation time is reached (Mishchenko et al., 2018).

4. Application to Toroidal Alfvén Eigenmodes and Internal Kink Modes

The pullback scheme has been verified in linear and nonlinear electromagnetic simulations of Toroidal Alfvén Eigenmodes (TAEs) and internal kink modes. In the ITPA-TAE benchmark, linear growth rates $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$5 and frequencies $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$6 exhibit the scaling:

$A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$7

where $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$8 is the fast-ion temperature and $A_\| = A_\|^{(s)} + A_\|^{(h)} \,.$9 is the fast-ion density.

A quantitative comparison for $A_\|^{(s)}$0 values of $A_\|^{(s)}$1 $A_\|^{(s)}$2 and $A_\|^{(s)}$3keV shows that results with pullback (PB) and the conventional control variate (CV) scheme agree within numerical error. The pullback allows time steps $A_\|^{(s)}$4 up to five times larger than control variate only.

$A_\|^{(s)}$5 [keV]	CV: $A_\|^{(s)}$6	PB: $A_\|^{(s)}$7	CV: $A_\|^{(s)}$8	PB: $A_\|^{(s)}$9
200	0.0120	0.0121	0.50	0.50
400	0.0168	0.0169	0.52	0.52
800	0.0205	0.0206	0.54	0.54

In nonlinear scenarios involving only wave-particle nonlinearity, the mode saturates by flattening the fast-ion distribution at resonance.

For the $A_\|^{(h)}$0 internal kink mode, control-variates-only runs at practical $A_\|^{(h)}$1 produce explosive, unphysical growth after a few steps, due to cancellation noise. In contrast, the pullback scheme produces clean physical mode evolution and radial mode structure matching analytic predictions. PB allows stable calculation at $A_\|^{(h)}$2, compared to the much smaller step size required with CV ($A_\|^{(h)}$3), resulting in up to a tenfold speedup with identical accuracy (Mishchenko et al., 2018).

5. Computational Performance and Robustness

Performance evaluations demonstrate that the pullback scheme enables accurate results for $A_\|^{(h)}$4, surpassing the stability threshold of CV ($A_\|^{(h)}$5). For fixed accuracy (e.g., $A_\|^{(h)}$6), PB needs five times fewer time steps and reduces CPU usage by a factor of 3–4, with the overhead for the pullback (scatter and weight update) constituting less than 5% of total computational cost.

The approach is most advantageous for high-$A_\|^{(h)}$7 electromagnetic simulations and when kinetic electrons are treated in the drift-kinetic approximation. For highly stiff regimes, such as very low electron-ion collision frequency ($A_\|^{(h)}$8) or fully kinetic electrons, a nonlinear pullback involving velocity shifts may be necessary (Mishchenko et al., 2018).

6. Algorithmic Extensions and Regime Limitations

Extensions to the pullback approach include incorporating collision operators into the pullback transformation, application to hybrid MHD-gyrokinetic models, and generalization to stellarator geometries. While highly effective in regimes with cancellation-dominated noise, the efficiency gain is most pronounced when conventional control variate techniques are insufficient to ensure stability at reasonable time steps. In exceedingly stiff configurations, further nonlinear algorithmic development is required to fully preserve accuracy and avoid residual noise.

7. Significance in Global Gyrokinetic Simulations

The excitation pullback scheme has enabled first-principles global gyrokinetic PIC simulations of challenging electromagnetic phenomena, such as realistic β internal kink modes, at computational cost and accuracy unattainable by previous methods. Its elimination of the Ampère equation’s large numerical cancellation, within the mixed-symplectic/Hamiltonian formulation, represents a robust methodological advance for plasma simulation codes such as ORB5. This suggests its applicability to a broader class of electromagnetic simulation challenges, particularly in regimes dominated by strong particle-field interaction and numerical instability due to cancellation effects (Mishchenko et al., 2018).