Peter Donnel
email: peter.donnel@cea.fr

J.-B Fontana

J Cazabonne

L Villard

S Brunner

S Coda

J Decker

Y Peysson

Electron-cyclotron resonance heating and current drive source for flux-driven gyrokinetic simulations of tokamaks

published or not. The documents may come

Introduction

Electromagnetic waves can transfer their momentum and energy to magnetized plasmas. Efficient coupling between the wave and the plasma is possible in the presence of a resonance between the wave oscillation and the cyclotron motion of particles constituting the plasma. Electromagnetic waves with frequencies matching the Doppler shifted electron cyclotron (EC) frequency, or one of its harmonics, can be used to heat electrons, so-called Electron-Cyclotron Resonant Heating (ECRH), and/or to drive a current, so-called Electron-Cyclotron Current Drive (ECCD).

One of the main advantages of EC waves is their ability to couple with the plasma in a narrow region of space. Hence, EC waves can be used either as a pure heating/current source or to control MHD instabilities by accurately driving current at a precise position in the plasma.

It has been experimentally observed that, under certain conditions, suprathermal electrons generated by the EC power deposition are found in a region much larger than the region where the power deposition is expected using quasilinear drift-kinetic Fokker-Planck simulations [START_REF] Coda | Electron cyclotron current drive and electron dynamics in the TCV tokamak[END_REF][START_REF] Nikkola | Modelling of the electron cyclotron current drive experiments in the TCV tokamak[END_REF], solving beam propagation and power deposition in quiet plasmas. Two complementary improvements of models have therefore been proposed to reconcile these predictions with experimental observations. The first improvement is to take into account fluctuations of plasma density, induced by turbulence, along the path of the beam. These fluctuations, indeed, lead to the scattering of the beam and a, generally broadened, deposition area fluctuating in time. This effect has been studied with different approaches, starting from raytracing and drift-kinetic modeling [START_REF] Peysson | RF current drive and plasma fluctuations[END_REF][START_REF] Decker | Effect of density fluctuations on eccd in iter and tcv[END_REF] to wave-kinetic modeling and full-wave studies [START_REF] Weber | Scattering of diffracting beams of electron cyclotron waves by random density fluctuations in inhomogeneous plasmas[END_REF][START_REF] Snicker | The effect of density fluctuations on electron cyclotron beam broadening and implications for ITER[END_REF][START_REF] Köhn | Microwave beam broadening due to turbulent plasma density fluctuations within the limit of the born approximation and beyond[END_REF][START_REF] Chellaï | Millimeter-wave beam scattering by edgeplasma density fluctuations in TCV[END_REF][START_REF] Chellaï | Millimeter-wave beam scattering and induced broadening by plasma turbulence in the TCV tokamak[END_REF], concluding that it can significantly broaden the power deposition profile, depending on the beam incidence onto the resonant layer and on the profile of density fluctuations.

The second improvement is to take into account the spatial transport of suprathermal electrons generated by the EC power deposition after the wave absorption. This has been first implemented as an ad-hoc additional transport term in drift-kinetic Fokker-Planck simulations [START_REF] Nikkola | Modelling of the electron cyclotron current drive experiments in the TCV tokamak[END_REF][START_REF] Choi | Study of suprathermal electron dynamics during electron cyclotron current drive using hard x-ray measurements in the TCV tokamak[END_REF]. Comparison with experiments suggests that this transport is directly proportional to the diffusion in phase space induced by EC wave absorption. This ad-hoc transport allows to reconcile simulation predictions with experimental observations, suggesting that turbulence can efficiently transport suprathermal electrons generated by the EC power deposition. More generally, the possible synergies between EC power deposition and turbulence have not been studied so far. Indeed, to study these potential synergies, self-consistent simulations of turbulence in presence of ECRH/ECCD are required. To the best of our knowledge, this has never been done before, especially because of the difficulty to self-consistently treat turbulence and the resonant plasma/wave interaction. In this paper we derive an analytical operator modeling ECRH and ECCD. This operator can be used as a source term in any global gyrokinetic code, thus enabling the self-consistent simulation of turbulence in presence ECRH/ECCD. This is therefore a major step forward to study the possibility of synergistic effects between turbulence and ECRH/ECCD deposition.

In a recent paper [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF], a quasilinear operator modeling a pure ECRH source in the specific case of a beam propagation in the equatorial plane of the tokamak was derived and implemented in the flux-driven gyrokinetic code ORB5 [START_REF] Lanti | Orb5: A global electromagnetic gyrokinetic code using the pic approach in toroidal geometry[END_REF]. In this paper the generalization of this source to model both ECRH and ECCD, still assuming a propagation of the beam in the equatorial plane, is reported. This generalization of the source to include an ECCD component is an important step toward the perspective of comparing numerical results with experimental observations. Indeed, pure heating is a very particular case of the use of EC waves, and the presence of a current drive component is more relevant for applications foreseen in future large devices, such as MHD mode mitigation. Moreover, observation of suprathermal electrons using Hard X-Ray spectrometry is easier in ECCD configuration, as the resonance takes place farther away from the bulk distribution, exciting electrons with higher velocities and increasing the photon count rate at higher energies.

The derivation of the improved source term is presented in section 2. The numerical treatment of this quasi-linear operator is described in section 3. In section 4, numerical tests of this new source term are reported.

Derivation of an ECRH/ECCD source for a beam propagating on the midplane

The gyrokinetic equation can be symbolically written

∂ Fs ∂t + Fs , Hs = s ′ C Fs , Fs ′ + Q EC Fs + S Fs ( 1
)
where Fs is the distribution function of gyro-centers of the species s. The second term of the left hand side is the Poisson bracket of Fs with the Hamiltonian Hs . It corresponds to particle drifts and contains nonlinear terms, leading to turbulence. The first term on the right hand side is the collision operator, which is a key ingredient for the saturation of the EC beam deposition at large power, as it smooths out small scale fluctuations in the velocity space that are generated by the plasma/wave coupling. The term Q EC Fs is the source term representing the ECRH/ECCD deposition which is derived in this paper. Note that it will be non zero only for electrons. Finally, the last term in Eq.1 represents other source and sink terms. Eq.1 needs to be solved for each species self-consistently with Poisson and Ampere equations to evolve the electromagnetic field. The Larmor radius of electrons is small compared with the ion Larmor radius. As Ion Temperature Gradient (ITG) and Trapped Electron Mode (TEM) turbulence develop at typical space scales of the order of the ion Larmor radius, the drift kinetic approximation is regularly used for electrons in gyrokinetic codes. This implies that no distinction is made between the gyro-center and the particle electron distributions for the electrons. Finite Larmor radius effects will however be kept in the realistic source term Q EC , as they are leading order terms of this operator.

The derivation of a quasilinear operator representing the coupling of an EC beam with a plasma has been presented in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]. We give here the main results of the derivation. The interested reader is invited to read [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF] for details of the computation. The operator describing the resonant coupling of an electromagnetic wave with the cyclotron motion of electrons reads

Q EC = ∞ n=-∞ Q n , (2)
where n stands for the harmonic number and

Q n = ∇ v • [D n • ∇ v F e ] , (3)
with F e the particle electron distribution function and ∇ v the standard nabla operator in velocity space. The diffusion matrix D n is defined as

D n = D n s n s T n . (4)
In Eq.4, the vector s n , which represents the direction of diffusion in the velocity space, is defined as

s n = s ⊥ n s ∥ n = -nΩe ω b γ λ v ∥ , v ⊥ v ⊥ c , (5)
where Ω e = eB/m e is the non-relativistic electron cyclotron frequency with e the elementary charge and m e the electron mass, ω b is the beam frequency, γ is the Lorentz factor, v ⊥ is the component of the velocity perpendicular to the magnetic field B, c is the light speed and the function λ v ∥ , v ⊥ is defined as

λ v ∥ , v ⊥ = 1 - nΩ e γω b c v ∥ , (6)
with v ∥ the component of the velocity parallel to the magnetic field. For electrons, the main interaction is with negative harmonics. Indeed, with the conventions given in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF], the resonance condition reads ω b = k ∥ v ∥ -nΩ e /γ, with k ∥ the parallel component of the wave vector. As ω b and Ω e are positive, n needs to be negative to fulfill the resonance condition. An approximate version of the resonant diffusion coefficient has been derived in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]. It reads

D -n ≃ D-n = πe 2 N (θ 0 ) 2m 2 e ω b c v ∥ |E| 2 sin (θ res ) exp -θres-θ0 σ 2 N 2 (θ res ) π 0 sin (θ) exp -θ-θ0 σ 2 N 2 (θ) dθ Θn k,res 2 . (7)
Eq.7 is rather general and relies on a few assumptions. The first one is the use of the cold plasma dispersion relation to express the plasma refractive index N . This assumption is used to express the direction of the electric field, which enters in the definition of Θn k,res . The derivation of Eq.7 also assumes that the electric field intensity possesses a Gaussian shape of width W 0 in the direction perpendicular to its propagation. This leads to a Gaussian shape of the Fourier transform of the electric field in the direction perpendicular to its mean wavevector k 0 . The consequence of this assumption is found in the exponential functions in Eq.7. Indeed, θ 0 is the angle between the magnetic field and the mean wavevector of the beam k 0 , and σ = 1/ (k 0 W 0 ) is the width in the distribution of angles resulting from the finite width of the beam W 0 . θ res is the resonant angle corresponding to the matching of the Doppler shifted beam frequency with the cyclotron frequency or one of its harmonics. In Eq.7, E is the electric field of the beam and Θn k,res is defined as

Θn k,res = [1 + C 1 (θ res )] J n+1 (ρ res ) + [1 -C 1 (θ res )] J n-1 (ρ res ) -2C 2 (θ res ) v ∥ v ⊥ J n (ρ res ) 2 1 + C 1 (θ res ) 2 + C 2 (θ res ) 2 , (8)
where J n is the n-th Bessel function of first kind of argument ρres = sin (θ res ) N (θ res ) ω b v ⊥ γ/ (Ω e c) is a normalized Larmor radius, and

C 1 (θ) = T S -N 2 (θ) , (9)
C 2 (θ) = N 2 (θ) cos (θ) sin (θ) P -N 2 (θ) sin 2 (θ) . ( 10
)
Note that if one neglects finite Larmor radius (ρ res = 0), all Bessel functions are equal to zero except J 0 . The consequence is that no absorption is possible except for the fundamental resonance. This is consistent with the fact that ECRH/ECCD deposition relies on a resonance between the electromagnetic wave and the cyclotron motion. The coefficients P, S and T are the "Stix coefficients" which are related to the plasma dispersion relation. They are defined as

P = 1 -(ω p /ω b ) 2 ,
where The quasilinear operator Eq.7, derived in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF], requires as an input some characteristic quantities of the beam (beam waist W 0 , mean wavevector of the beam k 0 as well as the power carried by the beam along its path, which is proportional to the square of the electric field of the beam integrated over a surface perpendicular to the beam propagation direction). These quantities can in principle be evaluated by coupling a code computing the beam quantities (i.e. a ray-tracing code) with a global gyrokinetic code. But the coupling of such codes is numerically challenging. To avoid the coupling of another code with the gyrokinetic code, the particular case of a propagation of the beam in the equatorial plane of an axisymmetric plasma is assumed. This assumption can be fulfilled or not depending on the experimental setup (position of the launcher and direction of injection). In [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF], a pure ECRH case is considered (θ 0 = π/2), leading to an one-dimensional beam propagation (along the major radius). However, the pure ECRH case is really specific and a generalization of the model to treat both ECRH and ECCD cases is required before envisaging a comparison with experiments. In this paper, the model is generalized by assuming an axisymmetric plasma. This makes the beam propagation 2D (in the equatorial plane), but still analytical, see section 2.2. The breaking of the axisymmetry, for instance by magnetic ripple or resonant magnetic perturbations is experimentally possible. But axisymmetry assumption, in addition to considering a beam of constant width W 0 and to neglecting the poloidal magnetic field (i.e. B θ ≪ B φ ), make it easier to have analytical models for the beam propagation. The possibility to include at the same time ECRH and ECCD requires a model, much more complex than the one for a pure ECRH case, to compute the power carried by the beam. Such a model has been derived by [START_REF] Albajar | Electron-cyclotron absorption in hightemperature plasmas: quasi-exact analytical evaluation and comparative numerical analysis[END_REF] and is presented in the subsection 2.4.

