Gyrokinetic equation in forms amenable to numerical computation

4 Gyrokinetic equation in forms amenable to numerical computation

4.1 Eliminate ∂⟨δϕ⟩_α∕∂t term on the right-hand side of Eq. (136)

Note that the coeﬃcient before ∂F₀∕∂𝜀 in Eq. (136) involves the time derivative of ⟨δϕ⟩_α, which is problematic if treated by using explicit ﬁnite diﬀerence in particle simulations (I test the algorithm that treats this term by implicit scheme, the result roughly agrees with the standard method discussed in Sec. 6). It turns out that ∂⟨δΦ⟩_α∕∂t can be eliminated by deﬁning another gyro-phase independent function δf by

q- ∂F0- δf = m ⟨δΦ ⟩α ∂𝜀 +δG0.

(139)

Substituting this into Eq. (136), we obtain the equation for δf:

[ ] -∂ + (ve + V + δV )⋅∇ δf ∂t ∥ ∥ D D X q-∂F0-[ ∂ ] − m ∂𝜀 ∂t + (v∥e∥ + VD + δVD ) ⋅∇X ⟨δΦ⟩α q [ ∂ ] ∂F − --⟨δΦ ⟩α -- + (v∥e∥ + VD + δVD )⋅∇X ---0 m ∂t ∂𝜀 = − δVD ⋅∇XF0 − q-∂⟨δL⟩α∂F0-. (140) m ∂t ∂𝜀

Noting that ∂F₀∕∂t = 0, e_∥⋅∇F₀ = 0, ∂F₀∕∂𝜀 ≈−F₀m∕T, we ﬁnd that the third line of the above equation is of order O(λ³) (after both sides being divided by F₀Ω). Therefore the third line can be dropped. Moving the second line to the right-hand side and noting that ⟨δL⟩_α = ⟨δϕ − v ⋅ δA⟩_α, the above equation is written as

[ ∂ ] --+ (v∥e∥ + VD + δVD )⋅∇X δf ∂t = − δ[VD ⋅∇XF0 ] −-q − ∂⟨v-⋅δA⟩α − (v e + VD + δVD )⋅∇X ⟨δΦ ⟩α ∂F0-, (141) m ∂t ∥ ∥ ∂𝜀

where two ∂⟨ϕ⟩_α∕∂t terms cancel each other. Note that the right-hand side of Eq. (141) contains a nonlinear term δV_D ⋅∇_X⟨δΦ⟩_α. This is diﬀerent from the original Frieman-Chen equation, where all nonlinear terms appear on the left-hand side. For the electrostatic limit, this term disappears because δV_D is perpendicular to ∇_X⟨δΦ⟩_α.

[Equation (141) corresponds to Eq. (A8) in Yang Chen’s paper[2] (where the ﬁrst minus on the right-hand side is wrong and should be replaced with q∕m; one q is missing before ∂(v ⋅ δA)∕∂t in A9).]

4.2 Eliminate ∂⟨δv ⋅ δA⟩_α∕∂t term on the right-hand side of GK equation

Similar to the method of eliminating ∂⟨δϕ⟩_α∕∂t, we deﬁne another gyro-phase independent function by

q ∂F δh = δf −--⟨v ⋅δA ⟩α--0-. m ∂𝜀

(142)

Most gyrokinetic simulations approximate the vector potential as δA ≈ δA_∥e_∥. Let us simplify Eq. (147) for this case. Then ⟨v ⋅ δA⟩_α is written as

⟨v ⋅δA⟩α ≈ ⟨v∥δA ∥⟩α.

(143)

Then expression (142) is written as

-q ∂F0- δh = δf − m ⟨v∥δA∥⟩α ∂𝜀 .

(144)

Then Eq. (141) is written in terms of δh as

[ ∂ ] ∂t + (v∥e∥ + VD +δVD ) ⋅∇X δh [( ) ] +-q ∂F0- ∂-+ v∥e∥ + VD + δVD ⋅∇X ⟨v∥δA ∥⟩α m ∂𝜀 ∂t ( ) -q ∂F0- +m ⟨v∥δA∥⟩α[(VD + δVD )⋅∇X ] ∂𝜀 = − δVD ⋅∇XF0 q [ ∂ ⟨v∥δA∥⟩α ] ∂F −-- −--------- − (v∥e∥ +VD + δVD )⋅∇X ⟨δΦ⟩α ---0, (145) m ∂t ∂𝜀

