5.4 Zero order term: the Fokker-Planck equation

5.4.1 Momentum dynamics

Fokker-Planck equation In this section, only the electron dynamics in momentum space is discussed, at fixed ψ_l+1∕2. Hence

l+1∕2,i+1∕2,j+1∕2(k+1)	= l+1∕2,i+1∕2,j+1∕2(k+1)
	-×
	l+1∕2,i+1∕2,j+1∕2(k+1)	(5.43)

where

(5.44)

Here, diffusion cross-terms between momentum and configuration spaces are neglected, i.e. D_pψ = D_ψp = D_ξψ = D_ψξ = 0. By definition, the discretization procedure must preserve the conservative nature of the momentum dynamics, as discussed in Sec. 3.5.1, and therefore, careful attention must be paid how to proceed from Eq. 5.4.1, in particular when diffusion or convection coefficients vary strongly on a grid step. Moreover, boundary conditions on flux grid must be as much as possible naturaly satisfied by the discrete form of the Fokker-Planck equation. Indeed, cross-derivatives never satisfy this condition simultaneously along the momentum and pitch-angle directions, and consequently additional external boundary conditions must be added to enforce this condition.

There are several ways to perform the discretization of the cross-derivatives. Since D_pξ and D_ξp terms do not contribute to the Maxwellian solution when no external perturbation is applied, it is not necessary to perform an interpolation between full and half grids for cross-derivatives in order to ensure a correct numerical flux balance in momentum space. In that case, the usual procedure is simply to center differences, as done in this section. However, there is a much more complex approach, but also more consistent one with respect to the use two grids, one for the fluxes and the other for the distribution, which is also presented in Appendix F. However, because of its complexity, it has not been yet implemented in the code even if this approach is in principle better designed for non-uniform grids with strong localized variations of D_pξ and D_ξp.

The starting point of the discretization procedure is the following expressions

	=
	+
	+ pDpξ	(5.45)

	=
	-
	-λDξp	(5.46)

where cross-derivatives terms have been developped directly from analytical formulas, and considered in separately. Here, it is assumed in an implicit manner that derivatives of pD_pξ and λD_ξp exists.

Hence

(5.47)

with

(5.48)

(5.49)

(5.50)

(5.51)

T	= -
	+	(5.52)

T	=
	+	(5.53)

(5.54)

(5.55)

Discrete expressions of the partial derivatives are,

(5.56)

(5.57)

(5.58)

(5.59)

(5.60)

(5.61)

and cross-derivatives

l+1∕2,i+1∕2,j+1∕2(k+1)	= l+1∕2,i+1∕2,j+1∕2(k+1)
	=
	-	(5.62)

while other derivatives become in discrete form

(5.63)

and

	l+1∕2,i+1∕2,j+1∕2
	=
	-	(5.64)

By definition, cross-derivatives are symmetric. Furthermore, the derivative

(5.65)

must be taken on the flux grid in order to fullfil naturally boundary condition at = 1 and avoid to specify D_ξp at this momentum space location. It is a crucial point, especially for the Lower Hybrid and the Ohmic current drive problems, since this domain of the momentum space plays a very important role for the determination of the plasma current level.

Since the distribution function is defined on the half grid, while fluxes on the full grid, it is necessary to interpolate, because in some derivatives, values of f₀ are taken on the full grid. In a general way, whatever the detailed value of the weighting factor δ which will be discussed in the Sec. 5.4.3, one may write for the terms proportional to D_pp and F_p,

f0,l+1∕2,i+1,j+1∕2(k+1)	= f0,l+1∕2,i+3∕2,j+1∕2(k+1)
	+ δp,l+1∕2,i+1,j+1∕2f 0,l+1∕2,i+1∕2,j+1∕2(k+1)	(5.66)

f0,l+1∕2,i,j+1∕2(k+1)	= f0,l+1∕2,i+1∕2,j+1∕2(k+1)
	+ δp,l+1∕2,i,j+1∕2f 0,l+1∕2,i-1∕2,j+1∕2(k+1)	(5.67)

and for terms proportional to D_ξξ and F_ξ

