A new quasilinear formulation for ICRF plasmas in a toroidal geometry

We present a new formulation for quasilinear velocity space diffusion for ICRF plasmas that considers two different aspects: (1) finite Larmor radius approximation and (2) includes the effect of toroidal geometry and constructs a positive definite form. In the first aspect, the Kennel-Engelmann (K-E) quasilinear diffusion coefficients are successfully approximated in a small Larmor radius limit and implemented for the numerical codes (TORIC-CQL3D). In the second aspect, the quasilinear diffusion is reformulated in a toroidal geometry in order to include the parallel dynamics in the inhomogeneous plasmas and magnetic fields. We use these two quasilinear formulations to simulate ITER plasmas with ICRF injection for minority fundamental heating and Tritium second harmonic cyclotron heating.


Introduction
The ion cyclotron range of frequency (ICRF) waves transfer their energy and momentum to plasmas, and as a result a significant amount of fast ions are produced.The kinetic description of the non-Maxwellian distribution plasmas is obtained using quasilinear diffusion coefficients due to the RF waves in a Fokker-Planck equation.The quasilinear diffusion coefficients are derived by Kennel and Engelmann (K-E) [1] with the assumptions of the homogenous plasmas and magnetic fields along the particle trajectory.In this proceeding, we present the theoretical modifications of the K-E diffusion coefficients in two different aspects.
In the first aspect, the coefficients are modified to be consistent with the dielectric tensor of the reduced model in the small Larmor radius approximation.In the reduced model, the dielectric tensor for the plasma current is expanded by a small parameter  ⊥   , in which  ⊥ can be replaced by the differential operator [2].Accordingly, the computation time of the reduced model can be significantly reduced by a factor of O(  2 ) compared to the full model [3] that uses  ⊥ explicitly without the small Larmor radius approximation.Here,   is the number of radial coordinate grid.In Section 2, we summarize the derivation of the equivalent quasilinear diffusion to the reduced model as in our recent publication [4], and present an additional example using the implementation in TORIC-CQL3D.
TORIC is already coupled with a sophisticated Fokker-Planck module, SSFPQL [5,6].It is worth noting the different features between our approach in TORIC-CQL3D [4] and that in TORIC-SSFPQL [5,6].CQL3D is the Fokker-Planck code, which has been validated against demonstrated by many RF wave experiments.It has fully nonlinear and relativistic collision operators and it can calculate the time evolution of distribution function for the tail formation and relaxation by transient RF waves.For the quasilinear diffusion, our formulation in TORIC-CQL3D guarantees the consistency with the dielectric tensor, while it cannot include the contribution from the higher order in  ⊥   that is retained in the K-E coefficients used in TORIC-SSFPQL.As the result of the self-consistency, the power absorptions by dielectric tensor and quasilinear diffusion theoretically match each other and the numerical iterations between the Maxwell's equation solver and the Fokker-Planck code likely converge without any normalization factor to the diffusion coefficients [4].
In the second aspect, we modify the formulation of K-E coefficients to be positive-definite in the toroidal geometry.The positive definite form can include some important features of the diffusion in the toroidal geometry which cannot be captured in the K-E coefficients, as will be explained in Section 3. Using this new positive-definite from, we can eliminate the unphysical growing mode violating H-theorem and the related numerical errors in the coupled code TORIC-CQL3D.To best of our knowledge, it is the first implementation of the positive-definite form for the quasilinear diffusion coefficient in a continuum Fokker-Planck solver for a toroidal geometry.

Kinetic energy change (𝑾𝑾
̇) The rate of the kinetic energy change due to the RF waves ( ̇) is used to define the quasiliner diffusion coefficients in a specific numerical formulation.For example, in the full model using a full Fourier spectral  EPJ Web of Conferences 157, 03028 (2017) space, the K-E quasilinear diffusion coefficients can be defined by where W  is the resonance Kernel (e.g.Eq. (11) of [3]).
For energy transfer, the quasilinear diffusion coefficients are gyroaveraged, and they can be described in the velocity space of speed-pitchangle coordinate ( ) by (2) The relation between the coefficients , , , and  are given to preserve the diffusion directions in velocity space [6].
In the reduced model, the kinetic energy can be defined to be consistent with the approximation of plasma currents using the expansion of the gyromotion around the gyrocenter.For example, when evaluating  ̇, the electric field at the gyrocenter is approximated by and it results in the expansion of  ̇ in terms of the small parameters where the superscript indicates the order of  ⊥   .The detailed derivations for Eq. ( 4) are in Section 2 of [4].

Quasilinear diffusion coefficient
The K-E quasilinear diffusion coefficient corresponding to  ̇ in the full model of Eq. ( 1) is where  is the integer to determine the harmonic number of ion cyclotron resonance and   is the effective potential [3].
To derive the quasilinear diffusion coefficients corresponding to  ̇ in the reduced model of Eq. ( 4), three important properties of the diffusion need to be considered : positive-definiteness, wave polarization, and, diffusion direction.We found a way to formulate the coefficients to preserve these three properties [4].The coefficient  is formulated from the lowest order of  ̇, and other coefficients are calculated from the original relations with .