where use has been made of that ∂F₀∕∂t = 0 and e_∥⋅∇F₀ = 0. Further noting that v_∥δA_∥∼ δΦ, qδΦ∕T ∼ O(λ¹), ∂F₀∕∂𝜀 ∼−F₀m∕T, V_D∕v_t ∼ O(λ¹), δV_D∕v_t ∼ O(λ¹), and ρ∇F₀ ∼ O(λ¹)F₀, we ﬁnd that the red term of the above equation (after divided by ΩF₀) is of O(λ³), hence can be dropped. Move the second line to the right-hand side, giving

[ ∂ ] ∂t +(v∥e∥ + VD + δVD )⋅∇X δh = − δV ⋅∇ F D X 0 + q-[(v∥e∥ + VD + δVD )⋅∇X (⟨δΦ − v∥δA ∥⟩α)]∂F0, (146) m ∂ 𝜀

where ∂⟨v_∥δA_∥⟩_α∕∂t terms cancel each other. Further note that δV_D (given by Eq. (138)) is perpendicular to ∇_X⟨δΦ − v_∥δA_∥⟩_α. Therefore the blue term in Eq. (146) is zero, then Eq. (146) simpliﬁes to

[ ] ∂-+ (v e + V + δV )⋅∇ δh ∂t ∥ ∥ D D X = − δVD ⋅∇XF0 -q ∂F0- +m [(v∥e∥ +VD ) ⋅∇X ⟨δΦ − v∥δA∥⟩α]∂𝜀 . (147)

Note that in terms of (X,𝜀,μ,α,σ) coordinates, v_∥ is written as

∘ --------- v∥ = σ 2𝜀− 2μB0,

(148)

where B₀(x) = B₀(X + ρ) with ρ = ρ(X,𝜀,μ,α). Since the scale length of B₀ is much larger than the thermal Larmor radius, B₀(x) ≈ B₀(X) and hence v_∥ can be approximated as a constant when gyro-angle α changes. Then v_∥ can be taken out of the gyro-averaging in expression (143), yielding

⟨v∥δA∥⟩α ≈ v∥⟨δA ∥⟩α.

(149)

Using this, the term related to δA_∥ in (147) can be further written as

(v∥e∥ + VD )⋅∇X ⟨v∥δA ∥⟩α = (v∥e∥ + VD )⋅∇X (v∥⟨δA∥⟩α) = ⟨δA ∥⟩α(v∥e∥ + VD ) ⋅∇X (v∥)+ v∥(v∥e∥ + VD )⋅∇X ⟨δA∥(⟩1α5.0)

Using expression (148), (v_∥e_∥ + V_D) ⋅∇_X(v_∥) is written as

⟨δA ⟩ (ve + V )⋅∇ (v ) ≈ ⟨δA ⟩ (v e )⋅∇ (v ) ∥ α ∥ ∥ D X ∥ ∥ α ∥ ∥ X (∥∘ ---------) = ⟨δA∥⟩α(v∥e∥)⋅∇X σ 2𝜀− 2μB0 (∘ ---------) = ⟨δA∥⟩ασ(v∥e∥)⋅∇X 2𝜀− 2μB0 − 2μe∥ ⋅∇XB0 = ⟨δA∥⟩ασv∥-2√2𝜀-−-2μB--- 0 = ⟨δA∥⟩αv∥−-2μe∥ ⋅∇XB0- 2v∥ = − ⟨δA ∥⟩α μe∥ ⋅∇XB0. (151)

(We can also obtain ∇_X(v_∥) = −μ(∇B₀)∕v_∥ by using Eq. (370)). Using the above results, equation (147) is written as

[ ] ∂- + (v∥e∥ + VD + δVD )⋅∇X δh ∂t -q ∂F0- = − δVD ⋅∇XF0 + m [(v∥e∥ + VD )⋅∇X ⟨δΦ ⟩α]∂𝜀 -q ∂F0- − m [v∥(v∥e∥ + VD )⋅∇X ⟨δA∥⟩α− ⟨δA ∥⟩α(μe∥ ⋅∇XB0 )]∂𝜀 , (152)

which agrees with the so-called p_∥ formulation given in the GEM manual (the ﬁrst line of Eq. 28). (It is said that p_∥ formulation uses p_∥ = v_∥ + q⟨δA_∥⟩_α∕m as an independent variable, but I can not relate this explanation to my derivation given above.)

4.3 Mixed-variable pullback method[4]

Deﬁne δf^(h) by

(h) q (h) ∂F0 δf = δf − m-⟨v∥δA ∥ ⟩α∂-𝜀 ,

(153)

where δA_∥^(h) is a part of δA_∥:

(s) (h) δA∥ = δA∥ + δA ∥ ,

(154)

with δA_∥^(s) determined by an evolution equation:

∂δA (s∥) ------ = − b ⋅∇δΦ. ∂t

(155)

