P Donnel

C Gheller

S Brunner

L Villard

E Lanti

N Ohana

M Murugappan

Moment approach of the multi-species nonlinear Coulomb collision operator adapted to particle-in-cell codes

In this article an approximated version of the multi-species, nonlinear Coulomb collision operator is derived via the use of a truncated moment expansion of the distribution function to compute the Rosenbluth potentials. The evolution of the distribution function due to the collision operator takes the form of a Fokker-Planck equation that can be efficiently solved by a Langevin approach in the particle-in-cell (PIC) framework. This kinetic collision operator is then implemented in the global gyrokinetic PIC code ORB5 while neglecting finite Larmor effects. This collision operator ensures exact conservation of density, total momentum and energy as expected. The neoclassical physics is also successfully retrieved by the code. Numerical aspects such as the noise control and the parallelisation of the code are presented in detail.

Introduction

In tokamaks, most of the transport is often attributed to turbulence. The most complete and appropriate approach to compute this turbulence is based on gyrokinetic theory. Even in weakly collisional plasmas, the inclusion of collisional proccesses is important for various reasons. Firstly, the evolution of the system produces fine scale structures in velocityspace. In a collisionless plasma, these structures are becoming finer in time until they reach the numerical limits of the code. Physically, these structures are smeared out by collisions. Secondly, the collisional physics referred to as neoclassical can be important for some key physics features. This is for instance the case for the bootstrap current, which is a key ingredient for steady state operation. Neoclassical transport can also be important in transport barriers such as the pedestal of H-mode plasmas where the turbulence level is reduced [START_REF] Battaglia | Kinetic neoclassical transport in the h-mode pedestal[END_REF]. Finally, turbulent and neoclassical physics can display important synergies. On the one hand, collisionnal process can modify turbulence either directly by collisionnal detrapping in the case of trapped electron modes or indirectly via the collisional damping of large scale flows [START_REF] Dif-Pradalier | Interplay between gyrokinetic turbulence, flows, and collisions: Perspectives on transport and poloidal rotation[END_REF]. On the other hand, turbulence can modify neoclassical transport of impurities [START_REF] Estève | Self-consistent gyrokinetic modeling of neoclassical and turbulent impurity transport[END_REF][START_REF] Donnel | Neoclassical impurity flux in presence of turbulent generated poloidal asymmetries and pressure anisotropy[END_REF] via the generation of relatively large poloidal asymmetries [START_REF] Donnel | Turbulent generation of poloidal asymmetries of the electric potential in a tokamak[END_REF].

Collision operators are derived in the kinetic theory and act on particles. In principle, these operators should be adapted to the gyrokinetic theory and act on gyrocenters. But this operation leads to a significant increase of the difficulty of the problem as the collision operator acts only in the velocity space in the kinetic framework whereas it is a 5D operator in the gyrokinetic framework. Moreover, conservation properties are expressed for particles. Their translation for gyrocenters is difficult. Fortunately, in the presence of long wavelength perturbations k ⊥ ρ < 1, Finite Larmor Radius (FLR) effects can be neglected [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]. This is called the drift-kinetic approximation. Neoclassical theory has been derived with this approximation. We will implicitly use the drift-kinetic approximation in all this paper. This implies that no distinction is done between the distribution functions of particles and gyrocenters. One should nevertheless keep in mind that this assumption could break down in certain regimes. This is for instance the case for energetic particles because of their large gyroradius. The development of gyrokinetic collision operator (i.e. including FLR effects) have made recent progress in the local gyrokinetic approach. Indeed, a nonlinear gyrokinetic collision operator has recently been derived by relying on a Fourier approach to compute FLR corrections [START_REF] Jorge | Nonlinear gyrokinetic coulomb collision operator[END_REF]. FLR effects have also been recently studied on a local version of the GENE code, by implementing the exact gyrokinetic linearized collision operator [START_REF] Pan | First implementation of gyrokinetic exact linearized landau collision operator and comparison with models[END_REF]. This collision operator is numerically expensive compared to standard collision operators and leads to modifications of results mostly in the regime k ⊥ ρ > 1 as expected. To the authors knowledge, collision operators with FLR corrections have never been implemented in global gyrokinetic codes. One of the reasons is that the Fourier approach used to compute the gyroaverage quantities in local codes is not possible in a global code.

Another difficulty encountered by global codes is that collisionality in tokamak plasmas is varying significantly across the plasma radius. Indeed, the core of the plasma is often in the low collisionality regime, while the edge of the plasma is much more collisional due to the lower temperatures. Furthermore, edge fluctuations are significant, making the linearization of the collision operator questionable. While an approximated collision operator can be used for the core, a complete collision operator is required to have a validity domain wide enough to cover all the collisionality regimes found in global simulations spanning from the core to the edge of tokamak plasmas.

The rest of this paper is organised as follow. Firstly, a derivation of the approximated nonlinear collision operator is presented. A numerical method to ensure exact conservation properties is then presented. In section 4, a linearization of this collision operator is derived to evaluate the effects of nonlinearities. In section 5 a noise control adapted to the presence of collisions is presented. The section 6 presents a serie of tests performed to verify the collision operator. The section 7 is dedicated to the high performance implementation of the collision operator which is needed for global gyrokinetic codes.

Derivation of the collision operator

Expression of the collision operator with Rosenbluth potentials

The Coulomb collision operator corresponding to a species a colliding on a species b can be written in the form of an advection diffusion operator

C ab (F a , F b ) = γ ab m a ∂ v • ∂ v • 1 2 ∂ v ∂ v G b F a -1 + m a m b ∂ v H b F a (1)
with γ ab =

(v) = dv F b (v ) |v -v | (2)
H b (v) = dv F b (v ) |v -v | (3)
The collisional evolution of the distribution function of the species a is simply given by the sum over the species of the collision operator

∂F a ∂t = b C ab (F a , F b )
This equation can be recast as a Fokker-Planck Equation

∂F a ∂t + ∂ v • [Γ a F a -∂ v • (D a F a )] = 0 (4)
with a drag term

Γ a = b γ ab m a 1 + m a m b ∂ v H b (5)
and a diffusion term

D a = b γ ab 2m a ∂ v ∂ v G b (6)
At this point, no approximation has been performed. In the PIC framework, this Fokker-Planck equation can be easily solved using a Langevin approach. The difficulty is hidden in the computation of the Rosenbluth potentials and their derivatives, which is detailed in the next section.

Computation of the Rosenbluth potentials

An exact expression of the Rosenbluth potentials based on an expansion of the distribution function on moments has been derived by Ji et al. [START_REF] Ji | Exact linearized coulomb collision operator in the moment expansion[END_REF]. The main steps of the derivation are described in this section. The idea is to expand the distribution function around an unshifted Maxwellian, which is in the kernel of the collision operator, using two sets of orthogonal polynomials. Irreducible harmonic polynomials P l are used for the direction dependence. Sonine (associated Laguerre) polynomials L (l+1/2) k are used for the energy dependence.

F b = ∞ l=0 ∞ k=0 F lk b σ l k (7) with σ l k = l! (l + k + 1/2)! 2 l k! (l + 1/2)! and F lk b = F M b L (l+1/2) k s 2 b P l (s b ) • M lk b In these expressions, F M b = N b π 3/2 v 3 T b e -
M lk b = 1 N b v l T b dvP l (v) L (l+1/2) k s 2 b F b (8)
Injecting the moment expansion of the distribution function Eq.( 7) in the Rosenbluth potentials Eq.( 2), (3), one gets

G b (v) = lk M lk b σ l k • dv F M b P l (s b ) L (l+1/2) k s 2 b |v -v | H b (v) = lk M lk b σ l k • dv F M b P l (s b ) L (l+1/2) k s 2 b 1 |v -v |
The Sonine polynomials can be written

L (l+1/2) k (x) = k m=0 c l km x m with c l km = (-1) m (l + k + 1/2)! (k -m)! (l + m + 1/2)!m!
Therefore the Rosenbluth potentials take the form

G b (v) = lk M lk b σ l k • k m=0 c l km G lm b H b (v) = lk M lk b σ l k • k m=0 c l km H lm b with G lm b = dv F M b P l (s b ) s 2m b |v -v | H lm b = dv F M b P l (s b ) s 2m b 1 |v -v |
One major step is the computation of the coefficients G lm b , H lm b that is detailed in [START_REF] Ji | Exact linearized coulomb collision operator in the moment expansion[END_REF]. The results are

G lm b = 2N b v T b P l (v) 2l + 1 1 2l + 3 I 2(m+l+2) b+ s l+1 b + s l+2 b I 2m+1 b- - 1 2l -1 I 2(m+l+1) b+ s l-1 b + s l b I 2m+3 b- (9)
H lm b = 2N b v T b P l (v) 2l + 1 I 2(m+l+1) b+ s l+1 b + s l b I 2m+1 b- (10)
with

I 2k b+ = (k -1/2)! 2   E b (1/2)! -E b k-1 j=0 s 2j+1 b (j + 1/2)!   (11)
and

I 2k+1 b- = k! E b 2 k j=0 s 2j b j! (12)
In these expressions, E b is the error function evaluated in s b ,

E b = 2 √ π s b 0 dte -t 2 , and E b its derivative, E b = 2 √ π e -s 2 b .

Truncation of the polynomial expansion

At this point no approximation was performed on top of the drift-kinetic approximation. In other words, the collision operator is still the full nonlinear Coulomb operator in the kinetic framework, and its drift-kinetic limit in the gyrokinetic framework.

