Energetic particles transport in constants of motion space due to collisions in tokamak plasmas

The spatio-temporal evolution of the energetic particles in the transport time scale in tokamak plasmas is a key issue of the plasmas confinement, especially in burning plasmas. In order to include sources and sinks and collisional slowing down processes, a new solver, ATEP-3D was implemented to simulate the evolution of the EP distribution in the three-dimensional constants of motion (CoM) space. The Fokker-Planck collision operator represented in the CoM space is derived and numerically calculated. The collision coefficients are averaged over the unperturbed orbits to capture the fundamental properties of EPs. ATEP-3D is fully embedded in ITER IMAS framework and combined with the LIGKA/HAGIS codes. The finite volume method and the implicit Crank-Nicholson scheme are adopted due to their optimal numerical properties for transport time scale studies. ATEP-3D allows the analysis of the particle and power balance with the source and sink during the transport process to evaluate the EP confinement properties.


Introduction
In tokamak plasmas, energetic particles (EPs) can be generated by fusion products, neutral beam injection, and other auxiliary heating with energy much greater than those of the background ions and electrons.The EPs can heat the plasma through Coulomb collisions with thermal particles, and improved confinements of EPs are needed for achieving high fusion gain.As the particle moves along its equilibrium trajectory in tokamak, its orbit can be described in the constants of motion (CoM) space [1].The CoM generally used in gyrokinetic simulations of fusion plasma are particle energy E, toroidal canonical momentum P ζ , and magnetic moment µ, which is determined by the magnetic field configuration.EPs' orbit can have a large radial excursion (∼ 0.1a) and the finite orbit width effect is important to understand the EP transport.The CoM variables provide an appropriate way to describe the EP dynamics.Specifically, in full f simulations the adoption of the CoM space is necessary to initialize a steady state distribution with the finite orbit width effects taken into account [2].A parametric equilibrium distribution function in the CoM space is applicable to gyrokinetic studies investigating the behavior of guiding centers with finite orbit width in axisymmetric tokamak plasmas, as discussed in [3].The characteristic velocities of EPs are close to the phase velocities of Alfvén eigenmodes (AEs), which EPs can destabilize through resonant interactions [4].The AE induced EP transport is more efficient and faster than the collisional transport, as summarized previously [5].However, collisions enhance the resonance width by moving particles in and out of the resonance islands [6] and affect the mode evolution [7].On a longer timescale, the distribution function is shaped by sources, sinks, wave-induced transport, and collisions, which need to be considered in the EP transport model.
The neoclassical transport and the bootstrap current induced by alpha particles and NBI were theoretically studied in Refs.[8,9], where was found that the bootstrap current induced by EP is about 10 ∼ 20% of the background bootstrap current.The ion poloidal gyroradius ρ p , which indicates the finite orbit width, can be comparable to or even larger than the characteristic length of the equilibrium L eq , which consequently violates the assumption in the traditional neoclassical theory based on the asymptotic solution of the drift equation using ρ p [10].It has been studied that the small poloidal gyroradius assumption of standard neoclassical theory breaks down due to the wide orbit width analytically [11] and numerically [12,13].Indeed, the neoclassical transport in the large ρ p /L eq regime has attracted attention due to its possible connection to the L-H transition [14].In addition, the collision effect is important in the reduced transport models to study the long-timescale steady-state profiles.In SOLPS [15] and ASTRA [16], the collision transport is simulated with coefficients predicted by neoclassical theory or matched to the experimental measurement.The transport of runaway electrons has been studied based on a fluid-kinetic code DREAM [17].For most codes using the Monte-Carlo collision operator, collisions represented in the velocity space (v ∥ , v ⊥ ) [18] or the pitch angle scattering [19] are adopted [20,21].Linearized collision operator in CoM coordinates E, µ is described and tested in [22].Besides the collisional process, a general theoretical framework for the transport of Phase Space Zonal Structures (PSZS) has been developed [1,23,24,25].PSZS is the long-lived toroidal symmetric (n = 0) structure that defines the nonlinear equilibrium in the presence of fluctuations such as Alfvénic instabilities.The PSZS theory is used to study the fluctuation-induced EP transport and in general includes collisional effects.Based on this framework, a reduced energetic particle (EP) transport code (named 'ATEP') is developed.The primary objective of ATEP is to systematically encompass the intricate physics embedded in the overarching equations, as detailed in [26,27].The general ATEP model addresses comprehensive physics, recovering critical gradient models [28], the 'kick' model [29], and quasi-linear resonance broadening models [30,31] within appropriate limits.
In this paper, we implement and investigate EP transport in the CoM space due to collisions.In order to include sources and sinks and collisional slowing down processes, a new solver, ATEP-3D is implemented for the energetic particle transport modeling on the transport time scale.The collision operator is represented in the three-dimensional (3D) CoM space and averaged over the unperturbed orbits.The collision frequencies are given by the sum of the contributions from different sources.The finite volume method (FVM) and the implicit scheme are adopted for optimized numerical properties.ATEP-3D allows the analysis of the particle and power balance with the source and sink during the transport process to evaluate the EP confinement properties.The following part of this article is organized as follows.The physics model and the numerical schemes are givens in Sections 2 and 3 respectively.The description of the experimental case and the simulations results are represented in Section 4, followed by Section 5 as the conclusion part.