Following the same procedure in Sec. 4.2, we obtain the following evoluation equation for δf^(h):

[ ∂ ] (h) ∂t + (v∥e∥ + VD +δVD ) ⋅∇X δf = − δVD ⋅∇XF0 + -q[(VD + δVD )⋅∇X ⟨δΦ⟩α]∂F0- m ∂𝜀 −-q[v∥(v∥e∥ + VD + δVD )⋅∇X ⟨δA(∥h)⟩α − ⟨δA(∥h)⟩α(μe∥ ⋅∇XB0 )]∂F0-. (156) m ∂𝜀

Next, consider the special case that F₀ is a local Maxwellian given by

( )3∕2 [ ] F0(X,𝜀) = n0(X ) --m----- exp − -m-𝜀-- , 2πT0(X ) T0(X)

(157)

Then ∂F₀∕∂𝜀 and ∇_XF₀ are written as

∂F ( m ) ---0= − -- F0, ∂𝜀 T

(158)

and

∂F0 ∂F0 ∇F0 = -∂n0∇n0 + ∂T-∇T0 ∇n ( mv2 3) ∇T = F0 --0-+ F0 ----− - ---0 n0[ ( 2T0 2) T]0 = − F κ + mv2- − 3 κ , (159) 0 n 2T0 2 T

where κ_n = − ∇nn00

, and κ_T = − ∇TT00

. Consider the case that F₀(X,𝜀) is independent of 𝜃 and α, then κ_n and κ_Tare written as κ_n = κ_n∇r and κ_T = κ_T∇r, where κ_n ≡− -1
n0

, κ_T ≡−

, and r is chosen to be the minor radius at the low ﬁeld side.

Using the above results, Eq. (156) is written as

[ ( ) ] dδf(h) -dr mv2- 3 q- dt = δVD ⋅∇ ψdψ κn + 2T0 − 2 κT F0 − T[(VD + δVD )⋅∇ ⟨δΦ ⟩]F0 q- (h) (h) + T [v∥(v∥e∥ + VD + δVD )⋅∇ ⟨δA∥ ⟩− ⟨δA∥ ⟩(μe∥ ⋅∇B0 )]F0. (160)

4.3.1 Weight evolution equation

The weight of the j th marker is deﬁned by

w = δf(h)(X ,v ,t)V , j j j pj0

(161)

where V _pj0 is the phase space volume occupied by the jth marker. Multiplying both sides of Eq. (160) by V _pj0, and noting that d(V _pj0)∕dt = 0, we obtain

[ ( 2 ) ] dwj- = δVD ⋅∇ψ dr- κn + mv--− 3 κT Vpj0F0 − q[(VD + δVD )⋅⟨∇δΦ ⟩]Vpj0F0. dt dψ 2T0 2 T + -q[v∥(v∥e∥ + VD + δVD )⋅⟨∇ δA(h)⟩ − ⟨δA (h)⟩(μe∥ ⋅∇B0 )]Vpj0F0, (162) T ∥ ∥

where use has been made of ∇⟨δΦ⟩ = ⟨∇δΦ⟩ and ∇⟨δA_∥^(h)⟩ = ⟨∇δA_∥^(h)⟩.

Note that ∇≡ ∂∕∂X. Using x = X + ρ, we obtain

∂δΦg(X,α,v⊥-) = ∂δΦp(X-+-ρ(X,α,v⊥-))- ∂X ∂X = ∂δΦp(x) = ∂δΦp(x) ⋅ ∂x- ∂X ( ∂x ) ∂X ∂δΦp(x) -∂ρ = ∂x ⋅ 1+ ∂X ∂δΦp(x) ≈ --∂x---, (163)

which is accurate up to O(λ) in gyrokinetic ordering. Note that ﬁeld values on a regular X grid are not available in our simulation because X is sampled by random markers rather than by a regular grid. This fact makes it diﬀcult to numerical compute ∂∕∂X. Then the relation ∂δΦ_g∕∂X ≈ ∂δΦ_p∕∂x, becomes useful because we can use δΦ deﬁned on x grid (which is available in simulations) to compute ∂δΦ_p∕∂x at x grid-points (using ﬁnite-diﬀerence), and then interpolate it to particle position. Similarly, we have ∂δA_∥g∕∂X ≈ ∂δA_∥p∕∂x.

Also note that ⟨δΦ_g⟩ is not known on X grid (it is known on random markers). So, to calculate ∂⟨δΦ_g⟩∕∂X, we need to exchange the operation order:

δ⟨δΦg-⟩ ∂δΦg- ∂δΦp- ∂X = ⟨ ∂X ⟩ ≈ ⟨ ∂x ⟩.

(164)