The evaluation of high order moments of the distribution function is sensitive to the numerical resolution. Furthermore, it can be shown that higher order moments are more effectively damped by collisions. It is then natural to truncate the expansion presented in the previous section. From now on only the components with l, k ∈ [0, 1] are kept. This truncation is somehow arbitrary, but it has been shown that this truncation is accurate for neoclassical physics (look table I of [START_REF] Hirshman | Neoclassical transport of impurities in tokamak plasmas[END_REF]). In principle, it would be relatively simple to extend the following derivation to include more moments. It would of course lead to a more accurate collision operator. But on top of the extra analytical work, one has to be careful to compute numerically higher order moments with a sufficient precision. To keep only the components with l, k ∈ [0, 1] is therefore a reasonable compromise between the physical accuracy and the numerical robustness and cost of the collision operator.

The truncation procedure is performed to compute the Rosenbluth potentials, and therefore modifies the collision operator. An index t is used for truncated quantities.

C t ab = -∂ v • Γ t ab F a -∂ v • ∆ t ab F a with an approximated drag term Γ t ab = γ ab m a 1 + m a m b ∂ v H t b
and an approximated diffusion term

∆ t ab = γ ab 2m a ∂ v ∂ v G t b
Note that no approximation is performed on the distribution function of the species a. The truncation removes the momentum and energy conservation of the collision operator. The momentum and energy conservation are also not fulfilled due to the finite numerical resolution of the problem. A correction term C cor ab is therefore added to reintroduce the conservation properties up to machine precision. It is detailed in the section 3.2. The model operator proposed in this article therefore takes the form

C mod ab = C t ab + C cor ab (13)
A comparison between the approximated version of the nonlinear collision operator proposed in this paper and some linear collision operator present in the literature is insightful at this stage. The full nonlinear collision operator reads

C ab (F a , F b ) = C ab (F M a , F M b ) + C ab (δF a , F M b ) + C ab (F M a , δF b ) + C ab (δF a , δF b )
In principle, the linear collision operator is simply obtained by dropping the last term. But the field part of the operator C ab (F M a , δF b ) is still difficult to handle and is in practice also approximated. Different approximations of the field operator have been proposed in the literature. Hirshman and Sigmar [START_REF] Hirshman | Approximate fokker-planck collision operator for transport theory applications[END_REF] proposed a similar truncation procedure to compute the approximated field operator. A naive truncation as proposed here leads to an easy computation of the approximated field operator but the operator loose some of its conservation properties. More elaborate variants of the truncature approach are then introduced to retrieved conservation properties. Lin [START_REF] Lin | Gyrokinetic particle simulation of neoclassical transport[END_REF] proposed to use the properties of the collision operator to construct a model field operator. The model collision operator proposed in this paper can be seen as an hybrid version of these two methods in the context of nonlinear collision. The truncation method is used to obtained a collision operator C t ab as close as possible to the original one and is then corrected C cor ab to reintroduce the required properties of the operator. To illustrate this point, it is shown in the section 4 that if only the components l = k = 0 are kept one retrieve exactly the Lin operator.

We now come back to the derivation of the truncated collision operator C t ab . Using the fact that c 0 00 = 1,

c 1 00 = 1, c 0 10 = 3 2 , c 1 10 = 5 2 , c 0 11 = -1, c 1 11 = -1 and σ 0 0 = 1, σ 1 0 = 1 2 , σ 0 1 = 3 2 , σ 1 1 = 5 4 , and
M 00 b = 1 M 10 b = V b v T b M 01 b = - V 2 b v 2 T b M 11 b = 2 m b N b v 3 T b 5 2 N b T b V b -Q b = - 2q b m b N b v 3 T b
where V b is the average velocity, Q b is the kinetic energy flux and q b the heat flux. It can be shown that the truncated Rosenbluth potentials are given by

G t b = G 00 b - V 2 b v 2 T b G 00 b - 2 3 G 01 b + 2 V b v T b • G 10 b - 4q b m b N b v 3 T b • G 10 b - 2 5 G 11 b ( 14
)
and

H t b = H 00 b - V 2 b v 2 T b H 00 b - 2 3 H 01 b + 2 V b v T b • H 10 b - 4q b m b N b v 3 T b • H 10 b - 2 5 H 11 b ( 15
)
The values of H lm b and G lm b with l, m ∈ [0, 1] are given in appendix A. After a bit of algebra, it is possible to compute the components of the Rosenbluth potentials

G t b N b v T b = 1 2s b + s b E b + E b 2 + 1 3 V 2 b v 2 T b E b s b + V b • v v 2 T b 1 2s 3 b - 1 s b E b - E b 2s 2 b + 2 5 q b • v m b N b v 4 T b E b s 3 b - E b s 2 b (16) v T b N b H t b = E b s b - V 2 b v 2 T b E b 3 + V b • v v 2 T b E b s 3 b - E b s 2 b - 4 5
q b • v m b N b v 4 T b E b (17)

Expressions of the drag and diffusion

At this point, the approximated Rosenbluth potentials have been computed Eq.( 16)- [START_REF] Brunner | Collisional delta-f scheme with evolving background for transport time scale simulations[END_REF]. The next step consists in expressing the drag force Eq.( 5) and the diffusion matrix Eq.( 6) associated. To do so, we use the fact that for any function of the norm of the velocity λ (v) and any constant vector A, one has

∂ v λ = 1 v dλ dv v ∂ v ∂ v λ = 1 v dλ dv I - vv vv + d 2 λ dv 2 vv vv ∂ v [λ (A • v)] = 1 v dλ dv (A • v) v + λA ∂ v ∂ v [λ (A • v)] = 1 v dλ dv I - vv vv + d 2 λ dv 2 vv vv (A • v) + 1 v dλ dv (Av + vA)
Using these relationships, it is possible to compute the derivatives of the Rosenbluth potentials

∂ v H t b = N b v 3 T b v E b s 2 b - E b s 3 b + 2 3 V 2 b v 2 T b E b + V b • v v 2 T b -3 E b s 5 b + 3 + 2s 2 b E b s 4 b + 8 5
q b • v m b N b v 4 T b E b + E b s 3 b - E b s 2 b N b V b v 3 T b - 4 5
q b m b v 5 T b E b ∂ v ∂ v G t b =    1 -1 2s 2 b E b s b + E b 2s 2 b + 1 3 V 2 b v 2 T b E b s 2 b -E b s 3 b + V b •v v 2 T b 1 s 3 b -3 2s 5 b E b + 3 2 E b s 4 b + 2 5 q b •v m b N b v 4 T b -3 E b s 5 b + 2 s 2 b + 3 s 4 b E b    N b v T b I - vv vv +    E b s 3 b - E b s 2 b + 2 3 V 2 b v 2 T b E b s 3 b -1 + 1 s 2 b E b + 2 V b •v v 2 T b 3 s 5 b -1 s 3 b E b -1 s 2 b + 3 s 4 b E b + 8 5 q b •v m b N b v 4 T b 3 E b s 5 b -1 + 2 s 2 b + 3 s 4 b E b    N b v T b vv vv + 1 s 3 b - 3 2s 5 b E b + 3 2 E b s 4 b N b V b v + vV b v 3 T b + -3 E b s 5 b + 2 s 2 b + 3 s 4 b E b 2 5 q b v + vq b m b v 5 T b

Expression of the collision operator in the Langevin approach

Until now, the effect of collisions on the distribution function was represented by a Fokker-Planck equation which is a partial differential equation (PDE). A Fokker-Planck equation is equivalent to a Langevin equation which is a stochastic differential equation (SDE) [START_REF] Tabar | Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems[END_REF]. SDE are easier to handle in the PIC framework. Itô and Stratonovich algebras are two equivalent algebras allowing to link one Fokker-Planck equation to its Langevin counterpart. In this article, Itô algebra is used. Then the Fokker-Planck equation representing the collision of a species a with a species b with the approximated collision operator

∂F ∂t = -∂ v • Γ t ab F a -∂ v • ∆ t ab F a is therefore equivalent to dv dt = Γ t ab + G ab N (18)
where G ab • G T ab = 2∆ t ab . In this expression the exponent T refers to the matrix transpose. It is different from the exponent t which stands for truncated quantities. The noise N is composed of random functions whose statistical characteristics are well-defined. Eq.( 18) is discretised using the Euler-Maruyama scheme.

∆v ab = Γ t ab ∆t + G ab √ ∆tR b (19)
where

R b =   R b x R b y R b z   with R b
x , R b y and R b z three independent random numbers chosen with a Gaussian distribution with zero mean and a variance of 1. To compute the matrix G ab , the diagonalisation of ∆ t ab can be performed. For this purpose, a basis needs to be chosen. In presence of a strong magnetic field, the mean flow and heat flux tend to align with the magnetic field. Perpendicular components of V b , q b are therefore neglected. For the basis, we choose e z the direction parallel to the marker velocity before collisions, e y the orthogonal direction in the plane containing the initial marker velocity and the parallel direction and finally e x to get an orthogonal basis. In this basis, any vector aligned with the magnetic field takes the form A = Ae = A (cos θe z + sin θe y ). Using also the fact that v cos θ = v and v sin θ = -v ⊥ , one gets A = A v v e z -v ⊥ e y . This allows to compute the three matrices

vv vv =   0 0 0 0 0 0 0 0 1   I - vv vv =   1 0 0 0 1 0 0 0 0   Av + vA = A   0 0 0 0 0 -v ⊥ 0 -v ⊥ 2v  
Then one can compute explicitly the collision drag and diffusion on the marker considered. The drag term takes the form