Physics models
The evolution equations for PSZS are derived in a recent work, which describe the slow evolution of macroscopic plasma profiles and the more rapid phase space corrugations on meso-scales due to fluctuations in toroidal geometry [23].The nonlinear gyrokinetic theory is utilized to define a neighboring nonlinear equilibrium consistent with toroidal symmetry-breaking perturbations.It introduces PSZS as a means to capture the fluctuation-induced corrugations of the reference state [23].The collisional transport and its synergistic interplay with fluctuation-induced transport are important to EP transport studies.Consequently, the focus of this work is to derive the collisional transport in the CoM space.By exploring this aspect, we aim to gain a better understanding of how collisions influence the transport of EPs and how they interact with the transport induced by fluctuations.This analysis will contribute to a comprehensive understanding of the complex mechanisms governing EP transport in fusion plasmas.Previously, in the CGM (critical gradient model) [28] a 1D radial diffusion coefficient was used with a recipe for the onset of transport (the growth rate of the toroidal Alfvén eigenmode γ TAE > E × B shearing rate) and an ad hoc saturation rule assuming a diffusion coefficient D = 1 ∼ 10m 2 /s [32].No phase space dependence is considered, with the assumption of diffusive transport, despite sufficient evidence of non-diffusive EP transport in theory and experiment.The kick model [29] resolves the phase space dynamics, but the self-consistent loop between the mode amplitude, the evolution of EP distribution, and background profiles is not taken into account.More fundamentally, the kick model is a probabilistic approach, while the PSZS theory solves transport equations in phase space.However, the comparison is an important step for model validation in the future [27].The RBQ (Resonance Broadened Quasilinear) model is built upon the foundation of the quasilinear and diffusive model, expanding to include the broadened resonance inherent in wave-particle interactions.The RBQ code is reduced into a 1D P ζ [30] and a 2D (E, P ζ ) [33] CoM space.Moreover, collisions are simplified into a diffusive collision operator.Our model aims to capture the 3D phase structure of EPs and the nonlinear interactions between EPs, the background, and the mode evolution.The distribution function of EPs in the CoM space is modeled in the form where f (X, t) is the EP distribution, X denotes the CoM phase space coordinates, the C(f ) is the collision operator, and γ(X, t) and S(X, t) are the sink/source terms.In tokamaks, an unperturbed particle trajectory is periodic in three canonical angles and can be described by the three CoM: particle energy E, toroidal canonical momentum P ζ and magnetic moment µ.We choose X = (P ζ , E, Λ), where Λ = µB 0 /E is also a constant of motion.

