1 Introduction

1.1 Gyrokinetic?

Electromagnetic perturbations of frequency lower than the ion cyclotron frequency are widely believed to be more important than high-frequency ones in transporting plasma in tokamaks. (This assumption can be verified numerically when we are able to do a full simulation including both low-frequency and high-frequency perturbations. This kind of verification is not possible at present due to computation costs.)

If we assume that only low-frequency perturbations are present, the Vlasov equation can be simplified. Specifically, symmetry of the particle distribution function in the phase space can be established if we choose suitable phase-space coordinates and split the distribution function in a proper way. The symmetry is along the so-called gyro-angle α in the guiding-center coordinates (X,μ,𝜀,α). In obtaining the equation for the gyro-angle independent part of the distribution function, we need to average the equation over the gyro-angle α and thus this model is called “gyrokinetic”.

In deriving the gyrokinetic equation, the perturbed electromagnetic field is assumed to be known and of low-frequency. To do a kinetic simulation, we need to solve the field equation to obtain the perturbed electromagnetic field. It is still possible that high frequency modes (e.g., compressional Alfven waves and Ω_H modes) appear in a gyrokinetic simulation. If the amplitude of high frequency modes is large, then the simulation is invalid since the gyrokinetic model is invalid in this case.

Our starting point is the Vlasov equation in terms of particle coordinates (x,v):

(1)

where f = f(t,x,v) is the particle distribution function, x and v are the location and velocity of particles. The distribution function f depends on 6 phase-space variables (x,v), besides the time t.

1.2 Guiding-center coordinates: a simple example

Choose Cartesian coordinates (x,y,z) for the configuration space x. Consider a simple case where the electromagnetic field is a time-independent field given by B = B₀(x,y,z) and E = 0. Let us examine the Vlasov equation in this case and see whether there is any coordinate system that can simplify the Vlasov equation.

Describe the velocity space using a right-handed cylindrical coordinates (v_⊥,α,v_∥), where v_∥ = v ⋅ e_∥, e_∥ = e_z is the unit vector along the magnetic field, α is the azimuthal angle of the perpendicular velocity (v_⊥ = v − v_∥e_z) relative to e_x. Note that these coordinates are defined relative to the local magnetic field, which, in more general cases, may vary in space.

Next, let us express the spatial gradient of f in terms of partial derivatives:

(3)

Using chain rule, expression (2) is written as

(6)

are generally nonzero since the defintion of (v_⊥,α,v_∥) depends on the local magnetic field direction. In our simple case, they are zero since the magnetic field direction is uniform. Similarly, we obtain

(7)

(8)

Then, in (x′,y′,z′,v_⊥,α,v_∥) coordinates, the spatial gradient of f is written as

(9)

Next, consider the gradient in velocity space

(12)

i.e.,

(13)

Define the following coordinates transform (guiding-center transform):

(14)

Next, we express the partial derivatives of f appearing in Eq. (13) in terms of the guiding-center coordinates. Using the chain rule, we obtain

(16)

and

(19)

Using the above results,  equation (13) in the guiding-center coordinates is written as

(20)

Here we find some terms cancel each other, giving

(21)

Equation (21) is the equation for f in the guiding-center coordinates (X,Y,Z,v_⊥′,α′,v_∥′).

Gyro-averaging

Define gyro-phase averaging operator

(22)

Then averaging Eq. (21) over α′, we get

(23)

As is assumed above, Ω is approximately independent of α′ and thus Ω in Eq. (23) can be moved outside of the gyro-angle integration, giving

(24)

Peforming the integration, we get

(25)

Since α′ = 2π and α′ = 0 correspond to the same phase space location, the corresponding values of the distribution function must be equal. Then the above equation reduces to

(26)

which is an equation for the gyro-angle independent part of the distribution function.

Next let us investigate whether the guiding-center transform can be made use of to simplify the kinetic equation for more general cases where we have a (weakly) non-uniform static magnetic field, plus electromagnetic perturbations of low frequency (and of small amplitude and k_∥ρ_i ≪ 1). And we will include the effect that B₀ depends on α′, i.e., grad-B and curvature drift.

2 Transform Vlasov equation from particle coordinates to guiding-center coordinates

Next, we define the guiding-center transform and then transform the Vlasov equation from the particle coordinates (x,v) to the guiding-center coordinates, i.e., express the gradient operators ∂∕∂x and ∂∕∂v in terms of the guiding-center coordinates.

2.1 Guiding-center transformation

In a magnetic field, given a particle location and velocity (x,v), we know how to calculate its guiding-center location X, i.e.,

(27)

where e_∥ = B₀∕B₀, Ω = qB₀∕m, B₀ = B₀(x) is the equilibrium (macroscopic) magnetic field at the particle position. We will consider Eq. (27) as a transform and call it guiding-center transform[1]. Note that the transform (27) involves both position and velocity of particles.

Given (x,v), it is straightforward to obtain X by using Eq. (27). On the other hand, the inverse transform (i.e., given (X,v), to find x) is in principle not easy because Ω and e_∥ depend on x, which usually requires solving a nonlinear equation for its root. Numerically, one can use

(28)

as an iteration scheme to compute x, with the initial guess chosen as x₀ = X. If we stop at the first iteration, then

(29)

In gyrokinetic PIC simulation, Eq. (29) is used as the inverse guiding-center transformation. This transformation needs to be performed numerically when we deposit markers to grids, or when we calculate the gyro-averaged field to be used in pushing guiding-centers.

[The equilibrium magnetic field we will consider has spatial scale length much larger than the thermal gyro radius ρ. In this case the difference between the values of e_∥(x)∕Ω(x) and e_∥(X)∕Ω(X) is negligible. The difference between equilibrium field values evaluated at X and x is usually neglected in gyrokinetic theory (except in deriving the gradient/curvature drift). Therefore it does not matter whether the above e_∥∕Ω is evaluated at x or X. What matters is where the perturbed fields are evaluated: at x or at X. The values of perturbed fields at x or at X are different and this is called the finite Larmor radius (FLR) effect.]

For later use, define ρ ≡−v × e_∥∕Ω, which is the vector gyro-radius pointing from the the guiding-center to the particle position.

2.2 Choosing velocity coordinates

The guiding-center transformations (27) and (29) involve the particle velocity v. It is the cross product between v and e_∥(x) or e_∥(X) that is actually used. Therefore, only the perpendicular velocity (which is defined by v_⊥ = v − v ⋅ e_∥) enters the transform. A natural choice of coordinates for the perpendicular velocity is (v_⊥,α), where v_⊥ = |v_⊥| and α is the azimuthal angle of the perpendicular velocity in the local perpendicular plane.

The parallel direction is fully determined by B₀(x), but there are degrees of freedom in choosing one of the two perpendicular basis vectors. In order to make the azimuthal angle α fully defined, we need to choose a way to define one of the two perpendicular directions. In GEM simulations, one of the perpendicular direction is chosen as the direction perpendicular to the magnetic surface, which is fully determined at each spatial point. (We need to define the perpendicular direction at each spatial location to make ∂α∕∂x|_v defined, which is needed in the Vlasov differential operators. However, it seems that terms related to ∂α∕∂x|_v are finally dropped due to that they are of higher order**check.)

In the following, α will be called the “gyro-angle” . [Note that, in the guiding-center coordinates (X,v_∥,v_⊥,α), α is a velocity coordinate rather than a spatial coordinate. When transformed back to particle coordinates, α affects both the velocity coordinate and the spatial coordinate. Consider a series of points in terms of guiding-center coordinates (X,v_∥,v_⊥,α) with (X,v_∥,v_⊥) fixed but with α changing. Using the inverse guiding-center transform (29), we know that the above points form a gyro-ring in space, i.e., α influences spatial location, in addtion to velocity.]