Γ t ab v T a = νab v 2 T a v 2 T b 1 + m a m b
Γy ab e y + Γz ab e z (20)

Γy ab = - v ⊥ v T b E b -s b E b s 4 b V b v T b - 4 5 q b m b N b v 3 T b E b s b Γz ab = - E b -s b E b s 2 b + 2V 2 b 3v 2 T b s b E b + 2 V b v v 2 T b - E b s 4 b + 1 + s 2 b E b s 3 b + 4q b v 5m b N b v 4 T b 2s 2 b -1 E b s b where νab = N b v 3 T a γ ab
ma . The diffusion term takes the form

D t ab v 2 T a = νab 2 v T a v T b Dab (21)
E b s 3 b + E b 2s 2 b - 1 3 V 2 b v 2 T b E b s 3 b - E b s 2 b + V b v v 2 T b s 2 b - 3 2 E b s 5 b + 3 2 E b s 4 b + 2 5 q b v m b N b v 4 T b -3 E b s 5 b + 2s 2 b + 3 E b s 4 b Dzz ab = E b s 3 b - E b s 2 b + 2 3 V 2 b v 2 T b E b s 3 b -s 2 b + 1 E b s 2 b + V b v v 2 T b 3 E b s 5 b -2s 2 b + 3 E b s 4 b + 4 5 q b v m b N b v 4 T b 3 E b s 5 b -2s 4 b + 2s 2 b + 3 E b s 4 b Dyz ab = - V b v ⊥ v 2 T b 1 s 3 b - 3 2s 5 b E b + 3 2 E b s 4 b + 2 5 q b v ⊥ m b N b v 4 T b 3 E b s 5 b -2s 2 b + 3 E b s 4 b
The diffusion matrix Dab is real and symmetric and can therefore be diagonalised. . This is the case if the distribution function is close enough to a Maxwellian with a subsonic mean velocity. This is the assumption that was made to truncate the polynomial expansion of the distribution function. The positivity of Dab is necessary to avoid anti-diffusion which is not manageable by the Langevin approach which is used.

It can be shown that the three eigenvalues of Dab are Dyy ab and λ ± ab =

Dyy ab + Dzz ab ± ( Dyy ab -Dzz ab ) 2 +4( Dyz ab ) 2 2
. Defining the orthogonal matrix P ab

P ab =      1 0 0 0 Dyz ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 Dyy ab -λ + ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 0 λ + ab -Dyy ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 Dyz ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2     
it can be shown that P T ab Dab P ab = ∆ ab ∆ ab , with

∆ ab =      Dyy ab 0 0 0 λ + ab 0 0 0 λ - ab     
Therefore G ab takes the form

G ab = νab v T a v T b v T a P ab ∆ ab P T ab = νab v T a v T b v T a        Dyy ab 0 0 0 ( Dyz ab ) 2 √ λ + ab +( Dyy ab -λ + ab ) 2 √ λ - ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 Dyz ab (λ + ab -Dyy ab ) √ λ + ab - √ λ - ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 0 Dyz ab (λ + ab -Dyy ab ) √ λ + ab - √ λ - ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2 (λ + ab -Dyy ab ) 2 √ λ + ab +( Dyz ab ) 2 √ λ - ab ( Dyz ab ) 2 +(λ + ab -Dyy ab ) 2        (22)
The combination of ( 19), ( 20) and ( 22) corresponds to the truncated collision operator C t ab in the Langevin approach. Its implementation is therefore straightforward. One should note that in the case Dyz ab = 0 there is an under determination of the eigenvalues λ ± ab . One should choose λ + ab = Dzz ab and λ - ab = Dyy ab for consistency.

3 Conservation and exchange of momentum and energy

Theoretical prediction

The parallel momentum evolution of species a due to collisions on species b during a time ∆t can be computed analytically if both distribution functions are truncated. In this case one gets

∆P th ab = ∆tC ab F t a , F t b m a v d 3 v = F t a m a ∆v ab d 3 v (23)
where the truncated distribution function reads

F t a = N a π 3/2 v 3 T a e -s 2 a 1 - V 2 a v 2 T a 1 - 2 3 s 2 a + 2ξs a V a v T a - 4q a m a N a v 3 T a ξs a 1 - 2 5 s 2 a ( 24
)
with ξ = v v and the statistical average of the evolution of the parallel velocity is

∆v ab = νab ∆tv T a P ab (25)
where

P ab = v T a v T b 2 1 + m a m b      ξ - E b -s b E b s 2 b + 2V 2 b 3v 2 T b s b E b + V b v T b 1 -3ξ 2 E b s 3 b + ξ 2 3 + 2s 2 b -1 E b s 2 b + 4q b 5m b N b v 3 T b 2s 2 b ξ 2 -1 E b      (26)
It can be shown (see appendix B) that the momentum exchange due to the collision of a species a on a species b takes the form

∆P th ab = ν ab ∆tN a m a         -V a 1 - V 2 b v 2 T b X 2 ba + V b 1 - V 2 a v 2 T a X 2 ab + 6 5 q a maNav 2 T a X 2 ab 1 -5 3 V 2 b v 2 T b X 2 ba -6 5 q b m b N b v 2 T b X 2 ba 1 -5 3 V 2 a v 2 T a X 2 ab         (27)
where

X ab = 1 + v T b v T a 2 -1 2 and ν ab = 4 3 √ π νab X 3 ab 1 + ma m b = 4 3 √ π N b v 3 T a e 2 a e 2 b m 2 a ln Λ 4π 2 0 X 3 ab 1 + ma m b
is the usual collision frequency. Similarly, the energy evolution of the species a due to collisions on species b during a time ∆t can be computed analytically if both distribution functions are truncated. In this case one gets

∆E th ab = C ab F t a , F t b 1 2 m a v 2 d 3 v = F t a 1 2 m a ∆v 2 ab d 3 v (28)
and the statistical average of the evolution of the energy

∆v 2 ab = νab ∆tv 2 T a E ab (29)
with at leading order in ∆t

E ab = 2 v T a v T b      E b s b -1 + ma m b E b -s b E b s b + V 2 b 3v 2 T b E b 2 1 + ma m b s 2 b -1 + V b v T b ξ E b -s b E b s 2 b + 2 1 + ma m b -E b s 2 b + 1 + s 2 b E b s b + 4 5 q b m b N b v 3 T b ξs b E b 1 + ma m b 2s 2 b -1 -1      (30)
It can be shown that (appendix C)

∆E th ab N a m a m a + m b ν ab ∆t                     3 (T b -T a ) + X ba V b v T b 2 2 m b ma + 3 T a -T b -X ab V a v T a 2 2 ma m b + 3 T b -T a + V a v T a V b v T b X ab X ba 2 2 ma m b + 3 T b -2 m b ma + 3 T a +2 V a v T a V b v T b X ab X ba ma m b + 2 T b -m b ma + 2 T a + 12 5 V b v T b q a maNav 3 T a X ba X 3 ab m b ma + 2 T a -4 + 3 ma m b T b -12 5 V a v T a q b m b N b v 3 T b X ab X 3 ba ma m b + 2 T b -4 + 3 m b ma T a + 24 5 q a maNav 3 T a q b m b N b v 3 T b (X ab X ba ) 3 4 + 3 ma m b T b -4 + 3 m b ma T a                     (31)
Using the fact that N a m a ν ab = N b m b ν ba , one gets the conservation of momentum ∆P th ab = -∆P th ba and energy ∆E th ab = -∆E th ba .

Numerical handling of conservation properties

The theoretical exchange rates are not exactly ensured numerically due to the finite number of markers, the finite time step, and the truncation of distribution function. This leads to inaccuracies in the conservation properties. To ensure conservations up to machine precision, the marker weights are modified after the Langevin kicks. One can imagine different ways to ensure conservations. But as pointed out by Lin [START_REF] Lin | Gyrokinetic particle simulation of neoclassical transport[END_REF], an enforcement of conservation without respecting the velocity dependence of the collision operator could lead to an unphysical distortion of the distribution function and ultimately discrepancies with respect to the neoclassical predictions. In the context of intra-species linearized collision operator, Satake [START_REF] Shinsuke | Development of a non-local neoclassical transport code for helical configurations[END_REF] implemented an approximated version of the field operator to correctly ensure conservations. This strategy was already used in the previous version of the collision operator (intra-species linearized collisions) in ORB5 [START_REF] Vernay | Neoclassical equilibria as starting point for global gyrokinetic microturbulence simulations[END_REF]. This strategy has been generalized in order to ensure the right exchange rates of parallel momentum and energy and particle conservation when applying the truncated nonlinear multi-species collision operator. It leads to the definition of a correction term of the form

C cor ab (v, ξ) = νab F M a (v) [a + bM 1,ab (v, ξ) + cM 2,ab (v, ξ)] (32)
where