For notation ease, the subscripts g and p will be dropped in the following.

4.3.2 In ﬁeld-aligned coordinates

In ﬁeld-aligned coordinates (x,y,z), ∂δΦ
∂x can be written as

∂δΦ- = ∂δΦ∇x + ∂δΦ∇y + ∂-δΦ-∇z. ∂x ∂x ∂y ∂z

(165)

Note that the direction of ∂r∕∂x at 𝜃 = −π is diﬀerent from that of ∂r∕∂x at 𝜃 = +π. This implies that ∂δΦ∕∂x is a non-periodic function of 𝜃. Similarly, ∇y is also a non-perioidic functon of 𝜃. Then both ∂δΦ∕∂x and ∇y are discontinuous across the 𝜃 cut (𝜃 = ±π).

The gyro-average of expression (165) is written as

⟨∂δΦ-⟩ ≈ ⟨∂δΦ-⟩∇x + ⟨∂δΦ-⟩∇y + ⟨∂δΦ⟩∇z, ∂x ∂x ∂y ∂z

(166)

where we approximate ∇x, ∇y, and ∇z as constants when performing the gyro-averge since they are determined by the equilibrim magnetic ﬁeld, which is nearly constant on the Larmor radius scale. As is mentioned above, ∇y is not continuous at the 𝜃 cut. Do we need to worry about this? We do not. This is because we must stick to the same branch when we perform gryo-average on it, and hence ∇y is always continous. Then, do we need to worry about the disconunity of ∂δΦ∕∂x across the 𝜃 cut when performing the gyro-averge on it? We do not either. The reason is the same: we must stick to a single branch. The disconunity is just irrelevant here. The discontinuity only manifest itself when we need to infer value on 𝜃 = −π from that on 𝜃 = +π (vice versa), i.e., when across branch communication is explicitly needed. In TEK, the ﬁeld equations are not solved at 𝜃 = +π and hence the ﬁeld values are not directly obtained. Instead, the ﬁeld values at 𝜃 = +π are infered from the ﬁeld values at 𝜃 = −π using the coninuity of ∇δΦ. At the 𝜃 cut, the continuity of ∇δΦ requires

⟨∂δΦ-⟩+∇x + ⟨∂δΦ⟩∇y+ = ⟨∂δΦ⟩− ∇x + ⟨∂δΦ-⟩∇y− , ∂x ∂y ∂x ∂y

(167)

where the superscript “+” and “−” refer to the location 𝜃 = +π and 𝜃 = −π, respectively. Dotting the above by ∇x, we obtain

∂δΦ- + ∂-δΦ- − ∂-δΦ- ∇y-− ⋅∇x-−-∇y+-⋅∇x-- ⟨ ∂x ⟩ = ⟨ ∂x ⟩ +⟨ ∂y ⟩ |∇x |2 ,

(168)

which is used in TEK to infer the value of ⟨ ∂δΦ
∂x ⟩⁺ from that of ⟨⟩⁻.

Similarly, for ∂δA_∥^(h)∕∂x, we obtain:

∂δA(h) ∂δA(h) ∂δA(h) ∂δA(h) ---∥--= ---∥-∇x + ---∥-∇y + ---∥-∇z. ∂x ∂x ∂y ∂z

(169)

(h) (h) (h) (h) ⟨∂δA∥--⟩ ≈ ⟨∂δA∥-⟩∇x + ⟨∂δA∥-⟩∇y + ⟨∂δA-∥-⟩∇z ∂x ∂x ∂y ∂z

(170)

Using these, Eq. (162) is written as

[ ( 2 ) ] dwj- = δVD ⋅∇ ψ dr-κn + mv--− 3 κT Vpj0F0 dt [ dψ 2T0 2 ] -q ∂δΦ- dX-′ ∂δΦ-dX-′ ∂δΦ- dX′- − T ⟨ ∂x ⟩ dt ⋅∇x + ⟨∂y ⟩ dt ⋅∇y + ⟨ ∂z ⟩dt ⋅∇z Vpj0F0 ( (h) (h) (h) ) -q ( ∂δA∥-- dX- ∂δA∥-- dX- ∂δA∥-- dX- ) + T v∥ ⟨ ∂x ⟩dt ⋅∇x + ⟨ ∂y ⟩dt ⋅∇y + ⟨ ∂z ⟩dt ⋅∇z Vpj0F0 q − --⟨δA(∥h)⟩(μb ⋅∇B0 )Vpj0F0, (171) T

where

≡ v_∥b + V_D + δV_D and dX′
dt

≡ V_D + δV_D. Deﬁne

q v δA(h) δΦ-= quδΦ-,δA(∥h) = -u-u--∥-, Tu Tu