ω p = N e
The rest of this section is organized as follows. In subsection 2.1, an analytical evaluation of the resonant angle θ res is performed. It allows an efficient implementation compared with the binary search used in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]. A similar computation, described in subsection 2.2, allows the computation of the beam path in the equatorial plane, and so of the mean angle θ 0 between the beam wave-vector and the plasma magnetic field. Once θ res and θ 0 are evaluated, the remaining difficulty consists of expressing the spatial distribution of the beam electric field E. To do so, the link between the beam electric field E and the power crossing a surface of constant major radius P (R) is established in subsection 2.3. An analytical model, described in subsection 2.4, is used to evaluate the beam power along its path.

Evaluation of the resonant angle

The resonant angle θ res , appearing in Eq.7, is the solution of the resonance condition ω b = k ∥ v ∥ -nΩ e /γ which can be rewritten as [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]

N (θ res ) cos (θ res ) = λ v ∥ , v ⊥ . ( 11
)
The cold plasma dispersion relation N (θ) reads

N 2 O/X (θ) = (RL + SP) tan 2 θ + 2SP ± G θ 2 S tan 2 θ + P , (12)
where O and X stand for the O-mode polarization and X-mode polarization respectively and

G θ = sign (ω b -Ω e ) (SP -RL) 2 tan 4 θ + P 2 (L -R) 2 tan 2 θ + 1 .( 13
)
Note that the definition of G θ has been generalized compared with the one in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF] to include the unusual case ω b < Ω e , corresponding to a configuration of injection from the high field side and a coupling at the fondamental resonance. Defining y = cos (2θ res ) and by squaring Eq.11, it is possible to show that y is the solution of (RL + SP) 1 -

y 2 + 2SP (1 + y) 2 -4λ 2 [S (1 -y) + P (1 + y)] = ± (1 + y) (SP -RL) 2 (1 -y) 2 + 8P 2 T 2 (1 + y). (14)
By squaring Eq.14 and using the definitions of the Stix coefficients, a long but straightforward computation allows one to compute the angle of resonance for the X-mode

θ X res = 1 2 arccos x + λ 2 if λ ≥ 0 = π - 1 2 arccos x + λ 2 if λ < 0 ( 15
)
and the O-mode

θ O res = 1 2 arccos x -λ 2 if λ ≥ 0 = π - 1 2 arccos x -λ 2 if λ < 0 (16) with x ± (y) = P 2 (2y -P) + Ω e ω b 2 P (1 -y) 2 -y 2 ± Ω e ω b (1 -P) y Ω e ω b 2 (1 -y) 2 + 4Py    / P P 2 - Ω e ω b 2 + y Ω e ω b 2 (P -1) . (17)
2.2. Model for the spatial dependence of θ 0

The spatial dependence of the mean angle θ 0 (R), appearing in Eq.7, also needs to be evaluated. To do so, we use the general expression which is implemented in the ray-tracing code C3PO [START_REF] Peysson | A versatile ray-tracing code for studying rf wave propagation in toroidal magnetized plasmas[END_REF]. The toroidal mode number of the beam is defined as ñ = Re φ • k 0 where e φ is the toroidal direction and k 0 is the mean wave vector associated with the beam. For an axisymmetric magnetic configuration, ñ is invariant along the ray trajectory [START_REF] Peysson | A versatile ray-tracing code for studying rf wave propagation in toroidal magnetized plasmas[END_REF]. By neglecting the poloidal component of the magnetic field, we obtain

N (θ 0 ) cos (θ 0 ) = R in R N (θ 0,in ) cos (θ 0,in ) , (18)
where θ 0,in = θ 0 (R in ) is the injection angle at the injection major radius R in , where the beam is entering the plasma. Eq.18 has the same structure as Eq.11. Therefore the mean angle is given for a X-mode beam by

θ X 0 = 1 2 arccos x + R in R N (θ 0,in ) cos (θ 0,in ) 2 if θ 0,in ≤ π 2 = π - 1 2 arccos x + R in R N (θ 0,in ) cos (θ 0,in ) 2 if θ 0,in > π 2 (19)
and for a O-mode beam by

θ O 0 = 1 2 arccos x - R in R N (θ 0,in ) cos (θ 0,in ) 2 if θ 0,in ≤ π 2 = π - 1 2 arccos x - R in R N (θ 0,in ) cos (θ 0,in ) 2 if θ 0,in > π 2 (20)
with x ± defined by Eq.17. It is noteworthy that, for a perpendicular injection θ 0,in = π/2, the mean angle is constant θ 0 (R) = π/2 , ∀R. This specific case was used in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF] for the pure ECRH case.

Spatial dependence of the electric field in the case of a beam propagating in the equatorial plane

The power crossing a surface with a major radius R reads

P (R) = R 2π 0 dφ ∞ -∞ dZ ϵ 0 |E| 2 2 v g • e R = 1 2 Rϵ 0 |v g | sin (θ 0 ) 2π 0 dφ ∞ -∞ dZ |E| 2 , (21)
Note that the electric field is implicitly taken at the major radius R. Consistently with the assumption used to derive Eq.7, we assume that the beam has a Gaussian distribution of width W 0 perpendicular to its propagation direction

|E| 2 = |E 00 (s)| 2 exp - X W 0 2 exp - Z W 0 2 , ( 22
)
where X is the coordinate in the direction perpendicular to the beam propagation in the equatorial plane and Z is the coordinate in the vertical direction. s stands for the arc length along the beam path from the injection location to the local position and |E 00 (s)| 2 is the maximum amplitude of the electric field at this arc length. When integrating the electric field on a surface of constant major radius R, s and X are such that the beam is evaluated on the surface of constant major radius R as illustrated by Fig. 1. Because we assume a propagation at the equatorial plane, and because the width of the beam is small compared to the size of the plasma, the integration in the vertical direction is trivial and gives

P (R) = √ π 2 RW 0 ϵ 0 |v g | sin (θ 0 ) 2π 0 dφ |E 00 (s)| 2 exp - X W 0 2 . ( 23
)
The integration in the direction φ is more challenging. In reality the EC beam is toroidally localised. However, assuming that the power deposited is quickly spread along the toroidal direction, we shall assume that the deposition is axisymmetric. This will have the advantage, in the context of a PIC code, to increase the sampling statistics. We then define a mean value for the electric field on the equatorial plane Ē (R) such that

Ē2 = 1 2π 2π 0 dφ |E 00 (s)| 2 exp - X W 0 2 , (24)
and so

P (R) = π 3/2 RW 0 ϵ 0 |v g | sin (θ 0 ) Ē2 . ( 25
)
We also define Ê such that

|E| 2 = Ê2 exp - Z W 0 2 . ( 26
)
We assume that W 0 ≪ R, and define a new coordinate X ′ which corresponds to the coordinate on the surface R (see Fig. 1). By noting that X ′ sin θ 0 = X and by imposing that Ê be independent of θ 0 (in the absence of absorption), one gets Ē = sin (θ 0 ) Ê and so

|E| 2 = P (R) π 3/2 RW 0 ϵ 0 |v g | sin 3 (θ 0 ) exp - Z W 0 2 . ( 27
)
Eq.27 links the amplitude of the electric field to the power crossing the surface of constant major radius P (R) and the group velocity v g , whose expression is derived in the appendix A.

Model of absorption for an oblique propagation

To evaluate the evolution of the beam power along its path, we use the quasi-exact analytical evaluation derived in [START_REF] Albajar | Electron-cyclotron absorption in hightemperature plasmas: quasi-exact analytical evaluation and comparative numerical analysis[END_REF]. The power is then related to the optical thickness of the plasma τ

P (R) = P in exp [-τ (R)] , (28)
where P in = P (R in ) is the injected power and τ is defined as the integral of the absorption coefficient α along the beam path τ

= L 0 αds = R Rin α |sin (θ 0 ) dR ′ |. The absorption coefficient is given by α = n≥n0 α n , (29)
with

α n = ω 2 p cΩ e π 2 µ 2 n 0 K 2 (µ) P n exp   - nµ n 0 1 -N 2 ∥   n n 0 2 -1, (30)
where

n 0 = ω b Ωe 1 -N 2 ∥ , µ = mec 2
Te , K 2 is the modified Bessel function of the second kind of order two, N ⊥ = N (θ) sin θ and N ∥ = N (θ) |cos θ|.

In practice, for all cases of interest the limit µ ≫ 1 holds. Within this limit, the sum in Eq.30 is very well approximated by its first non zero component (n = ⌈n 0 ⌉) which then corresponds to the harmonic of the beam considered. Moreover, one can use the fact that lim

x→∞ [K 2 (x) exp (x)] = π 2x to approximate the coefficient α n α n ≃ αn = ω 2 p cΩ e π 2 F n (µ) n 0 n n 0 2 -1, (31)
where

F n (µ) = µ 5/2 P n exp   µ   1 - n n 0 1 -N 2 ∥     . (32)
The function P n is defined as

P n = π (2n + 1)! (2 n n!) 2 nΩ e ω b N ⊥ 2 A J n+1/2 (z n ) 2 x n + B J n+3/2 (z n ) 2 x n , (33)
where

A (•) = |A xz | 2 + |ê y | 2 • +ℜ iA xz ê * y x n n ∂• ∂x n - x n n 2 n n + 1 |ê y | 2 • - ∂ 2 • ∂y 2 n +   x n n 1 -N 2 ∥   2 |ê z | 2 ∂ 2 • ∂y 2 n - x n 2ℜ (A xz ê * z ) + ℜ iê * y êz xn n ∂ ∂xn n 1 -N 2 ∥ ∂• ∂y n , (34)
and

B (•) = x n n 2 2n + 3 (n + 1) (n + 2) |ê y | 2 • - ∂ 2 • ∂y 2 n , (35)
In these expressions, A xz is defined as

A xz = êx + N ⊥ N ∥ 1 -N 2 ∥ êz , (36)
and the quantities x n , y n are respectively defined as

x n = ω b Ω e N ⊥ n n 0 2 -1, (37)
y n = µN ∥ 1 -N 2 ∥ n n 0 2 -1, (38)
and the complex argument of the Bessel functions z n is defined as

z n = 1 2 4x 2 n -y 2 n + iy n . ( 39
)
The derivatives of the function

|Jm+1/2(zn)| 2 xn
with respect to x n and y n , which are necessary to evaluate Eq.34 and Eq.35, are derived analytically in appendix B.

The use of the cold plasma limit allows one to link the different components of the normalized wave electric field ê = E k / |S| /ϵ 0 c. Indeed one has

i êy êx = i E k,y E k,x = T S -N 2 (40) and êz êx = E k,z E k,x = - N ∥ N ⊥ P -N 2 ⊥ (41)
In the cold plasma limit, the electromagnetic energy flux S is well approximated by the Poynting vector and leads us to express

|ê y | 2 = 1 N (θ) a 2 (θ) + b 2 (θ) (42)
with