M 1,ab (v, ξ) = F t a F M a P ab M 2,ab (v, ξ) = F t a F M a E ab
with P ab , E ab defined respectively by Eq.( 26) and Eq. [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF]. In practice, the ensemble average are applied on spatial bins. The coefficients a, b and c are adapted in order to have exactly the right exchange of momentum and energy in each spatial bin. In the case of full-F PIC codes, the application of this correction term is straightforward. On the other hand, ORB5 [START_REF] Lanti | Orb5: A global electromagnetic gyrokinetic code using the pic approach in toroidal geometry[END_REF] is a δF PIC code for which the full distribution function is split in two, a background noted F 0 and a a perturbation noted δF . Both distributions are represented with the same markers but different weights. Marker weights denoted with a p correspond to the background and the one noted with w correspond to the perturbation. Because of the presence of these two weights, there is a priori a freedom to ensure the right exchange rates of momentum and energy. We choose to modify the p-weights to handle the exchanges and to modify the w-weights to ensure a conservation of the three first moments of δF for each species. The idea behind this choice is that in principle the p-weights are much bigger than the w-weights. So the perturbations induced by the method should have a weaker relative effect compared to a case where corrections are performed on w-weights. Then one formally introduces two correction terms, one for p weights C cor,p ab and the other one for w weights C cor,w ab . The evolution of F 0 under this correction term is given by

dF 0 dt = C cor,p ab
In the ORB5 code, the background distribution function is represented by

F 0 = gP where g = N i=1 δ[z-zi(t)] J(z)
and P is the field of p-weights such that P (z r ) = p r . As the marker positions are untouched in this procedure, dg dt = 0. Therefore

dF 0 dt = g dP dt = F 0 P dP dt = C cor,p ab
Let us assume that the background distribution is equal to the local unshifted Maxwellian F 0 = F M a , then one has, for particle number r,

p a out,r -p a in,r = p a in,r ∆t C cor,p ab F M a = p a in,r a α p + b α p M 1,ab (v a out,r , ξ a out,r ) + c α p M 2,ab (v a out,r , ξ a out,r ) (33)
where the coefficients a α p , b α p , c α p have been redefined to include the factor νab ∆t. The numerical evolution of the three first moments of the background distribution function in a space bin α of volume V α after the application of the truncated collision operator C t ab and the correction term C cor,p ab reads

V α ∆N num p,a = r∈α p a out,r -p a in,r (34)
V α ∆P num p, ab = r∈α m a v a out,r p a out,r -v a in,r p a in,r (35)
V α ∆E num p,ab = r∈α m a 2 v a out,r 2 p a out,r -v a in,r 2 p a in,r (36)
One wants N num a = 0, ∆P num p, ab = ∆P th ab and ∆E num p,ab = ∆E th ab . Injecting Eq.(33), one gets the system of equations

M α ab   a α p b α p c α p   = V α   0 ∆P th ab ∆E th ab   -N α p,ab (37)
where

M α ab = r∈α p a in,r   1 m a v a out,r ma 2 v a out,r 2   1 M 1,ab (v a out,r , ξ a out,r ) M 2,ab (v a out,r , ξ a out,r ) N α p,ab = r∈α p a in,r     0 m a v a out,r -v a in,r ma 2 v a out,r 2 -v a in,r 2    
By a similar approach, one gets the correction term for the w-weights

w a out,r -w a in,r = p a in,r a α p + b α p M 1,ab (v a out,r , ξ a out,r ) + c α p M 2,ab (v a out,r , ξ a out,r )
where the coefficients are computed by solving the system

M α ab   a α w b α w c α w   = -N α w,ab
where

N α w,ab = r∈α w a in,r     0 m a v a out,r -v a in,r ma 2 v a out,r 2 -v a in,r 2    
The procedure described here ensures the conservation of density and the right exchange of momentum and energy after the collision of species a on the species b to machine precision in each spatial bin α.

Linearization of the collision operator

The model collision operator presented in this paper takes the form Eq.( 13)

C mod ab = C t ab + C cor ab
It is shown in this section that when only the component l = k = 0 is kept, C t ab corresponds to the test-particle operator and C cor ab corresponds to the field operator developed by Lin [START_REF] Lin | Gyrokinetic particle simulation of neoclassical transport[END_REF]. When only the component l = k = 0 is kept, the drag Eq.( 20) reduces to

Γy L,ab = 0 Γz L,ab = - E b -s b E b s 2 b
and the components of the diffusion matrix Eq.( 21) reduce to

Dyy L,ab = s 2 b - 1 2
E b s 3 b + E b 2s 2 b Dzz L,ab = E b s 3 b - E b s 2 b
Dyz L,ab = 0 Therefore, the Langevin kicks take the form

∆v L,ab v T a = νab ∆t v 2 T a v 2 T b 1 + m a m b Γz L,ab e z + νab ∆t v T a v T b      R b x Dyy L,ab R b y Dyy L,ab R b z Dzz L,ab      (38)
which is exactly the test-particle operator. The correction term developed in the part 3.2 is actually an extension of the procedure developed by Lin [START_REF] Lin | Gyrokinetic particle simulation of neoclassical transport[END_REF] to approximate the field operator in the context of like species collision. Indeed, when the correction term Eq.( 32) is applied with the simplifications due to linearization one gets

C cor,L ab = νab F M a (v) a + bM L 1,ab (v, ξ) + cM L 2,ab (v, ξ)
where

M L 1,ab = P L ab = - v v T b v T a v T b 2 1 + m a m b E b -s b E b s 3 b M L 2,ab = E L ab = 2 v T a v T b E b s b -1 + m a m b E b -s b E b s b
In the intra-species collision case, the approximated field-particle operator already implemented in several codes is retrieved [START_REF] Lin | Gyrokinetic particle simulation of neoclassical transport[END_REF][START_REF] Shinsuke | Development of a non-local neoclassical transport code for helical configurations[END_REF][START_REF] Vernay | Neoclassical equilibria as starting point for global gyrokinetic microturbulence simulations[END_REF]. In the multi-species case, the exchange rates of momentum Eq.( 27) and energy Eq.( 31) need to be adapted

∆P th,L ab = ν ab ∆tN a m a V b -V a + 6 5 q a m a N a v 2 T a X 2 ab - q b m b N b v 2 T b X 2 ba (39) ∆E th,L ab = N a m a m a + m b ν ab ∆t     3 (T b -T a ) + X ba V b v T b 2 2 m b ma + 3 T a -T b -X ab V a v T a 2 2 ma m b + 3 T b -T a     (40)

Moment conserving noise control

The Langevin approach used to implement the collision operator introduces weight spreading and therefore increases the numerical noise [START_REF] Brunner | Collisional delta-f scheme with evolving background for transport time scale simulations[END_REF]. The standard noise control procedure used in collisionless simulations in ORB5 [START_REF] Lanti | Orb5: A global electromagnetic gyrokinetic code using the pic approach in toroidal geometry[END_REF] relies on a modified Krook to ensure conservation of certain moments of the distribution function. This procedure is not compatible with the presence of collisions because it artificially modifies the distribution function with a rate similar or greater than typical collision frequencies. To reduce noise, a procedure adapted to the presence of collisions and described in [START_REF] Sonnendrücker | A split control variate scheme for pic simulations with collisions[END_REF] is therefore adopted. The idea is to smooth locally the weights using a gridless procedure in velocity space. The particles are first sorted in different space bins. Then, in each space bin, the velocity phase space is split into sub-domains using a quadtree procedure and having fixed the maximum number of particles inside each sub-domain. Once the sub-domains are defined, particles are grouped randomly by packs of four and their weights are modified in order to minimize the dispersion while maintaining conservations of mass, parallel momentum and energy. The procedure is identical for the p-weights and w-weights. It is presented only in the case of p-weights. The conservation properties read

i p cont i = i p old i = N p i v i p cont i = i v i p old i = P p i v 2 i p cont i = i v 2 i p old i = E p
The conservations give three constraints allowing to express the weights of the particles 2 to 4 as linear functions of the weight of the first particle.

p cont i = α i p cont 1 + β i with i ∈ [2, 4]
where

α 2 = v 2 4 -v 2 1 + v 1 -v 4 v 3 -v 4 v 2 3 -v 2 4 v 2 2 -v 2 4 - v 2 -v 4 v 3 -v 4 (v 2 3 -v 2 4 ) β 2,p = E p -N p v 2 4 - Pp-Npv 4 v 3 -v4 v 2 3 -v 2 4 v 2 2 -v 2 4 - v 2 -v 4 v 3 -v 4 (v 2 3 -v 2 4 ) α 3 = v 4 -v 1 + v 4 -v 2 α 2 v 3 -v 4 β 3,p = P p -N p v 4 + v 4 -v 2 β 2,p v 3 -v 4 α 4 = -1 -α 2 -α 3 β 4,p = N p -β 2,p -β 3,p p cont 1 is chosen to minimize the weights dispersion i<j p cont i -p cont j 2
. This is equivalent to the minimization of

i (p cont i )
2 . This minimization is obtained for

p cont 1 = - α 2 β 2,p + α 3 β 3,p + α 4 β 4,p 1 + α 2 2 + α 2 3 + α 2 4
The strength of the correction is chosen to be proportional to the mean distance between particles

p new i = Γp cont i + (1 -Γ) p old i with i ∈ [1, 4] with Γ = exp - (vi-vj ) 2 h 2 v where (v i -v j ) 2
is the mean square distance between particles and h v corresponds to the typical distance of smoothing in the velocity space. Note that this noise control can also be applied for collisionless simulations [START_REF] Sonnendrücker | A split control variate scheme for pic simulations with collisions[END_REF].

Verification of the collision operator

The implementation of the collision operators (linear and nonlinear) described previously has been verified with different physical tests described in this section. For all cases, conservation of particles, total momentum and energy is ensured up to machine precision in each real space bin thanks to the procedure described in section 3.2.

Isotropisation