(172)

Then, a typical term of the above is written as

(h) -(h) -- qv∥⟨∂δA-∥-⟩dX-⋅∇xtn = q--Tu-vnv∥⟨∂δA∥--⟩dX-⋅∇x T ∂x dt T quvu ∂x dt q∕q v ∂δA-(h) dX- -- = ---u--nv∥⟨----∥-⟩---⋅∇x, T ∕Tuvu ∂x dt

where

is the normlizing factor needed when coding. The mirror force term:

q (h) q Tu Bn μn -(h) -- --- T-⟨δA ∥ ⟩(μb ⋅∇B0 )tn = T-quvu-Ln--⟨δA∥ ⟩(μb ⋅∇B0 )tn (173) 2 -- --- = q∕qu-1--Bnvn⟨δA(∥h)⟩(μb ⋅∇B0 )tn T∕Tu vuLnBn = q∕qu-vn⟨δA(h)⟩(μb ⋅∇B- ) T∕Tu vu ∥ 0

where

is the normlizing factor needed when coding.

The pullback:

( ) -q⟨v∥δA(h)⟩∂F0-Vp = q-⟨v∥δA (h)⟩ − m- F0Vp m ∥ ∂𝜀 m ∥ T ≈ −-qv∥⟨δA(∥h)⟩F0Vp T = −-q∕quv∥⟨δA-(h∥)⟩F0Vp T ∕Tuvu

The δE × B₀ drift term is written as

∂δΦ ( ) − ⟨-∂x ⟩×2-B0 = − -12 ⟨∂δΦ-⟩∇x + ⟨∂δΦ⟩∇y + ⟨∂δΦ-⟩∇z × B0 B0 B0 ∂x ∂y ∂z = − ⟨∂δΦ⟩-1-∇x × B − ⟨∂δΦ⟩-1-∇y × B − ⟨∂δΦ⟩-1-∇z × B . ∂x B20 0 ∂y B20 0 ∂z B20 0

The ∇x, ∇y, and ∇z components of the above expression are written as

− ⟨∂δ∂Φx ⟩×-B0-⋅∇x = − ⟨∂δΦ⟩-1-∇y × B ⋅∇x − ⟨∂δΦ⟩-1-∇z × B ⋅∇x B20 ∂y B20 0 ∂z B20 0 ∂δΦ--1- ∂δΦ--1- = − ⟨∂y ⟩B2 ∇x × ∇y ⋅B0 − ⟨∂z ⟩B2 ∇x × ∇z ⋅B0. (174) 0 0

⟨∂δΦ-⟩× B ∂δΦ 1 ∂δΦ 1 − -∂x--2--0-⋅∇y = − ⟨----⟩-2∇x × B0 ⋅∇y − ⟨---⟩--2∇z × B0 ⋅∇y B 0 ∂x B0 ∂z B 0 = ⟨∂-δΦ-⟩-1-∇x × ∇y ⋅B0 − ⟨∂δΦ-⟩-1-∇y × ∇z ⋅B0 (175) ∂x B20 ∂z B20

∂δΦ- − ⟨∂x-⟩×-B0-⋅∇z = − ⟨∂δΦ-⟩ 1-∇x × B0 ⋅∇z − ⟨∂δΦ⟩-1-∇y × B0 ⋅∇z B20 ∂x B20 ∂y B20 ∂-δΦ- -1- ∂δΦ- -1- = ⟨ ∂x ⟩B20∇x × ∇z ⋅B0 + ⟨ ∂y ⟩B20 ∇y ×∇z ⋅B0 (176)

Using

dx= V ⋅∇x, dt G

(177)

i.e.,

dx-= t V ⋅∇x dt n G

(178)

To get the normalizing factor, consider only the contribution from E ×B drift. Then a typical term is written as

dx- = − tn⟨∂δΦ-⟩ 1-∇x × ∇y ⋅B0 dt ∂y B20 -- Tu--tn-- ∂δΦ- 1--- -- -- = − quBnL2n⟨ ∂y ⟩B2∇x × ∇y ⋅B0 --0 = − Tu---1---⟨∂δΦ-⟩ 1-∇x × ∇y ⋅B0, (179) quBnvnLn ∂y B20

where T_u∕(q_uv_nB_nL_n) is the normlizing factor we need when coding.

The magnetic ﬂuttering term is written as

⟨ ( ) ⟩ −-v∥⟨b× ∇ (δA )⟩ = −-v∥- b× ∂δA∥∇x + ∂δA∥∇y + ∂δA-∥∇z B0 x ∥ B0 ∂x ∂y ∂z -v∥- ( ∂δA∥- ∂δA∥- ∂δA-∥ ) ≈ −B0 b× ⟨ ∂x ⟩∇x + ⟨ ∂y ⟩∇y + ⟨ ∂z ⟩∇z

