4.3 Radio frequency waves

The volume-averaged quasilinear diffusion operator for radiofrequency waves in an infinite uniform plasma was first developped by Kennel and Engelmann [?]. The relativistic treatment was performed by Lerche [?], who proposes the following operator

Q	= - limV →∞∑ n=-∞∞d3k	(4.188)

with

P∥	= -	(4.189)
P⊥	= +	(4.190)
Ek,⊥	=	(4.191)

The electric field is assumed to be monochromatic at the frequency ω, and is decomposed into its Fourier components

(4.192)

which are projected on the rotating field frame

Ek,±	=	(4.193)
Ek,∥	= Ek,z	(4.194)

The wave vector is expressed in cylindrical coordinates as

and the argument of the Bessel functions is k_⊥v_⊥∕Ω. The relativistic cyclotron frequency is

(4.196)

4.3.1 Conservative formulation of the RF wave operator

In order to incorporate the RF wave physics in the Fokker-Planck or Drift Kinetic solvers, the operator (4.188) must be cast in a conservative form.

Spherical coordinates representation

The transformation from cylindrical to spherical momentum derivatives is given by

	= +	(4.197)
	= +	(4.198)

Since

p	=	(4.199)
ξ	=	(4.200)

one obtains

	= ξ	(4.201)
	=	(4.202)
	=	(4.203)
	= -	(4.204)

Hence rewritting (4.189-4.190) in cylindrical coordinates using (4.197-4.198)

P∥	= ξ +	(4.205)
P⊥	= -	(4.206)

In order to permute derivatives and the integral in expression (4.188) and obtain a conservative formulation, one must express both P_∥ and P_⊥ in terms of the derivatives associated with the divergence of a flux: 1∕p²p² and -1∕p.

Since

ξ	= p2ξ -	(4.207)
	= p2 -	(4.208)

and

	= - +	(4.209)
-	= -
	+ +	(4.210)

the operators in (4.188) can be expressed in a divergence form

P∥	= p2ξ -	(4.211)
P⊥-	= p2 -	(4.212)

so that (4.188) can be rewritten as

(4.213)

with

SpRF	= limV →∞∑ n=-∞∞d3k
		(4.214)

SξRF	= limV →∞∑ n=-∞∞d3k
		(4.215)

Considering the above expressions of the fluxes and the operators (4.205-4.206), RF-induced fluxes are purely diffusive, and therefore they are expressed as

(4.216)

with

(4.217)

Hence, using the expression of the gradient in spherical coordinates,

SpRF	= -DppRF + DpξRF	(4.218)
SξRF	= -DξpRF + DξξRF	(4.219)

with

DppRF	= -limV →∞∑ n=-∞∞d3k
		(4.220)
DpξRF	= limV →∞∑ n=-∞∞d3k
		(4.221)
DξpRF	= -limV →∞∑ n=-∞∞d3k
		(4.222)
DξξRF	= limV →∞∑ n=-∞∞d3k
		(4.223)

In the limit of a resonant diffusion,

(4.224)

and using the resonance condition

(4.225)

RF diffusion coefficients may be expressed in a simple form

DppRF	= ∑ n=-∞∞D nRF	(4.226)
DpξRF	= ∑ n=-∞∞-DnRF	(4.227)
DξpRF	= ∑ n=-∞∞-DnRF	(4.228)
DξξRF	= ∑ n=-∞∞2D nRF	(4.229)

where we define a diffusion coefficient

(4.230)

Using (4.191), we can define

Dk(n)	= e2
	2	(4.231)

so that (4.230) becomes

(4.232)

If rays are gathered in RF beams of different frequencies ω_b, diffusion coefficients are sums of each contribution ower all hamonics n , leading to the expression

DppRF	= ∑ n=-∞+∞∑ b(1 - ξ2)D b,nRF(p,ξ)	(4.233)
DpξRF	= ∑ n=-∞+∞∑ b -Db,nRF(p,ξ)	(4.234)
DξpRF	= ∑ n=-∞+∞∑ b -Db,nRF(p,ξ)	(4.235)
DξξRF	= ∑ n=-∞+∞∑ b2D b,nRF(p,ξ)	(4.236)