a (θ) =   1 + P N ∥ ω b ω 2 p -ω 2 b -Ω 2 e 1 -N 2 Ω e ω 2 p (P -N 2 ⊥ ) 2   sin θ (43) and b (θ) = 1 + P P -N 2 ⊥ ω b ω 2 p -ω 2 b -Ω 2 e 1 -N 2 Ω e ω 2 p 2 |cos θ| (44)
For the computation, one can assume without loss of generality that êy is real and positive. Then using Eq.40 and Eq.41, one obtains that êx and êz are purely imaginary quantities.

For the specific case of an O-mode close to perpendicular propagation, one has N ⊥ ≃ N ≃ P and then equations 41, 43 and 44 are not properly defined. For this specific case, one can approximate êx = êy = 0 and êz = P -1/2 .

Note that for oblique propagation, with low temperature, the numerical evaluation of Eq.32 is challenging as it involves the multiplication of very large numbers with very small ones. To circumvent this difficulty, one can use the logarithm to perform the multiplication and then exponentiate the result.

Numerical implementation of the source term

In the previous section, the model for the quasi-linear operator representing the ECRH/ECCD source has been derived. In this section, the numerical implementation of this source term in the global gyrokinetic code ORB5 is detailed. In subsection 3.1, the drag term is analytically derived, avoiding a difficult numerical evaluation of this term. The markers with a low parallel velocity are displaying a diverging drag term. To avoid this divergence, the numerical treatment of these markers is different from the one of the other markers. This treatment is detailed in subsection 3.2. Finally, the adaptation of the source term to the usage of "heavy electrons", often used in gyrokinetic code to lighten the numerical cost of simulations, is detailed in subsection 3.3.

Analytical derivation of the drag term

The quasi-linear operator Q n (Eq.3) can be written in the form of a Fokker-Planck operator

Q n = ∇ v • [∇ v • (D n F ) -Γ n F ] , (45)
where the drag force is defined as

Γ n = ∇ v • D n . (46)
Eq.45 can be efficiently implemented in a particle-in-cell code by finding the equivalent stochastic process thanks to Ito's lemma. This stochastic process is then discretized using a Euler-Maruyama scheme. More details about this can be found in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]. In this implementation, it is necessary to compute the drag term Eq.46 which takes the form

Γ n = ∇ v • D n = Γ ⊥ n ê⊥ + Γ ∥ n ê∥ , (47)
with

Γ ⊥ n = 1 v ⊥ ∂ v ⊥ Dn s ⊥ n s ⊥ n ∂v ⊥ + ∂ Dn s ⊥ n s ∥ n ∂v ∥ , (48)
Γ ∥ n = 1 v ⊥ ∂ v ⊥ Dn s ⊥ n s ∥ n ∂v ⊥ + ∂ Dn s ∥ n s ∥ n ∂v ∥ . (49)
The numerical evaluation of the derivatives in velocity space with finite differences is challenging as the resonant diffusion coefficient Dn (Eq.7) displays thin layers in velocity space, especially in the ECCD case. An alternative is to compute analytically the velocity derivatives appearing in Eq.48, 49. Without any simplification, the analytical derivation is challenging and would be numerically expensive. To alleviate this difficulty, some approximations are performed. The first one is to consider the limit of small wavelength as compared to the beam waist size σ ≪ 1 ⇔ N (θ 0 ) W 0 ω b ≫ c. For tokamaks with a relatively weak magnetic field, the small wavelength approximation is only marginally valid (i.e. for TCV W 0 ≥ 2cm, ω b = 82 GHz, N (θ 0 ) ∼ 1, σ ⪅ 0.2). For devices with stronger magnetic fields, this assumption is legitimately valid (i.e. for ITER

W 0 ≥ 2cm, ω b = 170 GHz, N (θ 0 ) ∼ 1, σ ⪅ 0.1).
The resonant diffusion coefficient Dn , Eq.7, possesses three dependencies with respect to θ res . The first one is in the dispersion relation N , the second in the function Θn k and the last one in the exponential term exp -(∆θ/σ) 2 . For N and Θn k , the difference between the function evaluated at θ res and at θ 0 is proportional to ∆θ ∼ σ ≪ 1. On the other hand, for the exponential term, one gets exp -(∆θ/σ)

2 -1 ∼ exp (-1) -1 which is a zero order term. Therefore, in the small wavelength limit, at leading order, the resonant diffusion coefficient reads

Dσ≪1 -n = √ πe 2 N (θ 0 ) 2m 2 e ω b σ c v ∥ |E 0 | 2 exp - ∆θ σ 2 Θn k,0 2 , ( 50
)
where

Θn k,0 = [1 + C 1 (θ 0 )] J n+1 (ρ 0 ) + [1 -C 1 (θ 0 )] J n-1 (ρ 0 ) -2C 2 (θ 0 ) v ∥ v ⊥ J n (ρ 0 ) 2 1 + C 1 (θ 0 ) 2 + C 2 (θ 0 ) 2 , ( 51
) with ρ0 = sin (θ 0 ) N (θ 0 ) ω b v ⊥ γ/ (Ω e c).
The velocity derivatives of the resonant diffusion coefficient in the limit σ ≪ 1 can then be computed

∂ Dσ≪1 -n ∂v ⊥ = 2 Dσ≪1 -n 1 Θ n k,0 ∂Θ n k,0 ∂v ⊥ - ∆θ σ 2 ∂∆θ ∂v ⊥ , (52)
∂ Dσ≪1 -n ∂v ∥ = Dσ≪1 -n 2 Θ n k,0 ∂Θ n k,0 ∂v ∥ - 1 v ∥ -2 ∆θ σ 2 ∂∆θ ∂v ∥ . ( 53
)
The derivatives of Θ n k,0 with respect to the velocity components are given in the appendix C. The velocity derivatives of ∆θ can be computed using

∂∆θ ∂v ∥,⊥ = - x ′ ± λ 2 |λ| 1 -x 2 ± (λ 2 ) ∂λ ∂v ∥,⊥ , (54)
where x ′ ± is the derivative of Eq.17, which is straightforward to compute, and the velocity derivatives of λ are given by

∂λ ∂v ∥ = nΩ e cω b γ - λ v ∥ , ( 55
)
∂λ ∂v ⊥ = nΩ e cω b v ⊥ v ∥ γ. (56)

Treatment of markers with a low parallel velocity

The drag coefficients diverge for markers with a low parallel velocity. This is due to the fact that the width of the resonant area in the perpendicular velocity direction goes to zero when the parallel velocity goes to zero. Therefore the drag is becoming singular for these low parallel velocities.

To circumvent this issue, markers with a low parallel velocity ( v ∥ < v max ) are treated differently. For these markers, the distribution function is assumed to be an unshifted Maxwellian F M e . This assumption is motivated by two considerations. First, in absence of a localized source, and even in presence of turbulence, distribution functions of electrons in the core are close to an unshifted Maxwellian. Second, the effective electron-electron collision frequency goes like

∼ (v T /v) 2 for v ≪ v T and ∼ (v T /v) 3 for v ≫ v T ,

whereas the effective electron-ion collision frequency goes like

∼ (v T /v) 3 [START_REF] Estève | A multi-species collisional operator for full-f gyrokinetics[END_REF]. For this reason, the low velocity part of the distribution function is often close to a Maxwellian distribution. The assumption that the electron distribution function is an unshifted Maxwellian allows to rewrite the quasi-linear operator Q n as

Q M n = ∇ v • [D n • ∇ v F M e ] = -∇ v • D n • v v 2 T F M e . (57)
By analogy with Eq.45 we impose for markers with low parallel velocities (

v ∥ < v max ) the changes D n → 0 and Γ n → Γ M n where Γ M n = D n • v v 2 T , (58)
where v T = T e /m e is the thermal velocity of electrons. This simple change of numerical treatment allows avoiding the singularity for low parallel velocities while keeping the same operator (provided that the distribution function is close to a Maxwellian). The choice for v max is a priori not settled. A numerical study led us to fix the choice of this parameter to v max = 0.1v T .

Usage of an artificial mass ratio between protons and electrons

Gyrokinetic simulations with a real mass ratio between protons and electrons are numerically challenging due to the time and spatial scales separation between the electron and ion dynamics. For this reason, simulations with an artificially reduced mass ratio between protons and electrons are often performed in the gyrokinetic community.

The electron mass appears in the plasma frequency ω p , in the cyclotron frequency Ω ce , in the prefactor of the resonant diffusion coefficient of Eq.7 and in the definition of the thermal velocity which appears for a term of the form

v ∥ /c = v ∥ /v T • (v T /c) or v ⊥ /c = (v ⊥ /v T ) • (v T /c).
In all cases, the real mass of electrons needs to be used in the definitions of the frequencies to ensure the right resonance condition and the proper dispersion TeeB between the real and artificial electron masses.

Numerical test of ECRH/ECCD source

In this section, testing of the source term described in the previous sections is detailed. In subsection 4.1, an analytical prediction of the minimum and the maximum major radii of absorption is given. The numerical evaluation of the power deposited in the plasma is detailed in subsection 4.2. In subsection 4.3, the evolution of the electron distribution function in the presence of the source term, but without other physical operators, is presented. Finally, subsection 4.4 presents a comparison of ORB5 with the code LUKE when activating all parts of the code, except the turbulence.

Analytical prediction of the deposition area in real space

The resonant condition

ω b = k ∥ v ∥ + n Ωe γ is equivalent to v ∥ -v ∥ ∆v ∥ 2 + v ⊥ ∆v ⊥ 2 = 1, (59)
which is the equation of an ellipse with an average parallel velocity

v ∥ c = ω b k ∥ c k ∥ c 2 + (nΩ e ) 2 , (60)
a semi-major axis in the perpendicular direction ∆v ⊥ given by

∆v ⊥ c 2 = k ∥ c 2 + (nΩ e ) 2 -ω 2 b k ∥ c 2 + (nΩ e ) 2 , (61)
and a semi-minor axis in the parallel direction ∆v ∥ given by

∆v ∥ c 2 = ∆v ⊥ c 2 (nΩ e ) 2 k ∥ c 2 + (nΩ e ) 2 . (62)
By neglecting the poloidal magnetic field, one has Ω e = Ω e,0 R0 R , where Ω e,0 = eB0 me is the electron cyclotron frequency on the magnetic axis. Eq.18 implies k ∥ R = k ∥,in R in . Using these relations allows one to express the characteristics of the ellipse as functions of the major radius R. From these expressions, one can deduce that the resonance between the wave and the plasma is possible only if R ≤ R res max with

R res max = k ∥,in cR in 2 + (nΩ e,0 R 0 ) 2 ω b . ( 63
)
For all major radii below R res max the resonance between the wave and the plasma is a priori possible. But for an efficient coupling between the wave and the plasma, particles fulfilling the resonance condition Eq.59 are needed. If one assumes the distribution function of electrons to be an unshifted Maxwellian of temperature T e , a reasonable coupling between the wave and the electrons is possible only if

-3v res T ≤ v ∥ -∆v ∥ ≤ 3v res T ( 64
)
where v res T is the thermal velocity of electrons at the position of maximum absorption. As v ∥ -∆v ∥ is an increasing function of the major radius, the condition Eq.64 allows determining a minimal major radius R eff min and maximal major radius R eff max for an efficient absorption, which read