The isotropisation of the distribution function due to collisions is a consequence of the H-theorem. For a single species simulation in the absence of drifts, the evolution of the isotropisation can be computed. Let us assume that perturbed distribution function takes the form

δF (s, ξ, t) = λ (s, t) 2 F M s 2 3ξ 2 -1 (41)
Due to the symmetry of the problem in the v direction, V and q are both equal to zero at all times. Therefore the truncated nonlinear collision operator C t ab reduces to the test-particle operator C T P ab which can be explicitly be expressed as a pitch angle scattering and a thermalization term

C T P = 1 2s 2 ∂ ∂s F M ν v s 4 ∂ ∂s δF F M + ν d 2 F M ∂ ∂ξ 1 -ξ 2 ∂ ∂ξ δF F M ( 42
)
where ν v is the velocity modulus diffusion rate

ν v = 2ν G (s) s 3 ( 43
)
ν d is the deflection frequency

ν d = ν [erf (s) -G (s)] s 3 (44)
and

G (s) = 1 2s 2 erf (s) -2 √ π se -s 2
is the Chandrasekhar function. The injection of Eq.(41) in Eq.( 42) leads to

∂λ ∂t = ν a ∂ 2 λ ∂s 2 + ν b ∂λ ∂s -ν T λ (45) with ν a = ν G s , ν b = ν erf(s)+(3-4s 2 )G(s) s 2
and ν T the characteristic rate of anisotropy relaxation

ν T = ν erf (s) + 8s 2 -3 G (s) s 3 (46)
Eq.( 45) implies that if the initial perturbed distribution function is of the form of Eq.( 41), it will keep this form when the distribution function evolves only due to collisions. Assuming that the initial perturbed distribution function is of the form of Eq.( 41) with an initial value λ (s, t = 0) = λ 0 independent of s, it is possible to show (appendix D) that

λ (s, t) λ 0 exp (-ν T t) 1 -ν a d 2 ν T ds 2 + ν b dν T ds t 2 2 + ν a dν T ds 2 t 3 3 (47)
The anisotropy with a distribution function of the form Eq.( 41) is given by

T -T ⊥ (t) = m N v 2 - v 2 ⊥ 2 δF d 3 v = 8T 5 √ π ∞ 0 λ (s, t) e -s 2 s 6 ds
Injecting Eq.( 47), one gets

T -T ⊥ (t) T -T ⊥ (t = 0) 16 15 √ π ∞ 0 e -ν T t 1 -ν a d 2 ν T ds 2 + ν b dν T ds t 2 2 + ν a dν T ds 2 t 3 3 e -s 2 s 6 ds (48)
In the general case, Eq.( 48) needs to be computed numerically. For early times, the isotropization rate can nevertheless be computed analytically by Taylor expansion at first order in time

T -T ⊥ (t) T -T ⊥ (t = 0) = 1 - 2 √ 2 5 νt ∞ 0 10s 2 -3 erf (s) + -8s 2 + 3 2 √ π se -s 2 e -s 2 sds + O ν 2 t 2 = 1 -1.6νt + O ν 2 t 2
where ν = 4ν 3 √ 2π . Note that there is a factor 2 missing in the isotropization rate computed in [START_REF] Donnel | A multi-species collisional operator for full-f global gyrokinetics codes: Numerical aspects and verification with the gysela code[END_REF]. Numerically, an initial perturbed distribution function of the form Eq.(41) with λ (v, t = 0) = K = 0.1 is chosen. Drifts are turned off, meaning that the distribution function evolves only due to intra-species collision d t F = C and noise control is activated. The test is performed with 10 7 markers and constant background profiles. The numerical evolution of the anisotropy is compared to the theoretical prediction Eq.( 48). Using the global feature of ORB5 and the local nature of this test, it is possible with one simulation to have the equivalent of several simulations with less markers by using spatial binning to group markers. It allows one to plot error bars as ensemble average estimates of the standard deviation for the evolution of the isotropisation. In this simulation, 16 radial bins are used. The comparison between numerical and theoretical solutions is shown in Fig. 1.

Exchange of momentum and energy between species

The theoretical exchange rates of momentum and energy between species have been computed in section 3.1. These exchange rates are numerically tested with simulations with two species. More specifically, hydrogen and deuterium ions are simulated. Initially, distribution functions are chosen to be Maxwellians with identical densities, temperatures but opposite mean velocities. As in the previous part, drifts are turned off, meaning that the distribution functions evolve only due to collisions d t F a = b C (F a , F b ) and noise control is activated. The test is performed with 10 7 markers for each species and constant background profiles. Initial Mach numbers are of order one to illustrate differences between linear and nonlinear collisions. Two simulations with the same initial distributions but different collision operators (linear or nonlinear) have been performed. In Fig. 2 exchange rates of momentum and energy obtained with the linear collision operator are compared with their theoretical predictions, Eq.( 39) and (40), respectively. In Fig. 3 exchange rates of momentum and energy obtained with the nonlinear collision operator are compared with their theoretical predictions, Eq.( 27) and [START_REF] Ohana | Towards the optimization of a gyrokinetic particle-in-cell (PIC) code on large-scale hybrid architectures[END_REF], respectively. Note that as predictions depend of instantaneous moments of the distribution function, the theoretical predictions (linear and nonlinear) are dependent of the collision operator used in the numerical simulation.

Neoclassical tests

In this section, neoclassical physics is tested in the simple case of a single ion species with an adiabatic electron response. An analytical axisymmetric ad hoc equilibrium model with circular concentric flux surfaces is used. It is described as

B = B 0 R 0 R r qR 0 e θ + e ϕ (49)
where θ, ϕ are respectively the geometrical poloidal and toroidal angles. q is the pseudo-safety factor, linked to the safety factor q = 1 2π 2π 0 B•∇ϕ B•∇θ dθ by the relation q = q √ 1 -2 . This simplified equilibrium is used for the numerical tests described in this section. Note that in this section the "safety factor" will refer to q that will be denoted q to simplify the notations. With this this simplified geometry, the flux surface average of a quantity A is defined as

A ψ = 2π 0 ARdθ 2π 0 Rdθ = 1 2π 2π 0 A (1 + cos θ) dθ
Numerically, the drift-kinetic approximation is used for charge deposition consistently with the derivation of neoclassical physics and the derivation of the collision operator. The non axisymmetric components of the electric potential are filtered out at each time step, to remove turbulence. The weak neoclassical transport allows to reach a quasi-steady state even in absence of sources to maintain fluxes.

Three quantities are compared with their theoretical predictions: the ion thermal conductivity χ, the poloidal rotation and the collisional damping of the zonal flow residual.

Heat diffusivity

A theoretical prediction for the heat diffusivity in all collisionality regimes in the presence of an impurity has been derived by Chang and Hinton [START_REF] Chang | Effect of impurity particles on the finite aspect ratio neoclassical ion thermal conductivity in a tokamak[END_REF]. In the case of single ion species plasma, the neoclassical ion thermal conductivity takes the form

χ th neo = χ BP + χ P S where χ BP ρ 2 Li ω ci = 2 √ 2q ρ Li R 0 ν 0.66 + 1.88 √ -1.54 1 + 1.03 √ ν + 0.31ν B 2 0 B 2 ψ χ P S ρ 2 Li ω ci = 2 √ 2q ρ Li R 0 ν 0.58ν 1 + 0.74ν 3/2 B 2 0 B 2 ψ - B 2 B 2 0 -1 ψ with ρ Li = Ti mi ω -1
ci the ion Larmor radius, ω ci = eiB mi the ion cyclotron frequency, = r R0 the local inverse aspect ratio and ν = qR0ν √ 2

3/2 v T the ion collisionality. In the case of a circular plasma without Grad-Shafranov shift, one has

B 2 0 B 2 ψ = 1 + 3 2 2 and B 2 B 2 0 -1 ψ = 1 -2
This theoretical prediction is compared with the numerical ion thermal conductivity computed as q H = -χ num N ∂T ∂r where q H = Q-5 2 T Γ is the heat flux, Q is the energy flux and Γ the particle flux. Collisionality scans have been performed with either the linear or the nonlinear collision operators. The results of these scans are depicted in Fig. 4, which shows an excellent agreement of ORB5 results with the theoretical prediction. The differences between the linear and the nonlinear collision operators are small in this case, as expected. Each simulation represented in Fig. 4 is performed with 3 • 10 8 markers. The simulations are performed for plasmas with ρ = 1 150 , inverse aspect ratio at the edge a R0 = 0.2, flat safety factor q = 1.6 and density profiles. Finally the temperature profile takes the form

d ln T d (r/a) = -κ T cosh -2 r -r ref ∆ T
with r ref = 0.5a, ∆ T = 0.2a and κ T = 1.2. This choice of profile for density and temperature is arbitrary as the neoclassical prediction is insensitive to the density gradient and proportional to the local gradient of temperature.

Poloidal rotation

In the neoclassical theory, the ion poloidal rotation is linked to the temperature gradient by the relation

V θ = k neo B ϕ e B 2 ψ

∂T ∂r

Different predictions of the coefficient k neo can be found in the literature. A first one derived by Kim et al. [START_REF] Kim | Neoclassical poloidal and toroidal rotation in tokamaks[END_REF] is based on approximate parallel viscosity coefficients which are obtained by a rational combination of the asymptotic viscosity coefficients in the various collisionality regimes. As the validity domain of the Plateau regime 1 ν 1 3/2 is never really reached with actual values of aspect ratio, this prediction is in fact accurate only in the Banana regime and the Pfirsch-Schlüter regime but only approximated in-between.