where

(4.237)

and D_k^b,(n) accounts for the polarization and the intensity of the RF wave. It is given by

Dkb,(n)	= e2
	2	(4.238)

Cylindrical coordinates representation

The Fokker-Planck equation is usually solved numerically in spherical coordinates, because of the spherical symmetry of the collisional operator. Therefore, in this coordinate system, the numerical problem is well-conditionned, ensuring stable convergence towards the solution. However, in many case the RF quasilinear operator has a more cylindrical symmetry, and it could be useful to derive its expression in this coordinate system, in order to get a more physical insight of the wave-particle interaction. Starting from the general expression of the flux divergence in cylindrical coordinates,

(4.239)

and taking into account of the diffusive nature of the wave-particle interaction,

(4.240)

where the cylindrical tensor elements are related to the spherical ones by

(4.241)

Applying this transformation, one obtains,

D⊥⊥RF	= ∑ n=-∞+∞∑ b2D b,nRF(p,ξ)	(4.242)
D⊥∥RF	= ∑ n=-∞+∞∑ bDb,nRF(p,ξ)	(4.243)
D∥⊥RF	= ∑ n=-∞+∞∑ bDb,nRF(p,ξ)	(4.244)
D∥∥RF	= ∑ n=-∞+∞∑ b2D b,nRF(p,ξ)	(4.245)

For n = 0 (which corresponds to the Lower Hybrid wave) the quasilinear diffusion is strictly along the parallel direction (i.e. magnetic field line). Also, at a cyclotron harmonic, where ω_b = nΩ, the diffusion is only perpendicular, when relativistic corrections are not considered.

4.3.2 RF Diffusion coefficient for a Plane Wave

4.3.3 Integration in k-space

The quasilinear diffusion coefficient (4.230), describing the interaction of the electrons with a given beam b at an harmonic n, is rewritten as

(4.246)

where,

Θkb,	= ebk,+e-iαJ n-1 + ebk,-e+iαJ n+1
	+ ebk,∥Jn	(4.247)

Here, the polarization vector in Fourier space

(4.248)

is introduced, whose components in the rotating field frame are

ebk,+	= =	(4.249)
ebk,-	= =
ebk,∥	= =

The electric field associated with a given ray is described by a plane wave, with given frequency ω_b and wave number k_b:

Eb	= Re = Re	(4.250)
	=	(4.251)

so that

Ebk	≡d3xe-ik⋅xE b
	=
	=	(4.252)

The electric field is the sum of two contributions with different frequencies

Ebk+	= Eb0δ with ω = ωb	(4.253)
Ebk-	= Eb0*δ with ω = -ωb	(4.254)

Therefore the QL diffusion tensor (4.217) is the sum of these two contributions:

(4.255)

associated with the diffusion coefficients

Db,nRF+(p,ξ)	= lim V →∞22δ	(4.256)
Db,nRF-(p,ξ)	= lim V →∞22δ

with, using (4.247)

Θkb+,	= eb+k,+e-iαJ n-1 + eb+k,-e+iαJ n+1
	+ eb+k,∥Jn	(4.257)
Θkb-,	= eb-k,+e-iαJ n-1 + eb-k,-e+iαJ n+1
	+ eb-k,∥Jn

The polarization components (4.249) are

eb+k,+	= = eb0,+	(4.258)
eb+k,-	= = eb0,-
eb+k,∥	= = eb0,∥

and

eb-k,+	= = eb0,-*	(4.259)
eb-k,-	= = eb0,+*
eb-k,∥	= = eb0,∥*

Expressing the condition k = ±k_b from the first delta function in (4.256), we find

Db,nRF+(p,ξ)	= lim V →∞2δ22δ	(4.260)
Db,nRF-(p,ξ)	= lim V →∞2δ22δ

with now (4.257) being