The ∇x, ∇y, and ∇z components of the above expression are written as

( ) -v∥ -v∥ ∂δA∥- ∂δA∥- −B0 ⟨b× ∇x (δA∥)⟩⋅∇x = −B0 b × ⟨ ∂y ⟩∇y + ⟨ ∂z ⟩∇z ⋅∇x v∥ ( ∂δA ∥ ∂ δA ∥ ) = B2- ⟨-∂y--⟩∇x × ∇y +⟨--∂z-⟩∇x × ∇z ⋅B0 (180) 0

v∥- -v∥ ( ∂δA∥- ∂δA-∥ ) − B0⟨b× ∇x (δA∥)⟩⋅∇y = −B0 b × ⟨ ∂x ⟩∇x + ⟨ ∂z ⟩∇z ⋅∇y v ( ∂δA ∂ δA ) = -∥2- ⟨----∥⟩∇y ×∇x +⟨----∥⟩∇y × ∇z ⋅B0 (181) B0 ∂x ∂z

( ) − v∥⟨b× ∇ (δA )⟩⋅∇z = −-v∥b × ⟨∂δA∥⟩∇x + ⟨∂δA-∥⟩∇y ⋅∇z B0 x ∥ B0 ∂x ∂y v∥-( ∂δA-∥ ∂δA-∥ ) = B2 ⟨ ∂x ⟩∇z × ∇x + ⟨ ∂y ⟩∇z × ∇y ⋅B0 (182) 0

Also to get the normalizing factor needed in the code, consider a typical term:

dx- v∥-∂δA-∥ dt = tn B20⟨ ∂y ⟩∇x × ∇y ⋅B0 T v v ∂δA- -- -- -- = --u-tn-2n---∥2⟨---∥⟩∇x × ∇y ⋅B0 quvu LnBn B 0 -∂y Tu 1 v∥ ∂δA∥ -- -- -- = quvuLnBn---2⟨-∂y--⟩∇x × ∇y ⋅B0 (183) B 0

where T_u∕(q_uv_uB_nL_n) is the normlizing factor I need when coding.

The basic terms involved in the above expressions can be written as

B2 (∇x × ∇y )⋅B0 = -0′, Ψ

(184)

( ∂x ∂z ∂x ∂z ) (∇x ×∇z )⋅B0 = ------− ------ Bϕ ≡ W 14 ∂Z ∂R ∂R ∂Z

(185)

-- ( 1 ∂α ∂𝜃) -- ( ∂α ∂𝜃 ∂α ∂𝜃 ) -- ( 1 ∂α ∂𝜃) (∇y × ∇z)⋅B0 = BR R-∂ϕ-∂Z- + B ϕ ∂Z-∂R-− ∂R-∂Z- + BZ − R-∂ϕ-∂R- ≡ W 15

(186)

In terms of δA_∥^(h), the Ampere’s law is written as

(∑ ω2 ) ∑ -p2s− ∇2⊥ δA(∥h) = μ0 δJ(||hs)+ ∇2⊥ δA(∥s), s c s

(187)

where δJ_||s^(h) is the parallel current carried by the distrbution funciton δf_s^(h), where the subscript s is species index (do not confuse it with the superscript s of δA_∥^(s)). Here δA_∥^(s) has been moved to the right-hand side since its value is already known when we solve the Ampere’s law for δA_∥^(h).

Compared with the p_∥ formulation, the improvement of the above procedure is that a part of δA_∥ is solved from an evolution equation (which is chosen by some physical consideration) and the remainder is solved from the Ampere’s law. We hope that the former is the dominant part of δA_∥ so that the remainder solved from the Ampere’s law is small and hence the cancellation problem (of the Ampere’s law) become less signiﬁcant. While a wise choice of the evolution equation for δA_∥^(s) may help ensure that δA_∥^(s) remain dominant over the entire simulation duration, we have another simple way to make δA_∥^(s) remain dominant: collect both δA_∥^(s) and δA_∥^(h) into δA_∥^(s) at the end of each time step:

(s) (s) (h) δA ∥new = δA∥old + δA∥old.

(188)

Then, to make δA_∥ untouched (so that electromagnetic ﬁeld remain unchanged), we set δA_∥^(h) to zero:

δA (h) = 0. ∥new

(189)