f0,i+1∕2,j+1(k+1)	= f0,l+1∕2,i+1∕2,j+3∕2(k+1)
	+ δξ,l+1∕2,i+1∕2,j+1f 0,l+1∕2,i+1∕2,j+1∕2(k+1)	(5.68)

f0,i+1∕2,j(k+1)	= f0,l+1∕2,i+1∕2,j+1∕2(k+1)	(5.69)
	+ δξ,l+1∕2,i+1∕2,jf 0,l+1∕2,i+1∕2,j-1∕2(k+1)	(5.70)

Gathering all terms, corresponding matrix coefficients for the zero order Fokker-Planck equation may be expressed as

(5.71)

where

Mp,l+1∕2,i+3∕2,j+3∕2	= Dpξ,l+1∕2,i+1∕2,j+1∕2
	+ Dξp,l+1∕2,i+1∕2,j+1∕2	(5.72)

Mp,l+1∕2,i+1∕2,j+3∕2	= Dpξ,l+1∕2,i+1,j+1∕2
	-Dpξ,l+1∕2,i,j+1∕2
	-Dξξ,l+1∕2,i+1∕2,j+1
	- pi+1∕2×
	Fξ,l+1∕2,i+1∕2,j+1	(5.73)

Mp,l+1∕2,i-1∕2,j+3∕2	= -Dpξ,l+1∕2,i+1∕2,j+1∕2
	-Dξp,l+1∕2,i+1∕2,j+1∕2	(5.74)

Mp,l+1∕2,i+3∕2,j+1∕2	= -Dpp,l+1∕2,i+1,j+1∕2
	+ ×
	Dξp,l+1∕2,i+1∕2,j+1
	-×
	Dξp,l+1∕2,i+1∕2,j
	+ Fp,l+1∕2,i+1,j+1∕2	(5.75)

Mp,i+1∕2,j+1∕2	= Dpp,l+1∕2,i+1,j+1∕2
	+ Fp,l+1∕2,i+1,j+1∕2δp,l+1∕2,i+1,j+1∕2
	+ Dpp,l+1∕2,i,j+1∕2
	-Fp,l+1∕2,i,j+1∕2
	+ Dξξ,l+1∕2,i+1∕2,j+1
	- pi+1∕2δξ,l+1∕2,i+1∕2,j+1×
	Fξ,l+1∕2,i+1∕2,j+1
	+ Dξξ,l+1∕2,i+1∕2,j
	+ pi+1∕2×
	Fξ,l+1∕2,i+1∕2,j	(5.76)

Mp,l+1∕2,i-1∕2,j+1∕2	= -Dpp,l+1∕2,i,j+1∕2
	-×
	Dξp,l+1∕2,i+1∕2,j+1
	+ ×
	Dξp,l+1∕2,i+1∕2,j
	-δp,i,j+1∕2F p,l+1∕2,i,j+1∕2	(5.77)

Mp,l+1∕2,i+3∕2,j-1∕2	= -Dpξ,l+1∕2,i+1∕2,j+1∕2
	-Dξp,l+1∕2,i+1∕2,j+1∕2	(5.78)

Mp,l+1∕2,i+1∕2,j-1∕2	= -Dpξ,l+1∕2,i+1,j+1∕2
	+ Dpξ,l+1∕2,i,j+1∕2
	-Dξξ,l+1∕2,i+1∕2,j
	+ pi+1∕2δξ,l+1∕2,i+1∕2,j×
	Fξ,l+1∕2,i+1∕2,j	(5.79)

Mp,l+1∕2,i-1∕2,j-1∕2	= Dpξ,l+1∕2,i+1∕2,j+1∕2
	+ Dξp,l+1∕2,i+1∕2,j+1∕2	(5.80)

For each coefficient, identical diffusion and convection terms appear, which must be evaluated for all physical processes that take place in the plasma (RF, Ohmic electric field, ...),

(5.81)

and

(5.82)

Knock-on source term The knock-on source term must be evaluated at each radial position ψ_l+1∕2, and also at each momentum and pitch-angle positions in momentum space. However, the poloidal position of the birth of the secondary electron plays a crucial role in the calculations, that significantly increases the difficulty of the problem.

As discussed in Sec. 7.2,

(5.83)