Θkb+,	= eb0,+e-iαb Jn-1 + eb0,-e+iαb Jn+1
	+ eb0,∥Jn	(4.261)
Θkb-,	= eb0,-*e-iJ n-1 + eb0,+*e+iJ n+1
	+ eb0,∥*J n
	= -eb0,-*e-iαb Jn-1-eb0,+*e+iαb Jn+1
	+ eb0,∥*J n	(4.262)

where we used the fact that k_b is transformed into -k_b by changing k_b∥ to -k_b∥ and α_b to α_b + π.

Using the properties of the Bessel functions, and defining n^′ = -n,

Jn+1	= n+1J n′-1	(4.263)
Jn-1	= n-1J n′+1
Jn	= nJ n′

so that (4.262) becomes

Θkb-,	= n
		(4.264)
	= n*	(4.265)

and the diffusion coefficients become

Db,nRF+(p,ξ)	= lim V →∞22δδ2	(4.266)
Db,nRF-(p,ξ)	= lim V →∞22δδ2

Using Parseval’s theorem,

limV →∞δ2	= limV →∞d3x = 1
limV →∞δ2	= limV →∞d3x = 1	(4.267)

one obtains

Db,nRF+(p,ξ)	= e2π2δ	(4.268)
Db,nRF-(p,ξ)	= e2π2δ	(4.269)

In the expression (4.234-4.236) for D^RF- tensor elements, we subsitute

(4.270)

Then, by redefining n^′→ n in the sum over all harmonics for D^RF-, we can combine the two contributions of D^RF+ and D^RF-, and finally we obtain an expression with one diffusion coefficient:

Db,nRF(p,ξ)	= e2π2δ
	= e2π2δ	(4.271)

with (4.261)

(4.272)

Nb∥	=	(4.273)
N∥res	=
	=	(4.274)

where

Ω	= = with ωce =	(4.275)
n′	= -n	(4.276)
βTe	= = =	(4.277)

4.3.4 Incident Energy Flow Density

Relation to the Electric Field

The time-averaged energy flow density in the beam is in general given by

(4.278)

where

(4.279)

is the flow of the electromagnetic energy or Poynting vector, and

(4.280)

is the kinetic energy flow density where K^ℍ is the Hermitian part of the dielectric tensor.

_b is the complex form of the electric field (4.250)

(4.281)

and the magnetic field _b is given by

(4.282)

using the relation N_b = k_bω_b∕c.

The energy flow density (4.278) can be formally rewritten as

(4.283)

where Φ_b is a non dimensional vector defined as

(4.284)

with

(4.285)

and

(4.286)

Using (4.282), we obtain

ΦbP	= Re
	= Re	(4.287)

In vacuum, Φ_b = , unit vector in the direction of propagation, and Φ_bT = 0. Here,

(4.288)

is the polarization vector. The incident power flux on a flux-surface is therefore given by

(4.289)

where is the local vector normal to the flux-surface.

Relation to the ray trajectories

The energy flows in the direction of the group velocity, noted

(4.290)

In ray-tracing calculations, this is the direction of the ray trajectories, which are parametrized by where s is the distance along the ray. Therefore _b can be determined from ray-tracing calculations or any other wave propagation model, and we have

(4.291)

In DKE calculations, flux-surfaces are discretized on a grid ψ_l+1∕2 and have a finite volume (). In consequence, the factor must be ”integrated” along the ray trajectory within this volume. We therefore define an incidence factor

(4.292)

where the ± sign accounts for the fact that the ray may be locally directed either inward or outward, and where Δl_l+1∕2^θ_b is the local thickness of the flux-surface.

(4.293)

The factor f_inc,b^l+1∕2 can be numerically calculated from the equilibrium and ray-tracing data. In the limit of an infinitely thin flux-surface, we have

(4.294)

We therefore substitute

(4.295)

4.3.5 Narrow Beam Approximation