The gyro-angle is an important variable we will stick to because we need to directly perform averaging over this variable (with X fixed) in deriving the gyrokinetic equation. We have multiple options for the remaining velocity coordinates , such as  (v,v_∥), or (v,v_⊥), or (v_∥,v_⊥). In Frieman-Chen’s paper, the velocity coordinates other than α are chosen to be (𝜀,μ) defined by

(30)

and

(31)

where Φ₀(x) is the equilibrium (macroscopic) electrical potential. Choosing μ as one of the phase space coordinates is nontrivial because it turns out a constant of motion. And this choice seems to be important in sucessfully getting the final gyrokinetic equation (I need to check this).

Note that (𝜀,μ,α) is not sufficient in uniquely determining a velocity vector. An additional parameter σ = sign(v_∥) is needed to determine the sign of v_∥ = v ⋅e_∥. In the following, the dependence of the distribution function on σ is often not explicitly shown in the variable list (i.e., σ is hidden/suppressed), which, however, does not mean that the distribution function is independent of σ.

Another frequently used velocity coordinates are (μ,v_∥,α). In the following, I will derive the gyrokinetic equation in (𝜀,μ,α) coordinates. After that, I transform it to (μ,v_∥,α) coordinates.

One important thing to note about the above velocity coordinates is that they are defined relative to the local magnetic field. If the field itself is spatially varying (such as in tokamaks), the above velocity coordinates are also spatially varying for a fixed velocity v. Specifically, the following derivatives are nonzero:

(32)

2.3 Summary of the phase-space coordinate transform

The transform from particle variables (x,v) to guiding-center variables (X,𝜀,μ,α,σ) is given by

(33)

As mentioned above, the dependence of the distribution function on σ will not be explicitly indicated in the following.

2.4 Distribution function in terms of guiding-center variables

Denote the particle distribution function expressed in particle coordinates (x,v) by f_p, and the same distribution expressed in the guiding-center variables (X,𝜀,μ,α) by f_g. Then

(34)

where (x,v) and (X,𝜀,μ,α) are related to each other by the guiding-center transform (29). Equation (34) along the guiding-center transform can be considered as the definition of f_g.

As is conventionally adopted in multi-variables calculus, both f_p and f_g are often denoted by the same symbol, say f. Which set of independent variables are actually assumed is inferred from the context. This is one subtle thing for gyrokinetic theory in particular and for multi-variables calculus in general. (Sometimes, it may be better to use subscript notation on f to identify which coordinates are assumed. One example where this distinguishing is important is encountered when we try to express the diamagnetic flow in terms of f_g, which is discussed in Appendix G.)

In practice, f_g is often called the guiding-center distribution function whereas f_p is called the particle distribution function. However, they are actually the same distribution function expressed in different variables. The name “guiding-center distribution function” is misleading because it may imply that we can count the number of guiding-centers to obtain this distribution function but this implication is wrong.

2.5 Spatial gradient operator in guiding-center coordinates

Using the chain-rule, the spatial gradient ∂f_p∕∂x is written

(35)

From the definition of X, Eq. (27), we obtain

(36)

where I is the unit dyad. From the definition of 𝜀, we obtain

(37)

where E₀ = −∂Φ₀∕∂x. Using the above results, equation (35) is written as

(38)

As mentioned above, the partial derivative ∂∕∂x is taken by holding v constant. Since B₀ is spatially varying, v_⊥ is spatially varying when holding v constant. Therefore and are generally nonzero. The explicit expressions of these two derivatives are needed later in the derivation of the gyrokinetic equation and is discussed in Appendix J. For notation ease, define

(39)

and

(40)

then expression (38) is written as

(41)

Note that ∂∕∂X is a shorthand for

|𝜀,μ,α

i.e., it is taken by holding (𝜀,μ,α) constant (rather than holding v constant). For notation ease, ∂∕∂X is sometimes denoted by ∇_X or simply ∇.

2.6 Velocity gradient operator in guiding-center coordinates

Next, let us express the velocity gradient ∂f∕∂v in terms of the guiding-center variables. Using the chain rule, ∂f∕∂v is written

(42)

From the definition of X, we obtain

(44)

From the definition of μ, we obtain

(45)

From the definition of α, we obtain

(46)

where e_α is defined by

(47)

Using the above results, expression (42) is written

(48)

2.7 Time derivatives in guiding-center coordinates

In terms of the guiding-center variables, the time partial derivative ∂f_p∕∂t appearing in Vlasov equation is written as

(49)

where V = (𝜀,μ,α). Here ∂X∕∂t and ∂V∕∂t are not necessarily zero because the equilibrium quantities involved in the definition of the guiding-center transformation are in general time dependent. This time dependence is assumed to be very slow in the gyrokinetic ordering discussed later. In the following, ∂X∕∂t and ∂V∕∂t will be dropped, i.e.,

(50)

2.8 Final form of Vlasov equation in guiding-center coordinates

Using the above results, the Vlasov equation in guiding-center coordinates is written

3 δf form of Vlasov equation in guiding-center variables

3.1 Electromagnetic field perturbation

Since the definition of the guiding-center variables (X,𝜀,μ,α) involves the equilibrium fields B₀ and E₀, to further simplify Eq. (52), we need to separate electromagnetic field into equilibrium and perturbation parts. Writing the electromagnetic field as

(53)

and

(54)

then substituting these expressions into equation (52) and moving all terms involving the field perturbations to the right-hand side, we obtain

(58)

where v_E0 is defined by v_E0 = cE₀ × e_∥∕B₀, which is the E₀ × B₀ drift. Further note that

Using Eqs. (58),  (59), and  (57), the left-hand side of equation (55) is written as

[Equation (60), corresponds to Eq. (7) in Frieman-Chen’s paper[3]. In Frieman-Chen’s equation (7), there is a term

(Emac − E0) ⋅ v

where E_mac is a given macroscopic electric field introduced when defining the guiding-center transformation. In my derivation  E_mac is chosen to be equal to the equilibrium electric field, and thus the above term is zero.]

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

(61)

i.e.

(65)

which agrees with Eq. (21) discussed in Sec. 1.2.

3.2 Distribution function perturbation

Expand the distribution function f_g as

(66)

where F_g is assumed to be an equilibrium distribution function, i.e.,

(67)

Using Eqs. (66) and (67) in Eq. (61), we obtain an equation for δF_g:

(68)

3.3 Gyrokinetic ordering

To facilitate the simplification of the Vlasov equation in the low-frequency regime, we assume the following orderings (some of which are roughly based on experiment measure of fluctuations responsible for tokamak plasma transport, some of which can be invalid in some interesting cases.) These ordering are often called the standard gyrokinetic orderings.

Assumptions for equilibrium quantities

Define the spatial scale length L₀ of equilibrium quantities by L₀ ≈ F_g∕|∇_XF_g|. Assume that L₀ is much larger than the thermal gyro-radius ρ_i ≡ v_t∕Ω, i.e., λ ≡ ρ_i∕L₀ is a small parameter, where v_t = is the thermal velocity. That is

(69)

and

(70)

The equilibrium E₀ × B₀ flow, i.e.,

(71)

is assumed to be weak with

(72)

Assumptions for perturbations

We assume that the amplitudes of perturbations are small. Specifically, we assume

(73)

where δΦ is the perturbed scalar potential defined later in Eq. (79).

We consider low frequency perturbations with ω∕Ω ∼ O(λ¹), then

(74)

The perturbation is assumed to have a long wavelength (much longer than ρ_i) in the parallel direction

(75)

and have a short wavelength comparable to the thermal gyro-radius in the perpendicular direction

(76)

[Combining Eq. (75) and (76), we obtain

(77)

i.e., the parallel wave number is one order smaller than the perpendicular wave-number.]