Collision operator in
The linearized collision operator appeared in early literature as the approximate Fokker-Planck collision operator for the transport theory applications [34].The CoM orbit representation is especially useful for treating EPs that have large orbits [35].The neoclassical transport studies considering large orbits width is derived [11].Transport fluxes formulated using CoM first appeared in [36] and later in [37].In this work, we derive the linearized collision operator in the (P ζ , E, Λ) CoM coordinates.We assume that the energetic particles generated by NBI or fusion processes are immediately ionized, and we take into account the Coulomb collisions between charged particles.These EPs exhibit azimuthal symmetry with respect to the magnetic field.The collision can be derived from the Rosenbluth expression [38].Inter-species EP collisions are neglected due to their high energy and low density.The collisions experienced by EPs primarily involve interactions with the background thermal ions and electrons, which can be described by a Maxwellian distribution.In our work, a linear collision operator in the CoM space is adopted.The collision operator of the guiding center distribution function in the velocity space coordinate (v ∥ , v 2 ⊥ ) is formulated by the Eq. ( 15) of [18], where the first two terms of the Right Hand Side (RHS) are the slowing down terms in (v ∥ , v 2 ⊥ ) directions respectively, and the last three terms are the diffusion.The collision coefficients are Functions F , G and H are defined by respectively.Here x = v 2 /v 2 thβ and ϕ(x) is the Maxwellian integral defined by Here, m β and v thβ = 2k B T β /m β represent the mass and thermal velocity of the background Maxwellian-distributed particles.The basic collision frequency of EP colliding with background β particle here is defined by where ln Λ EP,β is the Coulomb logarithm, ϵ 0 is the vacuum permittivity.Due to the constraint that the particle's perpendicular energy must always be smaller than the total energy if utilizing coordinates such as E and µ, the CoM space tends to have a significant portion where f vanishes, namely f (µ > E/B) = 0.This issue necessitates specific treatment, such as employing triangular meshes, to appropriately handle the loss cone boundaries [22].In order to simplify the numerical simulation process, an alternative coordinate, Λ, is chosen instead of µ.Here, Λ is defined as Λ = µB 0 /E, which is a dimensionless quantity that provides a more convenient way to describe the equations.The three CoM are [39,40] with ψ is the poloidal flux coordinate, g is an equilibrium function.In HAGIS code, Boozer coordinates are used.The covariant form of the magnetic field is B = g(ψ)∇ζ + I(ψ)∇θ + δ(ψ, θ)∇ψ, where g is the covariant field component in ζ direction and g is a function of ψ [39,40].Then we can derive where σ = sign(v ∥ ).Substitute Eq. 7 to Eq. 2 we can derive the collision operator in CoM space.The collision operator is formulated The collision coefficients used in Eq. 8 are averaged over unperturbed orbits.The leading terms are kept in the equations and the off-diagonal diffusion terms are omitted.Additionally, the derivation in the Appendix substantiates the absence of off-diagonal diffusion terms.The diffusion coefficient in Λ direction D ΛΛ is zero when Λ = 0. To include the finite orbit width effect, an average over the unperturbed orbit is adopted as follows, where τ b is the bounce time, during which particles complete their poloidal orbits once, i.e. τ b = 1 θ dθ with the range of integration for a trapped particle being (−θ b , θ b ) and for a passing particle being (0, 2π), θ b is the bounce angle.For simplicity, the average ⟨⟩ is omitted.The numerical calculation of collision coefficients in the CoM space is detailed in Section 4.

Code structures
The transport equation Eq. 1 is solved in the ATEP-3D.The spatial differential is discretized using the finite volume method (FVM), a widely used technique in fluid simulation due to its favorable conservation properties.Specifically, a fully implicit Crank-Nicholson scheme is employed in order to ensure the accurate representation of the long-term behavior of EP transport.In the boundary, the finite difference boundary condition is used with different choices.The ATEP-3D numerical solver has been implemented using MATLAB, and a Fortran version is currently under development.To facilitate efficient testing of conservation properties and steady state studies, 1D and 2D solvers are utilized.The ATEP-3D code follows an object-oriented programming paradigm, consisting of two classes: atep 3d which serves as the core solver, and dist exp 3d cls an Experimental Data Interface.

Finite volume method for spatial discretization
The Finite Volume Method operates by subdividing the simulation domain into discrete cells and solving the equation at the centroid of each cell.The simulation domain is divided into N P × N E × N Λ cells with center values denoted using integer indices (i, j, k).The center values within the cells are denoted using integer indices (i, j, k), while the points located on the cell faces are indicated by indices with half integers.Since the advection terms in the equations can be negative and positive, the centered discretization scheme is used.The phase-space derivatives of the partial differential equation is evaluated at the cell face, ∂f ∂x m = , where x m is the mth coordinate and i m is the index in x m direction of this cell.Since the advection terms in the equations can be negative and positive, the centered discretization scheme is used.The phase-space derivatives are discretized using a linear centered scheme, i.e., the value of a quantity X at the cell face is taken as the average of the adjacent cells, The transport coefficients are dependent on the phase-space coordinates but remain constant over time.