Let consider a ray which crosses the flux surface at the respective poloidal locations θ_b. Here, it is assumed that all electrons interact with the beam, whatever its toroidal location, because of the axisymmetry. Therefore, calculations are only valid for circulating electrons, except for those localized on low n order rational q-surfaces, since their trajectories are exactly periodic in configuration space. Consequently, either they may never crosses the beam or if they cross it, the quasilinear assumption might likely fail, since some correlations could remain between two momentum kicks. However, the total weight of these electrons is marginal, and their contribution to RF power absorption and current generation is simply neglected. It is important to recal that trapped particle interaction with RF wave is not addressed in this code for similar arguments. For all particles whose trajectory is periodic, a specific treatment of the quasilinear interaction is necessary, which is at this stage beyond the scope of this development.

If we assume that the extension of the beam power on the flux-surface is very small in the poloidal directions, we can approximate

(4.296)

This assumption is in contradication with the plane-wave decomposition of the RF beam, since plane waves are essentially present in the entire configuration space. However, it is assumed that plane waves interfere destructively everywhere except in the beam region, which can be limited to a narrow poloidal and toroidal extention. By definition, _b is determined so that P_b,inc is the total incident power on the flux-surface

(4.297)

which gives in the system , using (A.199)

(4.298)

Hence,

(4.299)

and the incident power flux is simply

(4.300)

Using (4.295) and (4.300),

(4.301)

For concentric circular flux surfaces, = , then = 1. We obtain

(4.302)

where Ψ_p = Ψ_p = 1 + ϵ.

4.3.6 Normalized Diffusion Coefficient

The diffusion coefficient is usually normalized to the collisional diffusion coefficient ν_ep_Te², thus defining

(4.303)

with

(4.304)

Using (4.271) and (4.301), the diffusion coefficient (4.303) becomes

(4.305)

with

Db,n,0RF	= Pb,inc
	= Pb,inc	(4.306)

where ω_pe is the electron plasma frequency

(4.307)

We also recall the expressions (4.272), (4.274) and (4.284-4.287)

(4.308)

(4.309)

Φb	= ΦbP + ΦbT	(4.310)
ΦbP	= Re	(4.311)
ΦbT	= -eb*⋅ ⋅ eb	(4.312)

For concentric circular flux surfaces,

(4.313)

4.3.7 Bounce Averaged Fokker-Planck Equation

In the Fokker-Planck equation, the diffusion and convection elements are bounce-averaged according to the expressions (3.189) and (3.194), which gives, using (4.233-4.236)

DppRF(0)	=	(4.314)
DpξRF(0)	= σ
DξpRF(0)	= σ
DξξRF(0)	=

where ()

(4.315)

Since 1 - ξ² = Ψ,

DppRF(0)	= ∑ n=-∞+∞∑ b(1 - ξ02)D b,nRF(0)(p,ξ 0)	(4.316)
DpξRF(0)	= ∑ n=-∞+∞∑ b -Db,nRF(0)(p,ξ 0)
DξpRF(0)	= ∑ n=-∞+∞∑ b -Db,nRF(0)(p,ξ 0)
DξξRF(0)	= ∑ n=-∞+∞∑ b2D b,nRF(0)(p,ξ 0)

where

(4.317)

and Ω₀ is the cyclotron frequency taken at the minimum B value, according to the relation

(4.318)

with

(4.319)

Therefore, only the bounce averaged diffusion coefficient D_b,n^RF(0)(p,ξ₀) has to be calculated for the different types of RF waves, keeping the same formalism. It is determined by inserting expression (4.305) into (4.317),

(4.320)

which gives

(4.321)

We perform the poloidal integration

	= T ∫ θminθmax dθΨDb,n,0RF2δδ	(4.322)
	= ΨθbDb,n,0RF,θb HHT δ2	(4.323)

with

θb	= (ψ,θb)	(4.324)
rθb	= r	(4.325)
Rθb	= R	(4.326)
Bθb	= B	(4.327)
BP θb	= BP	(4.328)
Ψθb	= =	(4.329)
ξθb	= ξ = σ	(4.330)