k Kim neo = √ 2µ K 01 µ K 00 µ K 11 + √ 2 -µ K 01 2
where, for i = 0, 1

µ K ii = g μB ii 1 + 2.92ν Kim μB ii μP ii 1 + qR0ν √ 2 6v T μP ii μP S ii with g = ft
fc where f c is the circulating fraction, f t = 1 -f c is the trapped fraction and

ν Kim = g 1.46 ν √ 2 qR0v T B 2 ψ 2 (e •∇B) 2 ψ .
In the large aspect ratio limit, one has f t 1.46 √ -0.46 3/2 and ν Kim

g 1.46 qR0ν √ 2 2 v T . The coefficient µ K 01 is computed via the relation µ K 01 = 2.5µ K 00 -K K 01 where K K 01 = g KB 01 1 + 2.92ν Kim KB 01 KP 01 1 + qR0ν √ 2 6v T KP 01 KP S 01
The different coefficients for the Kim prediction are summed up in the table Tab.1.

A second theoretical prediction has been proposed by Schaing [START_REF] Shaing | An approximate analytic expression for plasma viscosity in finite aspect ratio tokamaks and its applications[END_REF]. As Kim's prediction, it relies on a moment approach of the neoclassical physics with approximate parallel viscosity coefficients. But the approximation is more accurate than the one made Kim's, especially in the Plateau regime. One should, however, note an important issue concerning the validity of Shaing's prediction in the Pfirsch-Schlüter regime. Indeed, even if the parallel viscosity coefficients derived by Shaing Eq.( 51) are valid in the Pfirsch-Schlüter regime, their link with k neo , Eq.( 50), has been derived by neglecting the parallel friction. This approximation is not valid for the Pfirsch-Schlüter regime [START_REF] Helander | Collisional Transport in Magnetized Plasmas[END_REF]. Having this remark in mind, Shaing's prediction takes the form

k Schaing neo = - µ S 2 √ 2ν µ S 1 µ S 3 -µ S 2 2 + µ S 1 √ 2ν (50) with   µ S 1 µ S 2 µ S 3   = 8 3 √ π ∞ 0   1 s 2 -5 2 s 2 -5 2 2   K B K P S K B + K P S s 4 e -s 2 ds (51)
The Banana contribution is given by

K B (s) = f t f c ν D (s)
where ν D is the deflection frequency Eq.( 44). The Plateau and Pfirsch-Schlüter contributions are given by

K P S = 3 2 v 2 T s 2 ∞ m=1 F m ν T I m R ν T (52) F m = 2 B 2 ψ B • ∇θ ψ sin (mθ) e • ∇B ψ × sin (mθ) (B • ∇θ) e • ∇B ψ + cos (mθ) e • ∇B ψ × cos (mθ) (B • ∇θ) e • ∇B ψ and ν T I m R = - 3 2 ν T ω m 2 - 9 2
ν T ω m 4 + 1 2 + 3 + 9 2 ν T ω m 2 ν T ω m 2 ν T ω m tan -1 ω m ν T
where ω m = sv T me • ∇θ and ν T is the characteristic rate of anisotropy relaxation Eq.( 46). Note that these expressions are derived in a set of coordinates such that the Jacobian takes the form √ g = F (ψ) B•∇θ where F is an arbitrary function of the poloidal flux and e • ∇θ is a flux function. These hypothesis are fulfilled with the simplified geometry of circular concentric flux surfaces used in this section Eq.( 49). Indeed, in this case the Jacobian takes the form √ g = rR and

e • ∇θ = 1 qR0 1+( q )
2 which is a flux function. Then one has

B 2 ψ = B 2 0 1 + 2 1 2 + 1 q 2 + O 4 B • ∇θ = B 0 qR B • ∇θ ψ = B 0 qR 0 e • ∇B = B 0 R 0 qR 2 sin θ
and so

sin (mθ) e • ∇B ψ = 1 2π B 0 qR 0 2π 0 sin (mθ) sin θ (1 + cos θ) -1 dθ = B0 2qR0 + O 3 if m = 1 O 2 otherwise sin (mθ) (B • ∇θ) e • ∇B ψ = 1 2π B 0 qR 0 2 2π 0 sin (mθ) sin θ (1 + cos θ) -2 dθ =    2 B0 qR0 2 + O 3 if m = 1 O 2 otherwise cos (mθ) e • ∇B ψ = 1 2π B 0 qR 0 2π 0 cos (mθ) sin θ (1 + cos θ) -1 dθ = 0
At leading order in , the sum in Eq.( 52) is therefore dominated by the term m = 1

K P S = 3 8 2 ω 1 -3 ν T ω 1 -9 ν T ω 1 3 + 1 + 6 + 9 ν T ω 1 2 ν T ω 1 2 tan -1 ω 1 ν T + O 4
where

ω 1 = sv T qR 0
Finally a third prediction, obtained by Sauter et al. [START_REF] Sauter | Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime[END_REF] is based on a numerical fitting of multiple neoclassical simulations performed with the code CQL3D. One should be aware that there is an errata [START_REF] Sauter | Erratum: Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime[END_REF] for one sign on this prediction that finally reads

k Sau neo = - α 0 + 0.25 1 -f 2 t √ ν 1 + 0.5 √ ν + 0.315ν 2 f 6 t 1 1 + 0.15ν 2 f 6 t where α 0 = - 1.17fc 1-0.22ft-0.19f 2 t .
The coefficient k neo is numerically computed using the usual relation of the neoclassical theory

k num neo = 1 ∂ r T T N ∂ r N + ∂ r T + e∂ r φ -e q BV ψ 1 ∂ r T T N ∂ r N + ∂ r T + e∂ r φ -e q B 0 V
The results are shown in Fig. 5. An excellent agreement of ORB5 results with Shaing's prediction is observed. The differences between the linear and the nonlinear collision operator are small in this case as expected. The simulations used to obtain Fig. 5 are exactly the same as the ones described in the section 6.3.1.

Collisional damping of the zonal flow residual

Zonal flows are key players in the saturation of turbulence [START_REF] Diamond | Zonal flows in plasmas -a review[END_REF]. It is then critical for turbulent codes to obtain the correct level of zonal flows. Collisionless plasmas are known to exhibit a linearly undamped axisymmetric component of the potential φ 00 [START_REF] Rosenbluth | Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks[END_REF]. Hinton and Rosenbluth [START_REF] Hinton | Dynamics of axisymmetric exb and poloidal flows in tokamaks[END_REF] have predicted the collisional damping of this residual flow.

φ 00 (t) φ 00 (t = 0) = exp β 2 α 2 t erfc β α √ t + erf (ν 0 t) α B 2 θ B 2    1.8 (ν 0 t) 5/9 exp -(3ν 0 t) 2/3 +1 + 1.4 exp (-ν d t)    (53)
Where α = 1 + 1.6 q 2 √ is the collisionless residual, β = 3πq 2 ν1/2 Λ 3/2 , Λ = λ ln 16 νt 1/2 with λ 1. Different frequencies associated with collisions are defined as ν = 0.61ν ii , ν 0 = 1.9ν ii , ν d = νii 0.64 √ . In this formula, the initial potential φ 00 (t = 0) corresponds to a case for which all GAM oscillations have been Landau damped.

The first term of Eq. ( 53) corresponds to the fast collisional smoothing of the distribution function in the trapping boundary layer. The last term is the contribution from high-energy ions which have small collision frequencies and thus dominate the long time behaviour. An extra term erf (ν 0 t) has been added in front of the last term of the right hand side to ensure the relation also holds for initial time [START_REF] Donnel | A multi-species collisional operator for full-f global gyrokinetics codes: Numerical aspects and verification with the gysela code[END_REF].

The collisional damping of the residual can be tested by first launching a collisionless simulation, waiting for most of the GAM to be damped then adding collisions. When collisions are activated, a new burst of GAM is observed as can be seen in Fig. 6. A good agreement is anyway found between the predicted and computed φ 00 decay once collisions are activated.

High Performance Implementation

Gyrokinetic codes like ORB5 need outstanding computational resources to solve increasingly complex problems. The effective exploitation of state-of-the-art High Performance Computing (HPC) architectures represents therefore a mandatory task to tackle the most challenging questions in plasma physics. Current HPC performance is being achieved through massively parallel multi-CPU, many-core, heterogeneous systems equipped with accelerators like GPUs. The ability to efficiently exploit such architectures is therefore critical.

ORB5 addresses this challenge implementing two levels of parallelism (see [START_REF] Ohana | Gyrokinetic simulations on many-and multicore architectures with the global electromagnetic particle-in-cell code orb5[END_REF] for details). The first level is based on the distribution of data among different processes, using the Message Passing Interface (MPI) [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF]. Two different MPI based solutions are implemented. The first adopts the domain decomposition along the toroidal dimension, where each subdomain works with field and particle data of a toroidal slice. This slicing is uniform in order to guarantee a balanced workload. Hence, the number of subdomains is constrained to be a divider of the number of toroidal grid cells, which is frequently a power of 2 for best fast Fourier transform performance. The second MPI-based solution is domain cloning. It consists in replicating field data subdomains on different tasks, splitting the particle workload. MPI communications among clones are reductions of the grid data before the field solvers, and broadcasts afterwards. Domain cloning allows one to run with more MPI tasks than the number of toroidal grid points. In addition, it proved to be useful on multi-socket nodes, where the optimum configuration is to use one clone per socket.

