Polarization density matrix

6 Polarization density matrix

We need to perform the gyro-average and velocity integration in the polarization density (expression (230)) to obtain a matrix form that relates the polarization density at gridpoints to δΦ at gridpoints, so that the corresponding matrix form of the left-hand side of Poisson’s equation (231) can be inverted to solve for δΦ.

The dependence of the ﬁrst term in expression (230) (often called adiabatic term) on δΦ is trivial (because δΦ is independent of the velocity in particle coordinates):

∫ -q ∂F0- δnad = m (δΦ) ∂𝜀 dv q ∫ ∂F0 = m-δΦ -∂𝜀 dv. (232)

If we assume F₀ is Maxwellian, then the velocity integration can be analytically performed:

q ∫ ( m ) δnad = --δΦ −--fM dv. m T = − qδΦ-n0. (233) T

Next, let us consider the more complicated term, the second term in expression (230), i.e.,

q ∫ ( ∂F0) A (x) ≡ − m dv ⟨δΦ⟩-∂𝜀- .

(234)

At particle location x, we construct local (velocity) cylindrical coordinates (v_⊥,v_∥,α′). Then the velocity element is given by dv = v_⊥dv_⊥dv_∥dα. Then the above integration is written as

∫ ∞ ∫ ∞ ∫ 2π A(x) = −-q dv dv ∂F0v dα⟨δΦ⟩. m − ∞ ∥ 0 ⊥ ∂𝜀 ⊥ 0

(235)

Note that ∂F₀∕∂𝜀 (in terms of particle coordinates) dpends on α. But this dependence is usually weak and neglected. So it is moved outside of the α integration. Note that the above is a velcoity space integration for a ﬁxed particle position.

Using the deﬁnition of the gyro-averaging ⟨…⟩ = (2π)⁻¹ ∫ ₀^2π(…)dα′, expression (235) is written as

( ) -q ∫ ∞ ∫ ∞ ∂F0- ∫ 2π -1-∫ 2π ′ A (x) = −m −∞ dv∥ 0 dv⊥ ∂𝜀 v⊥ 0 dα 2π 0 δΦ(x)dα .

(236)

Here the gyro-averaging can be understood as time integration of the gyro-orbit that has the point x on its trajectory. The gyro-phase when the particle is passing x is α. So its guiding-center location X can be calculated as

e∥ X = x + v⊥(v⊥,α)× --. Ω

(237)

Using X, the full gyro-ring is written as (parameterized by α′)

x′ = X − v⊥(v⊥,α′)× e∥, Ω

i.e.,

′ e∥ ′ e∥ x = x + v⊥(v⊥,α) × Ω − v⊥ (v⊥,α )× Ω ,

(238)

Then expression (236) is written as

∫ ∞ ∫ ∞ A (x) = − q dv∥ dv⊥∂F0-v⊥ m −∞ 0 ∂𝜀 ∫ 2π [ 1 ∫ 2π ( e∥ ′ e∥) ′] × dα 2π- δΦ x + v⊥(v⊥,α)× Ω- − v⊥(v⊥,α )× Ω- dα . (239) 0 0

To proceed, we have two options. One is to assume no analytical expression of δΦ(x) and use pure numerical interpolation to handle the above integration, and ﬁnally expression it as combination of δΦ values on gridpoints. Another option is to assume basis function expansion of δΦ and perform analytically the integration for each basis function. The ﬁrst option is further discussed in Sec. 6.2. Next, we focus on the second option and use Fourier expansion.

6.1 In terms of Fourier basis functions

We Fourier expand δΦ(x) as

∫ -dk-- δΦ(x) = δΦk exp(ik ⋅x)(2π)3,

(240)

where δΦ_k ≡∫ δΦ(x)exp(−ik ⋅ x)dx is the expansion coeﬃcients. Then expression (239) is written as