and where H is the usual Heaviside step function

Therefore,

Db,nRF(0)(p,ξ 0)	= ΨθbDb,n,0RF,θb ×
	HHT δ2	(4.331)

where we define, using (4.308), (4.309) and (4.306),

(4.332)

(4.333)

(4.334)

with

zbθb	=
	= -Nb⊥
	= -Nb⊥	(4.335)

4.3.8 Bounce Averaged Drift-Kinetic Equation

In the Drift Kinetic equation, the diffusion and convection elements are bounce-averaged according to the expressions (3.218) and (3.223), which gives, using (4.233-4.236)

ppRF(0)	= σ	(4.336)
pξRF(0)	=
ξpRF(0)	=
ξξRF(0)	= σ

and

pRF(0)	=	(4.337)
ξRF(0)	= σ

where

(4.338)

Since 1 - ξ² = Ψ and Ω = ΨΩ₀,

ppRF(0)	= ∑ n=-∞+∞∑ b(1 - ξ02) b,nRF(0)D(p,ξ 0)	(4.339)
pξRF(0)	= ∑ n=-∞+∞∑ b -b,nRF(0)D(p,ξ 0)	(4.340)
ξpRF(0)	= ∑ n=-∞+∞∑ b -b,nRF(0)D(p,ξ 0)	(4.341)
ξξRF(0)	= ∑ n=-∞+∞∑ b2 b,nRF(0)D(p,ξ 0)	(4.342)

pRF(0)	= ∑ n=-∞+∞∑ b -b,nRF(0)F(p,ξ 0)	(4.343)
ξRF(0)	= ∑ n=-∞+∞∑ b2 b,nRF(0)F(p,ξ 0)

where

b,nRF(0)D(p,ξ 0)	= σ	(4.344)
b,nRF(0)F(p,ξ 0)	= σ

and Ω₀ is the cyclotron frequency taken at the minumum B value, as defined for the Fokker-Planck equation. Therefore, only the two bounce averaged diffusion coefficient _b,n^RF(0))D(p,ξ₀) and _b,n^RF(0)F(p,ξ₀) have to be calculated for the different types of RF waves, keeping also the same formalism. It is important to recall that only = σξ₀ can be put outside the corresponding of the bounce-averaging integral.

These coefficient are determined by inserting expression (4.305) into (4.317),

b,nRF(0)D(p,ξ 0)	= σ	(4.345)
b,nRF(0)F(p,ξ 0)	= σ

which gives

b,nRF(0)D(p,ξ 0)	= σ	(4.346)
b,nRF(0)F(p,ξ 0)	= σ

We perform the poloidal integration

	σ
	= σT ∫ θminθmax dθσDb,n,0RF2δδ
	= Db,n,0RF,θb HHσT σδ2	(4.347)

	σ
	= σT ∫ θminθmax dθσDb,n,0RF
	×2δδ
	= Db,n,0RF,θb HH
	× σT σδ2	(4.348)

where H is the usual Heaviside step function and D_b,n,0^RF,θ_b is given by (4.331).

Therefore, we find

b,nRF(0)D(p,ξ 0)	= Db,n,0RF,θb ×	(4.349)
	HHσT σδ2

b,nRF(0)F(p,ξ 0)	= Db,n,0RF,θb ×	(4.350)
	HHσT σδ2

4.3.9 Modeling of RF Waves

The quasilinear diffusion coefficients (4.331), (4.349) and (4.350) describe the interaction between electrons and a discrete set (b) of RF plane waves (4.250) with frequencies ω_b. In addition, we can linearly superpose a discrete set of plane waves with the same frequency but different parallel wave number k_∥ because of the linearity of the QL operator (4.188) in k. This set of plane waves is typically represented by rays. Each ray is characterized by: the wave frequency ω_b, the poloidal location θ_b, the index of refraction N_b, the polarization vector e_b and the power flow Φ_b. In general, the ray propagation path and the evolution of the index of refraction N_b are determined by ray-tracing (RT) calculations. When the wave properties are determined by RT calculations, the launched power spectrum in k_∥ is typically decomposed in an array of rays with different poloidal launching angle and/or k_∥ power spectrum. This distribution depends on the antenna and other launching parameters. The evolution of each ray is determined separately, by the RT code, and the contributions of each ray to wave-particle interaction are separately accounted for in the sum over b in (4.233-4.236). The separation of these contributions is justified if their k_∥ power spectrums do not overlap.