R eff max/min = ±3v res T k ∥,in R in + nΩ e,0 R 0 1 -3 v res T c 2 ω b . (65)
Note that the computation of R eff min and R eff max is valid only for an unshifted Maxwellian. For a shifted Maxwellian, or a distribution function with a non negligible deviation to a Maxwellian, the absorption boundaries can be shifted. For most of the numerical tests presented in this article, unshifted Maxwellians are a good proxy of the actual distribution function. Note that this is no longer true for the cases with a high power deposition as presented in subsection 4.4. Furthermore, it is noteworthy that the inner bound R eff min is a theoretical limit for efficient absorption of the beam. For cases with large density and temperature, most of the beam power is absorbed before reaching this inner major radius.

Computation of the power deposition

The source is numerically implemented in ORB5 [START_REF] Lanti | Orb5: A global electromagnetic gyrokinetic code using the pic approach in toroidal geometry[END_REF] via a Langevin equation as in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF]. The power deposited on a given marker i reads

∆ v 2 i ∆t = 2 n v ∥,i Γ ∥ n,i + v ⊥,i Γ ⊥ n,i + D n,i s ∥ n,i 2 + s ⊥ n,i 2 + ∆t n Γ ∥ n,i 2 + Γ ⊥ n,i 2 , ( 66
)
where the brackets ⟨.⟩ stand for the ensemble average. By multiplying this quantity by the weight of the marker times m e /2, and by summing over the markers in a spatial bin, one gains access to the ensemble average power deposition on this bin. Note that Eq.66 allows putting a quantitative limit on the maximum time step which is allowed for a given source amplitude, as the term proportional to ∆t in Eq.66 should always be small compared with the term independent of the time step.

Test with only the source

In this section the electron distribution function is evolved only by the source operator dF e /dt = Q n , where Q n is the quasi-linear operator described above in Eq.3. In all the tests presented in this section and the next one, a beam with an X polarization and interacting with the second harmonic of the cyclotron frequency ('X2' in short) is simulated and an analytical axisymmetric magnetic geometry with B = B0R0 R ϵ q e χ + e φ is used. In this expression, χ is the poloidal angle, φ is the toroidal angle, ϵ = r/R 0 is the inverse aspect ratio and q = q √ 1 -ϵ 2 , where q is the safety factor. In all the tests presented in this section and the next one, typical TCV values have been used (B 0 = 1.4 T, R 0 = 0.89 m, a = 0.25 m).

In this subsection only, to simplify the test of the source, the safety factor q = 10 is constant and the plasma density N e and temperature T e are chosen to have no spatial dependence. The density and temperature are scanned around reference values

N ref e = 2•10 18 m -3 , T ref e = 1
.17 keV, to test the source behavior. The reference density has been chosen quite low to have approximately half the power absorbed. This case is the most stringent test for the model. Indeed, with a larger density, the beam would be fully absorbed and the sensitivity to the numerical parameters would be weaker than in the case of half absorption. On the other hand, a really weak absorption would not numerically test the model of absorption described in section 2.4.

The frequency of the beam f beam = 78 GHz is chosen to have the maximum of absorption near the magnetic axis. It allows one to simulate only a restricted volume of the plasma r < 0.2a and so to have a large density of markers with a limited cost for the simulation. For the results shown in Fig. 2, N p = 2 • 10 8 markers were used. With this choice, the simulation results are numerically converged as discussed in Fig. 3. The power of the beam is chosen to be small enough so that the distribution function stays close to a Maxwellian, allowing a direct comparison between ORB5 results and Albajar's prediction, which assumes a Maxwellian distribution for electrons.

The profiles of absorption for a pure ECRH and a mixed ECRH/ECCD case with different densities and temperatures are shown in Fig. 2. In all cases, excellent agreement is found between the numerical results obtained with ORB5 and the analytical prediction [START_REF] Albajar | Electron-cyclotron absorption in hightemperature plasmas: quasi-exact analytical evaluation and comparative numerical analysis[END_REF]. It validates the derivation of the model and its implementation in ORB5.

A convergence test is performed in Fig. 3. Fast convergence is obtained both in a pure ECRH and a mixed ECRH/ECCD case. More markers are needed to converge in the pure ECRH case as the deposition is narrower in real space, but, in both cases, a good quantitative agreement with the analytical prediction is found, even with a low resolution (which would be under-resolved for a simulation in the presence of turbulence).

Benchmark against drift-kinetic simulations

In this section, the code ORB5 is compared with the code LUKE, which is a 3D relativistic bounce-average drift-kinetic Fokker-Planck solver for the electron distribution function F e (ρ, v ∥ , v ⊥ ) [START_REF] Peysson | Numerical simulations of the radio-frequency-driven toroidal current in tokamaks[END_REF], coupled to the ray-tracing code C3PO [START_REF] Peysson | A versatile ray-tracing code for studying rf wave propagation in toroidal magnetized plasmas[END_REF], and is mainly used to study suprathermal electron physics, especially Electron-Cyclotron and Lower-Hybrid current-drive problems. Physics operators included in LUKE are collision, Ohmic electric field, quasilinear wave-plasma interaction and additional adhoc radial transport of fast electrons, in the form of bounce-average Fokker-Planck operators. The basis of the quasilinear treatment of the wave-plasma interaction is similar in LUKE and ORB5 [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF].

In this subsection, all parts of the ORB5 code are activated (source, collisions, Vlasov and quasi-neutrality). The non axisymmetric modes of the electric potential are filtered out at each time step to remove the turbulence as LUKE cannot simulate it. As a consequence, the ad-hoc radial transport of fast electrons is not activated in LUKE, as it is supposed to model turbulent transport. Within this framework, both codes solve similar problems. Some differences remain between the two codes. Firstly, LUKE solves a bounce-average problem whereas ORB5 solves the 5D problem. The collision operators are also different as ORB5 possesses a nonlinear collision operator, but without relativistic corrections [START_REF] Donnel | Moment approach of the multi-species nonlinear coulomb collision operator adapted to particle-in-cell codes[END_REF], whereas LUKE uses a linearized collision operator with a relativistic model for electron-electron collisions [START_REF] Braams | Conductivity of a relativistic plasma[END_REF], and non-Maxwellian ion background model for electron-ion collisions, as described in [START_REF] Karney | Fokker-planck and quasilinear codes[END_REF]. collisions are not treated in LUKE, as ions are not simulated and a Maxwellian distribution is assumed for their distribution. Finally, the computation of the power along the beam is different between the two models. Indeed, ORB5 uses Albajar's linear model to compute the local power of the beam and therefore disregards some quasi-linear effects. On the other hand, the calculation scheme of LUKE/C3PO, detailed in [START_REF] Peysson | Numerical simulations of the radio-frequency-driven toroidal current in tokamaks[END_REF], retains all quasi-linear effects. In LUKE, a first guess of the linear power deposited by the EC waves is found by the ray-tracing, assuming a Maxwellian plasma. The total wave power is carried by a total of 24 independent rays, and the different ray contributions are integrated together over real and velocity space. Then, the associated diffusion coefficient is evaluated and used in the Fokker-Planck solver to find the new distribution function. The deposited power is re-evaluated, and ray power content is corrected by wave power balance so that a new diffusion coefficient is found. A convergence loop is made such that the calculations of the absorbed power density using the diffusion coefficient and the electron distribution function are consistent with a relative accuracy of ≤ 1%, dealing with the quasilinear nature of the absorbed power (interdependence of absorbed power and distribution function).

The chosen equilibrium for the benchmark is an analytic TCV-like equilibrium generated with the method explained in the appendix of [START_REF] Tronko | Intercode comparison of gyrokinetic global electromagnetic modes[END_REF], with minor radius a = 0.25 m, major radius R 0 = 0.89 m, magnetic field on axis B 0 = 1.4 T, a parabolic safety factor profile q(r = R -R 0 ) = q 0 + (q a -q 0 )(r/a) 2 (with q 0 = 1 and q a = 3), constant density N ref e and a temperature profile corresponding to realistic temperature profiles obtained in TCV, as described in [START_REF] Sauter | On the non-stiffness of edge transport in l-mode tokamak plasmas[END_REF] T e (ρ P )

T edge =        Taxis + λ T ρ 2 P if ρ P ≤ ρ P,1 Tped exp [-κ T (ρ P -ρ P,2 )] if ρ P,1 ≤ ρ P ≤ ρ P,2 1 -µ T (ρ P -ρ P,3 ) if ρ P,2 ≤ ρ P ≤ ρ P,3 1 if ρ P ≥ ρ P,3 (67)
where ρ P = ψ P /ψ P,a (and ψ P is the poloidal magnetic flux), ρ P,1 = 0.27, ρ P,2 = 0.72, ρ P,3 = 0.90 and

   Taxis = Tcore -λ T ρ 2 P,1 Tcore = Tped exp (-κ T (ρ P,1 -ρ P,2 )) Tped = 1 -µ T (ρ P,2 -ρ P,3 ) (68)
with λ T = -κ T Tcore /(2ρ P,1 ), κ T = 3, µ T = 12 and T edge = 140 eV. A pure hydrogen plasma is assumed with N e = N i and T e = T i , and no Ohmic electric field is considered (this option is not available in ORB5). In that case, the system is invariant by the change of coordinates v ∥ → -v ∥ and θ 0 → π -θ 0 . Then, for the energy absorption, there is no difference between θ 0 ≤ π/2 and θ 0 ≥ π/2.

For all ORB5 simulations, 2 • 10 8 and 10 9 markers have been used for representing the ions and electrons respectively. With these resolutions, the results are converged with respect to the number of markers. In LUKE, the size of the phase space grid (p, ξ = p ∥ /p, ρ) is 151×121×30. The points in the spatial grid are not linearly spaced, they are concentrated around the power absorption location, guessed from the linear ray-tracing calculation.

ORB5 simulations have been performed using two different proton-electron mass ratios (m i /m e = 100 or 1836). The advantage of using an artificially reduced mass ratio is to increase the time step of the simulation (∆t = 1.0 Ω -1 ci the simulations with i /m e = 100 and ∆t = 0.1 Ω -1 ci for the simulations with m i /m e = 1836), leading to lighter numerical cost. In order to have the collision physics as close as possible between the two types of simulations, a modification of the collisionality is applied for the simulations with the artificial mass ratio. Indeed, the collision frequency of a species a on a species b reads [START_REF] Donnel | Moment approach of the multi-species nonlinear coulomb collision operator adapted to particle-in-cell codes[END_REF]

ν ab = ln Λ 3π 3/2 ϵ 2 0 N b T 3/2 a e 2 a e 2 b √ m a 1 + ma m b 1 + maT b m b Ta 3/2 . (69)
By considering the fact that electrons are much lighter than ions, one can see that the electron-electron ν ee and electron-ion ν ei collision frequencies are proportional to m

-1/2 e
. For simulations with heavy electrons, the correct ratio between the collision frequency of electrons and the electron cyclotron frequency, which is the most relevant dimensionless parameter to study the impact of the source in the presence of collisions, can be retrieved by simply multiplying the collision frequency of reference in ORB5 by m real e /m arti e . This choice is made to retain the correct physics as much as possible when using an artificial mass ratio between protons and electrons. Note that, with this choice, the electron collisionality ν ⋆ e = (νee+νei)qR ϵv T e is lower than the real one. This choice therefore lowers the neoclassical transport of electrons. Furthermore, this hypothesis could have an impact in the presence of trapped electron mode turbulence, which is sensitive to the electron collisionality.

