Implementation and verification of the avalanche source in a 3D full-f particle-in-cell model of relativistic electrons for studies of tokamak disruptions

Technical Summary: Implementation and Verification of the Avalanche Source in a 3D Full-f PiC Model

Problem Statement
Tokamak disruptions pose severe risks to reactor-scale devices, not only due to thermal and electromagnetic loads but also because of the acceleration of electrons to relativistic energies, known as runaway electrons (REs). These REs can multiply exponentially via knock-on collisions with thermal electrons, a process known as the RE avalanche. This phenomenon is exponentially sensitive to the pre-disruption plasma current, potentially leading to multi-megaampere RE beams in future devices like ITER, which could cause catastrophic localized wall damage.

Existing simulation tools often lack the capability to self-consistently model RE avalanche dynamics within full 3D magnetohydrodynamic (MHD) active fields. While fluid models and lower-dimensional kinetic codes exist, they often fail to resolve the full 3D MHD dynamics crucial for understanding RE transport and phase-space distributions. Furthermore, common source terms (e.g., Rosenbluth-Putvinski) rely on simplifying assumptions, such as infinite incoming energy or parallel motion, which may not accurately capture the complex pitch-angle distributions and trapped particle populations observed in realistic scenarios. There is a critical need for a high-fidelity model that treats REs kinetically in 3D fields while conserving energy and momentum during collisional processes.

Methodology
The authors implemented a novel avalanche source within JOREK, a 3D nonlinear MHD code that utilizes a hybrid fluid-kinetic framework. In this framework, the bulk plasma is treated as a conducting fluid, while REs are modeled kinetically using a full-f particle-in-cell (PiC) approach. The REs can be evolved using either full-orbit (FO) or guiding-center (GC) pushers.

The implementation consists of two primary components:

1. Energy and Momentum Conserving Knock-on Collision Operator:

   • The avalanche is modeled using Møller scattering, the relativistic description of electron-electron scattering.
   • Unlike previous source terms that assume infinite energy or parallel motion, this operator treats the incoming RE with finite energy and pitch angle.
   • The differential cross-section is integrated with a momentum cut-off ($\gamma_{min}$) to distinguish between small-angle and knock-on collisions.
   • The operator strictly conserves energy and momentum. Upon collision, a marker splits into three: the original marker (reduced weight), a scattered primary marker, and a newly generated secondary marker.
   • This approach allows for the generation of secondary REs with high pitch angles, including those that may become trapped or even move opposite to the electric field, which are subsequently removed only when their energy drops below a specific threshold.

2. Resampling Algorithm:

   • Because knock-on collisions cause an exponential increase in the number of markers (tripling them per collision event), a resampling technique was implemented to maintain computational feasibility.
   • The resampling is based on the 5D guiding-center phase space ($R, Z, \phi, p_\parallel, \mu$).
   • Spatial Binning: Markers are binned into volume elements based on the JOREK finite element grid (refined if necessary) and toroidal bins.
   • Momentum Binning: Within each spatial bin, markers are sorted into momentum bins based on parallel momentum ($p_\parallel$) and magnetic moment ($\mu$).
   • New markers are generated with weights and momenta derived from the distribution of the old markers within each bin, ensuring the total weight (density) and the momentum space distribution are preserved.
   • The algorithm supports both pseudorandom (PCG) and quasirandom (Sobol) number generators, with the latter offering better error scaling.

Key Contributions

• First Kinetic Avalanche in 3D MHD: To the authors' knowledge, this is the first model capable of simulating the RE avalanche kinetically while self-consistently evolving 3D MHD active fields.
• Conservative Collision Operator: The implementation moves beyond approximate source terms (like Rosenbluth-Putvinski or Chiu-Harvey) to a fully conservative Møller scattering operator. This enables accurate modeling of growth rates and pitch-angle distributions, including the generation of trapped REs.
• Resampling for Long-Term Simulation: The development of a 5D resampling algorithm allows for the restriction of marker counts without significantly degrading the representation of the phase-space distribution, making long-timescale simulations of re-avalanching computationally viable.
• Verification: The model was verified against analytical expressions for growth rates and compared against non-simplified source terms, demonstrating convergence to theoretical limits when time steps are sufficiently small.

Results
The model was applied to a JET-like termination scenario (based on shot #95135) to demonstrate its applicability to realistic 3D MHD active scenarios.

• Scenario Setup: The simulation began near the end of a termination event where flux surfaces were reforming. The initial RE population was derived from previous fluid simulations.
• Re-avalanching: As the magnetic flux surfaces reformed and stochasticity decreased, the model successfully captured the onset of re-avalanching.
• Spatial Distribution: The results showed that while the total RE current density was highest in the core, the secondary REs (generated by the avalanche) were primarily localized in a ring around the core. This was attributed to the higher parallel electric field in this region, which drove higher growth rates and energy gains for secondary electrons.
• Field Consistency: The study highlighted the importance of self-consistency; when the RE current was subtracted from the total current to calculate the parallel electric field ($E_\parallel$), the avalanche gain was reduced compared to cases where the field was artificially high. However, re-avalanching was still observed.
• Resampling Performance: Tests demonstrated that resampling could reduce the marker count to less than a quarter of the original while preserving 3D structures (such as magnetic islands) and maintaining current density profiles with relative errors on the order of 0.01%.

Significance and Claims
The paper positions this work as a significant milestone toward a more complete understanding of RE interactions with MHD modes during disruptions. By overcoming the limitations of fluid models and non-conservative source terms, this model provides a tool to:

• Accurately assess risks associated with RE beams in large devices like ITER.
• Investigate mitigation scenarios involving re-avalanching, a process that may be more efficient in high-current future devices than in current ones.
• Study the role of trapped REs as seeds for future avalanches.

The authors note that while the current implementation allows for uncoupled particle tracing to reduce costs, the model is designed for self-consistent coupling. Future work is required to port the code to accelerated high-performance computing systems to address the long time scales involved in periodic termination and re-avalanching, and to implement nuclear sources for studies in activated devices.