The second level of parallelism corresponds to multi-threading. For multi-core CPUs, it has been implemented by using the OpenMP1 application programming interface (API). OpenMP represents a consolidated and stable solution to exploit shared memory devices, ensuring good scalability up to many cores and high portability. It is based on a set of directives which instruct the compiler on how to split the work among the available threads, with minimum impact on the source code, although specific customization has to be done in order to ensure good performance and scalability. In ORB5, OpenMP is used exclusively in routines looping over the particles. Such loop iterations are in general independent from one another, which make them trivially parallelizable over OpenMP threads.

An alternative to the use of OpenMP has been recently introduced in ORB5 [START_REF] Ohana | Gyrokinetic simulations on many-and multicore architectures with the global electromagnetic particle-in-cell code orb5[END_REF] in order to effectively exploit modern GPU accelerated systems. The GPU implementation adopts the OpenACC2 programming model to enable to the GPU the most time consuming parts of the code. Likewise OpenMP, OpenACC is a directive-based approach, which allows one to keep the source code almost unchanged, leaving to the compiler the task of creating the GPU version. OpenACC provides two main classes of directives: those to implement work sharing and those to control data transfer. Both are designed with performance as the main objective. In ORB5, all the computing steps involving particles have been offloaded to the GPU, with only the field operations performed on the host. Loops over particles are particularly well suited for GPUs since millions of independent iterations can be executed by thousands of threads. Communication between GPU and CPU is minimized by keeping particle data on the GPU only.

The collision operator, acting on the particles, supports all three levels of parallelism. Particles are split among MPI tasks, each calculating collisions on their set of points. Communication is limited to the handling of conservation properties, detailed in section 3.2, as spatial bins are chosen to lay on 2D poloidal plane. Note that the computation of local moments of the distribution function and the noise control presented in the section 5 are not specific to collisions and are therefore not accounted for in the performance analysis presented below.

All main loops within the collision operator are multi-threaded using OpenMP on the CPU and OpenACC on the GPU. For most loops this has a straightforward implementation, each particle being evolved independently from the others. Only in two cases, both corresponding to the handling of conservation properties, specific care has to be taken in order to manage race conditions, with particles properties being binned on a regular spatial mesh. This is accomplished adopting the atomic directive, that was demonstrated in [START_REF] Ohana | Towards the optimization of a gyrokinetic particle-in-cell (PIC) code on large-scale hybrid architectures[END_REF], to provide an effective solution for thread-safe accumulation operations.

The GPU implementation adopts the cuRAND3 library for random number generation. The cuRAND library focuses on the simple and efficient generation of high-quality pseudorandom numbers directly on the GPU, avoiding any hostdevice data movement. Several algorithms are available, highly optimized for the NVIDIA GPU architectures. In case cuRAND is not available, or non NVIDIA GPUs are used, the random numbers generator of CPU can be adopted. In this case, the generated data have to be moved from the host to the device at each call, causing a performance drop due to the copy overhead. The cuRAND functions can be used both on the GPU and on the CPU, generating exactly the same random number sequence for the same initial seed. This facilitates the effective comparison of CPU and GPU runs.