where γ = , ξ^*2 = and the poloidal angle θ^* is defined by the relation

(5.84)

.On the discrete mesh, the source term _R at grid point becomes

(5.86)

on each flux surface ψ_l+1∕2. Here

(5.87)

For numerical convenience and for avoiding singularities at boundaries of the integration domain, the source term _{R,l+1∕2,i+1∕2,j+1∕2} must be multiplied by a part of the Jacobian J, i.e.J_ψJ_p, and matrix coefficients in the Fokker-Planck equations are

5.4.2 Spatial dynamics

In this section, only the electron dynamics in configuration space is discussed, at fixed p_i+1∕2 and ξ_0,j+1∕2,

(5.90)

since cross-terms between momentum and configuration spaces are neglected, i.e. D_pψ = D_ψp = D_ξψ = D_ψξ = 0.

Therefore,

(5.91)

with

(5.96)

(5.97)

As for the dynamics in momentum space, interpolation between full and half grids must be performed, since some values of the distribution function f₀ must be determined on the flux grid. Therefore

Gathering all terms, matrix coefficients for the zero order Fokker-Planck equation, according to the relation,

(5.104)

(5.105)

As expected, spatial transport corresponds to a simple tridiagonal matrix.

5.4.3 Grid interpolation

Momentum grid interpolation The usual method for determining interpolation coefficients δ_{p,l+1∕2,i+1,j+1∕2}, δ_{p,l+1∕2,i,j+1∕2}, δ_{ξ,l+1∕2,i+1∕2,j+1} and δ_{ξ,l+1∕2,i+1∕2,j} so that the distribution function f₀ may be estimated on the full grid (flux grid) from its value on the half grid (distribution grid) is to satisfy the constraint that numerically, the stationary solution lim_t→+∞f₀ = f_M is an exact Maxwellian whithout any external accelerating mechanism (RF wave or Ohmic electric field). In this case, under the influence of collisional slowing-down and pitch-angle only, coefficients D_pξ = D_ξp = 0, and

Since this constraint corresponds to the relation

(5.114)

one can deduce,

Hence, the ratios are

(5.131)

(5.132)

(5.133)

(5.134)

with

(5.135)

(5.136)

(5.137)

(5.138)

Furthermore, since ∇_pS_p = 0, one deduces naturally that

(5.139)

(5.140)

and since F_ξ = 0 for collisions, because of the isotropic nature of pitch-angle scattering

(5.141)

(5.142)

which indicates that when collisions is the single physical process considered, f₀^(∞) is independent of ξ₀. Therefore, ∂∕∂p → d∕dp, and integrating the first equation of (5.141) leads to the relation

(5.143)

Values of f₀^(∞) on the respective half grid positions and are then

(5.144)

(5.145)

and ratios are finally easily determined

(5.146)

(5.147)

while, since f₀^(∞) does not depend of ξ₀,

(5.148)

It is then straightforward to evaluate coefficients δ_{p,l+1∕2,i+1,j+1∕2} et δ_{p,l+1∕2,i,j+1∕2} by the trapezoidal method,

Since

Much as the same way, using index conversion p_i+3∕2 → p_i+1∕2, p_i+1 → p_i and p_i+1∕2 → p_i-1∕2 in the previous relation, one obtains for the second integral,

From the two formulations of the ratios

(5.157)

one obtains easily

Much in the same way,

(5.160)

with ß

For a uniform grid, w_{p,l+1∕2,i,j+1∕2}^(nu) = 0, and the well known expression of the Chang and Cooper function is recovered,

(5.162)

with x = w_p^(u). In the limit where x ≪ 1, an expansion of the function gives g = 0.5 - x∕12 + x³∕720 + O with lim_x→0g = , as shwon in Fig.5.2. For a non-uniform grid, a generalized Chang and Cooper function must be introduced, with two arguments

(5.163)

where the term y = w_p^(nu) is zero for the uniform case, as shown in Ref.[?].

With this definition,

(5.164)

where

(5.166)

and

(5.167)

whilewith

(5.168)

and

(5.169)

This expression simplifies and becomes

(5.170)

Furthermore, when plasma relaxation by collisions is the single process considered, the relation

(5.171)