The beam frequency is slightly decreased compared with section 4.3 (75GHz instead of 78GHz) to have an absorption of the beam power on the low field side. Indeed LUKE uses polar coordinates, which makes the deposition of the EC power and the beam propagation near the magnetic axis difficult to treat numerically. Note that there is no difficulty for ORB5 close to the magnetic axis due to the usage of (R, Z) coordinates for the EC module. The two angles of injection θ 0,in used in the benchmark are illustrated in figure 4, which shows C3PO ray-tracing results (poloidal and top views) for equatorial propagation of the beam. In the pure ECRH case (θ 0,in = 0.5π) the ray goes straight in the plasma as expected, whereas for the case with an ECCD component (θ 0,in = 0.4π), the beam is slightly curved, accordingly to predictions in the subsection 2.2. Comparisons between ORB5 and LUKE/C3PO for the two angles of injection and for two different densities (N ref In each of these plots, the linear predictions of Albajar and C3PO are given as reference points. Error bars on LUKE absorbed power come from the relative gap on absorbed power density between the last iteration and the one before (basically less than 1% for completely converged simulations). In LUKE, error bars on the Full-Width at Half Maximum (FWHM) of the deposition profile and on absorption location come from the spatial grid resolution. In ORB5, error bars on all quantities correspond to the standard deviation with respect to time.

At low power, an excellent agreement is found between all the predictions for the pure ECRH cases Fig. 5,7. The agreement for the ECCD cases Fig. 6, 8 at low power is still good, but a mismatch of approximately 10% is found for the total power absorbed predicted by ORB5 compared with the other predictions. As the total power absorbed is not constrained in the model implemented in ORB5, the fact that the total power absorbed is above 100% is not a problem. If this point becomes critical for later applications of this source term, a convergence loop between the energy deposited on the distribution function and the one included in the ray could be implemented in a similar way as the one used in LUKE. This possibility has not been added in ORB5 yet, as it would significantly increase the numerical cost of the source term. The overall agreement at low power for all cases validates the implementation of the model in ORB5.

For high-power beams, quasi-linear effects become important, especially for the total absorbed power. For the pure ECRH cases Fig. 5,7, quasi-linear effects lead to a significant decrease of the total power deposited both in LUKE and ORB5. These quasi-linear effects are stronger for LUKE than for ORB5. The quantitative difference between the results from the two codes can come from the fundamental differences between the two codes, i.e. 3D bounce-average for LUKE compared with 5D driftkinetic for ORB5, but also from the difference between the collision operators (linear with relativistic correction for LUKE, nonlinear but without relativistic correction for ORB5), or the feedback loop on the power absorbed by the plasma and the power of the beam, which is included in LUKE but not in ORB5. It is also noteworthy that ORB5 simulations using the real mass ratio between electrons and protons display larger quasi-linear effects than the one using heavy electrons. This difference is probably due to the difference in collisionality between the two versions of the code. It implies that simulations using the true electron mass will have to be used for the simulations dedicated to the study of quasi-linear effects. Simulations with a non zero current drive component, Fig. 6, 8, display smaller quasi-linear effects in both codes.

Quasilinear dependencies at high injected power can be interpreted by writing the equation for the deposited power density in the momentum-pitch-angle coordinates with D n the diffusion tensor resulting from wave-particle interaction, as defined in [START_REF] Donnel | Quasilinear treatment of wave-particle interactions in the electron cyclotron range and its implementation in a gyrokinetic code[END_REF], and D pp , D pξ , D ξp and D ξξ its in (p, ξ) coordinates. The dominant diffusion term is D pp (at least by an order of magnitude in all cases considered in paper), so the quasilinear behavior can be partly understood by looking the product pp ∂p at the absorption location in phase space. The quasilinear evolution of the total absorbed power results from a balance between the relative evolution of D pp with injected power due to the modification of the beam power along its path (quasilinear effect due to the loop between the energy deposited on the plasma and the power carried by the beam, which is not included in ORB5) and the deviation F e from a Maxwellian, resulting the competition between the source, collisions and the neoclassical physics is not inluded in these simulations).

In the pure ECRH case, the flattening of the distribution by the source, illustrated on Fig. 9 is the dominant effect and leads to a drop of the fraction of power absorbed as seen in Fig. 5 and7. The flattening of the distribution function increase with the absorbed power, explaining the dependence of quasi-linear effects with respect to the injected power seen on Fig. 5 and7. On the opposite, an increase of the collisionality, for instance by increasing the density, reduces the deviations from the Maxwellian and explains the difference of quasilinear effects between Fig. 5 and 7 at large power. For the case with ECCD, represented in Fig. 6 and8, the combined effect of the source and the collision pitch-angle scattering can result in a large pitch-angledependent population of suprathermal electrons, illustrated in Fig. 10, which shows the electron distribution in (p, ξ) space at the absorption location. The presence of a suprathermal population increases the power absorbed. In the case of ECCD, the two mechanisms described above, i.e. flattening of the distribution function by the source which leads to a reduction of the absorbed power on the one hand, and the increase of the suprathermal population by the combination of the source and collisions leading to an increase of the absorbed on the other hand, are direct competition. The overall effect is therefore case dependent and sensitive to the collision operator.

explains the difference between ORB5, which predicts a decrease of the power in the ECCD case studied here, and LUKE, which predicts a small increase of the absorbed power as can be seen in Fig. 6 and8.

The combination of low density, high temperature and high power leads to difficulty in convergence of the quasilinear power in LUKE, as the collision frequency is low and the quasilinear diffusion is high. This translates the fact that the power density becomes leading to strong quasilinear distortions of the distribution, with less collisions to the hole wave in the distribution. LUKE simulations convergence is ensured in figures 5, 6 and 7, but it could not be reached for P 0 > 100 kW at lower density in ECCD (N ref e = 5 • 10 18 m -3 , Fig. 8). The quantitative agreement between ORB5 and LUKE for low power deposition and the qualitative agreement at high power validate the derivation of the source term and its implementation in ORB5. The ORB5 code is therefore ready to study the turbulence in the presence of the ECRH/CD source. This study is left for a future paper.

Figure 10. 2-D distribution function of electron obtained with LUKE, at maximum of absorption ρ G,abs for each angle, for high injected power. The dominant term of the RF quasilinear bounce-averaged diffusion coefficient Dpp is overlaid, as well as a theoretical prediction for the resonance location (see [START_REF] Prater | Heating and current drive by electron cyclotron waves[END_REF]). In ECCD, the resonance is pitch-angle dependent, leading to wave interaction with particles at higher p. At low power, the distribution remains close to Maxwellian.

Conclusion

In this paper, a quasilinear operator representing the electron cyclotron resonance heating and current drive is derived in the specific case of a beam propagating in the equatorial plane of an axisymmetric tokamak. This assumption allows one to treat the beam propagation in the plasma analytically, removing the need for the coupling with a dedicated beam propagation code, such as a ray-or beam-tracing code.

This operator can be used as a source term in any global gyrokinetic code. It has been implemented in the particle-in-cell code ORB5. The numerical treatment specific to particle-in-cell code is detailed. The numerical implementation is validated against an analytical model for the linear absorption. The quasilinear corrections are also tested by benchmarking ORB5 against the bounce-average drift-kinetic code LUKE, which is dedicated to the simulation of the electron distribution function in the presence of plasma/wave interactions. A good agreement is found between the two codes despite some differences in the limit of large power deposition, which are discussed in detail. A further development of the model may include the implementation of a convergence loop on deposited power (as it is done in LUKE), but this will dramatically increase the cost of already heavy simulations, especially if a realistic electron-to-ion mass ratio is required.

The implementation of source term described in this paper in any global gyrokinetic code enables the self-consistent study of turbulence in presence of ECRH/ECCD by this code. The next step of this project is to study interplay between this realistic source of momentum and energy on the one hand, and the turbulence on the other hand. This study will help to elucidate a long standing question about the possibility for turbulence to radially transport suprathermal electrons generated by the resonant wave/particle interaction outside the resonant area, a mechanism has been proposed to reconcile theoretical predictions with experimental observations made in some cases [START_REF] Coda | Electron cyclotron current drive and electron dynamics in the TCV tokamak[END_REF][START_REF] Nikkola | Modelling of the electron cyclotron current drive experiments in the TCV tokamak[END_REF].

Acknowledgment

This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project s1067. Some of the simulations presented here were obtained on the EUROfusion partition of the MARCONI100 supercomputer. This work was supported in part by the Swiss National Science Foundation. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

Appendix A: Group velocity for an oblique propagation

For the evaluation of Eq.27, it is necessary to compute the group velocity in the general case. It can be shown that the group velocity is given by

v g,O/X = cN O/X N 2 O/X + ω b 2 ∂N 2 O/X ∂ω b (72)
The cold plasma dispersion relation Eq.12 can be written differently

N 2 O/X = f 1 ω 2 b ± Ωeω 2 p ω b f 2 (ω 2 b ) f 3 (ω 2 b ) (73)
where The analytical expression of the group velocity Eq.77 can be used for an efficient implementation linking the electric field of the beam to its power Eq.27 where z n = 1 2 4x 2 n -y 2 n + iy n is derived with respect to x n and y n . These derivatives are useful to evaluate Eq.33 in the theoretical prediction.

f 1 (x) = 2 tan 2 θ + 1 x -ω 2
Case ℜ (z n ) > 0 When ℜ (z n ) > 0 , i.e. when 4x 2 n > y 2 n , one can use the spherical Bessel function of first kind j m which is related to J m+ 1 2 via the relation (Eq. 10.1.1 of [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF])

J m+ 1 2 (z) = 2z π j m (z) (81)
which is valid for any integer m and any complex argument z. In practice, only the low m integer are useful. One can then use the fact that (Eq.10.1.10 of [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF])

j
The first values of f m (z) are given in Tab 1. Using the fact that |z n | = x n , it is possible to show that

m f m (z) m f m (z) -5 105z -4 -10z -2 0 z -1 -4 -15z -3 + z -1 1 z -2 -3 3z -2 2 3z -3 -z -1 -2 -z -1 3 15z -4 -6z -2 -1 0 4 105z -5 -45z -3 + z -1
J m+1/2 (z n ) 2 x n = 2 π |j m (z n )| 2 ( 84
)
where zn is the complex conjugate of z n . Then using the relationship linking j m to its derivative (Eq. 10.1.21 of [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF])

zj ′ m (z) = -(m + 1) j m (z) + zj m-1 (z) , (85)
It is possible to show that

∂ ∂x n    m+ 1 2 (z n ) 2 x n    = 4x n [G ′ m (z - n ) G m (z + n ) -G ′ (z + n ) G m (z - n )] π z + n -z - n (96) ∂ ∂y n    J m+ 1 2 n ) 2 x n    = 2 [z + n G ′ m (z + n ) G m (z - n ) -z - n G ′ m (z - n ) G m (z + n )] π z + n -z - n (97) ∂ 2 ∂y 2 n    J m+ 1 2 (z n ) 2 x n    = 2 π z + n -z - n 2 z + n 2 G ′′ m z + n G m z - n -z - n 2 G ′′ m z - n G m z + n -2x 2 n G ′ m z + n G ′ m z - n (98) ∂ 2 ∂x n ∂y n    J m+ 1 2 (z n ) 2 x n    = 4x n π z + n -z - n 2 G ′ m z + n G m z - n -G ′ m z - n G m z + n +y n G ′ m z + n G ′ m z - n -z + n G ′′ m z + n G m z - n -z - n G ′′ m z - n G m z + n (99) m g m (x) m g m (x) -5 -105x -4 -10x -2 0 x -1 -4 15x -3 + x -1 1 -x -2 -3
-3x -2 2 3x -3 + x -1 -2

x -1 3 -15x -4 -6x -2 -1 0 4 105x -5 + 45x -3 + x -1