In terms of the scalar and vector potentials δΦ and δA, the perturbed electromagnetic field is written as

(78)

and

(79)

Then

(80)

(81)

Using the above orderings, it is ready see that δE_∥ is one order smaller than δE_⊥, i.e.,

(82)

Most gyrokinetic simulations approximate the vector potential as δA ≈ δA_∥e_∥.

(83)

We will assume v_tδA_∥∼ δΦ.

3.4 Equation for macroscopic distribution function F_g

The evolution of the macroscopic quantity F_g is governed by Eq. (67), i.e.,

(84)

where the left-hand side is written as

(85)

i.e., F_g0 is independent of the gyro-angle α. The balance on O(λ¹) gives

(86)

Performing averaging over α, ∫ ₀^2π(…)dα, on the above equation and noting that F_g0 is independent of α, we obtain

(87)

Note that a quantity A = A(x) that is independent of v will depend on v when transformed to guiding-center coordinates, i.e., A(x) = A_g(X,v). Therefore A_g depends on gyro-angle α. However, since ρ_i∕L ≪ 1 for equilibrium quantities, the gyro-angle dependence of the equilibrium quantities can be neglected. Specifically, e_∥, B₀ and Ω can be considered to be independent of α. As to v_∥, we have v_∥ = ±. Since B₀ is considered independent of α, so does v_∥. Using these results, equation (87) is written

(88)

Using E₀ = −∇Φ₀, the above equation is written as

(89)

Note that

(90)

where the error is of O(λ²)Φ₀, and thus, accurate to O(λ), the last term of equation (89) is zero. Then equation (89) is written as

(91)

which implies that F_g0 is constant along a magnetic field line.

3.5 Equation for δF_g

Using F_g ≈ F_g0, equation (68) is written as

(92)

where δRδF_g is a nonlinear term which is of order O(λ²) or higher, L_gδF_g and δRF_g0 are linear terms which are of order O(λ¹) or higher. The linear term δRF_g0 is given by

(93)

In obtaining (93), use has been made of ∂F_g0∕∂α = 0. Another linear term L_gδF_g is written as

Next, to reduce the complexity of algebra, we consider the easier case in which ∂F_g0∕∂μ = 0.

Balance on order O(λ¹): adiabatic response

The balance between the leading terms (terms of O(λ)) in Eq. (92) requires that

(95)

where δF_a is a unknown distribution function to be solved from the above equation. It is ready to verify that

(96)