holds², since f₀^(∞) must be Maxwellian f_M, and one finds

(5.172)

Using identities introduced in Sec. 5.2,

(5.173)

it is easy to demonstrate that

(5.174)

and

(5.175)

Finally, since

(5.177)

which is exactly corresponding to a numerical Maxwellian distribution function. Hence, numerical errors do not propagate with the weighting here considered.

Finally, for the non-uniform p grid, the final results for interpolation rules are

(5.178)

with

(5.179)

and

(5.180)

and where

(5.181)

with

(5.182)

and

(5.183)

Much in the same way,

(5.184)

where

(5.185)

and

(5.186)

It is important to note that since only collisions are considered for evaluating interpolating coefficients δ_p, momentum and pitch-angle dynamics are naturally decoupled. Consequently, the ratios F_p∕D_pp are only functions of p, and it is the reason why they haves similar values at different pitch-angle grid points ξ₀.

By definition, δ_p must be lower than unity, a condition that is naturaly satisfied for the uniform momentum grid, as shown in Fig.5.2. However, special care must be taken for the non-uniform case, since there is no exact cancellation between 1∕x and 1∕. In the limit y < x ≪ 1, it can be easily shown that

(5.187)

and the condition g ≤ 1 is equivalent to the relation

(5.188)

For a uniform grid, this condition is always fulfilled, since it corresponds to y = 0. For the non-uniform case, it is clear that the range of validity of the extended Chang and Cooper function is much more restricted, since y∕x is finite. Consequently, if x ≃ pΔp is small, y which is proportional to the variation of the grid step as function of p must be significantly smaller which means that the momentum grid is nearly uniform in this region of the phase space. For larger values of x, the condition is easier to satisfy, which indicates that the non-uniformity must be a growing function of the momentum value p, and nearly flat close to p = 0.

There are however additional limitations on the choice of the momentum for large p values. Indeed, in this domain, other mechanismes are at play, and the usual technique is to extrapolate calculations carried out for collisions to other accelaration mechanismes (Ohmic electric field, RF waves,...), using the generalized weighting factor based on the simple rule,

(5.189)

where ∑ _m is the sum over all the physical processes m that take place in the plasma. Since the same syntax may be kept, one obtains simply

(5.190)

and

(5.191)

for the uniform contribution, while for the non-uniform one

(5.192)

and

(5.193)

In that case, the interpolating weights exhibit a dependence with ξ₀. This has however a weak importance for the accuracy of the calculations, since in the domain where fluxes are strongly modified since usually _effective → 0 in presence of RF quasilinear diffusion. The Maxwellian distribution function is consequently no more the solution of the Fokker-Planck equation.

The reason why an effective expression of F_p∕D_pp is used for evaluating δ_p where collisions are not the single process is based on the fact that the Chang and Cooper function tends towards δ_p ≃ for a uniform mesh when _effective increases. In that case, it corresponds exactly to the standard linear interpolation, i.e. the well known arithmetic mean. The interpolation procedure is therefore consistent with the grid, an important characteristic for reducing the rate of convergence.

Unfortunately, such a property is not valid for a non-uniform grid, and the interpolation may be wrong , leading to possible an anomalous behaviour of the code. Consequently, even at larger values of p, only a uniform mesh may provide a consistent solution with the numerical grid in presence of external sources of acceleration. From this analysis, it turns out that the only purpose for using a non-uniform momentum mesh is to establish a link between the fine mesh in the vicinity of p = 0, and a coarse grid in the region where other physical mechanismes are at play. Consequently, the momentum mesh is build from the relation

(5.194)

with the recurrence relation

(5.195)

where p₀ = 0, Δp_np-1 = p_max∕ and Δp₀ is an arbitrary value. Here i_ref is the index value corresponding to the transition between the fine and the coarse grids. When Δp₀ = Δp_np-1, the case of an uniform mesh is well recovered, whatever i_ref. Here Δp₀ and i_ref must be chosen so that the non-uniform part of the momentum grid is far enough from p = 0 in order the relation δ_p > 1 is satisfied, but also far from the region where external forces play a role³.

Finally, w_p^(nu) require evaluations of the quantities at the respective grid points and . The usual method to deal with this problem is shifting indexes, according to the relation