∂v ∥ = [1 + C 1 (θ 0 )] J ′ n+1 (ρ 0 ) + [1 -C 1 (θ 0 )] J ′ n-1 (ρ 0 ) -2C 2 (θ 0 ) v ∥ v ⊥ J ′ n (ρ 0 ) 2 1 + C 1 (θ 0 ) 2 + C 2 (θ 0 ) 2 ∂ ρ0 ∂v ∥ - C 2 (θ 0 ) J n (ρ 0 ) v ⊥ 1 + C 1 (θ 0 ) 2 + C 2 (θ 0 ) 2 (100) ∂ Θn k,0 ∂v ⊥ = [1 + C 1 (θ 0 )] J ′ n+1 (ρ 0 ) + [1 -C 1 (θ 0 )] J ′ n-1 (ρ 0 ) -2C 2 (θ 0 ) v ∥ v ⊥ J ′ n (ρ 0 ) 2 1 + C 1 (θ 0 ) 2 + C 2 (θ 0 ) 2 ∂ ρ0 ∂v ⊥ + C 2 (θ 0 ) v ∥ J n (ρ 0 ) v 2 ⊥ 1 + C 1 (θ 0 ) 2 + C 2 (θ 0 ) 2 (

e 2 /

 2 (ϵ 0 m e ) is the plasma frequency, S = (R + L) /2, T = (R -L) /2, where R = (P -Ω e /ω b ) / (1 -Ω e /ω b ) and L = (P + Ω e /ω b ) / (1 + Ω e /ω b ).

Figure 1 .

 1 Figure 1. Link between X and X ′ coordinates.

2 T e Ωce = Dnm 2 e

 22 relation. One needs to use an artificial speed of light c arti = c real m real e /m arti e in all terms of the form v T /c to ensure that v arti T /c arti = v real T /c real . Concerning the electron mass appearing in the prefactor of the resonant diffusion coefficient, one needs to use the artificial electron mass. The purpose of this choice is to keep the same normalized diffusion coefficient Dn v

Figure 2 .= 2 •

 22 Figure 2. Fraction of power absorbed as a function of the major radius for a pure ECRH case (θ in = 0.5π) on the left and a mixed case ECRH/ECCD (θ in = 0.4π) on the right. The density and temperature, are scanned around reference values N ref e = 2 • 10 18 m -3 , T ref e = 1.17 keV.

Figure 3 .

 3 Figure 3. Convergence test for a pure ECRH case (θ in = 0.5π) on the left and a mixed ECRH/ECCD case (θ in = 0.4π) on the right. The electron density Ne = 2 • 10 18 m -3 and temperature Te = 1.17 keV are constant.

e = 5 •= 1 •

 51 10 18 m -3 and N ref e 10 19 m -3 ) are shown in figures 5, 6, 7 and 8.

Figure 4 .

 4 Figure 4. Results of C3PO ray tracing for equatorial midplane ECCD, launched from the Low Field Side, with an EC beam frequency of 75 GHz, and the angle θ 0,in between the beam wave vector k and the magnetic field B at the beam entrance in the plasma ranging from π/2 (pure ECRH) to 0.4π.

Figure 5 .

 5 Figure5. Results of injected power P 0 scan, for θ 0,in = π/2, for Albajar's theoretical prediction, C3PO (linear calculation), quasilinear LUKE, and ORB5, showing total power fraction (left), power absorption location ρ G,abs (right) and full width at half maximum (∆ρ G ) abs of the power deposition profile dP abs /dρ G , where ρ G = (R -R 0 )/a, (center ).

Figure 6 .

 6 Figure 6. Results of injected power P 0 scan, for θ 0,in = 0.4π, for Albajar's theoretical prediction, C3PO (linear calculation), quasilinear LUKE, and ORB5, showing total absorbed power fraction (left), power absorption location ρ G,abs (right) and full width at half maximum (∆ρ G ) abs of the power deposition profile dP abs /dρ G , where ρ G = (R -R 0 )/a, (center ).

1 -ξ 2 S

 12 ξ dξ, (70) where the flux S in (p, ξ) is given byS = -D n • ∇ p F e = S p

Figure 7 .= 5 •

 75 Figure 7. Results of injected power P 0 scan, for θ 0,in = 0.5π at lower density (N ref e = 5 • 10 18 m -3 ), for Albajar's theoretical prediction, C3PO (linear calculation), quasilinear LUKE, and ORB5, showing total absorbed power fraction (left), power absorption location ρ G,abs (right) and full width at half maximum (∆ρ G ) abs of the power deposition profile dP abs /dρ G , where ρ G = (R -R 0 )/a, (center ).

Figure 8 .= 5 •

 85 Figure 8. Results of injected power P 0 scan, for θ 0,in = 0.4π at lower density (N ref e = 5 • 10 18 m -3 ), for Albajar's theoretical prediction, C3PO (linear calculation), quasilinear LUKE, and ORB5, showing total absorbed power fraction (left), power absorption location ρ G,abs (right) and full width at half maximum (∆ρ G ) abs of the power deposition profile dP abs /dρ G , where ρ G = (R -R 0 )/a, (center ).

Figure 9 .

 9 Figure 9. Pitch-angle integrated distribution function of electrons obtained with LUKE (left) and its derivative with respect to p (right), for the pure ECRH, at the maximum of absorption ρ G,abs , compared with the initial distribution. The distribution function is for the high injected power case. The shaded area shows a rough approximation of the resonance location. The flattening of the distribution function leads to a decrease of the local power deposited by ECRH

p 2 -xΩ 2 e + tan 2 θ + 2 ω 2 p Ω 2 ee tan 4 θx + 4 x -ω 2 p 2 2 e ( 76 ) 2 b g 1 1 Ω 2 eg 2 4 pg 3 2 e

 222227621122432 tan 2 θ + 1 (75)f 3 (x) = 2 x -ω 2 p -Ω 2 e x tan 2 θ + 1 + ω 2 p ΩThis formulation allows one to express the group velocity asv g,O/X = cN O/X f 3 ω (x) = xf 1 (x) = 2 tan 2 θ + 1 x 2 x -ω 2 p -(x) = xf ′ 2 (x) -f (x) = 4 tan 2 θ + 1 x 2 -ω (x) = xf ′ 3 (x) -f 3 (x) = 2 2 θ + 1 x 2 -ω 2p Ω

  Appendix B: Evaluation of the Bessel functions and their derivatives for the theoretical prediction In this section, the function |Jm+1/2(zn)| 2 xn

  m (z) = f m (z) sin (z) + (-1) m+1 f -m-1 (z) cos (z) (82)where f 0 (z) = z -1 , f 1 (z) = z -2 and (Eq. 10.1.19 of[START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF])f m-1 (z) + f m+1 (z) = (2m + 1) f m (z) /z.

Table 1 .

 1 Values of fm (z) for low m integers

Table 2 .

 2 Values of gm (x) for low m integers Appendix C: Velocity derivatives of Θn

	k,0
	∂ Θn k,0

one can show that

When ℜ (z n ) = 0 , i.e. when 4x 2 n < y 2 n , it is convenient to replace the Bessel functions of first kind occurring in Eq.33 by modified Bessel functions of first kind of same order and real arguments by using the property

where

One can use the property (Eq.10.2.12. of [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF])

with

which is valid for any integer m and real argument x. The functions g m are defined by recurrence as g 0

The first values of g m (x) are given in Tab 2. Using the fact that z + n z - n = x 2 n , one gets


Ryuki Hyodo
email: hyodo@elsi.jp

Ida Shigeru

Tristan Guillot

A "no-drift" runaway pile-up of pebbles in protoplanetary disks in which midplane turbulence increases with radius

Keywords:

Context. A notable challenge of planet formation is to find a path to directly form planetesimals from small particles. Aims. We aim to understand how drifting pebbles pile up in a protoplanetary disk with a nonuniform turbulence structure. Methods. We consider a disk structure in which the midplane turbulence viscosity increases with the radius in protoplanetary disks, such as in the outer region of a dead zone. We perform 1D diffusion-advection simulations of pebbles that include back-reaction (the inertia) to the radial drift and the vertical and radial diffusions of pebbles for a given pebble-to-gas mass flux. Results. We report a new mechanism, the "no-drift" runaway pile-up, that leads to a runaway accumulation of pebbles in disks, thus favoring the formation of planetesimals by streaming and/or gravitational instabilities. This occurs when pebbles drifting in from the outer disk and entering a dead zone experience a decrease in vertical turbulence. The scale height of the pebble subdisk then decreases, and, for small enough values of the turbulence in the dead zone and high values of the pebble-to-gas flux ratio, the back-reaction of pebbles on gas leads to a significant decrease in their drift velocity and thus their progressive accumulation. This occurs when the ratio of the flux of pebbles to that of the gas is large enough that the effect dominates over any Kelvin-Helmholtz shear instability. This process is independent of the existence of a pressure bump.

Introduction

Forming planetesimals directly from small particles in an evolving protoplanetary disk is a major challenge in planet formation due to the "growth barrier" (Blum & Wurm 2000;[START_REF] Zsom | [END_REF]) and the "drift barrier" [START_REF] Whipple | From Plasma to Planet[END_REF]Weidenschilling 1977). Streaming instability (SI) may be a promising mechanism for forming planetesimals directly from pebbles, although it requires special conditions: the Stokes number τ s 0.01 and a local elevated solid-to-gas ratio (Z ≡ ρ p /ρ g 1, where ρ p and ρ g are the spatial densities of pebbles and gas, respectively; [START_REF] Youdin | Planetary Systems in the Universe[END_REF]Youdin & Goodman 2005;Carrera et al. 2015). How, when, and whether such conditions are met in an evolving protoplanetary disk are still a matter of intense debate.

Different mechanisms have been proposed for planetesimal formation via the local enhancement of materials. They include the diffusive redistribution and recondensation of water vapor outside the snow line (e.g., Stevenson & Lunine 1988;Ciesla & Cuzzi 2006;Ros & Johansen 2013), and detailed numerical simulations have been performed (e.g., [START_REF] Drżkowska | A92 L9[END_REF]Schoonenberg & Ormel 2017;Hyodo et al. 2019Hyodo et al. , 2021;;Gárate et al. 2020). Other mechanisms include pile-ups at pressure maxima. These pressure maxima can be created by the planet's gravity (e.g., Dipierro & Laibe 2017;Kanagawa et al. 2018), by a sharp change in the local magneto-rotational instability-(MRI-) driven turbulence structure at an evaporation front (e.g., Kretke & Lin 2007;Brauer et al. 2008), by a change in the gas accretion velocity at the inner (e.g., Chatterjee & Tan 2014;Ueda et al. 2019;Charnoz et al. 2019) and outer (e.g., Pinilla et al. 2016) edges of a dead zone, and/or by the change in the gas profile due to gas-pebble friction (back-reaction) combined with pebble growth (Gonzalez et al. 2017). However, a positive pressure gradient may not be likely because the non-magnetohydrodynamic (non-MHD) effects erase a clear transition between dead and active zones (e.g., Mori et al. 2017). The back-reaction of solids onto the gas would smooth out the bump (e.g., Taki et al. 2016;Kanagawa et al. 2018). Further studies are required to assess these likelihoods in an evolving disk.