Examples
The quasilinear diffusion coeffcients for the reduced model are evaluated in the wave code by the reduced model, TORIC [2].Using the diffusion coefficients, TORIC is coupled with a Fokker-Planck code CQL3D [6], to have the self-consistent solutions of distribution and wave fields.In the recent paper [4], we simulate the scenario for ITER 10MW He3 minority species damping (n=1) by the reduced model TORIC-CQL3D, and compared the results with the full model by AORSA-CQL3D [3].The comparisons show that the reduced model is acceptably valid when the considerable amount of fast ions satisfies  ⊥   ≾ 1.0 and the wave power density is less than 1MW/m 3 in the scenario.
In this proceeding, we introduce the results of the second harmonic damping (n=2) scenarios of ITER 5.3T.We simulate two ion species with the ratio of (D, T)= (50,50)%, and the ICRF wave frequency is 52.5MHz.The dominant wave power is absorbed by the second harmonic damping of the major species Tritium around the magnetic axis because ω = 2Ω at R = R 0 +0.17m.The power decompositions are reasonably similar between TORIC and AORSA: for Maxwellian plasmas, 49% of T second harmonic damping and 48% electron damping in TORIC, while 48% of T second harmonic damping, and 51% electron damping in AORSA.For self-consistent non-Maxwellian plasmas, 48% of T second harmonic damping, and 46% electron damping in TORIC, while 61% of T second harmonic damping, and 39% electron damping in AORSA.
Figure 1 and 2 show the distribution functions of TORIC-CQL3D and AORSA-CQL3D, which show the similar patterns but some difference in the tails around  ∥ ⋍ 0. The main reason of difference could be the error of the approximation on the Bessel function factor  1 ( ⊥   ) ⋍  ⊥   /2 for the high energetic Tritium of several MeV.The error of the small Larmor radius approximation can be larger for n=2 than that of n=1 in [4].

Positive definite form
In a toroidal geometry, the magnetic field is inhomogeneous along the particle trajectory and the assumption of K-E coefficients is violated.The violation is typically ignorable because of the small correlation length of the wave-plasma interaction compared to the variation of the magnetic fields.However, we found that it does not result in the positive-definiteness when the quasilinear diffusion coefficients are bounce-averaged for the bounce-averaged Fokker-Planck code, CQL3D [6].As a consequence, the negative values in the bounce-averaged diffusion coefficient 〈〉  = ∫   ∥ cause unphysical growing modes and numerical errors when enforcing to eliminate them.
The positive-definiteness of the bounce-averaged coefficients are recovered when treating the trajectory integral and the bounce integral as the same way in the toroidal geometry [7].For example, the bounce-averaged coefficient is reformulated by where 〈… 〉  is the gyroaverage and the symmetric operator (  ,   ) =  * (  ,   ) is obtained by the change of variables  1 =    and  2 = , giving ( 2 , ∥,2 )−  ( 1 , ∥,1 ) , ( 9) and it guarantees the positive-definiteness in the toroidal geometry.In Eq. ( 5), K-E coefficients have the Diracdelta function depending on only  ∥,1 because of the uniform phase only in the trajectory integral, and it breaks the symmetry shown in Eq. ( 9).

Trajectory integral
The evaluation of the trajectory integral in Eq. ( 9) is not trivial because of the phase variation in the toroidal geometry, The stationary approximation around the resonance point   / = 0 is useful for the short correlation length.By expanding   in terms of the small temporal distance from the resonance     up to third order, it results in the Airy functions (), as many previous studies investigated [8,9], where  = 0 3   /  3 | ,  + =  2   /  2 / 3   /  3 , and  1 =  + 2 at  =   , result in the negative argument of the Airy function.
When  2   / 2 = 0 is satisfied, the correlation length becomes large so that the variation of the phase could be important even for the non-resonance locations with   / ≠ 0. This condition typically occurs at the inner-midplane and outer-midplane for passing particles and trapping tips and outer-midplane for trapped particles [9].We can include this non-resonant contribution to the diffusion by the positive arguments of the Airy function, ∫  ⊥ ()   ( ∥ ) where  2 = 2  / / 3   /  3  at the specific conditions  =   .The resonant contribution in Eq. ( 11) and the nonresonant contribution in Eq. ( 13) are smoothly connected by the Airy function, and it results in the continuous diffusion coefficients in velocity space.

Implementation
This positive definite form of the bounce-averaged coefficients in Eq. ( 8) is implemented in TORIC.The resonant contribution is evaluated when the spectral mode has at least one resonance location on a flux surface, and the non-resonant contribution is included when the spectral mode has no resonant location.
The phase integral   is not poloidally periodic on a general flux surface with the parallel spectral mode number (m+nq) not being an integer.Thus, the evaluation of the diffusion coefficient depends on the initial position of the integral range, and the average diffusion changes depending on the number of poloidal periods in the evaluations.We will show the importance of the initial position and evaluation range in a future publication.
Additionally, the positive definite form can show the correlation between the consecutive resonance crossings, while the kicks by the Dirac-delta function of the K-E diffusion coefficients are likely randomized.This correlation between resonances are important in the two consecutive resonance around the location of  =   (e.g.outer-midplane), where the toroidal geometry effect needs to be considered within the relatively long correlation length.
This new quasilinear diffusion is expected to be beneficial for ICRF modelling in a low aspect ratio tokamak that has the significantly varying magnetic field along particle trajectories.In spite of the good features to capture the toroidal effects in the quasilinear diffusion, this positive definite form is computationally expensive because it requires about         ~10 6 times more number of floating-point operations than evaluating K-E quasilinear diffusion in Eq. ( 5).Here,   ~100,   ~100   ~10, and   ~10 are the number of velocity space grid in each dimension, the number of polodal spectral grid, the order of interpolation, the required number of interpolation, respectively.This work was supported by US DoE Contract No. DE-FC02-01ER54648 under a Scientific Discovery through Advanced Computing Initiative.

Fig. 1 .Fig. 2 .
Fig. 1.Reduced model results by TORIC-CQL3D : (a) 2-D contour plots of the distribution function in velocity space at r/a=0.1 (b) 1-D distribution functions in terms of speed for several pitch-angles.Here, u norm is the momentum corresponding to the energy of Tritium 4 MeV.(a) (  )   (   ∥ ) (  )