(5.196)

and

(5.197)

Hence, all quantities have the same size , which is crucial for the 3 - D matrix representation. Furthermore, boundary conditions are naturally satisfied with this technique, since arbitrary values may be attributed to the variable X at these grid points.

Pitch-angle grid interpolation For the pitch-angle terms, calculation is very simple, and from previous relations, one deduces directly that w_{ξ,l+1∕2,i+1∕2,j+1} = w_{ξ,l+1∕2,i+1∕2,j} = 0 since F_ξ = 0. Therefore, δ_{ξ,l+1:2,i+1∕2,j+1} may be defined in an arbitrary manner. It is important to note that this choice is different from the one described in Ref. [?], where the pitch-angle weighting factors are defined in the same way as for the momentum p. However, it turns out that in presence of RF waves, this leads sometimes to spurious numerical evolutions, while the simple approach here described avoid them.

The most natural way is therefore to perform a linear interpolation, between f_{0,l+1∕2,i+1∕2,j+3∕2}^(∞) and f_{0,l+1∕2,i+1∕2,j+3∕2}^(∞), and also between f_{0,l+1:2,i+1∕2,j+1∕2}^(∞) and f_{0,l+1∕2,i+1∕2,j-1∕2}^(∞) according to the relations

Taking into account of the non-uniform grid steps

(5.200)

(5.201)

and for a uniform grid, the relation δ_{ξ,l+1∕2,i+1∕2,j+1} = δ_{ξ,l+1∕2,i+1∕2,j} = 1∕2 are exactly recovered.

Finally, for the non-uniform pitch-angle grid ξ, the interpolation rules are

(5.202)

Furthermore, δ_ξ require evaluations of the quantities at the respective grid points and . As for the momentum grid interpolation, indexes are shifted,

(5.203)

Spatial grid interpolation The determination of δ_ψ is based on the same method used for δ_ξ, since no specific condition is required for interpolating f₀ on the flux grid. A linear interpolation is considered between f_{0,l+3∕2,i+1∕2,j+1∕2}^(∞) and f_{0,l+1∕2,i+1∕2,j+1∕2}^(∞) and also between f_{0,l+1∕2,i+1∕2,j+1∕2}^(∞) and f_{0,l-1∕2,i+1∕2,j+1∕2}^(∞),

(5.206)

(5.207)

For a uniform grid the relation δ_ψ,l+1 = δ_ψ,l = 1∕2 is well recovered.

Furthermore, δ_ψ require evaluations of the quantities at the respective grid points and . As for the momentum grid interpolation, indexes are shifted according to the rule

(5.208)

5.4.4 Discrete description of physical processes

5.4.5 Collisions

According to the bounce averaged expression,

(5.209)

and

(5.210)

where each coefficient is the sum of the electron-electron and electron-ion interactions.

Electron-electron collision operator

Belaiev-Budker relativistic collision model Here, coefficients corresponding to the Belaiev-Budker collision operator that ranges from non-relativistic to fully relativistic regimes may be expressed as, according to the notation used in Ref.[?],

(5.211)

(5.212)

with

(5.213)

(5.214)

Expressions of coefficients F_{11,l+1∕2,i+1∕2}, F_{12,l+1∕2,i+1∕2} and F_{21,l+1∕2,i+1∕2} are

(5.215)

(5.216)

(5.217)

Coefficients A^ee and F^ee must be also calculated on the flux grids and therefore, according to the previous definition,

(5.218)

and

(5.219)

where

(5.220)

(5.221)

and

(5.222)

(5.223)

(5.224)

for the grid points i, while for the grid points i + 1, one has just to replace i by i + 1 in the set of above expressions.

Furthermore, the expression of coefficient B_t^ee is

(5.225)

with

(5.226)

and

(5.227)

where

(5.228)

(5.229)

(5.230)

(5.231)

(5.232)

and

(5.233)

(5.234)

(5.235)

(5.236)

(5.237)

Here,

(5.238)

(5.239)

with

(5.240)

(5.241)

and

(5.242)

Coefficients A, B_t and F for the Beliaev-Budker relativistic collision operator need to calculate accurately integrals of type

(5.243)