We do not address the RT problem in this work. Therefore, we assume that the path θ_b and the parallel index of refraction N_b∥ are determined either by coupling of DKE to a RT code, or by any other propagation model. RT calculations require to calculate the dispersion tensor and solve the dispersion relation, and therefore they can provide input for the remaining wave properties . However, these properties can also be determined by solving locally the wave equation, once ψ,θ_b and N_∥ are known. Keeping the calculation of within DKE allows us to use a more advanced model - such as the fully relativistic dispersion solver R2D2 - for the evaluation of these important wave properties. A numerical code is required to solve the wave equation, and evaluate the diffusion coefficient (4.331) in the general case. However, it is possible to obtain analytic expressions in some cases, assuming that the waves can be treated in the cold plasma limit. The calculation of D_b,n^RF(0)(p,ξ₀) in the cold plasma model is performed in Appendix D. In addition, simplified expression and the comparison with operators used in the litterature are given for lower-hybrid (LH) and electron-cyclotron (EC) waves in appendices D.2 and D.3 respectively.

Note that the cold plasma model, if used in the RT calculations, does not account for the wave-particle interaction and the power absorption from the rays. The power absorbed from each ray and deposited on electrons can be evaluated if a hot plasma model is used, but even then, the calculation is not consistent with DKE because quasilinear effects are not accounted for. In order to have a consistency between RT and DKE calculations, the power absorption calculated by DKE must be inserted back in the RT calculations. In the case where rays turn back toward the edge (as it can happen in LHCD), iterations between RT and DKE calculations are necessary to ensure self-consistency.

Modeling of RF k_∥ spectrum

We consider two different models for the k_∥ spectrum of a given ray.

Square power spectrum in k_∥ Here, the power spectrum of a given ray is taken as being constant in N_∥ between two limits N_{b∥ min} and N_{b∥ max}. Therefore, we operate the following transformation

(4.351)

where

(4.352)

In this case, the RT calculations and wave equation are solved for the central value.

(4.353)

which is a good approximation only when

(4.354)

For simplification and benchmarking purposes, in LHCD, one can consider only one ray with a square spectrum. Then, the limits N_{b∥ min} and N_{b∥ max} are respectively given by the accessibility condition and the condition for strong linear Laudau damping.

Gaussian power spectrum in k_∥ The power spectrum of a given ray can be assumed to have a Gaussian dependence on N_∥, centered around a value N_b∥,0 and with a Gaussian width ΔN_b∥. Then, we operate the transformation

(4.355)

In this case too, the RT calculations and wave equation are solved for the central value N_b∥,0, which is a good approximation only when

(4.356)

For simplification and benchmarking purposes, in ECCD, one can consider only one ray with a Gaussian spectrum. Then, the wave properties N_b∥,0 and ΔN_b∥ are determined by the beam characteristics, as described in 4.3.9.

Beam size and power spectrum

If we consider a Gaussian beam of waist d, of negligible diffraction angle, and of negligible dispersion, such that the wave vector is in the direction of the beam, and N ≃ 1, we get that by the properties of Fourier transform, the width ΔN_∥ is related to the beam waist d as

(4.357)

This conditions are approximately satisfied for EC beamsif the density is low (ω_pe ≪ ω).

Then, and we can verify that for a beam waist of a few centimeters (d = 5 cm, f = 110 GHz, N_∥≃ 0.2 → ΔN_∥≃ 0.01), the condition (4.356) is satisfied, and there is no nead to decompose the beam spectrum into many rays.