Over the last few decades, MHD simulations have shown that a region where ionization is too low for the MRI to operate (i.e., a "dead zone"; Gammie 1996) might ubiquitously exist in the inner part of the disk midplane, and that only surface layers are magnetically active, supporting accretion (Fig. 1; Gressel et al. 2015;Simon et al. 2015;Bai & Stone 2013;Bai et al. 2016;Mori et al. 2017). The transition between active and dead zones might be smooth (i.e., the α-parameter in the disk midplane α mid increases as a function of the distance to the star r), and α dead = 10 -5 -10 -3 is reported for the α-parameter within a dead zone (e.g., Gressel et al. 2015;Simon et al. 2015;Mori et al. 2019). We note that the sound waves propagated from MRI-active surface layers (e.g., Okuzumi & Hirose 2011;[START_REF] Yang | [END_REF]) and the vertical shear instability (VSI; e.g., Fig. 1. Schematic illustration of pebble drift and its pile-up within a protoplanetary disk with a dead zone. The disk gas accretion is characterized by α-parameter a acc , while the midplane diffusivity, being a dead zone (α dead α acc ) in the inner region, is characterized by α mid . During the inward drift of pebbles, the pebble scale height H p decreases as H p ∝ α 1/2 mid until a KH instability prevents it from becoming smaller. A smaller H p leads to an elevated local midplane concentration of pebbles within a thinner midplane layer. The elevated midplane pebble-to-gas ratio causes the back-reaction to be more effective in reducing the radial drift velocity of pebbles. Such a physical interplay with a sufficiently large pebble-to-gas mass flux results in a progressive accumulation of pebbles in a runaway fashion (i.e., the "no-drift" runaway pile-up). Stoll & Kley 2016;Flock et al. 2020) could induce vertical mixing with an equivalent α for vertical diffusivity up to ∼10 -3 , even within a dead zone. These results imply nonuniform disk turbulence structures.

In this Letter, we consider a toy model of a protoplanetary disk with a dead zone located at the inner region of the disk. We consider different α-parameters for the gas accretion and for the pebble motion. Including the effects of the back-reaction of pebbles onto gas, which slows the drift velocity of pebbles, we study how drifting pebbles within a dead zone pile up. We demonstrate that a runaway pile-up of pebbles occurs for a sufficiently large pebble-to-gas mass flux F p/g when pebbles reach a critical low-level of turbulence; we call this the "no-drift" (ND) runaway pile-up (the summary of this mechanism is shown by a schematic illustration in Fig. 1). It is worth mentioning that the newly reported ND mechanism does not require the snow line or the pressure maxima.

In Sect. 2, we describe the disk models and settings of our 1D simulations that solve the diffusion-advection of drifting pebbles. In Sect. 3, we show the results of our 1D simulations as well as analytical arguments. In Sect. 4, we summarize this Letter.

Models and settings

Gas structure

Here, we adopt a toy model where the gas accretion toward the central star is characterized by a non-dimensional α-parameter α acc (Fig. 1). Using the prescription of the classical α-accretion disk model (Shakura & Sunyaev 1973;Lynden-Bell & Pringle 1974), the surface density of the gas is given as

Σ g = Ṁg 3πν acc = Ṁg 3πα acc c 2 s Ω -1 K , (1)
where Ṁg and ν acc = α acc c 2 s Ω -1 K are the gas mass accretion rate and the effective viscosity, respectively (c s is the gas sound velocity and Ω K is the Keplerian orbital frequency).

The gas rotates at sub-Keplerian speed. The degree of the deviation of the gas rotation frequency from that of Keplerian η is given by

η ≡ Ω K -Ω Ω K = - 1 2 ∂ ln P g ∂ ln r H g r 2 = C η H g r 2 , (2)
where Ω, P g , and H g are the gas orbital frequency, the gas pressure, and the gas scale height, respectively, and C η is defined as

C η ≡ - 1 2 ∂ ln P g ∂ ln r , ( 3
)
which depends on the temperature profile (e.g., C η = 11/8 for T ∝ r -1/2 ).

Using the α-disk prescription, the gas accretion velocity v g is written as

v g = - 3ν acc 2r = - 3α acc H 2 g Ω K 2r = - 3α acc 2 H g r 2 v K , (4)
- 3 2 α acc ηv K - 1 2 d ln P g d ln r -1 = - 3 2 α acc ηv K C -1 η , (5)
where a negative sign indicates accretion toward the central star.

Pebbles in the disk midplane

The radial drift velocity of pebbles, including the effects of gas-solid frictiondrift back-reaction (hereafter Drift-BKR)is given as (Ida & Guillot 2016;Schoonenberg & Ormel 2017;Hyodo et al. 2019)

v p = - Λ 1 + Λ 2 τ 2 s 2τ s Ληv K -v g , (6)
where τ s is the Stokes number of pebbles. Λ ≡ ρ g /(ρ g + ρ p ) = 1/(1 + Z) characterizes the strength of the back-reaction due to the pile-up of pebbles, where Z ≡ ρ p /ρ g is the midplane pebbleto-gas density ratio (ρ g = Σ g / √ 2πH g and ρ p = Σ p / √ 2πH p , where Σ p and H p are the surface density and scale height of pebbles, respectively).

The disk midplane could be too weakly ionized for the MRI to operate (the dead zone with extremely weak turbulence; Gammie 1996). Thus, the disk midplane could have very different turbulence from that of gas accretion, which controls the disk surface density (i.e., α acc ). In this Letter, we use an effective viscous parameter α mid to characterize the radial and vertical diffusion processes of pebbles near the disk midplane (Fig. 1).

The scale height of pebbles characterizes the degree of pebble concentration in the disk midplane (ρ p ∝ H -1 p ). In the steady state, the scale height of pebbles is regulated by the vertical turbulent stirring (Dubrulle et al. 1995;Youdin & Lithwick 2007;Okuzumi et al. 2012;Hyodo et al. 2019) as

H p,tur = 1 + τ s α mid (1 + Z) -K -1/2 H g , (7)
where a coefficient K characterizes the strength of the backreaction onto the diffusivity (Hyodo et al. 2019;Ida et al. 2021) and K = 1 is used (the choice of the K value does not alter the conclusion -see Appendix A for the K = 0 case). Fig. 2. α mid /α acc (top panels) and Z = ρ p /ρ g (bottom panels) as a function of the distance to a star in the 1D numerical simulations. Included are the cases where α acc = 10 -2 and τ p = 0.1. From left to right: cases where (F p/g , α dead , r * , ∆r tra ) = (0.11, 10 -4 , 9 au, 5 au), (0.17, 10 -4 , 9 au, 5 au), (0.47, 10 -4 , 9 au, 5 au), and (0.17, 10 -4 , 5 au, 0.5 au). The gray lines in the top panels represent the analytical critical α mid /α acc below which the ND runaway pile-up is expected to occur for a given F p/g (Eq. ( 23) with F p/g > F p/g,crit2 ). The gray lines in the bottom panels show the critical Z (i.e., Z crit = 1) above which the ND runaway pile-up occurs (Sect. 3.2). The blue, green, and red lines in the bottom panels are those at t = 3 × 10 4 yrs, t = 1 × 10 5 yrs, and t = 3 × 10 5 yrs from the beginning of the calculations, respectively. Panel a: system that reaches a steady-state, panels b-d: ND runaway pile-ups.

For a small α mid (i.e., for a small H p,tur ), a vertical shear Kelvin-Helmholtz (KH) instability prevents a further decrease in the pebble scale height. This minimum scale height, H p,KH , is (Hyodo et al. 2021)

H p,KH Ri 1/2 Z 1/2 (1 + Z) 3/2 C η H g r H g = Ri 1/2 Z 1/2 (1 + Z) 3/2 ηr, (8)
which is valid for Z 1, and Z 1 indicates a gravitational collapse. Thus, the scale height of pebbles H p is given as

H p = max H p,tur , H p,KH .
(9)

Numerical settings

We performed 1D diffusion-advection simulations that included the back-reaction (the inertia) to radial drift of pebbles that slows the pebble drift velocity as pile-up proceeds (see details in Hyodo et al. 2019Hyodo et al. , 2021)). The drift velocity of pebbles is given by Eq. ( 6), while the diffusivity is given as

D p = α mid c 2 s Ω -1 K Λ K /(1 + τ 2 s )
. We included back-reaction onto the diffusivity of pebbles (with K = 1). We assumed that pebbles drift outside the snow line, and the sublimation of pebbles was neglected (i.e., the dead zone exists beyond the snow line). The Stokes number of pebbles was set to be constant: τ s = 0.1 (Okuzumi et al. 2016;Ida & Guillot 2016). The disk midplane temperature was set to T (r) = 150 K × (r/3 au) -β (β = 1/2). The surface density of the gas (molecular weight of µ g = 2.34) is described by Σ g = Ṁg /3πν g , where ν g = α acc c 2 s Ω -1 K . We used α acc = 10 -2 and Ṁg = 10 -8 M yr -1 . These led to C η = 11/8 for Σ g ∝ r -1 and T ∝ r -1/2 .

We considered a dead zone in the inner region of the disk midplane, and we used a non-dimensional turbulence parameter in the midplane α mid , which differs from the one that characterizes gas accretion, α acc (Eq. ( 1)). The midplane α mid is modeled as

α mid (r) = α dead + α acc -α dead 2 erf 3 + 6 (r -r * ) ∆r tra + 1 , (10)
where: α acc and α dead are those outside and inside a dead zone; r * is the innermost radial distance, where α mid = α acc ; and ∆r tra is the radial width of the transition from α mid = α acc to α mid = α dead . As shown below, the ND mode occurs for an arbitrary choice of the dead zone structure, that is, irrespective of a sharp or a smooth change between the active and dead zones (i.e., an arbitrary choice of α dead , α acc , r * , and ∆r tra ), as long as it satisfies that α mid is smaller than a threshold value (Sect. 3.2.2). At the beginning of the 1D simulations, we set the pebbleto-gas mass flux F p/g at the outer boundary (r out = 15 au), and F p/g at r out was fixed throughout the simulations. We note that Elbakyan et al. (2020) performed a time-dependent simulation of disk formation to find a high F p/g variation between O(10 -4 ) to O(1) in a single disk evolution (see also Ida et al. 2021).

Results

Numerical results

Figure 2 shows the results of our 1D simulations for different combinations of F p/g and the midplane turbulence structures (i.e., different choices of α dead , r * , and ∆r tra in Eq. ( 10)) for the case of α acc = 10 -2 . The top panels show given turbulence structures (i.e., α mid (r)/α acc ) and the bottom panels show the resultant midplane pebble-to-gas ratio in their spatial densities (i.e., ρ p /ρ g ). The gray lines in the top panels are the analytically derived critical α mid /α acc below which a runaway pile-up of pebbles is expected to occur for a given F p/g (Eq. ( 23) with τ s = 0.1).

For a small F p/g (e.g., F p/g = 0.11; panel a in Fig. 2), the system reaches a steady state. For this given F p/g , analytical arguments (the gray line) predict the requirement of a smaller α mid /α acc for a runaway pile-up of pebbles to occur (i.e., the black line in Fig. 2 needs to be below the gray line).

For a larger F p/g (e.g., F p/g = 0.17; panel b in Fig. 2) that has the same α mid /α acc structure as that of panel a, pebbles progressively pile up over time. Continuous pile-up occurs as the drift velocity of pebbles progressively decreases due to strong Drift-BKR and as the pile-up efficiently proceeds. For the F p/g given here, the gray line overlaps with the black line (i.e., the analytical requirement for the runaway pile-up is met) and the analytical arguments (Sect. 3.2) show a very good accordance with 1D numerical simulations. For an even larger F p/g with the Fig. 3. Midplane pebble-to-gas ratio Z ≡ ρ p /ρ g in the α mid /α acc -F p/g space. Included are the cases where r = 5 au, C η = 11/8 for T ∝ r -1/2 , τ s = 0.1, and Ri = 0.5. Left and right panels: cases where α acc = 10 -2 and α acc = 10 -3 , respectively. The color contours are obtained by directly solving Eq. ( 11). The red regions indicate the no-drift runaway pile-up. The diagonal dashed black line (Eq. ( 16)) and the horizontal dashed black line (Eq. ( 20)) are the analytically derived critical F p/g,crit above which the ND runaway pile-up takes place. The critical Z above which the ND runaway pile-up occurs (i.e., Z crit = 1) is shown by the white lines. The parameter map is very weakly dependent on the radial distance to the star (see Eqs. ( 16) and ( 20)), and the ND mode can occur irrespective of the shape and position of a dead zone. same α mid (r)/α acc as those in panels a and b, a runaway pileup of pebbles occurs more quickly and significantly (panel c in Fig. 2).