∫ ∫ ∞ ∫ ∞ A(x) = − q --dk3δΦk exp(ik ⋅x) dv∥ dv⊥∂F0-v⊥ ∫ m (2π) −∞ 0∫ ∂𝜀 2π ( ( e∥) -1- 2π ′ ( ( ′ e∥) ) × 0 dα exp ik⋅ v ⊥(v⊥, α)× Ω 2π 0 dα exp − ik ⋅ v⊥(v⊥,α )× Ω

Let k_⊥ deﬁne one of the perpendicular direction ₁, i.e., k_⊥ = k_⊥₁. Then another perpendicular basis vector is deﬁned by ₂ = e_∥×₁. Then v_⊥ is written as v_⊥ = v_⊥(₁ cosα′ + ₂ sinα′), which deﬁnes the gyro-angle α′. Then the blue term is written as

e∥ ′ ′ e∥ − ik ⋅v⊥ × Ω-= − ik⋅v⊥(ˆe1cosα + ˆe2sinα )× Ω- v⊥- ′ ′ = − ik⋅ Ω (− ˆe2cosα + ˆe1sinα ) k⊥v⊥- ′ = − i Ω sinα . (241)

Then the α′ integration is written as

∫ ( ) -1- 2π dα′exp − ik⊥v⊥-sinα ′ 2π 0 Ω ( k⊥v⊥) = J0 -Ω--- . (242)

where use has been made of the deﬁnition of the zeroth Bessel function of the ﬁrst kind.

Similarly, the α integration is also reduced to J₀ (k⊥v⊥)
Ω . Then A(x) is written as

q∫ dk [∫ ∞ ∫ ∞ ∂F0 (k⊥v ⊥) ] A(x) = − m (2π)3-δΦk exp(ik ⋅x) dv∥ dv⊥ ∂𝜀-v⊥2πJ20 --Ω-- . (243) −∞ 0

The remaining velocity integration can be performed analytically if F₀ is Maxwellian

If F₀ is a Maxwellian distribution, then the remaining velocity integration in Eq. (243)can be analytically performed. Maxwelliand distribution is given by

( ) ----n0(X-)----- −-mv2- F0 = fM = (2πT (X)∕m )3∕2 exp 2T (X ) (244) n ( − v2) = ----30∕23 exp --2- , (245) (2π) vt 2vt

where v_t = ∘ ----
T∕m , then

∂F0- m- ∂𝜀 = − T fM .

(246)

As is mentioned above, we ignore the weak dependence of n₀(X) and T(X) on v (introduced by X = x + v × e_∥∕Ω) when transformed back to particle coordinates. (For suﬃciently large velocity, the corresponding Larmor radius will be large enough to make the equilibrium undergo substantial variation. Since the velocity integration limit is to inﬁnite, this will deﬁnitely occur. However, F₀ is exponentially decreasing with velocity, making those particles with velocity much larger than the thermal velocity negligibly few and thus can be neglected.)

Parallel integration Using Eq. (246), the expression in the square brackets of Eq. (243) is written as

∫ ∫ ( k⊥v⊥) ∂F0 2π dv⊥dv ∥J20 -Ω--- -∂𝜀-v⊥ ∫ ∫ ( ) ( 2 2 ) m---n0-- 2 k⊥v-⊥- -1 v∥ +-v⊥ = − T(2π)1∕2 dv⊥dv ∥J 0 Ω v3t exp − 2v2t v⊥ (247) ∫ ∫ ( ) ( -2 -2) = − m--n0-- dv⊥dv-∥J2 k⊥v⊥- exp − v∥ +-v⊥ v⊥, (248) T(2π)1∕2 0 Ω 2

where v_∥ = v_∥∕v_t, v_⊥ = v_⊥∕v_t. Using

∫ ∞ ( 2) √ --- exp − x- dx = 2π, −∞ 2

(249)

the integration over v_∥ in expression (248) can be performed, yielding

m ∫ ∞ - ( k⊥v⊥) ( v2 )- − T-n0 dv⊥J20 -Ω--- exp −-⊥2- v⊥ 0

(250)

Perpendicular integration Using (I veriﬁed this by using Sympy)

∫ ( ) ∞ 2 x2 2 2 0 J0(ax)exp − 2 xdx = exp(− a )I0(a ),

(251)

where I₀(a) is the zeroth modiﬁed Bessel function of the ﬁrst kind, expression (250) is written

− m-n0exp(− b)I0(b) T

(252)

where b = k_⊥²v_t²∕Ω² = k_⊥²ρ_t². Then the corresponding density (234) is written as

∫ A(x) = qn0 δΦk exp(ik ⋅x)exp(− b)I0(b)-dk-. T (2π)3

(253)

In Fourier space, the adiabatic term in expression (233) is written as

∫ ∫ q- ∂F0- qn0 qn0 -dk-- m(δΦ) ∂𝜀 dv = − T δΦ = − T δΦk exp(ik ⋅x)(2π)3 .

(254)

Plugging expression (253) and (254) into expression (230), the polarization density n_p is written as

qn0∫ -dk-- np = − T δΦk exp(ik ⋅x)[1− exp(− b)I0(b)](2π )3. (255)

Deﬁne

Γ 0 = exp(− b)I0(b),

(256)

then Eq. (255) is written as

∫ n = − qn0 δΦ exp(ik ⋅x)[1 − Γ ]-dk-, p T k 0(2π)3

(257)

Expression (257) agrees with the result given in Yang Chen’s notes. Note that the dependence on species mass enters the formula through the Larmor radius ρ_t in Γ₀.

Pade approximation

Γ₀ deﬁned in Eq. (256) can be approximated by the Pade approximation as

Γ ≈ -1--. 0 1+ b

(258)

The comparison between the exact value of Γ₀ and the above Pade approximation is shown in Fig. 1.

pict

Fig. 1: Comparison between the exact value of Γ₀ = exp(−(k_⊥ρ)²)I₀((k_⊥ρ)²) and the corresponding Pade approximation 1∕(1 + (k_⊥ρ)²).

Using the Pade approximation (258), the polarization density n_p in expression (257) can be written as

∫ n ≈ − qn0 δΦ exp(ik⋅x)--k2⊥ρ2---dk--. (259) p T k 1 + k2⊥ρ2(2π)3

Long wavelength approximation of the polarization density

In the long wavelength limit, k_⊥ρ ≪ 1, expression (259) can be further approximated as

∫ np ≈ − qn0 δΦk exp(ik ⋅x)k2⊥ρ2-dk-3, T (2π) = qn0 ρ2∇2⊥δΦ. (260) T

Then the corresponding term in the Poisson equation is written as

q-n = q2n0-ρ2∇2δΦ 𝜀0 p 𝜀0T ⊥ ρ2- 2 = λ2∇ ⊥δΦ, (261) D

where λ_D is the Debye length deﬁned by λ_D² = T𝜀₀∕(n₀q²). For typical tokamak plasmas, the thermal ion gyroradius ρ_i is much larger than λ_D. Therefore the term in expression (261) for ions is much larger than the space charge term ∇²δΦ ≡∇_⊥²δΦ + ∇_∥²δΦ ≈∇_⊥²δΦ in the Poisson equation. Therefore the space charge term can be neglected in the long wavelength limit.

Equation (261) also shows that electron polarization density is smaller than the ion polarization density by a factor of ρ_e∕ρ_i ≈ 1∕60. Note that this conclusion is drawn in the long wavelength limit. For short wavelength, the electron polarization and ion polarization density can be of similar magnitude (to be discussed later).

Polarization density expressed in terms of Laplacian operator

The polarization density expression (260) is for the long wavelength limit, which partially neglects FLR eﬀect. Let us go back to the more general expression (259). The Poisson equation is written

2 − 𝜀0∇⊥ δΦ = qiδni + qeδne.

(262)

Write δn_i = n_pi + δn_i′, where δn_pi is the ion polarization density, then the above expression is written

− 𝜀0∇2⊥δΦ − qinpi = qiδn′i + qeδne.

(263)

Fourier transforming in space, the above equation is written

− 𝜀0k2⊥δˆΦ − qiˆnpi = qiδˆn′i + qeδˆne,

(264)

where _pi is the Fourier transformation (in space) of the polarization density n_pi and similar meanings for δ , δ_i′, and δ_e. Expression (259) implies that _pi is given by

qn k2ρ2 ˆnpi = − i-i0δˆΦ--⊥-i2-2. Ti 1 + k⊥ρi

(265)

Using this, equation (264) is written

( ) 2 ˆ qini0--k2⊥-ρ2i-- ˆ ′ − 𝜀0k⊥ δΦ− qi − Ti 1+ k2⊥ρ2i δΦ = qiδˆni + qeδnˆe,

(266)

Multiplying both sides by (1 + k_⊥²ρ_i²)∕𝜀₀, the above equation is written

( ) − (1 + k2⊥ρ2i)k2⊥δˆΦ− qi − qini0(k2⊥ρ2i)δˆΦ = -1(1+ k2⊥ρ2i)(qiδˆn′i + qeδˆne). 𝜀0 Ti 𝜀0

(267)

Next, transforming the above equation back to the real space, we obtain

( ) 2 2 2 qi qini0-2 2 -1 2 2 ′ − (1− ρi∇⊥ )∇ ⊥δΦ − 𝜀0 Ti ρi∇ ⊥δΦ = 𝜀0(1− ρi∇ ⊥)(qiδni + qeδne).

(268)

Neglecting the Debye shielding term, the above equation is written

( 2 ) − ρi2-∇2⊥δΦ = -1(1− ρ2i∇2⊥)(qiδn′i + qeδne), λDi 𝜀0

(269)

which is the equation actually solved in many gyrokinetic codes, where λ_Di² = 𝜀₀T_i∕(q_i²n_i0).

6.2 Polarization density obtained by numerical integration

In Sec. 6, to evaluate the polarization density, the potential δΦ is Fourier expanded in space, and then the double gyro-angle integration of each harmonic is expressed as the Bessel function. (This method is often used in numerical codes that use specrum expansion in both radial and toroidal direction, where the local perpendicular wave number can be easily calculated numerically. The GEM code uses this method.) In this section, we avoid using the local Fourier expansion, and directly express the double gyro-angle integral as linear combination of values of δΦ at spatial grid-points.

6.3 Direct evaluation of the double gyrophase integration

Then A(x) is written as

A(x) q-∫ ∞ ∫ ∞ ∂F0- ≈ − m −∞ dv∥ 0 dv⊥ ∂𝜀 v⊥ N1 N2 ( ) × 2π-∑ -1-∑ δΦ x + v⊥(v⊥,αi)× e∥ − v⊥(v⊥,α′)× e∥ , (270) N1 i=1 N2 j=1 Ω j Ω

where N₁ is the number of points chosen for α discretization, N₂ is the number of points chosen for α′ discretization. For notation ease, deﬁne

′ e∥ Δ ρij = [v⊥ (v⊥,αi)− v⊥ (v⊥, αj)]× -Ω ,

(271)

which is a function of (x,v_⊥,α_i,α_j′). Then Eq. (270) is written as

∫ ∞ ∫ ∞ N1 N2 A (x ) = − q dv∥ dv⊥∂F0-v⊥2π-∑ -1-∑ δΦ (x + Δρij). (272) m −∞ 0 ∂𝜀 N1 i=1N2 j=1

The guiding-center transform and its inverse used in the above are illustrated in Fig. 2.

Fig. 2: The guiding-center transform and its inverse used in Eq. (270). Here x and v_⊥ are given.

Let us summarize the above algorithm. Given v_⊥, the double gyro-angle integration at a particle location x is evaluated by the following steps:

Use the guiding center transform to get guiding-center locations (four locations are shown in Fig. 2, namely X_i with i = 1,2,3,4, corresponding gyro-angle α_i with i = 1,2,3,4;
For each guiding-center location, use the inverse guiding-center transform to calculate points on the corresponding gyro-ring (four points are shown for each guiding-center X_i in this case, namely x_ij with j = 1,2,3,4, corresponding gyro angle α_j′ with j = 1,2,3,4.)
Then the double gyro-angle integral appearing in the polarization density is approximated as

$∫ 2π ( ∫ 2π ) ∑N1∑N2 dα -1- δΦ(x)dα′ ≈ -2π-- δΦ(xij), 0 2π 0 N1N2 i=1j=1$ (273)

where N₁ = 4,N₂ = 4 for the case shown in Fig. 2.

The above 3 steps specify how to approximate the double gyro-angle integration as the sum of values of δΦ at discrete spatial points x_ij. The spatial points x_ij are not necessarily grid points. Interpolations are used to further express δΦ(x_ij) as weighted combination of δΦ values at gridpoints.

6.4 Toroidal DFT of polarization density in ﬁeld-aligned coordinates

In the ﬁeld-aligned coordinates (x,y,z), Fourier expand δΦ along y:

Ny∕2 ( ) δΦ (x) ≈ ∑ exp ιn 2πy δΦn(x,z), n=−Ny∕2 Ly

(274)

where ι = √ ---
− 1 , N_y is an even integer specifying how many harmonics are included, δΦ_n is the expansion coeﬃcient. Using expression (274) in Eq. (272), we obtain

q ∫ ∞ ∫ ∞ ∂F 2π ∑N1 1 ∑N2 A (x) = −-- dv∥ dv⊥---0v⊥--- --- m −∞ 0 ∂𝜀 N1 i=1 N2 j=1 N∑y∕2 ( ) exp ιn 2π(y+ Δ ρijy) δΦn (x+ Δ ρijx,z + Δρijz) n=−Ny∕2 Ly Ny∕2 ( ) ( )∫ ∞ ∫ ∞ N1 N2 = ∑ exp ιn 2πy − q- dv dv ∂F0v 2π-∑ -1-∑ n=−N ∕2 Ly m −∞ ∥ 0 ⊥ ∂𝜀 ⊥ N1 i=1 N2 j=1 (y ) exp ιn 2πΔ ρijy δΦn (x+ Δ ρijx,z + Δρijz). (275) Ly

From Eq. (275), the Fourier expansion coeﬃcient of A(x) can be identiﬁed:

( q-)∫ ∞ ∫ ∞ ∂F0- 2π∑N1 1--N∑2 An (x,z) = − m −∞ dv∥ 0 dv⊥ ∂𝜀 v⊥ N1 N2 ( ) i=1 j=1 exp ιn 2πΔ ρ δΦ (x + Δρ ,z + Δρ ). (276) Ly ijy n ijx ijz

Assume each toroidal harmonic amplitude δΦ_n(x,z) is weakly varying along ﬁeld line, then we can make the following approximation:

δΦn(x+ Δ ρijx,z + Δ ρijz) ≈ δΦn(x + Δρijx,z),

(277)

because the variation along a ﬁeld line over a distance of a Larmor radius is small. Then expression (276) is written as

( )∫ ∞ ∫ ∞ ∑N1 N∑2 An (x,z) ≈ − q- dv∥ dv⊥ ∂F0v⊥ 2π- 1-- m −∞ 0 ∂𝜀 N1 i=1 N2 j=1 ( 2π ) exp ιn L-Δ ρijy δΦn (x + Δρijx,z). (278) y

The approximation (277) makes the dependence of A_n(x,z) on δΦ_n(x,z) become local in the z direction, i.e., A_n(x,z) at z = z₀ only involves values of δΦ_n(x,z) at z = z₀. Meanwhile A_n(x,z) at x = x₀ involves values of δΦ_n(x,z) at x≠x₀ besides x = x₀.

For Maxwelliand distribution

Assume that F₀ is Maxwellian, then

( ) ∂F0- m- m-----n0---- − mv2∥ ( − mv2⊥) ∂ 𝜀 = − T fM = −T (2πT ∕m)3∕2 exp 2T exp 2T .

(279)

Note that Δρ_ij is independent of v_∥. Then the integration over v_∥ in Eq. (278) can be analytically performed:

( )∫ ∞ (− mv2 ) ∫ ∞ ( ( 2) ) An (x,z) = q- exp ----∥- dv∥ dv⊥ ----n0-3∕2 exp − mv-⊥ v⊥ T − ∞ 2T 0 (2πT ∕m) 2T 2π N∑1 1 ∑N2 ( 2π ) × --- --- exp ιn--Δ ρijy δΦn(x+ Δ ρijx,z) N1 i=1N2 j=1 Ly q ∫ ∞ ( − v2 ) 1 ∑N1 1 N∑2 = --n0 dv⊥v⊥ exp --⊥- --- --- T 0 2 N1 i=1 N2 j=1 ( 2π ) exp ιn LyΔ ρijy δΦn (x + Δρijx,z).

where v_⊥ = v_⊥∕ ∘ ----
T∕m , and use has been made of

∫ ( ) ∞ exp − x2 dx = √2-π. −∞ 2

(280)

Using sine expansion for the radial dependence of δΦ

Expand δΦ_n(x,z) in terms of sine basis functions sin(mπx∕L_x):

N∑x−1 (m π ) δΦn(x,z) = δΦmn (z)sin -Lx x , m=1

(281)