Implicit scheme for time evolution
The time derivatives are discretized using the Crank-Nicholson method, l represents the time step index.Currently, a uniform grid and a fixed time step size are employed.
For simplicity, the source is assumed S l ≈ (S l + S l+1 )/2.An explicit source term is utilized.The transport coefficients D x , D xx are time-invariant.The derivatives appearing in the diffusion terms are also conveniently discretized using a central difference approximation.Move all terms at time step l + 1 in to the left hand side and those at step l to the right hand side, we obtain a system of linear equations.
where F and S are arrays and Fi+(j−1)N P +(k−1)N P N E = f ijk .The degree of freedom of this system is N dof = N P × N E × N Λ .The matrix representation of a single time step in the process involves the differential operator matrix M * , which is a square matrix with size N dof × N dof .It should be noted that the matrices M * are constructed based on selected boundary conditions and are independent of time.

Numerical Verification with Varied Boundary Conditions
Manufactured solutions are employed to assess the solver's performance, such as the pure advection and diffusion equations.Throughout our investigation, we investigate the mass conservation and identify the crucial role of boundary conditions in ensuring conservation.This section presents the results and discussions regarding the advection in the transport equation.The advection terms, such as the slowing down of energy due to collisions, induce a shift in the distribution function f while preserving its shape.A positive advective coefficient results in a negative flow.The EPs experience deceleration due to interaction with the background.In the transport model, EPs at the low energy tail is considered to merge to background ions, leading to EP "losses" at the boundaries.To address this, an "open" boundary condition is implemented, allowing the flux to pass through the boundary without distorting the distribution function or affecting the interior solution.The open boundary condition is represented by the equation: ∂ 2 f /∂x 2 = 0.In this section, a 1D test case is utilized to elucidate the numerical calculation of advection and the significance of open boundary conditions, The advective coefficient is assumed to be constant, and the initial distribution follows a Gaussian function.When D x = 1, D xx = 0 and the initial distribution f 0 = exp(−x 2 /0.025), the analytical solution of Eq. 13 is f = exp[−(x + D x t) 2 /0.025] with free boundaries.Different boundary conditions are implemented and tested.Since the simulation domain is not infinite, when the distribution function passes through the boundary the numerical error will arise rapidly due to the Dirichlet and Neumann boundary conditions.With using the open boundary condition, the conservation is satisfactory and the flux out of the system is properly calculated.As shown in Fig. 1(a), the left figure shows the result using zero boundary condition.When the time is longer than slowing down time, the numerical error increases at the lower boundary.The Neumann boundary condition is better but the conservation of total particle number is not ensured as shown in Fig. 1(b).With the open boundary condition, the numerical error is negligible and the particle number is conserved to machine precision as shown in Fig. 1(c2).The Courant number C = Dx∆t ∆x = 1×0.020.01 = 2 in these cases, and the numerical results have converged.

Orbit averaged collision in the CoM space
The collision coefficients are numerically calculated by the HAGIS code.The neoclassical transport of particles with wide orbits has been studied and validated [13].
To obtain collision coefficients in the CoM space, the implementation involves following test particles in the equilibrium and then averaging the coefficients over one poloidal orbit.Calculation of collision coefficient is added in the IMAS [41] versions of HAGIS and related wrapper module (FINDER).And the averaged diffusion coefficients are stored in distribution IDS.

ITER pre-fusion discharge
In this work, we used the ITER pre-fusion power operation (PFPO) discharge (shot number 100015, run number 1) [42].The PFPO phases are featured by low magnetic field, plasma current and density (1.8 T and 5 MA for a density between 40% and 50% of the Greenwald density).In the following, we adopt the equilibrium and other variables from the ITER PFPO discharge for the calculation of the EP collisions and the transport.The collision coefficients vary along particle orbits and the effects of the finite obit width and variation of the background profiles are important.In this work, uniform background density and temperature profiles of thermal ion and electron are assumed.The averaged collision coefficients in the CoM space for this case are calculated.In Fig. 2, the collision coefficient D P of co-going particles with E = 298 keV is shown for example.For the co-passing particles, the D P is positive, it is due to the velocity drag.The D P of passing particle is larger since it is related to the parallel velocity.For trapped particles, D P is small.This illustrates that the collision coefficients depend on the particles orbit types, which can be described in the CoM space.All collision coefficients in the CoM space are constructed in this way.