which require a special attention in the vicinity of p = 0, and also of type,

(5.244)

The goal is to ensure an acceptable numerical accuracy which preserves the conservative nature of the equations to be solved by numerical method, without any use of ad-hoc boundary conditions to compensate spurious particle leak arisinf from an improper flux balance at each grid point. For this purpose, a new momentum grid called “super-grid” is introduced, correponding to a refined mesh. Since no specific condition is required as compared to the links between distribution and flux grids, the momentum super-grid is simply defined as a sum of two different meshes for integrals of type ∫ ₀^p_i+1∕2 Xdp^′

(5.245)

one for very fine calculations below p_1∕2, and a second, less refined up to p_i+1∕2. For integrals of type ∫ ₀^p_iXdp^′ and ∫₀^p_i+1Xdp^′, the corresponding super-grids are defined in the same way,

(5.246)

and

(5.247)

It is important to note that p₁ correponds to the second point of the grid p_i while it is the first one for the grid p_i+1, so that numerical integration starts at the same momentum value. Much in the same way, integrals of type ∫ _pi^∞Xdp^′ , ∫ _pi+1∕2^∞Xdp^′ and ∫ _pi+1^∞Xdp^′ are simply defined on the following super-grids

(5.248)

All integrals are performed by the trapezoidal method, even if any more accurate technique like the Simpson method may be used in that case. A crucial point is that integrals exactly end or start at points p_i, p_i+1∕2 or p_i+1 so that no overlap between ∫ ₀^p_i+1∕2 Xdp^′ and ∫ _pi+1∕2^∞Xdp^′ can take place. This is especially important to avoid spurious numerical particle leak, that could break the conservative nature of the equations to be solved.

Concerning the first order Legendre correction that ensures momentum conservation, one must calculate

(5.249)

on the distribution function grid, where

(5.250)

as shown in Sec. 4.1.6.

Therefore, by definition,

(5.251)

where for the Belaiev-Budker relativistic collision operator,

(5.253)

and

(5.254)

The set of coefficients ₁ is

(5.255)

(5.256)

(5.257)

(5.258)

(5.259)

(5.261)

(5.262)

(5.263)

(5.264)

and the coefficients ₂ are

(5.265)

(5.266)

(5.267)

(5.268)

(5.272)

(5.273)

and

(5.274)

For this case, integrals are simply calculated according to the rule

(5.275)

so that f_{0,l+1∕2,i+1∕2,j+1∕2} must not be interpolated between 0 and p_i+1∕2 for refined calculations as done for coefficients A^ee, B_t^ee and F^ee. Even if by this techniqe, the numerical accuracy is poor in the vicinity of p = 0, consequences are fairly negligible for the momentum conservation and the current level, since first-order Legendre corrections are weighted by p.

Relativistic Maxwellian background The relativistic Maxwellian background corresponds to that case where the first order Legendre correction for momentum conservation is neglected. Matrix coeffients A^ee, B_t^ee and F^ee determined in the previous section are used.

Non-relativistic Maxwellian background For this case, matrix coefficients are

(5.276)

(5.277)

and

	Bt,l+1∕2,i+1∕2ee =
		(5.278)

where

(5.279)

For values of the matrix coefficients A^ee, B_t^ee and F^ee on the momentum flux grids, one just has to replace i + 1∕2 by i or i + 1.

High velocity limit Though the high velocity limit corresponds to a restricted range of application regarding the full collision operator, it can contribute to useful comparisons with some theoretical calculations. Therefore, this possibility has been implemented in the code. In that case, expressions of the coefficients are greatly simplified,

(5.280)

and

(5.281)

while

(5.282)

Since no integrals appear in the coefficients, the expressions at momentum flux grid points i and i + 1 can be obtained in a straightforward manner, by just replacing i + 1∕2 by the corresponding index values. In that limit, the first-order Legendre correction is neglected. In the calculations, both electron-electron and electron-ion collision model in the high velocity limit are used for a consistent description of the collisions.

Non-relativistic Lorentz model This very simple case corresponds to

(5.283)

and also on the momentum flux grid points i and i + 1. Obviously, the first-order Legendre correction is also neglected in that case.

Electron-ion collision operator