The performance and the scalability of the collision operator has been analyzed on the Piz Daint HPC platform, operated by the Swiss National Supercomputing Center (CSCS -https://www.cscs.ch). Piz Daint consists of 5704 nodes, each combining an Intel Xeon E5-2690 v3 (Haswell) CPU and an NVIDIA Tesla P100 GPU. The nodes are interconnected by a Dragonfly Aries network. The adopted Fortran compiler is PGI 19.7.0, supporting OpenACC 2.7 and CUDA toolkit 10.1.

Single node performance

The single node performance of the collision operator has been tested adopting an overall configuration representative of a typical physical problem and fitting the GPU memory, with 2 • 10 7 markers (adiabatic electrons are used) and a computational mesh of 192 × 192 × 96 cells, in the radial, poloidal and toroidal directions respectively. The grid for spatial binning in the collision operator is 32 × 32 cells. Fig. 7 show the results of both the linear (L) and the nonlinear (NL) solver for the four main algorithmic components of the operator: the random number generator, the Langevin kick, the moment evolution due to Langevin kicks and the moment correction. The tests on the CPU have been run using all the 12 cores available on the node. The results of the MPI implementation are shown. Results obtained with the OpenMP implementation are similar, and not shown here for brevity. The GPU tests are run using only one MPI task. Ten simulation time steps are performed for each run. The random number generator takes approximately the same time on the GPU and on the CPU, where the work is split, each core computing 1/12 of the entire random set. In the GPU implementation, the adoption of the cuRAND library, calculating the random numbers directly on the GPU, avoids any overhead due to data movement. In case the cuRAND library cannot be used, numbers can be generated on the host, using the CPU algorithm. However, as already stated above, this leads to a data transfer overhead to move data from the processor to the device.

For the remaining three algorithmic components, the GPU gives a significant performance boost, since each of them consists in a large parallel loop on particles, with no dependencies between loop iterations, high arithmetic intensity and effective memory access (as detailed in [START_REF] Ohana | Towards the optimization of a gyrokinetic particle-in-cell (PIC) code on large-scale hybrid architectures[END_REF] and [START_REF] Hariri | A portable platform for accelerated pic codes and its application to gpus using openacc[END_REF]).

On the CPU the calculation of the NL kicks is about 3 times slower than the L version. For the GPU implementation the NL and L algorithms take the same time, slightly less than that of the L kick operator running on the CPU. The computation of the moments evolution is the same for the L and NL cases, hence related timings are comparable on both the CPU and on the GPU. The GPU algorithm results to be approximately 5.5 times faster than the CPU one. Finally, the higher performance boost on the GPU is given by the moment correction algorithm. On the CPU this represents both in its linear and non linear implementation the most time consuming part. The GPU version is approximately 5.5 (L) and 5.9 (NL) times faster than the CPU.

Overall, the NL operator takes around 32 sec on the 12 cores of the CPU, while it takes ∼7.8 sec on the GPU, to run 10 time steps. Hence, the GPU implementation results to be around 4 times faster than the CPU. The L operator runs in 20.8 sec on the CPU and 7.6 sec on the GPU, the GPU being 2.7 times faster than the CPU. We note that the higher speed-up is achieved for the NL operator, which is the most computational intensive of the two operators.

Scalability

The parallel scalability of the collision operator has been checked on two configurations. The first is adopting a grid for spatial binning of 32×32 cells and using the same number of markers and the same mesh as in Section 7.1 scaling from 1 to 8 nodes (small configuration). For the second configuration, the number of markers and the mesh size are both increased by a factor of 8 for the tests from 8 to 64 nodes (large configuration). Again, we used 12 MPI tasks per node on the CPU and 1 MPI task per node for the GPU tests.

Fig. 8 shows the speed-up on the two different set-ups with the number of nodes. The speed-up is calculated with respect to the time measured using the NL collision operator on one CPU (taken as baseline). The results for the large setup are further divided by 8, for direct comparison with the small one.

On the CPU an almost linear scalability is achieved from 1 to 64 nodes both for the NL and for the L case (the latter being slightly faster). It is worth reminding that in terms of MPI tasks, this corresponds to a range between 12 and 768 processes. Using 8 nodes the small and the large cases can be compared directly, providing a measure of weak scaling, which results to be almost ideal (i.e. ≈ 1), thanks to the limited impact of the MPI communication, which introduces a negligible overhead in the CPU version and in the range of considered MPI tasks.

On the GPU, the L and NL operators perform similarly, as already expected from the single node analysis, but here extended to multiple nodes. Differently from the CPU case, however, scalability tends to lose the linear trend already at 4 and at 32 nodes (and MPI tasks/GPUs) for the small and the large setups, respectively. This has a twofold origin. First, increasing the number of MPI tasks each GPU has a smaller and smaller dataset to work on, hence a progressively smaller workload. As a result, the efficiency of the accelerator tends to decrease, resulting in a loss of performance. This is confirmed from the weak scaling results at 8 MPI tasks: when the data size is increased, the workload per GPU increases and the computational efficiency and the speed-up grow considerably. Second, the smaller computing time compared to the CPU leads to a higher impact of the MPI communication overhead. As shown in Fig. 8, the collision operators scale well on both CPUs and GPUs. The overall cost of the NL operator on CPUs on a production case, like the one presented in section 6.3, is of the order of 20%. On GPUs, the relative cost of the collisions is smaller as other parts of the code have not been ported on GPUs yet. Once the code fully ported on GPUs, similar relative cost is expected on GPUs.

Conclusions

A nonlinear multispecies collision operator has been derived. Only two assumptions are performed. Firstly, the drift kinetic approximation is done because of the difficulty to treat FLR effects in the context of global gyrokinetic simulations. The second approximation consists in the truncation of the moment expansion used to compute the Rosenbluth potentials. This approximation is the consequence of the numerical cost required to compute accurately high order moments of the distribution function. The truncation proposed in this paper is consistent with the one performed in the standard neoclassical physics. The derivation presented in this article could be easily extended to include more moments to study highly nonlinear physics.

The Langevin formalism allows for an efficient implementation of this collision operator using a PIC approach eventough it introduces weight spreading. A weight spreading reduction scheme combined with a conservation procedure is implemented to solve this issue while conserving particle, total momentum and energy up to machine precision. An extensive verification of the collision operator has been performed with success. The collision operator implemented on ORB5 can be used both on CPUs and GPUs thanks to the use of MPI combined with OpenMP and OpenACC directives.

1 3 I 4 b+ s b + s 2 b I 1 b-+ s b I 2 b+ + I 3 b- G 01 b = 2N b v T b 1 3 I 6 b+ s b + s 2 b I 3 b-+ s b I 4 b+ + I 5 b- G 10 b = 2N b v T b v 3 1 5 I 6 b+ s 2 b + s 3 b I 1 b--I 4 b+ + s b I 3 b- G 11 b = 2N b v T b v 3 1 5 I 8 b+ s 2 b + s 3 b I 3 b--I 6 b+ + s b I 5 b-
The values of H lm b with l, m ∈ [0, 1] are deduced from Eq.( 10)

H 00 b = 2N b v T b I 2 b+ s b + I 1 b- H 01 b = 2N b v T b I 4 b+ s b + I 3 b- H 10 b = 2N b v T b v 3 I 4 b+ s 2 b + s b I 1 b- H 11 b = 2N b v T b v 3 I 6 b+ s 2 b + s b I 3 b-
In these expressions, the values of the terms I 2k b+ are computed using Eq.( 11) are extracted from Eq.( 12)

I 2 b+ = 1 2 [E b -s b E b ] I 4 b+ = 1 2 3 2 E b -E b 3 2 s b + s 3 b I 6 b+ = 1
I 1 b-= E b 2
I 3 b-= E b 2 1 + s 2 b I 5 b-= E b 1 + s 2 b + s 4 b 2
Appendix B: Exchange of momentum

In this part, the theoretical prediction of collisionnal exchange of momentum for a species a colliding on a species b is computed. Putting together Eq.( 23), ( 24), ( 25), ( 26

×         2 3 v T b v T a V a v T a - 2q a maNav 3 T a 1 -2 5 v T b v T a 2 s 2 b -s b E b + s 2 b E b + 2V 2 b 3v 2 T b s 4 b E b + 2 3 V b v T b 1 - V 2 a v 2 T a 1 -2 3 v T b v T a 2 s 2 b s 2 b E b + 4q b 5m b N b v 3 T b 1 - V 2 a v 2 T a 1 -2 3 v T b v T a
        2 3 v T b v T a V a v T a -A 0 + B 1 + 2V 2 b 3v 2 T b B 2 + 2 3 V b v T b B 1 - V 2 a v 2 T a B 1 -2 3 v T b v T a 2 B 2 -4 3 v T b v T a q a maNav 3 T a -A 0 + B 1 + 2V 2 b 3v 2 T b B 2 -2 5 v T b v T a 2 -A 1 + B 2 + 2V 2 b 3v 2 T b B 3 + 4q b 5m b N b v 3 T b 2 3 B 2 -B 1 - V 2 a v 2 T a 2 3 B 2 -B 1 -2 3 v T b v T a 2 2 3 B 3 -B 2        
Integrating by part, it is possible to show that

A 0 = 1 2 v T a v T b 2 B 0 A n = 1 2 v T a v T b 2 [2nA n-1 + B n ] B 0 = X ab B n = X 2n+1 ab n k=1 k - 1 2 ∀n > 0
Using these results, and the fact that 1 -

X 2 ab = v T b v T a 2
X 2 ab , one gets

∆P th ab 4 √ π νab ∆tN a m a v T b 1 + m a m b X 3 ab         1 3 v T b v T a V a v T a -v T a v T b 2 + V 2 b v 2 T b X 2 ab + 1 3 V b v T b 1 - V 2 a v 2 T a X 2 ab -2 5 v T b v T a q a maNav 3 T a X 2 ab -v T a v T b 2 + 5 3 V 2 b v 2 T b X 2 ab - 2q b 5m b N b v 3 T b v T b v T a 2 X 2 ab 1 -5 3 V 2 a v 2 T a X 2 ab        
This leads to Eq. [START_REF] Rosenbluth | Poloidal flow driven by ion-temperature-gradient turbulence in tokamaks[END_REF].

Appendix C: Exchange of energy

In this part, the theoretical prediction of collisionnal exchange of energy for a species a colliding on a species b is computed. Putting together Eq.( 28), [START_REF] Sauter | Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime[END_REF], (29) one gets

∆E th ab m a 2N a √ π νab ∆tv 2 T b ∞ 0 s 2 b ds b e -s 2 a × 1 -1 dξ 1 - V 2 a v 2 T a 1 - 2 3 s 2 a + 2ξs a V a v T a - 4q a m a N a v 3
                    -ma m b A 0 + 1 + ma m b B 1 + 1 + 2 ma m b (B 1 -A 0 ) + 2 1 + ma m b B 2 + 8 15 V b v T b v T b v T a q a maNav 3 T a v T b v T a 2 1 + 2 ma m b (B 2 -A 1 ) + 2 1 + ma m b B 3 + 8 15 q b m b N b v 3 T b v T b v T a V a v T a -
                   
Using the expressions for A n and B n , one gets

∆E th ab m a 4N a √ π νab ∆tv 2 T b X 3 ab ×                             1 2 1 -ma m b v T a v T b 2 + 1 2 V 2 b v 2 T b X 2 ab 2 3 -1 3 v T b v T a 2 + ma m b -1 2 V 2 a v 2 T a X 2 ab ma m b 2 3 -1 3 v T a v T b 2 + 1 + 1 3 V 2 a v 2 T a V 2 b v 2 T b X 4 ab 3 2 + ma m b v T b v T a 2 -1 + 3 2 ma m b + 1 3 V a v T a V b v T b v T b v T a X 2 ab 2 + ma m b -1 + 2 ma m b v T a v T b 2 + 2 5 q a maNav 3 T a V b v T b v T b v T a X 4
v T b v T a 2                            
which is equivalent to Eq.( 31).

e 2 a e 2 b 4π 2 0

 2 ma ln Λ, where ln Λ is the Coulomb logarithm. G b and H b are the Rosenbluth potentials associated with the species b G b

s 2 b

 2 is the local unshifted Maxwellian and s b = v v T b is the velocity normalized to the thermal speed v T b = 2T b m b . The normalized moments of the distribution function are given by

Figure 1 :Figure 2 :

 12 Figure 1: Isotropisation of the distribution function due to intra-species collisions. Comparison between the numerical results (red) and the theoretical prediction (black)

Figure 3 :

 3 Figure 3: Exchange rates of momentum (left) and energy (right) using the nonlinear collision operator. The numerical results (red) should be compared with the nonlinear theory (blue). The linear prediction (black) is shown for information.

Figure 4 :

 4 Figure 4: Comparison of the numerical heat diffusivity obtained with linear (black cross) and the nonlinear (red cross) collision operator with the theoretical prediction of Chang-Hinton (blue curve).

Figure 5 :

 5 Figure 5: Comparison of the numerical k neo obtained with the linear (black crosses) and the nonlinear (red crosses) collision operator with the theoretical prediction of Kim (green curve), Sauter (magenta curve) and Shaing (blue curve).

Figure 6 :

 6 Figure 6: Comparison of the numerical evolution of the zonal component of the electric potential φ 00 obtained with the nonlinear collision operator with the Hinton-Rosenbluth prediction of the collisional damping of the zonal flow residual.

Figure 7 :

 7 Figure 7: Wall clock time for 10 iterations of the collision operator running on one node of the Piz Daint system, with or without the GPU. The CPU only tests exploit all 12 cores available on the host. The cost of the four main routines inside the collision operator is shown, using both the Linear and the Non Linear solvers

Figure 8 :

 8 Figure8: Relative speed-up of the collision operator in two different setups: small (20M) and large (160M). Speed-ups are calculated with respect to the time measured using the NL collision operator on one CPU (taken as baseline). The numbers for the large setup are further divided by 8 for direct comparison with the small case. Ideal linear scaling (dashed blue line) is shown for reference.

+ s 7 b

 7 and the values of I 2k+1 b-

) one finds ∆P th ab 2 √ π νab ∆tN a m a v T b 1 + m a m b 1 - 1 dξ ∞ 0 ds b s 2 b e -s 2 a 1 - 1 ξth ab 4 √ π νab ∆tN a m a v T b 1

 211021141 n dξ = 1+(-1) n n+1 , one gets ∆P

2 b-

 2 and B n = ∞ 0 ds b s 2n b E b e a m a v T b 1 + m a m b ×

Using the fact that 1 - 1

 11 ξ n dξ = 1+(-1) n n+1 , one gets

3 - 2 +

 32 ma m b B 3 + 2 1 + ma m b B 4

  m b N b v 3 T b X 6 ab -3 + 4 ma m b + 4 + 3 ma m b

Table 1 :

 1 Coefficients for Kim prediction

	Collisionality regime	B	P	PS
	μ00	0.53 3.54 1.35
	K01	0.71 10.63 5.57
	μ11	1.39 11.52 6.90
	Kim's prediction gives, in the case of a single ion species plasma,

https://www.openmp.org

https://www.openacc.org

https://docs.nvidia.com/cuda/curand

Acknowledgements

The authors would like to thank B. Frei, X. Garbet and Y. Sarazin for fruitful discussions. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project s909. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was partly supported by the Swiss National Science Foundation.

Appendix A: Computation of the approximated Rosenbluth potentials

In this part, the values of H lm b and G lm b with l, m ∈ [0, 1] are given. Using Eq.( 9) and the fact that

Appendix D: Approximated solution for the isotropisation test

The goal of this appendix is to derive an approximated solution for Eq.( 45) in the special case of an initial condition independent of the velocity λ (s, t = 0) = λ 0 . We assume a solution of the form