Collisional EP transport
In the following, we show the results of the NBI slowing down process simulated by ATEP-3D.The simulation regime of energy is 0.1 M eV to 1 M eV .In this work, we are more concerned with the slowing down process of EPs and thus we use 0.1 M eV as the lower boundary along E. In the future we will use a similar rule as TRANSP [43], that treat all particles below 1.5v thermal as the main thermal species.In addition, at the even lower energy range, as the EPs merge to the background ions, the evolution of the background ion also needs to be modeled, which merits more effort in our future work.In simulations, the energy is normalized to 1 M eV .1D integrated normalized collision coefficients are shown in Fig. 4, giving the structure and magnitude of the advection and diffusion coefficients in all directions.In order to model the process of the EP evolution with source, the NBI source is modeled by a Gaussian function, where A s = 2×10 18 , δP ζ = 0.2, δE = 0.05, δΛ = 0.4.The initial distribution is f 0 = 0. Steady-state of EP distribution due to collisions and source is reached at about 1s.
Clearly, the main processes are energy slowing down and pitch angle scattering.We compare the NBI case with NEMO/SPOT.NEMO is to simulate the NBI deposition.SPOT (simulation of particle orbits in a tokamak) is an orbit following Monte Carlo code with Fokker-Planck collision operator [21].The birth energy is 745 KeV for both ATEP-3D and NEMO/SPOT, without the consideration of the half or third energy of EPs since ITER will employ a neutral beam injector with negative ion source.The SPOT result is shown by the black line in Fig. 5.The steady-state solution of ATEP-3D agrees reasonably well with SPOT result.The discrepancies may come from that SPOT using the realistic density and temperature profiles.Additionally, the physics model and simplifications in NEMO/SPOT are not identical to those used in ATEP-3D, such as the parametrization of source and sink terms.Further detailed comparisons with other transport codes will be carried out in the future.The backmapping of EP distribution in CoM space to the real space will be performed as described in [26].

Conclusions
This study focuses on the transport of energetic particles in tokamak plasmas.Due to the relatively large orbit width resulting from the high energy of the energetic particles, the collision operator is expressed in constants of motion space and averaged along the unperturbed particle trajectory, which yields the diffusion and drag coefficients in  constants of motion space.The numerical tool ATEP-3D has been developed, using the finite volume method and implicit Crank-Nicolson scheme.Favorable numerical behaviors are demonstrated in terms of conservation properties, large time step size and long time scale numerical behavior.The open boundary condition is shown to be a key ingredient for the modelling of the transport problem in the presence of source and sink as well as the collision terms.As the numerical verification, the numerical results in one dimensional case show good agreement with the analytical results.The ITER prefusion discharge serves as a test case, demonstrating the ATEP-3D code's capability in three-dimensional scenarios.The collisional coefficients in the constants of motion space are calculated and the slowing-down process of the EP transport is simulated.Adding the transport process due to Alfvénic perturbations [26], the ATEP-3D code will be applied to the studies of the evolution of the EP distributions with the phase space zonal structure taken into account in a consistent way.

Figure 2 .
Figure 2. Collisional transport coefficient D P of E = 298 KeV particles in the constants of motion space.
The normalization of P ζ and Λ is chosen so that the normalized values are from 0 to 1.The time unit is second.The collision coefficients are normalized accordingly.The three-dimensional structure of the normalized collision coefficients are shown in Fig. 3.The advection along energy direction (D E ) is the dominant one indicated by its magnitude.The advection along P ζ is contributed by both co-passing and counter-passing particles, corresponding different signs in D P .The diffusion terms are dominated by D ΛΛ .The

Figure 3 .
Figure 3. Collision coefficients in the three dimensional constants of motion space.The collision coefficients are normalized according to the normalization of energy, P ζ and time.

Figure 5 .
Figure 5.The figure shows the time evolution of the EP distribution simulated by ATEP-3D, represented by the color lines with circle markers.The black dashed line with triangle markers represents the SPOT result.