Here “old” and “new” refers to before and after the re-spliting, respectively. Note that this re-spliting is made at the end of each time step and does not correspond to any time evolution. This re-splitting keeps the value of δA_∥ untouched and hence does not inﬂuence the electromagnetic ﬁeld. Meanwhile, the deﬁnition of Eq. (153) indicates that, for a given δf, the re-splitting will make the value of δf^(h) increase by Δ = q-
m ⟨v_∥δA_∥old^(h)⟩_α ∂F0
∂𝜀 since δA^(h) is reduced by δA_old^(h). In other words, to keep δf unchanged, which is what we need to ensure, we must add Δ to δf^(h):

δf(h) = δf(h) + q-⟨v δA (h)⟩ ∂F0-. new old m ∥ ∥old α ∂𝜀

(190)

After the operations in (188) (189) and (190), both the electromagnetic ﬁeld and the distribution function δf remain unchanged. Therefore we do not change any physical process of the system. Hence the above operations can only have eﬀects on numerical round-oﬀ errors. As mentioned above, we hope the numerical eﬀects will reduce the cancellation errors.

Deﬁne

-- quδΦ -(h) quvuδA(∥h) -(s) quvuδA (s∥) -(h) δJ(∥ss) δΦ = -T--,δA∥ = ---T-----,δA∥ = ---T----,δJ||s = n-q-v--, u u u u u ns

(191)

where q_u,n_u,T_u,v_u are chosen units, v_ns is a velocity unit for a species, whose value is diﬀerent among diﬀerent species [v_ns = L_n∕(2π∕(B_n|q_s|∕m_s))]. In TEK, v_u is chosen as v_ns of the ﬁrst ion species.

Then the Ampere’s law (187) is written as

( 2 ) 2 2 ∑ ωps− ∇2 δA(h)= -ω-pu--vu ∑ vnsδJ(h)+ ∇2δA-(s), s c2 ⊥ ∥ Tu∕mu c2 s vu ||s ⊥ ∥

(192)

where ω_pu² = n_uq_u²∕(m_u𝜀₀), m_u is a mass unit, and use has been made of μ₀ = 1∕(c²𝜀₀).

δf pullback in Eq. (190) is now written as (for the case of F₀ being Maxwellian):

q ( m) δf(nhe)w = δf(ohld)+ --⟨v∥δA (h∥o)ld⟩α − -- F0 (h) qm (h) T = fold − T-⟨v∥δA∥old⟩αF0 (h) q - vn -(h) = fold − T-v∥vu⟨δA∥old⟩αF0 (193)

In terms of units used in TEK, Eq. (155) is written as

= −

b ⋅∇δΦ,

where t = t∕t_u, t_u = L_n∕v_u. Then the above equation is written as

--(s) ∂δA-∥- = − b ⋅∇δΦ, ∂t

(194)

i.e.,

∂δA-(s∥) -- ∂δΦ- ------ = − (b⋅∇z )-∂z-. ∂t

(195)

4.4 Discretizing Laplacian operator

The Laplacian operator is approximated as

2 2 2 ∇2 δA∥ ≈ ∂-δA∥|∇x|2 + ∂-δA-∥|∇y |2 + ∂-δA∥2∇x ⋅∇y. ⊥ ∂2x ∂2y ∂x∂y

(196)

To approximate δA_∥ of zero boundary condition A_∥(x = 0) = A_∥(x = L_x) = 0, the sine expansion can be adopted. In the y direction, full Fourier expansion is needed. Then, at each value of z, A_∥(x,y,z) is approximated by the following two-dimensional expansion:

N∑y∕2 ( 2π )Nx∑−1 ( m π ) δA ∥(x,y,z) ≈ exp in--y Am,n(z)sin ---(x− x0) , n=− Ny∕2 Ly m=1 Lx

(197)

Using this expression, Eq. (196) is written as

( ) 2 Ny∑∕2 -2π N∑x−1 − ∇ ⊥δA∥ = exp inLy y δAm,n(z) n{=[−Ny∕2 m=1 ] } ( m-π)2 2 ( 2π-)2 2 (m-π ) 2π-( mπ-) ( mπ- ) × Lx |∇x | + nLy |∇y| sin Lx (x − x0) − inLy Lx 2∇x ⋅∇y cos Lx (x − x0)

At x = x_j, with x_j − x₀ = jΔ and L_x∕Δ = N_x, the above expression is written as

2 Ny∑∕2 ( 2π ) N∑x−1 − ∇ ⊥δA∥ = exp inL--y δAm,n(z) n{=[−Ny∕2 y m=1 ] } ( m π)2 ( 2π )2 (jm π) 2π ( mπ ) ( jmπ ) × L-- |∇x |2 + nL-- |∇y|2 sin -N-- − inL-- L-- 2∇x ⋅∇y cos N--(198,) x y x y x x