Comparing panels b and d in Fig. 2 (the same F p/g , α dead , and α acc but different shapes and locations of the transition between active and dead zones), 1D simulations demonstrate that a runaway pile-up of pebbles occurs irrespective of the shape and location of the dead zone as long as F p/g is sufficiently large. This shows the robustness of this new mechanism of runaway pebble pile-up.

Numerical results show that when the change in the midplane turbulence on r is smooth (i.e., α mid ), the pile-up of pebbles efficiently propagates toward the outer region in the disk (compare panels b and d in Fig. 2). This indicates that, when the ND runaway pile-up occurs at a given radial distance, a progressive pile-up of pebbles could occur even at a greater radial distance where the ND conditions are not met, leading to a global pile-up of pebbles for the SI to operate.

When α acc = 10 -3 , 1D simulations only show a runaway pile-up of pebbles for a very large value of F p/g (F p/g 0.8). This is because the surface density of the gas increases (Eq. ( 1)) for a smaller α acc , and ρ p /ρ g correspondingly decreases.

Analytical arguments

Midplane concentration of pebbles

Below, we discuss the concentration of pebbles in the disk midplane considering Drift-BKR. Our arguments are made using the analytical equation presented below (Eq. ( 11)). We solve this equation using two distinct approaches: a direct numerical approach and an analytical approach with some approximation. In both ways, we identify a parameter space (F p/g -α acc -α mid space) where the ND runaway pile-up of drifting pebbles could take place.

The concentration of pebbles at the midplane is written as

Z ≡ ρ p ρ g = Σ p Σ g h -1 p/g = v g v p h -1 p/g F p/g , (11)
where h p/g ≡ H p /H g and v p are functions of Λ(Z). We directly find a solution (i.e., solve for Z) for Eq. ( 11) for a given F p/g and α mid . The color contours in Fig. 3 show directly obtained ρ p /ρ g and the red region is where no steady-state solution is found (i.e., the radial drift of pebbles is halted due to Drift-BKR). The dashed lines in Fig. 3 are analytically derived critical F p/g,crit , above which no steady-state solution is found for ρ p /ρ g (i.e., the ND runaway pile-up). The analytical estimations reproduce the direct solutions well. The F p/g,crit can be divided into two different regimes. In the first regime, H p,tur > H p,KH,max , where H p,KH,max is the maximum scale height of pebbles regulated by a KH instability. In this case, using Eqs. ( 5) and ( 6), the vertically averaged metallicity of pebbles is given as

Z Σ ≡ Σ p Σ g = F p/g × v g v p , (12)
F p/g × -3 2 α acc ηv K C -1 η - Λ 1+Λ 2 τ 2 s 2Λτ s + 3α acc 2C η ηv K , ( 13
)
3α acc F p/g 4τ s C η Λ 2 = 3α acc F p/g 4τ s C η 1 + h -1 p/g Z Σ 2 , (14)
where an approximation is made for α acc τ s 1, Z 1, Λ 1. Solving this equation gives

Z Σ = (1 -2ab) ± √ 1 -4ab 2ab 2 , ( 15
)
where a = 3α acc F p/g 4τ s C η and b = h -1 p/g . The Z Σ has a real solution when 1 -4ab > 0 and a critical F p/g,crit1 for the first regime is given as

F p/g,crit1 = C η 3 τ s α acc h p/g (α mid τ s ) 1/2 C η 3α acc , (16) 0
.15 × α acc 10 -2 -1 α mid 10 -4 1/2 τ s 0.1 1/2 C η 11/8 (17)
for H p,tur > H p,KH,max , where h p/g τ s / α mid (1 + Z) -K -1/2 (τ s /α mid ) -1/2 for Z 1 and α mid 1 + τ s (Eq. ( 9)). The above criteria for F p/g,crit1 are independent of r, and thus the ND mechanism occurs irrespective of the position and shape of a dead zone as long as it satisfies its criteria (i.e., Eq. ( 16)). The diagonal dashed black line in Fig. 3 shows Eq. ( 16) and is in accordance with the direct solutions of Eq. ( 11) (color contours in Fig. 3). The critical Z above which the ND occurs, Z crit , is obtained by inserting Eq. ( 16) into Eq. ( 11), and Z crit ∼ 1 for the first regime. Analytical Z crit = 1 agrees with the direct solution (color contours in Fig. 3). We note that because F p/g,crit1 is proportional to C η , it increases with a steeper Σ g gradient at an outer boundary of a dead zone. However, the ND mode does not need a sharp boundary. If the boundary is very gradual, the change in C η is negligible.

The second regime appears when a KH instability plays a role, which corresponds to H p,tur < H p,KH,max . Thus, in this case, the scale height of pebbles is described by H p,KH , which is independent of α mid . As ρ p /ρ g increases from zero, H p,KH initially increases and reaches its maximum H p,KH,max at Z = ρ p /ρ g = 1/2. As ρ p /ρ g further increases, H p,KH decreases1 . By smoothly connecting with the critical value of Z crit = 1 in the first regime (i.e., Z crit = 1 as well in the second regime, which agrees with the direct solution; see the color contours in Fig. 3), the critical F p/g in the second regime F p/g,crit2 adopts H Z=1 p,KH (H p,KH with Z = 1), and Z crit is given as

Z crit = ρ p ρ g | H Z=1 p,KH = Σ p Σ g h Z=1 p/g,KH -1 , ( 18
)
3α acc F p/g,crit2 4τ s C η (Λ Z=1 KH ) 2 h Z=1 p/g,KH -1 , (19)
where h Z=1 p/g,KH ≡ H Z=1 p,KH /H g and an approximation is made for α acc τ s 1 and Λ 1. The Λ Z=1 KH is Λ using H Z=1 p,KH . Here, Z = 1 and Λ Z=1 KH = 1/2. Thus, F p/g,crit2 with C η = 11/8 is given as

F p/g,crit2 = C η 3 τ s α acc h Z=1 p/g,KH = 11 24 τ s α acc h Z=1 p/g,KH , (20) 0
.06 × α acc 10 -2 -1 τ s 0.1 Ri 0.5 1/2 H g /r 0.04 , (21)
for H p,tur < H p,KH,max which is independent of α mid . We note that the above criterion for F p/g,crit2 has a very weak dependence on r through H g /r. Equation ( 20) is plotted by a horizontal dashed black line in Fig. 3, and it shows a good consistency with the direct solution of Eq. ( 11). Including a KH instability (i.e., the second regime and Eq. ( 20)) prevents the scale height of pebbles from becoming increasingly smaller as α mid decreases, reducing the parameter range of the ND runaway pile-up2 .

Criteria for the ND runaway pile-up

Summarizing the analytical arguments above, the scale height of pebbles decreases with decreasing α mid until a KH instability plays a role. As pebbles drift inward from the active to the dead zones in the disk midplane, α mid becomes smaller with decreasing distance to the star r, while keeping a constant F p/g . This implies that the evolutionary path in the α mid /α acc -F p/g space (Fig. 3) is a horizontal shift from right to left for a given fixed F p/g as pebbles drift inward. Therefore, for the ND mode to take place, the following two conditions need to be met.

First, F p/g needs to be sufficiently large to satisfy F p/g ≥ F p/g,crit2 . As F p/g,crit2 ∝ 1/α acc (Eq. ( 20)), the ND mode is limited to relatively large F p/g values for the α acc = 10 -3 case (F p/g 0.8 at 5 au), while only moderate F p/g 0.08 is required in the case of α acc = 10 -2 .

Second, a sufficiently small α mid /α acc is required in order to have a sufficiently small scale height of pebbles. Rewriting Eq. ( 16) gives a critical α mid /α acc below which the ND occurs for a given F p/g as

α mid,crit α acc ≡ 3F p/g C η 2 α acc τ -1 s , (22) 4
.76 × 10 -3 F p/g 0.1 2 C η 11/8 -2 α acc 10 -2 τ s 0.1 -1 , (23)
which is independent of r. These two criteria are in very good accordance with the 1D simulations (Sect. 3.1).

Conclusions

In this Letter, we studied how drifting pebbles pile up in a protoplanetary disk with a nonuniform turbulence structure, including the back-reaction of pebbles onto the gas that slows the drift velocity of pebbles as pile-up proceeds. We considered that gas accretion is regulated by an α-parameter that is distinct from that of midplane turbulence. In the disk midplane, the turbulence strength can decrease with decreasing radial distance (e.g., the case where the inner region of the disk is an MRI-inactive deadzone and the outer region is MRI-active). Thus, drifting pebbles are further concentrated in the midplane as the scale height decreases.

We demonstrated a new mechanism for a runaway pile-up of pebbles to occur for a moderate F p/g when pebbles reach a critical level of low turbulence in the dead zone via a continuous slowing of their drift velocity due to the back-reaction (the "nodrift" (ND) runaway pile-up). The ND runaway pile-up occurs irrespective of the shape of the dead zone. Our results imply that SI could also occur as a natural consequence of pebble drift in a protoplanetary disk with a dead zone.

The minimum scale height of pebbles in a dead zone is a critical parameter that characterizes the ND mode. Further studies that include detailed physical processes -such as a KH instability, VSI, vertical shear streaming instability, and the velocity fluctuation by the sound waves propagated from MRI-active surface layers -are required.

Appendix A: The case of K = 0

Here, we show the case of K = 0, where the diffusivity of pebbles is a fixed constant: D p = α mid c 2 s Ω -1 K /(1 + τ 2 s ) (i.e., the backreaction onto the diffusivity is neglected). Figure A.1 is the same as Fig. 2 but for the case of K = 0.

As seen in the case of K = 1 (Fig. 2), a runaway pile-up of pebbles within a dead zone occurs as long as F p/g meets the conditions discussed in Sect. 3.2.2 (the analytical arguments were made for K = 0). Compared to the case of K = 1, the radial width and the timescale of the pile-up are in general narrower and longer in the case of K = 0 because the diffusivity is a constant for K = 0, while it becomes smaller for K = 1 as pile-up proceeds. These additional simulations with K = 0 further show the robustness of the physical mechanism presented in this Letter and demonstrate that Drift-BKR, which slows the drift velocity of pebbles, is a fundamental ingredient for the ND mechanism.

  α mid,crit /α acc K=0, α acc =10 -2 , F p/g =0.17

  α mid,crit /α acc K=0, α acc =10 -2 , F p/g =0.47(d) Analytical α mid,crit /α acc K=0, α acc =10 -2 , F p/g =0.17

  Fig. A.1. Same as Fig. 2 but for the case of K = 0.

We note that H p,KH is a mean square root of the vertical height of particles but not the truncated height of the distribution.

The vertical shear streaming instability (VSSI; Lin 2020) might affect the pebble scale height when the midplane diffusivity is very small. If the VSSI dominated over a KH instability for the scale height of pebbles, the boundary of the ND regime would be regulated by the VSSI.

Acknowledgements. We thank Satoshi Okuzumi and Shoji Mori for helpful comments and fruitful discussion. R.H. acknowledges the financial support of JSPS Grants-in-Aid (JP17J01269, 18K13600). also acknowledges JAXA's International Top Young program. S.I. acknowledges the financial support (JSPS Kakenhi 15H02065, MEXT Kakenhi 18H05438). T.G. was partially supported by a JSPS Long Term Fellowship at the University of Tokyo.