Non-relativistic Maxwellian background In that case matrix coefficients are

	Fl+1∕2,i+1∕2ei = ∑ s ∑ s′
	Zss′2	(5.285)

and

In the case T_{ss^′,l+1∕2} = 0, to avoid a singularity,

(5.287)

while

(5.288)

since erf = 1 and erf^′ = 0.

Finally, for values of the matrix coefficients A^ei,F^ei and B_t^ei on the momentum flux grids, one just has to replace i + 1∕2 by i or i + 1.

High-velocity limit In this limit, matrix coefficients are

(5.289)

(5.290)

The double sum ∑ _s ∑ _s^′ takes into account of all ions species s in ionization state s^′. Here, n_{ss^′,l+1∕2} is the normalized ion density at ψ_l+1∕2, as introduced in Sec. 6.3.1, and m_s is the ion rest mass normalized to the electron rest mass m_e.

Coefficients A^ei and F^ei must be also calculated on the flux grids and therefore, according to the previous definition,

(5.291)

and

(5.292)

for the grid points i, while for the grid points i + 1, one has just to replace i by i + 1 in the set of above expressions.

Furthermore, the expression of coefficient B_t^ei is

(5.293)

Non-relativistic Lorentz model In that limit

(5.294)

while

(5.295)

A straightforward extrapolation may be done for the momentum flux grid i by i + 1.

5.4.6 Ohmic electric field

According to the bounce averaged expression,

(5.296)

and

(5.297)

where E_∥0,l+1∕2 is the parallel component of the Ohmic electric field along the magnetic field line direction normalized to the Dreicer field taken at the poloidal position where the magnetic field B is minimum, as explained in Sec.4.2.

It is important to recall that the Ohmic electric field operator has a cylindrical symmetry, while the description of the electron dynamics in momentum space is based on the spherical symmetry of the collision operator. As a consequence, there is a fundamental contradiction for the internal boundary at p = 0 where p²S_p = 0. Indeed, for large values of E_∥0, the maximum of the distribution function f₀ is no more at p = 0, but may be significantly shifted along the axis p_⊥ = 0. Since in that extrem case p²S_p≠0 at p = 0 while it is naturally enforced to be null by construction with the grid definition, the distribution function has a wrong shape close to p = 0 and the conservative scheme is no more preserved. It is also important to note that the external boundary ∂f₀∕∂ξ₀ = 0 at ξ₀ = 0 is also no more consistent with initial bounce averaged assumptions, which represents also an important failure of the use of the code.

Consequently, in order to avoid a misuse of the code, the range of validity of E_∥0 is restricted so that condition

(5.298)

is fullfiled at p∕p_th = 1. Using the high-velicity limit of the collision operator, this corresponds roughly to

(5.299)

since n_e ≃ 1 in normalized units, as defined in Sec.6.3.1. Consequently, flux surface averaged value of the Ohmic electric field should be restricted to 0.05, in Dreicer units. A warning is indicated when this value is exceeded.

Radio-frequency waves

From the expressions given in Sec. 4.3.7, the components of the tensor D_p^RF are

(5.300)

and

(5.301)

where Ω_{0,l+1∕2,i+1∕2} is the electron cyclotron frequency taken at the minimum B value

(5.302)

Here the dependence of Ω_{0,l+1∕2,i+1∕2} with the index i arises from relativistic down-shift. The indexes n and b correspond respectively to the wave harmonics and the narrow beam label (for ray-tracing calculations). The discrete expression of the quasilinear diffusion coefficient D_{b,n,l+1∕2,i+1∕2,j+1∕2}^RF(0) is

(5.305)

(5.306)

with

(5.307)

and

θb,l+1∕2	= (ψl+1∕2,θb)	(5.308)
rθb,l+1∕2	= r	(5.309)
Rθb,l+1∕2	= R	(5.310)
Bl+1∕2θb	= B	(5.311)
BP,l+1∕2θb	= BP	(5.312)
Ψθb,l+1∕2	= =	(5.313)
ξθb,l+1∕2,j+1∕2	= ξ = σ	(5.314)

All coefficients corresponding to different indexes may be obtained readily by performing the adequate index transformation, i + 1∕2 → and j + 1∕2 →