where ∇x and ∇y are evaluated at x_j. Plugging the DST coeﬃcients

N − 1 ( ) δA (z) =-2- ∑x A (x ′,z)sin j′m-π- , m,n Nx j′=1 n j Nx

(199)

into expression (198) gives

Ny∑∕2 ( 2π ) N∑x−1 2 Nx∑−1 (j′mπ ) − ∇2⊥δA∥ = exp in---y --- An(xj′,z)sin ----- n=−Ny∕2 Ly m=1 Nx j′=1 Nx {[( )2 ( )2 ] ( ) ( ) ( )} × m-π |∇x |2 + n2π- |∇y|2 sin jm-π − in2π- mπ- 2∇x ⋅∇y cos jmπ- Lx Ly Nx Ly Lx Nx

For a single toroidal harmonic, the above expression reduces to

N∑x−1 Nx∑−1 ( ′ ) -2- An(xj′,z)sin j-mπ- m=1 Nx j′=1 Nx {[( )2 ( )2 ] ( ) ( ) ( )} × m-π |∇x |2 + n2π- |∇y|2 sin jm-π − in2π- mπ- 2∇x ⋅∇y cos jmπ(200) Lx Ly Nx Ly Lx Nx

Deﬁne

( ) Nx∑− 1-2- j′mπ- Mjj′,n ≡ Nx sin Nx {m=[1( )2 ( )2 ] ( ) ( ) ( ) } × mπ- |∇x|2 + n2π- |∇y |2 sin jm-π − in2π- m-π 2∇x ⋅∇y cos jm-π Lx Ly Nx Ly Lx Nx

then expression (200) is written as

Nx∑ −1 Mjj′,nAn (xj′,z). j′=1

(201)

4.5 Summary of distribution function split

In simulations, the seemingly trivial thing on how to split the distribution function is often considered to be a big deal. Separating the perturbation from the total distribution function gets the famous name “δf particle method”, in contrast to the conventional particle method which is now called full-f particle method.

In the above, the perturbed part of the distribution function, δF, is split at least three times in order to (1) simplify the gyrokinetic equation by splitting out the adiabatic response; (2) eliminate the time derivatives, ∂δϕ∕∂t, and (3) ∂δA∕∂t, on the right-hand. To avoid confusion, I summarize the split of the distribution function here. The total distribution function F is split as

F = F0 +δF,

(202)

where F₀ is the equilibrium distribution function and δF is the perturbed part of the total distribution function. δF is further split as

q- ∂F0- q- ∂F0- δF = δh + m (δΦ− ⟨δΦ⟩α) ∂𝜀 + m ⟨v⋅δA ⟩α∂ 𝜀 ,

(203)

where δh satisﬁes the gyrokinetic equation (147) or (152). In PIC simulations, δh is evolved by using markers and its moment in the phase-space is evaluated via Monte-Carlo integration.

The integration of the red and blue terms can be approximated analytically in terms of δΦ and δA. But this can also be evaluated using markers, i.e., using Monte-Carlo method, to avoid the inaccurate cancellation between the integration of these parts and the integration of δh.

Another issue is how to accurately solve for δΦ from the Poisson equation and δA from the Ampere’s law. To make the solution accurate, one need to have a large term that is explictly dependent on δΦ and δA to appear on the left-hand side of the equations and small terms appear on the right-hand side.

The red term in expression (203) gives “the polarization density” when integrated in the velocity space (discussed in Sec. 6). The reason for the name “polarization” is that (δΦ −⟨δΦ⟩_α) is the diﬀerence between the local value and the averaged value on a gyro-ring, expressing a kind of “separation”.

are also needed to be evaluated using Monte-Carlo markers.

This term is moved to the left-hand side of the Poisson equation and is utilized in solving the Poisson equation. The blue term also has an explicit dependence on δA, which, however, will cause numerical problems in particle simulations if it is moved to the left-hand side of the Ampere equation. This is the “cancellation problem” in gyrokinetic simulations discussed in Sec. 4.5.

The blue and red terms in the above expression explicitly depends on the perturbed ﬁeld. The velocity integrations of these two terms can be performed analytically. However, in some cases, the phase space integration of the blue terms must be evaluated using markers, i.e., using Monte-Carlo method, to avoid the inaccurate cancellation between the integration of these parts and the integration of δh (the latter is computed using Monte-Carlo method). When will the inaccurate cancellation is signiﬁcant depends on the problem being investigated and thus can only be determined by actual numerical experiments. Electromagnetic particle simulation experiments indicate that the parallel current carried by the blue term must be evaluated via Monte-Carlo method, otherwise inaccurate cancellation between this term and δh will give rise to numerical instabilities.