is a solution to the above equation, accurate to O(λ). [Proof: Substitute expression (96) into the left-hand side of Eq. (95), we obtain

Separate δF_g into adiabatic and non-adiabatic part

As is discussed above, the terms of O(λ) can be eliminated by splitting a so-called adiabatic term form δF_g. Specifically, write δF_g as

(101)

where δF_a is given by (96), i.e.,

(102)

which depends on the gyro-angle via δΦ and this term is often called adiabatic term. Plugging expression (101) into equation (92), we obtain

(103)

Next, let us simplify the linear term on the right-hand side, i.e, δRF_g0 −L_gfδF_a, (which should be of O(λ²) or higher because Ω∂δF_a∕∂α cancels all the O(λ¹) terms in δRF_g0).

L_gδF_a is written

The consequence of this is that, as we will see in Sec. 3.6, δG is independent of the gyro-angle, accurate to order O(λ¹). Therefore, separating δF into adiabatic and non-adiabatic parts also corresponds to separating δF into gyro-angle dependent and gyro-angle independent parts.

Linear term expressed in terms of δΦ and δA

Let us rewrite the linear term (108) in terms of δΦ and δA. The δE + v ×δB term in expression (108) is written as

(109)

Note that this term needs to be accurate to only O(λ). Then

(110)

where the error is of O(λ²). Using the vector identity v ×∇_x ×δA = (∇δA) ⋅v − (v ⋅∇)δA and noting v is constant for ∇_x operator, the above equation is written

(111)

Note that Eq. (41) indicates that ∇_xδΦ ≈∇_XδΦ, where the error is of O(λ²), then the above equation is written

(112)

Further note that the parallel gradients in the above equation are of O(λ²) and thus can be dropped. Then expression (112) is written

(114)

Using expression (113), equation (108) is written

(115)

where all terms are of O(λ²).

3.6 Equation for the non-adiabatic part δG

Plugging expression (115) into Eq. (103), we obtain

(116)

where L_g is given by Eq. (94), i.e.,

Expansion of δG

Expand δG as

where δG_i ∼ O(λⁱ⁺¹)F_g0, and note that the right-hand side of Eq. (116) is of O(λ²), then, the balance on order O(λ¹) requires

(118)

i.e., δG₀ is gyro-phase independent.

The balance on order O(λ²) requires (for the special case of E₀ = 0):

Gyro-averaging

Define the gyro-average operator ⟨…⟩_α by

(120)

where h = h(X,α,𝜀,μ) is an arbitrary function of guiding-center variables. The gyro-averaging is an integration in the velocity space. [For a field quantity, which is independent of the velocity in particle coordinates, i.e., h = h(x), it is ready to see that the above averaging is a spatial averaging over a gyro-ring.]

Gyro-averaging Eq. (119), we obtain

(122)

The second gyro-averaging on the left-hand side of Eq. (121) can be written as

(123)

where V_D is the magnetic curvature and gradient drift (Eq. (123) is derived in Appendix xx, to do later). Then Eq. (121) is written

Simplification of the nonlinear term

Next, we try to simplify the nonlinear term ⟨δRδF_g⟩_α appearing in Eq. (124), which is written as

(127)

Accurate to O(λ²),the first term on the right-hand side of the above is zero. [Proof:

(129)

Using the expression of δR given by Eq. (56), the above expression is written as

(131)

The other two terms in Eq. (130) can be proved to be zero. [Proof:

(134)

Using this in Eq. (129), we obtain

(135)

which is of O(λ²).

Final equation for the non-adiabatic part of perturbed distribution function

Using the above results, the gyro-averaged kinetic equation for δG₀ is finally written as

(137)

where the first term is called “adiabatic term”, which depends on the gyro-phase α via δΦ. Equation (136) is the special case (∂F_g0∕∂μ|_𝜀 = 0) of the Frieman-Chen nonlinear gyrokinetic equation given in Ref. [3]. Note that the nonlinear terms only appear on the left-hand side of Eq. (136) and all the terms on the right-hand side are linear. The term

(138)

consists of the δE × B₀ drift and magnetic fluttering term (refer to  expression (330) in Sec. C.3). For notaiton ease, this term is denoted by δV_D in the following:

(139)

Note that ∇≡ ∂∕∂X, which is taken in the guiding-center coordinates (X,𝜀,μ,α) while holding (𝜀,μ,α) constant. How do we numerically calculate ∂Φ_g∕∂X in PIC simulations? The difficulty is that X grid is not available, which makes it difficult to use the finite difference method. We note x grid is available in PIC simulations. It turns out we can make use of x grid to approximatly calculate the derivatives in X space. Using x = X + ρ and the definition δΦ_g(X,α,v_⊥) = δΦ_p(x) we obtain

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:

(141)

For notation ease, the subscripts g and p are usually dropped. Which one is intended should be obvious from the context.

4 Gyrokinetic equation suitable for numerical simulation

The gyrokinetic equation (136) contains time derivatives on the right-hand side, which is problematic if treated by using explicit finite difference in PIC simulations. Next, we discuss some methods that can eliminate these terms, making the gyrokinetic equation more amenable to PIC simulations.

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

Note that the coefficient before ∂F₀∕∂𝜀 in Eq. (136) involves the time derivative of ⟨δϕ⟩_α, which is problematic if treated by using explicit finite difference (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 defining another gyro-phase independent function δf by

(142)

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

Using

δVD = −∇⟨δΦ − v ⋅ δA⟩×

in the right-hand side of Eq. (144) yields

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

Similar to the method of eliminating ∂⟨δϕ⟩∕∂t, we define another gyro-phase independent function by

(146)

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

(147)

Then expression (146) is written as

(148)

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

(152)

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 (147), yielding

(153)

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

4.3 Mixed-variable pullback method[4]

Define δf^(h) by

(157)

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

(158)

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

(159)

Begining with Eq. (144) and following the same procedure in Sec. 4.2, we obtain the following evoluation equation for δf^(h):

(161)

where δJ_||j^(h) is the parallel current carried by the distrbution funciton δf_j^(h), where the subscript j is species index. 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 A_∥^(s) 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 significant. While a clever 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:

(162)

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

(163)

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 influence the electromagnetic field. Meanwhile, the definition of Eq. (157) indicates that, for a given δf, the re-splitting will make the value of δf^(h) increase by Δ = ⟨v_∥δA_∥old^(h)⟩ _α since δA^(h) is reduced by δA_old^(h). To keep δf unchanged, which is what we need to ensure, we must add Δ to δf^(h):

(164)

After the operations in (162) (163) and (164), both the electromagnetic field and the distribution function δf remain unchanged. Therefore we do not change any physical state of the system. Hence the above operations can only have effects on numerical round-off errors. As mentioned above, we hope the numerical effects will reduce the cancellation errors.

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

(165)

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

(166)

and

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

Weight evolution equation

The weight of the j th marker is defined by

(169)

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

In field-aligned coordinates

In field-aligned coordinates (x,y,z),   can be written as

(171)

Note that the direction of ∂r∕∂x at 𝜃 = −π is different 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 (171) is written as

(172)

where we approximate ∇x, ∇y, and ∇z as constants when performing the gyro-averge since they are determined by the equilibrim magnetic field, 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 field equations are not solved at 𝜃 = +π and hence the field values are not directly obtained. Instead, the field values at 𝜃 = +π are infered from the field values at 𝜃 = −π using the coninuity of ∇δΦ. At the 𝜃 cut, the continuity of ∇δΦ requires

(173)

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

(174)

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

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

(175)

(176)

Using these, Eq. (170) is written as

(178)

Then, a typical term of the above is written as

The pullback:

Using

(183)

i.e.,

(184)

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

The magnetic fluttering term is written as

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

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

(190)

(191)

Define

(193)

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

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

(194)

where use has been made of μ₀ = 1∕(c²𝜀₀).

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

= −b ⋅∇δΦ,

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

(196)

i.e.,

(197)

4.4 Summary of distribution function split

In simulations, the seemingly trivial thing on how to split the distribution function is often a big deal. Separating the perturbation from the total distribution function gives rise to the the famous “δ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 ∂δ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

(198)

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

(199)

where δh satisfies the gyrokinetic equation (151) or (156). In PIC simulations, δh is evolved by using markers and its moment in the phase-space is evaluated via Monte-Carlo integration. The blue and red terms explicitly depends on the unknown perturbed field. After being integrated in the velocity space, these two terms give  the polarization density and the skin current, respectively. The polarization density is discussed in Sec. 6. The skin current is discussed in Sec. (4.5) and the so-called “cancellation problem” is discussed in Sec. 5.

4.5 Skin current

Let us calculate the velocity space moment of ⟨v ⋅δA⟩. Consider the approximation δA ≈ δA_∥e_∥, then the blue term in Eq. (199) is written as

(200)

Notice that v_∥ can be taken out of the gyro-averaging. Then the above equation is written

(201)

If we neglect the FLR effect, then the above expression is written

(202)

The zeroth order moment (number density) is then written as

(203)

which is zero if F₀ is Maxwellian. Next, consider the parallel current carried by distribution (202), which is written

(204)

If F₀ is a Maxwellian distribution given by

(205)

then

(206)

Then expression (204) is written

(207)

Working in the spherical coordinates, then v_∥ = v cos𝜃 and dv = v² sin𝜃dvd𝜃dϕ. Then expression (207) is written

(209)

where the c∕ω_ps is often called the skin depth of species “s”. This current δj_∥ is often called “skin current”, which can also be written as

(210)

where ρ = ∕Ω is the thermal Larmor radius.

5 Parallel Ampere’s Law

(211)

where the parallel currents are given by

(212)

(213)

where δJ_∥i′ and δJ_∥e′ is the parallel current carried by the distribution function δh in Eq. (199), which are updated from the value at the nth time step using an explicit scheme and therefore does not depends on the field at the (n + 1)th step. The blue terms  in Eqs. (212) and (213) are called “skin current”, which depend on the unknown field at the (n + 1)th step and thus need to be moved to the left-hand side of Ampere’s law (211) if we want to solve this equation by direct methods. In this case, equation (211) is written as

Therefore we go back to Eq. (211) and try to solve it using iterative methods. However, it is found numerically that directly using Eq. (211) as an iterative scheme is usually divergent. To obtain a convergent iterative scheme, we need to have an approximate form for the blue terms, which is independent of markers and so that it is easy to construct its matrix, and then subtract this approximate form from both sides. After doing this, the iterative scheme has better chance to be convergent (partially due to that the right-hand side becomes smaller). An approximate form is that derived by neglecting the FLR effect given in Sec. 4.5. Using this, the iterative scheme for solving Eq. (211) is written as

Because the ion skin current is less than its electron counterpart by a factor of m_e∕m_i, its accuracy is not important. The cancellation error is not a problem and hence can be neglected. In this case, equation (215) is simplified as

5.1 Discretizing Laplacian operator

In TEK, the Laplacian operator is approximated as

(217)

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:

(218)

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

(220)

into expression (219) gives

(222)

6 Poisson’s equation and polarization density

Poisson’s equation is written as

(223)

where −𝜀₀∇²δΦ is called the space-charge term. Since we consider modes with k_∥≪ k_⊥, the space-charge term is approximated as ∇²δΦ ≡∇_⊥²δΦ + ∇_∥²δΦ ≈∇_⊥²δΦ. Then Eq. (223) is written as

(224)

This approximation eliminates the parallel plasma oscillation from the system. The perpendicular plasma oscillations seem to be only partially eliminated in the system consisting of gyrokinetic ions and drift-kinetic electrons. There are the so-called Ω_H modes (also called electrostatic shear Alfven wave) that appear in the gyrokinetic system which have some similarity with the plasma oscillations but with a much smaller frequency, Ω_H ∼ (k_∥∕k_⊥)(λ_D∕ρ_s)ω_pe.

Using expression (199), the perturbed number density δn is written as

(226)

which has an explicit dependence on δΦ and is usually moved to the left hand of Poisson’s equation when constructing the numerical solver of the Poisson equation, i.e., equation (224) is written as

(227)

where δn_i′ = δn_i −δn_pi = ∫ δh_idv, which is evaluated by using Monte-Carlo markers. Since some parts depending on δΦ are moved from the right-hand side to the left-hand side of the field equation, numerical solvers (for δΦ) based on the left-hand side of Eq. (227) probably behaves better than the one that is based on the left-hand side of Eq. (224), i.e., −𝜀₀∇_⊥²δΦ.

6.1 Discussion on cancellation scheme

The polarization density is part of the perturbed density that is extracted from the source term and moved to the left-hand side of the Poisson equation. The polarization density will be evaluated without using Monte-Carlo markers, whereas the remained density on the right-hand side will be evaluated using Mote-Carlo markers. The two different methods of evaluating two parts of the total perturbed density can possibly introduce significant errors if the two terms are expected to cancel each other and give a small quantity that is much smaller than either of the two terms. This is one pitfall for PIC simulations that extract some parts from the source term and move them to the left-hand side. To remedy this, rather than directly moving a part of the distribution function to the left-hand side, we subtract an (approximate) analytic expression from both sides of Eq. (224). The analytical expressions on both sides are evaluated based on grid values of perturbed electromagnetic fields and are independent of makers. All the original parts of the distribution functions are kept on the right-hand side and are still evaluated by using markers, which hopefully avoids the possible cancellation problem. This strategy is often called a cancellation scheme. Since unknown perturbed electromagnetic fields appear on the right-hand side, iteration is needed to solve the field equation.

Note that two things appear here: What motivates us to move parts of the distribution function to the left? It is the goal of hopefully making the left-hand side matrix more well-behaved (such as good condition number, etc.) Why do we need the cancellation scheme? Because we want to avoid the numerical inaccuracy that appears when large terms cancels each other. Note that iteration is needed when the cancellation scheme is used because the right-hand side explicitly contains unknown electromagnetic fields.

It turns out that the cancellation scheme is not necessary for Eq. (227), but for the field solver for Ampere’s equation (discussed later), this cancellation scheme is necessary in order to obtain stable results.

7 Polarization density expressed as local Fourier expansion

Since δΦ is independent of the velocity in the particle coaordinates, the first term (adiabatic term) in expression (226) is trivial and the velocity integration can be readily performed (assume F₀ is Maxwellian), giving

(229)

7.1 Gyro-averaging of δΦ in guiding-center coordinates

In order to perform the gyro-averaging of δΦ, we Fourier expand δΦ in space as

(230)

and then express x in terms of the guiding center variables (X,v) since the gyro-averaging is taken by holding X rather than x constant. The guiding-center transformation gives

(231)

Using expressions (230) and (231), the gyro-average of δΦ is written as

(235)

7.2 Gryo-angle integration in particle coordinates

Next, we need to perform the integration in velocity space, which is done by holding x (rather than X) constant. Therefore, it is convenient to transform back to particle coordinates. Using X = x + v ×, expression (235) is written as

(236)

Then the velocity integration is written as

(238)

Since this is at x rather than X, k_⊥, v_⊥, and Ω are different from those appearing in expression (233). However, since this difference is due to the variation of the equilibrium quantity e_∥∕Ω in a Larmor radius, and thus is small and is ignored in the following.

Plugging expression (238) into expression (237) and using dv = v_⊥dv_⊥dv_∥dα, we get

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

In order to perform the remaining velocity integration in expression (240), we assume that F₀ is a Maxwellian distribution given by

(243)

Again we will 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 sufficiently large velocity, the corresponding Larmor radius will be large enough to make the equilibrium undergo substantial variation. Since the velocity integration limit is to infinite, this will definitely 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. (243), the expression in the square brackets of Eq. (240) is written as

(246)

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

(247)

Perpendicular integration Using (I verified this by using Sympy)

(248)

where I₀(a) is the zeroth modified Bessel function of the first kind, expression (247) is written

(249)

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

(250)

Final form of polarization density

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

(251)

Plugging expression (250) and (251) into expression (226), the polarization density n_p is written as

(253)

then Eq. (252) is written as

(254)

Expression (254) 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 Γ₀.

7.3 Pade approximation

Γ₀ defined in Eq. (253) can be approximated by the Pade approximation as

(255)

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

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

Long wavelength approximation of the polarization density

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

Equation (258) 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 (257) is for the long wavelength limit, which partially neglects FLR effect. Let us go back to the more general expression (256). The Poisson equation is written

(259)

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

(260)

Fourier transforming in space, the above equation is written

(261)

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

(262)

Using this, equation (261) is written

(263)

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

(264)

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

(265)

Neglecting the Debye shielding term, the above equation is written

(266)

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

8 Polarization density matrix obtained by numerically integrating in phase space using grid

In Sec. 7, to evaluate the polarization density, the potential δΦ is Fourier expanded in space using local Cartesian coordinates, and then the double gyro-angle integration of each harmonic is expressed as the Bessel function. (It seems that the original motivation of using the Fourier expansion here is to facilitate analytical treatment and is not designed for numerical use. GEM code does make use of the local Fourier expansion in its numerical implementation, where the local perpendicular wave number needs to be estimated numerically, which seems awkward.)

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. The polarization density is given by Eq. (226), i.e.,

(267)

Our goal is to write n_p in the following real-space discrete from:

(268)

8.1 Direct evaluation of the double gyrophase integration

Define

(269)

Using dv = v_⊥dv_⊥dv_∥dα, the above integration is written as

(270)

where use has been made of the assumption that ∂F₀∕∂𝜀 is uniform in α in (x,v_⊥,v_∥,α) coordinates, and thus is moved outside of the α integration. Using the definition of gyro-averaging (2π)⁻¹ ∫ ₀^2π(…)dα′, the above integration is written as

(271)

Note that the gyro-averaging is performed in the guiding-center space, i.e., performed by varying the gyroangle α′ while keeping guiding-center position X, v_⊥, and v_∥ constant. Using the definition of δΦ_g (i.e., its relation with δΦ):

(272)

where the particle location x is computed from the guiding-center location by

(273)

then A(x) is written as

Next, in order to perform the remaining velocity space integration, transform back to the particle coordinates (because the velocity integration is performed in the particle coordinates, i.e., it is performed by keeping the particle coordinate x constant):

(275)

where Δρ_ij is a function of (x,v_⊥,α_i,α_j′). Then Eq. (274) is written as

(277)

where N₁ = 4,N₂ = 4 for the case shown in Fig. 2. This gives how to approximate the double gyro-angle integration using the discrete values of δΦ_ij. The spatial points x_ij appearing in Eq. (277) are not necessarily grid points. Linear interpolations are used to express δΦ(x_ij) as linear combination of values of δΦ at grid-points.

8.2 Performing the parallel velocity integration

Assume that F₀ is Maxwellian, then

(278)

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

(280)

8.3 Toroidal DFT of polarization density in field-aligned coordinates

In field-aligned coordinates (x,y,z), Fourier expansion of δΦ along y is written

(281)

where ι = . Use this in Eq. (279), yielding

(284)

which should be a good approximation since the variation of δΦ along a field line over a distance of a Larmor radius is small. Then expression (283) is written as

8.4 Using MC integration, to be continued, not necessary

The integration we try to express is given by

A Adiabatic electron response

Assume that the total electron density satisfies the Boltzmann distribution on each magnetic surface, i.e.,

(287)

where N_e is a radial function. Note that this does not imply that the equilibrium density is N_e (it just implies that the total density is N_e at the location where δΦ = 0, which can still be different from the equilibrium density).

Further assume that the magnetic surface average of δn_e = n_e − n_e0 is zero, i.e.,

(288)

where n_e0 is the equilibrium electron density, ⟨…⟩_f is the magntic surface averaging operator. Using Eq. (287) in the above condition, we get

(289)

Then expression (287) is written as

(290)

Then the perturbation δn_e = n_e − n_e0 is written as

(291)

This model of electron response is often called adiabatic electron.

A.1 Poisson’s equation with adiabatic electron response

Pluging expression (291) into the Poisson equation (227), we get

(292)

When solving the Poisson equation, the equation is Fourier expanded in toroidal harmonics and each harmonic is independent of each other, so that they can be solved independently. For n≠0 harmonics, the ⟨δΦ⟩_f terms is zero and thus it is trivial to treat the electron term. For the n = 0 harmonic, the ⟨δΦ⟩_f term is nonzero and needs special treatment. I use the following method to obtain ⟨δΦ⟩_f. First slove the n = 0 harmonic of the following equation:

(293)

(i.e., Eq. (292) with electron term dropped), and then take the magnetic surface average of the solution δΦ′ to get ⟨δΦ′⟩_f. It can be proved that ⟨δΦ′⟩_f is equal to ⟨δΦ⟩_f. Then solving Eq. (292) becomes easier since ⟨δΦ⟩_f term is known and can be moved to the right-hand side as a source term.

B Characteristic curves of Frieman-Chen nonlinear gyrokinetic equation

Examining the left-hand side of Eq. (), it is ready to find that the characteristic curves of this equation are given by the following equations:

(294)

(295)

(296)

(It is instructive to notice that the kinetic energy 𝜀 is conserved along the characteristic curves while the real kinetic energy of a particle is usually not conserved in a perturbed electromagnetic field. This may be an indication that Frieman-Chen equation neglects the velocity space nonlinearity.) For notation ease, we denote the total guiding-center velocity by V_G, i.e.,

(297)

where

(298)

In the above, v_∥ can be obtained by using v_∥ = σ. However σ is not easy to determine since it changes with time for trapped particles. In practice, we get v_∥ from an evolution equation. This evolution equation for v_∥ can be obtained by combining the equations for 𝜀 and μ. Next, we derive this equation.

B.1 Time evolution equation for v_∥

Using the definition μ = v_⊥²∕(2B₀), equation (296), i.e., dμ∕dt = 0, is written as

(299)

which is written as

(300)

which can be further written as

(301)

Using the definition of the characteristics, the right-hand side of the above equation can be expanded, giving

(302)

where dX∕dt, d𝜀∕dt, and dμ∕dt are given by Eq. (294), (295), and (296), respectively. Using Eqs. (294)-(296) and ∂B₀∕∂t = 0, equation (302) is reduced to

(303)

On the other hand, equation (295), i.e., d𝜀∕dt = 0, is written as

(304)

which can be further written as

(305)

Using Eq. (303), the above equation is written as

(306)

which is the equation for the time evolution of v_∥. This equation involves dX∕dt, i,e., the guiding-center drift, which is given by Eq. (294). Equation (306) for v_∥ can be simplified by noting that the Frieman-Chen equation is correct only to the second order, O(λ²), and thus the characteristics need to be correct only to the first order O(λ) and higher order terms can be dropped. Note that, in the guiding-center drift dX∕dt given by Eq. (294), only the v_∥e_∥ term in  is of order O(λ⁰), all the other terms are of O(λ¹). Using this, accurate to order O(λ¹), equation (306) is written as

(307)

which is the time evolution equation ready to be used for numerically advancing v_∥. Note that only the mirror force −μe_∥⋅∇B appears in Eq. (307) and there is no parallel acceleration term qv_∥δE_∥∕m in Eq. (307). This is because δE_∥ = −b ⋅∇δΦ − ∂δA_∥∕∂t is of order O(λ²) and (**check** the terms involving E_∥ are of O(λ³) or higher and thus have been dropped in deriving Frieman-Chen equation.)

C From (δΦ,δA) to (δE,δB)

C.1 Expression of δB_⊥ in terms of δA

Note that

(311)

Note that the parallel gradient operator ∇_∥≡ e_∥⋅∇ = ∂∕∂z acting on the the perturbed quantities will result in quantities of order O(λ²). Retaining terms of order up to O(λ), equation (311) is written as

(312)

i.e., only the parallel component survive, which exactly cancels the last term in Eq. (310), i.e., equation (310) is reduced to

(313)

Using basis vectors in field-aligned coordinates

In terms of δB_xL and δB_yL,  δB_⊥ is written as

(314)

Dotting the above equation by ∇x and ∇y, respectively, we obtain

(315)

(316)

Equations (315) and (316) can be further written as

(317)

and

(318)

The solution of this 2 × 2 system is expressed by Cramer’s rule in the code.

Use B₀ = Ψ′∇x ×∇y

b = Ψ′∇x ×∇y∕B₀

C.2 Expression of δB_∥ in terms of δA

(321)

Note that the parallel gradient operator ∇_∥≡ e_∥⋅∇ = ∂∕∂z acting on the the perturbed quantities will result in quantities of order O(λ²). Retaining terms of order up to O(λ), equation (311) is written as

(322)

Using this, equation (320) is written as

(323)

However, this expression is not useful for GEM because GEM does not use the local coordinates (r,ϕ,z).]

C.3 Expressing the perturbed drift in terms of δE and δB

The perturbed drift δV_D is given by Eq. (139), i.e.,

(324)

Using δL = δΦ − v ⋅ δA, the above expression can be further written as

Accurate to O(λ), the last term on the right-hand side of expression (325) is written

Using Eqs. (326), (328), and (329), expression (325) is finally written as

(330)

Using this, the first equation of the characteristics, equation (294), is written as

C.4 Expressing the coefficient before ∂F₀∕∂𝜀 in terms of δE and δB

[Note that

(333)

where ∂δA_⊥∕∂t is of O(λ²). This means that δE_⊥ + ∇_⊥δϕ is of O(λ²) although both δE_⊥ and δϕ are of O(λ).]

Note that

D Split-weight scheme for electrons in GEM code

The perturbed distribution function is decomposed as given by Eq. (199), i.e.,

(337)

where the term in blue is the so-called adiabatic response, which depends on the gyro-angle in guiding-center coordinates. Recall that the red term ⟨δΦ⟩_α, which is independent of the gyro-angle, is introduced in order to eliminate the time derivative ∂⟨δΦ⟩_α∕∂t term on the right-hand side of the original Frieman-Chen gyrokinetic equation.

The so-called generalized split-weight scheme corresponds to going back to the original Frieman-Chen gyrokinetic equation by introducing another ⟨δΦ⟩_α term with a free small parameter 𝜖_g. Specifically, δh in the above is split as

(338)

(If 𝜀_g = 1, then the two ⟨δΦ⟩_α terms in Eq. (337) and (338) cancel each other.) Substituting this expression into Eq. (), we obtain the following equation for δh_s:

special case of 𝜖_g = 1

For the special case of 𝜖_g = 1 (the default and most used case in GEM code, Yang Chen said 𝜖_g < 1 cases are sometimes not accurate, so he gave up using it since 2009), equation (340) can be simplified as:

For 𝜖_g = 1, δF is written as

Equation (341) actually goes back to the original Frieman-Chen equation. The only difference is that ⟨v ⋅δA⟩_α is further split from the perturbed distribution function. Considering this, equation (341) can also be obtained from the original Frieman-Chen equation (136) by writing δG₀ as

(343)

In this case, δF is written as

(344)

Substituting expression (343) into equation (136), we obtain the following equation for δh_s:

In GEM, the split weight method is used only for electrons, and ∂⟨δΦ⟩_α∕∂t is approximated by ∂δΦ∕∂t, which is obtained from the vorticity equation (rather than from time-difference scheme).

When using the split weight scheme, a ∂δϕ∕∂t terms appear in the right-hand-side of the weight evolution equation. GEM makes use of the vorticity equation (time derivative of the Poissson equation) to evaluate ∂δϕ∕∂t.

E Coordinate system and grid in TEK code

In (ψ,𝜃,ϕ) coordinates:

• ψ is the normalized poloidal magnetic flux, ψ = (Ψ − Ψ_axis)∕(Ψ_lcfs − Ψ_axis).
• 𝜃 ∈ [−π,π) is increasing along the anti-clockwise direction when viewed along ∇ϕ direction with 𝜃 = −π at the high-field-side midplane.
• ϕ is the toroidal angle of the right-handed cylindrical coordinates (R,ϕ,Z).

In this convention, the Jacobian of the (ψ,𝜃,ϕ) coordinate system, 𝒥 = (∇ψ ⋅∇𝜃 ×∇ϕ) is negative, i.e., (ψ,𝜃,ϕ) is a left-handed system. The field-line-following coordinate system (ψ,𝜃,α) is also a left-handed system. The coordinate system (x = ψ,y = α,z = 𝜃) is a right-handed system.

E.1 Poloidal grid

I do not need to make connection with GEM’s equilibrium poloidal array because there is no coupling of equilibrium quantities between the code written by me and the original code in GEM. The coupling happens for the perturbed quanties, whose poloidal grids need to be consistent.

• The gridpoint number along 𝜃 for perturbations is denoted by mpol2. Its value is not chosen by users. Instead, it is equal to numproc/ntube, i.e., mpol2 is the total number of cells along 𝜃 for MPI parallelation. (numproc and ntube must be chosen in a way that makes  numproc/ntube an integer.)
• Poloidal grid-points for perturbations are indexed as 0:mpol2-1, with 0 corresponding to 𝜃 = −π and mpol2 corresponding to 𝜃 = π. Field equations are solved at 0 : mpol2 − 1. (The field at mpol2, i.e., 𝜃 = +π, is not identical to that at 𝜃 = −π due to toroidal shift of a magnetic field line when it finish one poloidal loop, and is obtained by toroidally inerpolating the fields at 𝜃 = −π.) No array for this is actually used in the code. The indexes just correspond to gclr of MPI processors.
• The equilibrium poloidal gridpoint number, mpol, must be chosen in a way that makes (mpol-1)/mpol2 an integer, which can be larger than 1, i.e., one simulation cell can span several equilibrium cells.

E.2 Toroidal grid

E.3 Radial grid

GEM array index system is better than TEK’s because my system is not consistent: sometimes I use 0-based index and sometimes I use 1-based index, sometimes the index ends at n and sometimes ends at n+1. It is important to know accurately the transformation between the two systems.

F Implementation of gyrokinetics in particle-in-cell (PIC) codes

F.1 Monte-Carlo evaluation of distribution function moment at grid-points

Suppose that the 6D guiding-center phase-space (X,v) is described by (ψ,𝜃,ϕ,v_∥,α,v_⊥) coordinates, where α is the gyro-angle. Denote the Jacobian of the coordinate system by 𝒥_rv.

Suppose that we sample the 6D phase-space by using a probability function P(ψ,𝜃,ϕ,v_⊥,α,v_∥). (Then the effecitive probability function used in rejection method is P𝒥_r𝒥_v.) Note that we will sample a distribution function δf_g(X,v_⊥,v_∥,α) that happens to be independent of the gyro-angle α. Since δf_g is uniform in α, it is natural to smaple α using a uniform probability function, i.e., P(X,v_⊥,α,v_∥) can be chosen to be independent of α. Then the weight is also independent of α. In terms of sampling δf_g, there is no need to specify α because both f_g and marker distribution g = N_pP are indepedent of α. However, sampling α is still necessary in simulations. From the perspective of programming, it is ready to understand why α needs to be specified: when computing density at grid-points, we need to compute particle locations from the guiding-center locations which requires us to specify the gyro-angle α for each marker. Specifically, markers in guiding-center space (X,v_⊥,α,v_∥) with fixed (X,v_⊥,v_∥) but different values of α correspond to markers in particle space with different locations (and different α′). In the PIC method of computing δn or δj_∥ at fixed x, the contribution of each particle marker to the density is affected by the distance of the particle marker to the fixed x. Therefore different values of α contributes differently to δn or δj_∥, and thus the resolution of α matters.

For a marker with coordinators (X,v_⊥,v_∥,α), the corresponding particle position can be calculated by using the inverse guiding-center transformation (29). Then we can deposit the marker weight to grid-points in the same way that we do in conventional PIC simulations. Looping over all markers, we build the particle density at grids.

Compared with conventional PIC methods, where particle positions are directly sampled, what is the benefit of sampling guiding-center positions and then transform them back to the particle positions? More computations are involved since we need to numerically perform the inverse guiding-center transformation. The answer lies in the important fact that the distribution function f_g(X,v_⊥,v_∥,α) that needs to be numerically evaluated in gyrokinetic simulation is actually independent of the gyro-angle α. Furthermore, we use a probability density function P that is independent of α to sample the 6D phase space (X,v_⊥,v_∥,α). Then the marker weight w ≡ f_g∕(N_pP) is independent of α, where N_p is number of markers loaded. This fact enable resolution in particle coordinates can be increased in a trivial way, as described below.

Suppose we have a concrete sampling of f_g in the 6D phase-space (X,v_⊥,v_∥,α), i.e., (X_j,v_⊥j,v_∥j,α_j,w_j) with j = 1,…,N_p, then we can do the inverse guiding-center transform and then deposit particles to the grids to obtain a moment (e.g. density).

Since P𝒥_rv is independent of α, the sampling α_j with j = 0,1,…,N_p are uniform distributed random numbers. Therefore we can generate another sampling of α (denoted by α_j′ with j = 1,…,N_p) using random number generators and combine α_j′ with the old sampling (X_j,v_⊥j,v_∥j) to obtain a new sampling (X_j,v_⊥j,v_∥j,α_j′,w_j′).  The values of w at the new sampling points are equal to the original values, i.e., w_j′ = w_j since the particle weight w = f_g∕(N_pP) is independent of α.  Using the new sampling and following the same procedures given above, we can estimate values of the moment at grid-points again, which will differ from the estimation obtained using the old sampling because the resulting sampling of x and α′ for f_p is different. Taking the average of the two estimations will give a more accurate estimation because the resolution in x and α′ for f_p is increased.

[Note that even if P𝒥_rv is dependent on α due to the possible dependence of 𝒥_rv on α, we still can easily generate a new set of sampling which differs from the old sampling only in α. Specifically, in the rejection method, we use the old sampling (X_j,v_⊥j,v_∥j) for each reject step and only adjust α to satisfy the acceptance criteria.]

We can also construct a new sampling by replacing α_j with α_j + Δ, where Δ is a constant. Then the new set of sampling (X_j,v_⊥j,v_∥j,α_j + Δ,w_j) with j = 1,…,N_p is still a consistent set of sampling.

In doing the deposition, each marker has a single gyro-angle. Due to the independence of the weight (w = f_g∕(N_pP)) of the gyro-angle, the particle phase space resolution can be increased in a way that there can be several gyro-angles for a single (X_j,v_⊥j,v_∥j). This gives the wrong impression that the gyro-angle of a guiding-center marker is arbitrary. The correct understanding is that given above, i.e., we do several (say 4) separate sampling of the phase space and then average the results.

In the code, the gyro-angle is defined relative to the direction ∇ψ∕|∇ψ| at the guiding-center position. Specifically, the gyro-angle is defined as the included angle between ∇ψ∕|∇ψ| and −v × e_∥(x)∕Ω(x), where e_∥(x)∕Ω(x) is approximated by the value at the guiding-center location.

Then (X_j,v_⊥j,v_∥j) can be evolved by using the guiding-center motion equation. (It is obvious how to evalute the gyro-averaging of the electromagnetic fields needed in pusing markers.)

F.2 Monte-Carlo sampling of 6D guiding-center phase-space

Suppose that the 6D guiding-center phase-space (X,v) is described by (ψ,𝜃,ϕ,v_∥,v_⊥,α) coordinates. The Jacobian of the coordinate system is given by 𝒥 = 𝒥_rv_⊥, where 𝒥_r = 𝒥 (ψ,𝜃) is the Jacobian of the coordinates (ψ,𝜃,ϕ). Suppose that we sample the 6D phase-space by using the following probability density function:

(347)

where V _r is the volume of the spatial simulation box, T is a constant temperature. P given above is independent of ψ,𝜃,ϕ and α. It is ready to verify that the above P satisfies the following normalization condition:

(348)

I use the rejection method to numerically generate N_p markers that satisfy the above probability density function. [The effective probability density function actually used in the rejection method is P′, which is related to P by

(349)

Note that P′ does not depend on the gyro-angle α.]

Then the weight of a marker sampling δf_g is written as

(350)

Since both f_g and P are independent of the gyro-angle α, w is also independent of α.

The numerical representation of δf_g is written

(351)

Although the distribution function δf_g to be sampled is independent of the gyro-angle α, we still need to specify the gyro-angle because we need to use the inverse guiding-center transformation, which needs the gyro-angle. Each marker needs to have a specific gyro-angle value α_j so that we know how to transform its X_j to x_j and then do the charge deposition in x space.

To increase the resolution over the gyro-angle (the quantity of intest to us, e.g., density, is the integration of the particle distribution function at fixed spatial location, and the distribution at the fixed location is not uniform in the gyro-angle. So the resolution over α matters), we need to load more markers. However, thanks to the fact that both sampling probability density function P and δf_g are independent of α, the resolution over the gyro-angle can be increased in a simple way.

(352)

This corresponds to sampling the 6D phase-space 4 separate times (each time with identical sampling points in (X,v_∥,v_⊥) but different sampling points in α) and then using the averaging of the 4 Monte-Carlo integrals to estimate the exact value. This estimation can also be (roughly) considered as a Monte-Carlo estimation using 4 times larger number of markers as that is originally used (the Monte-Carlo estimation using truly 4 time larger number of markers is more accurate than the result we obtained above because the former also increase the resolution of (X,v_∥,v_⊥) while the latter keeps the resolution of (X,v_∥,v_⊥) unchanged.)

In numerical implementation, we choose N sampling points that are evenly distributed on the gyro-ring (N is usually 4 as a compromise between efficiency and accuracy). And the weight is evenly split by the N sub-markers on the gyro-ring. Therefore each sub-marker have a Monte-Carlo weight w_j∕N, where w_j is the weight of the j th marker. Then calculating the integration (358) at a grid corresponds to depositing all the N sub-markers associated with each guiding-center marker to the grid, as is illustrated in Fig. 6. However, interpreting in this way is confusing to me because: why can you split a single sampling into N different samplings? I prefer the above interpretation that the 6D phase space is sampled 4 separate times and thus we get 4 estimations and finally we take the averaging of these 4 estimations. It took a long time for me to finally find this way of understanding.

In summary, the phase-space to be sampled in gyrokinetic simulations are still 6D rather than 5D. In this sense, the statement that gyrokinetic simulation works in a 5D phase space is misleading. We are still working in the 6D phase-space. The only subtle thing is that the sixth dimension, i.e., gyro-angle, can be sampled in an easy way that is independent of the other 5 variables.

In numerical implementation, the gyro-angle may not be explicitly used. We just try to find 4 arbitrary points on the gyro-ring that are easy to calculate. Some codes (e.g. ORB5) introduces a random variable to rotate these 4 points for different markers so that the gyro-angle can be sampled less biased.

From the view of particle simulations, the gyrokinetic model can be considered as a noise reduction method, where the averaging over the gyro-angle is equivalent to a time averaging over a gyro-period, which reduces the fluctuation level (in both time and space) associated with evaluating the Monte-Carlo phase integration. Here the averaging in gyrokinetic particle simulation refers to taking several points on a gyro-ring when depositing markers to spatial grids to obtain the density and current on the grids. (Another gyro-averaging appears in evaluating the guiding-center drift.) In gyrokinetic particle simulation, even a step size smaller than a gyro-period is taken, the quantities used in the model is still the ones averaged over one gyro-period. In this sense, a gyrokinetic simulation is only meaningful when the time step size is larger than one gyro-period. [**Some authors may disagree with that the gyro-averaging is a time-averaging. They may consider the gyro-averaging as the phase-space integration over the gyro-angle coordinate. This view seems to be right in Euler simulations but seems to be wrong in particle simulations. The reason is as follows. For each marker, choose a random gyro-phase and then do the inverse transformation to obtain particle position, and sum over all markers (this corresponds to phase-space Monte-Carlo integration, which include the gyro-angle integration, so no further gyro-angle integration is needed); choose another random gyro-phase and repeat the above procedure (this can be interpreted as do the phase-space Monte-Carlo integration at another time), choose further random gyro-phase for each marker and repeat. Finally averaging all the above values to obtain the final estimation of the phase-space integration. This amounts to a time-averaging over a gyro-motion. In summary, sampling several times with different gyro-phases for each marker and taking the average amounts to the time averaging over gyro-motion**]

When doing the time-average over the ion cyclotron motion, the time variation of the low-frequency mode is negligible and only the spatial variation of the modes is important. For the gyro-motion, only the gyro-angle is changing and all the other variables, (X,v_⊥,v_∥),  are approximately constant. As a result, this time averaging finally reduces to a gyro-averaging.

I prefer to reason in terms of particle position and velocity, and consider the guiding-center location as an image of the particle position. When working in the guiding-center coordinates, I prefer to reason by transforming back to the particle position. This reasoning is clear and help me avoid some confusions I used to have.

In coding, the initial sampling of the phase space is performed for particles ((v,x) of particles) and then is transformed to guiding-center coordinates.

F.3 Distribution function transform**check

In the above, we assume that X and x are related to each other by the guiding-center transformation (27) or (29) , i.e., x and X are not independent. For some cases, it may be convenient to treat x and X as independent variables and express the guiding-center transformation via an integral of the Dirac delta function. For example,

(353)

where x and X are considered as independent variables, ρ is the gyroradius evaluated at x, δ³(x − X − ρ) is the three-dimensional Dirac delta function. [In terms of general coordinates (x₁,x₂,x₃), the three-dimensional Dirac delta function is defined via the 1D Dirac delta function as follows:

(354)

where 𝒥 is the the Jacobian of the general coordinate system. The Jacobian is included in order to make δ³(x) satisfy the normalization condition ∫ δ³(x)dx = ∫ δ³(x)|𝒥|dx₁dx₂dx₃ = 1.]

Expression (353) can be considered as a transformation that transforms an arbitrary function from the guiding-center coordinates to the particle coordinates. Similarly,

(355)

is a transformation that transforms an arbitrary function from the the particle coordinates to the guiding-center coordinates.

F.4 Moments of distribution function expressed as integration over guiding-center variables

In terms of particle variables (x,v), it is straightforward to calculate the moment of the distribution function. For example, the number density n(x) is given by

(356)

However, it is a little difficult to calculate n(x) at real space location x by using the guiding-center variables (X,v). This is because holding x constant and changing v, which is required by the integration in Eq. (356), means the guiding-center variable X is changing according to Eq. (27). Using Eq. (34), expression (356) is written as

(357)

As is mentioned above, the dv integration in Eq. (357) should be performed by holding x constant and changing v, which means the guiding-center variable X = X(x,v) is changing. This means that, in (X,v) space, the above integration is a (generalized) curve integral along the the curve  X(v) = x −ρ(x,v) with x being constant. Treating X and x as independent variables and using the Dirac delta function δ, this curve integral can be written as the following double integration over the independent variables X and v:

(358)

Another perspective of interpreting Eq. (358) is that we are first using the transformation (353) to transform f_g to f_p and then integrating f_p in the velocity space.

G Diamagnetic flow **check**

The perturbed distribution function δF given in Eq. () contains two terms. The first term is gyro-phase dependent while the second term is gyro-phase independent. The perpendicular velocity moment of the second term will give rise to the so-called diamagnetic flow. For this case, it is crucial to distinguish between the distribution function in terms of the guiding-center variables, f_g(X,v), and that in terms of the particle variables, f_p(x,v). In terms of these denotations,  equation () is written as