The Eulerian variational formulation is presented to obtain governing equations of the electromagnetic turbulent gyrokinetic system. A local momentum balance in the system is derived from the invariance of the Lagrangian of the system under an arbitrary spatial coordinate transformation by extending the previous work [H. Sugama et al., Phys. Plasmas 28, 022312 (2021)]. Polarization and magnetization due to finite gyroradii and electromagnetic microturbulence are correctly described by the gyrokinetic Poisson equation and Ampère's law which are derived from the variational principle. Also shown is how the momentum balance is influenced by including collisions and external sources. Momentum transport due to collisions and turbulence is represented by a symmetric pressure tensor, which originates in a variational derivative of the Lagrangian with respect to the metric tensor. The relations of the axisymmetry and quasi-axisymmetry of the toroidal background magnetic field to a conservation form of the local momentum balance equation are clarified. In addition, an ensemble-averaged total momentum balance equation is shown to take the conservation form even in the background field with no symmetry when a constraint condition representing the macroscopic Ampère's law is imposed on the background field. Using the WKB representation, the ensemble-averaged pressure tensor due to the microturbulence is expressed in detail and it is verified to reproduce the toroidal momentum transport derived in previous works for axisymmetric systems. The local momentum balance equation and the pressure tensor obtained in this work present a useful reference for elaborate gyrokinetic simulation studies of momentum transport processes.

Gyrokinetics1–7 is a powerful theoretical framework based on which a large number of analytical and numerical studies on microinstabilities and turbulent processes in magnetized plasmas8 have been done. The original (or classical) gyrokinetic theory1–3,9 adopts the WKB approximation (or ballooning representation)10 and treats the perturbed parts of particle distribution functions and electromagnetic fields with gyroradius-scale perpendicular wavelengths. This type of gyrokinetic theory is widely employed as the basic model for local flux-tube gyrokinetic simulations11–15 to evaluate turbulent particle and heat fluxes. The other type of (or modern) gyrokinetic theory uses the Lie transformation method16 to obtain gyrocenter coordinates which obey the Lagrangian and/or Hamiltonian dynamics derived from the variational formulation.4,5 The modern theory guarantees favorable conservation properties17–22 of gyrokinetic equations for total distribution functions (including both background and fluctuation parts), which are generally used for long-time global gyrokinetic simulations.23–30 It is also noted here that classical gyrokinetic equations are shown to be consistently derived from modern ones by properly taking account of different phase-space coordinate systems used in the two type of theories.31 

Over the years, momentum transport processes have been attracting much attention because they determine profiles of plasma flows such as background plasma rotations and E × B zonal flows, which are regarded as important factors for stabilizing or regulating instabilities and improving plasma confinement.32 Also, there are much activities in designing advanced magnetic configurations such as toroidal systems with quasi-symmetry,33–36 in which reduction of neoclassical transport and increase in plasma flows are expected. In this paper, local momentum balance equations, which describe the momentum transport processes in electromagnetic gyrokinetic turbulence, are derived by extending the previous work37 on the momentum balance in electrostatic turbulence based on the Eulerian variational formulation,38 which is also called the Euler–Poincaré reduction procedure.39–44 

In conventional studies, momentum balance equations are obtained by taking the first-order velocity-space moment of a kinetic equation or from a variation in an action integral of the Lagrangian under infinitesimal translation or rotation. In the former derivation, it is unclear how the momentum balance in the direction perpendicular to the background magnetic field can be obtained from the gyrokinetic equation in the gyrophase-averaged form. In the latter, Noether's theorem can be applied to connect the symmetry condition of the system directly with the canonical momentum conservation equation,45,46 in which, however, local momentum transport is represented by the asymmetric canonical pressure tensor because of the vector potential included in the canonical momentum. In this work as well as in the previous works,37,44 the invariance of the Lagrangian under arbitrary infinitesimal transformations of general spatial coordinates is used to derive the local momentum balance equation which contains the symmetric pressure tensor obtained by taking the variational derivatives of the Lagrangian with respect to the 3 × 3 metric tensor components. This is analogous to the derivation of energy-momentum conservation laws from the invariance of an action integral under arbitrary transformations of spatiotemporal coordinates in the theory of general relativity.47 The relations of the symmetry and quasi-symmetry properties of the background magnetic field to the momentum balance equation are investigated with the help of the symmetric pressure tensor. In addition, the effects of collisions and/or external momentum sources can be easily included in the local momentum balance equation, by which both collisional48–50 and turbulent transport processes are described.

In extending the previous study for electrostatic turbulence37 to the case of electromagnetic turbulence, one need to consider the average and fluctuating parts of the magnetic field and accordingly those parts of the magnetic potential. Then, as shown in Ref. 31 and this paper, the variational derivative of the gyrokinetic Lagrangian with respect to the fluctuating part of the vector potential is used to correctly represent both the average and fluctuating parts of the local particle flux and the current density which appears in Ampère's law. On the other hand, the variational derivative of the gyrokinetic Lagrangian with respect to the average part of the vector potential also takes a form similar to the particle flux and it appears in the momentum balance equation derived using the variational technique in the present study. Comparison between the above-mentioned two types of particle fluxes shows that their average parts coincide with each other to the leading order in the small gyroradius expansion although their fluctuating parts do not. It is also shown in the present work that, using the ensemble average of the latter particle flux to self-consistently determine the average part of the magnetic field from Ampère's law, the conservation form of the ensemble-averaged local momentum balance equation can be derived.

Currently, large-scale gyrokinetic simulations such as global ones solving phenomena from a device size to an ion gyroradius-scale and cross-scale flux-tube simulations treating interactions between both ion and electron gyroradius scales are actively conducted. Huge simulations including all scales from the machine size to the electron gyroradius remain a challenging future task, for which global simulations need to treat full finite ion gyroradius effects at least as done in flux-tube simulations. In principle, the gyrokinetic model and the momentum balance equation presented in the present work contain all scales ranging from macroscopic equilibrium gradient lengths to microscopic turbulence wavelengths of the order of the electron gyroradius. The macroscopic behaviors of the momentum transport processes are described by the ensemble-averaged momentum balance equation, which is shown to take the conservation form under a condition to adjust the background field to the macroscopic Ampère's law. Furthermore, the WKB representation is used to explicitly express the full gyroradius effects of the electromagnetic turbulence on the symmetric pressure tensor, the ijth component of which represents the turbulent transport of the ith momentum component in the jth direction. These expressions can be applied to evaluation of the local momentum transport by the flux-tube simulations.

The rest of this paper is organized as follows. In Sec. II, equations of the gyrocenter motion in turbulent electromagnetic fields are derived as the Euler–Lagrange equations from the Lagrangian given as a function of the gyrocenter coordinates. In Sec. III, the Lagrangian for the whole gyrokinetic system consisting of particles of all species and electromagnetic fields is presented to derive gyrokinetic Vlasov equations for gyrocenter distribution functions and the gyrokinetic Poisson and Ampère equations for electrostatic and vector potentials based on the Eulerian (or Euler–Poincaré) variational formulation. Then, the gyrokinetic and field parts of the Lagrangian are all represented in terms of general spatial coordinates in Sec. IV and the invariance of the Lagrangian under an arbitrary infinitesimal transformation of spatial coordinates is used to derive the momentum balance equations for a single-particle-species system and for a system including all particle species and electromagnetic fields in Sec. V. In Sec. VI, axisymmetric, non-axisymmetric, and quasi-axisymmetric toroidal systems are investigated from the viewpoint of momentum balance, and Sec. VII presents the ensemble-averaged momentum balance equation, which is shown to take the conservation form when the background field is determined by the condition representing the macroscopic Ampère's law. The ensemble-averaged pressure tensor caused by the electromagnetic turbulence is expressed in detail using the WKB representation in Sec. VIII. Finally, conclusions are given in Sec. IX. In  Appendix A, the potential field included in the gyrocenter Hamiltonian is represented by gyroradius expansion around the gyrocenter, which is used in  Appendix B to expand the electromagnetic interaction part of Lagrangian density in terms of the electrostatic and vector potentials and their derivatives. In the same way as in  Appendix B, charge and current densities are expanded in  Appendixes C and  D, respectively, where the polarization and magnetization parts are identified. Energy balance equations in electromagnetic gyrokinetic turbulence are presented in  Appendix E.

The Lagrangian for describing the gyrocenter motion of the charged particle is given by4,5,31
L GYa ( Z , Z ̇ , t ) e a c A a * ( X , U , t ) · X ̇ + m a c e a μ ϑ ̇ H GYa ( Z , t ) ,
(1)
where the modified vector potential A a * is defined by A a * ( X , U , t ) A ( X , t ) + ( m a c / e a ) U b ( X , t ), the subscript a represents the particle species with mass ma and charge ea, and ̇ d / d t represents the time derivative along the motion of the particle in phase space. The gyrocenter phase-space coordinates X, U, μ m v 2 / ( 2 B ), and ϑ denote the gyrocenter position, the velocity component parallel to the magnetic field, the magnetic moment, and the gyrophase angle, respectively. The vector potential and the unit vector parallel to the background magnetic field B are written by A and b B / B, respectively. Here, it is supposed that A can weakly depend on time t and accordingly the background magnetic field B = × A is allowed to slowly vary in time. Thus, we can treat the inductive electric field E T c 1 A / t, which drives the Ohmic current in tokamaks.
The gyrocenter Hamiltonian which appears on the right-hand side of Eq. (1) is given by
H GYa ( Z , t ) 1 2 m a U 2 + μ B ( X , t ) + e a Ψ a ( Z , t ) ,
(2)
and the potential Ψ a including effects of the turbulent electromagnetic fields is defined by31,
e a Ψ a ( Z , t ) e a ψ a ( Z , t ) ϑ e a c v B a · A ̂ ( X + ρ a , t ) ϑ + e a 2 2 m a c 2 | A ̂ ( X + ρ a , t ) | 2 ϑ e a 2 2 B μ ( ψ ̃ a ) 2 ϑ ,
(3)
where ρ a is the gyroradius vector given by ρ a b × v / Ω a, Ω a e a B / ( m a c ) is the gyrofrequency, and
ψ a ( Z , t ) ϕ ( X + ρ a , t ) v c · A ̂ ( X + ρ a , t ) .
(4)
The gyrophase average Q ϑ and the gyrophase-dependent part Q ̃ of an arbitrary function Q of the gyrocenter phase-space coordinates Z ( X , U , μ , ϑ ) are represented by
Q ϑ 1 2 π Q d ϑ and Q ̃ Q Q ϑ ,
(5)
respectively. The particle's velocity v is written as
v U b ( X , t ) W [ sin ϑ e 1 ( X , t ) + cos ϑ e 2 ( X , t ) ] ,
(6)
where ( e 1 , e 2 , b ) are unit vectors which form a right-handed orthogonal triad and are regarded as functions of ( X , t ). The magnetic moment is given by μ m a W 2 / 2 B, and
v B a c e a B b × ( m a U 2 b · b + μ B )
(7)
is the first-order drift velocity consisting of curvature drift and B drift. Second-order terms retained in Eq. (3) are necessary for correctly deriving the gyrokinetic Poisson and Ampère equations from variational derivatives with ϕ and A ̂, respectively, as shown in Sec. III. Especially, it is shown in Ref. 31 that the second-order term ( e a / c ) v B a · A ̂ ϑ, which is often neglected in conventional studies, should be kept for the variational derivative to obtain the gyrokinetic Ampère's law including both equilibrium and turbulent parts accurately. In  Appendix A, Ψ a is expanded in the gyroradius and decomposed into several parts which have different dependences on electrostatic and magnetic fluctuations. It is noted in Refs. 51–53 that basic equations including terms of higher order, which are not considered here, are required for accurately describing the flux-surface-averaged momentum balance along the symmetry direction in up-down symmetric tokamaks and stellarator-symmetric quasisymmetric stellarators where the low-flow ordering is assumed as in the present work.
Using the gyrocenter Lagrangian in Eq. (1), the Euler–Lagrangian equations are given by
d d t ( L GYa Z ̇ ) L GYa Z = 0 ,
(8)
from which the gyrocenter motion equations are obtained as
d Z d t = { Z , H GYa } + { Z , X } · e a c A a * t ,
(9)
with the Poisson brackets defined by
{ X , X } = c e a B a * b × I , { X , U } = B a * m a B a * , { X , ϑ } = 0 , { U , ϑ } = 0 , { ϑ , μ } = e a m a c .
(10)
Equation (9) is rewritten as
d X d t = 1 B a * [ ( U + e a m a Ψ a U ) B a * + c b × ( μ e a B + Ψ a + 1 c A a * t ) ] ,
(11)
d U d t = B a * m a B a * · ( μ B + e a Ψ a + e a c A a * t ) ,
(12)
d μ d t = 0 ,
(13)
and
d ϑ d t = Ω a + e a 2 m a c Ψ a μ ,
(14)
where B a * and B a * are defined by
B a * × A a * and B a * B a * · b ,
(15)
respectively. The gyrocenter motion given by Eqs. (11)–(14) satisfies Liouville's theorem, which is expressed as
D a ( Z , t ) t + Z · ( D a ( Z , t ) d Z d t ) = 0 ,
(16)
where the Jacobian D a ( Z , t ) is given by D a ( Z , t ) = B a * / m a.
The action integral for the gyrokinetic Vlasov–Poisson–Ampère system is given by
I t 1 t 2 d t L GKF t 1 t 2 d t ( L G K + L F ) ,
(17)
where the gyrokinetic Lagrangian LGK is defined by the phase-space integral of the gyrocenter distribution function Fa multiplied by the gyrocenter Lagrangian LGYa [see Eq. (1)] as
L G K a L GKa a d 6 Z F a ( Z , t ) L GYa ( Z , u a Z ( Z , t ) , t ) .
(18)
Here, based on the Eulerian picture, the temporal change rates of the gyrocenter coordinates are regarded as functions of ( Z , t ) and they are represented by
u a Z ( Z , t ) = [ u a X ( Z , t ) , u a U ( Z , t ) , u a μ ( Z , t ) , u a ϑ ( Z , t ) ] ,
(19)
which are used in Eq. (18) to evaluate L GYa ( Z , u a Z ( Z , t ) , t ). Then, the distribution function Fa satisfies
F a ( Z , t ) t + Z · [ F a ( Z , t ) u a Z ( Z , t ) ] = 0 ,
(20)
where the functional form of u a Z ( Z , t ) is determined later by the Euler–Poincaré variational principle.
Here, the Darwin approximation is made to remove electromagnetic waves propagating at light speed, and the Lagrangian LF is defined by5 
L F 1 8 π V d 3 x [ | E L ( x , t ) | 2 | B ( x , t ) + B ̂ ( x , t ) | 2 + 2 c λ ( x , t ) · A ̂ ( x , t ) ] ,
(21)
where V represents the spatial domain of the system, λ plays the role of a Lagrange undetermined multiplier to derive the Coulomb gauge condition
· A ̂ = 0
(22)
from the variational condition δ I / δ λ = 0 (or δ L GKF / δ λ = δ L F / δ λ = 0), and E L is the longitudinal (or irrotational) part of the electric field written in terms of electrostatic potential ϕ as
E L ϕ .
(23)
Now, the trajectories of particles in the phase space as well as the electrostatic potential and perturbed vector potential are virtually varied to derive the governing equations of the collisionless electromagnetic gyrokinetic turbulent system from the Eulerian variational principle. Following the same procedure as in Refs. 31 and 44 the infinitesimal variation in the phase-space trajectory is represented in the Eulerian picture by
δ Z a E ( Z , t ) = [ δ X a E ( Z , t ) , δ U a E ( Z , t ) , δ μ a E ( Z , t ) , δ ζ a E ( Z , t ) ] .
(24)
Then, the variations in the functional forms of u a Z = ( u a X , u a U , u a μ , u a ϑ ) and F a ( Z , t ) are written in the Eulerian picture as
δ u a Z ( Z , t ) = ( t + u a Z ( Z , t ) · Z ) δ Z a E ( Z , t ) δ Z a E ( Z , t ) · Z u a Z ( Z , t )
(25)
and
δ F a ( Z , t ) = Z · [ F a ( Z , t ) δ Z a E ( Z , t ) ] ,
(26)
respectively. The variations in the electrostatic potential and the perturbed vector potential are denoted by δ ϕ and δ A ̂, respectively. Using Eqs. (17), (18), (25), and (26), the variation in the action integral IGKF is given by
δ I GKF = a t 1 t 2 d t d 6 Z F a × [ ( L GYa Z ) u ( d d t ) a ( L GYa u a Z ) ] · δ Z a E + t 1 t 2 d t V d 3 x [ δ L GKF δ ϕ δ ϕ + δ L GKF δ A ̂ · δ A ̂ ] + B . T .
(27)
Here, B.T. represents boundary terms that appear due to partial integrals and ( d / d t ) a denotes the time derivative along the phase-space trajectory defined by
( d d t ) a t + u a Z · Z ,
(28)
from which one obtains ( d / d t ) a Z = u a Z. The variational derivative δ F [ f ] / δ f of any functional F [ f ] of a function f in three-dimensional space is defined as a function in the space that satisfies
d 3 x δ F [ f ] δ f ( x ) φ ( x ) = lim ε 0 F [ f + ε φ ] F [ f ] ε ,
(29)
from which one can also write
δ F [ f ( y ) ] δ f ( x ) d d ε F [ f ( y ) + ε δ 3 ( y x ) ] | ε = 0 ,
(30)
where δ 3 ( y x ) δ ( y 1 x 1 ) δ ( y 2 x 2 ) δ ( y 3 x 3 ). From Eq. (30), one has
δ ϕ ( X + ρ a ) δ ϕ ( x ) = δ ( X + ρ a x ) .
(31)
Now, it is required that δ I GKF = 0 holds for any variations δ Z a E, δ ϕ, and δ A ̂, which vanish on the boundaries of the integral region. Then, it is found from Eq. (27) that
( d d t ) a ( L GYa u a Z ) ( L GYa Z ) u = 0 ,
(32)
δ L GKF / δ ϕ = 0, and δ L GKF / δ A ̂ = 0 need to be satisfied. Here, since Eq. (32) is equivalent to Eq. (8), one finds that u a Z ( Z , t ) should be given by the right-hand side of Eq. (9). Thus, the gyrokinetic Vlasov equation is given by Eq. (20) with Eq. (9) substituted into u a Z ( Z , t ).
The gyrokinetic Poisson equation is derived from δ L GKF / δ ϕ = 0 and written as
· E L = 4 π ρ c ,
(33)
where the charge density ρc is given by
ρ c = δ L G K δ ϕ = a δ L GKa δ ϕ = a e a d 6 Z δ 3 ( X + ρ a x ) ( F a + e a ψ ̃ a B F a μ ) ,
(34)
and δ L F / δ ϕ = ( 1 / 4 π ) · E L is used.
The gyrokinetic Ampère's law is derived from δ L GKF / δ A ̂ = 0 as
× ( B + B ̂ ) = 4 π c j 1 c λ ,
(35)
where the electric current density j is given by
j = c δ L G K δ A ̂ = c a δ L GKa δ A ̂ = a e a d 6 Z δ 3 ( X + ρ a x ) × [ F a ( Z , t ) ( v e a m a c A ̂ + v B a ) + e a ψ ̃ a B F a μ v ] ,
(36)
and δ L F / δ A ̂ = ( 1 / 4 π ) × ( B + B ̂ ) ( 1 / 4 π c ) λ is used. It is noted here that an arbitrary vector field a is written as a = a L + a T where the longitudinal (or irrotational) and transverse (or solenoidal) parts of a are given by a L ( x ) = ( 4 π ) 1 d 3 x ( · a ( x ) ) / | x x | and a T ( x ) = ( 4 π ) 1 × ( × d 3 x a ( x ) / | x x | ), respectively.54 Then, the longitudinal part of Eq. (35) gives
λ = 4 π j L .
(37)
From the transverse part of Eq. (35), the gyrokinetic Ampère's law is written as
× ( B + B ̂ ) = 4 π c j T .
(38)
In Eqs. (37) and (38), j L and j T represent the longitudinal and transverse parts of j, respectively.

In this section, general spatial coordinates are used to represent the Lagrangian of the electromagnetic turbulent gyrokinetic system defined in Sec. III. The Lagrangian is given as the integral of the Lagrangian density with respect to the general spatial coordinates, and it is invariant under an arbitrary spatial coordinate transformation.

The action integral IGKF in Eq. (17) is written here as
I GKF t 1 t 2 d t L GKF t 1 t 2 d t V d 3 x L GKF ,
(39)
where the Lagrangian density L GKF is given by
L GKF L G K + L F L G K a d 3 v F a ( x , v , t ) L GYa ( x , v , t ) L F g ( x ) 8 π [ g i j ( x ) ( E L ) i ( x , t ) ( E L ) j ( x , t ) g i j ( x ) { B i ( x , t ) + B ̂ i ( x , t ) } { B j ( x , t ) + B ̂ j ( x , t ) } + 2 c λ ( x , t ) g i j ( x ) i A ̂ j ( x , t ) ] .
(40)
Here, x ( x i ) i = 1 , 2 , 3 , v ( U , μ , ϑ ), d 3 x d x 1 d x 2 d x 3, and d 3 v dUd μ d ϑ are used, and j represents a covariant derivative. In the equation for L G K shown in Eq. (40), x ( x i ) i = 1 , 2 , 3 represents the coordinates not of the position of the particle but that of the gyrocenter (denoted by X in Sec. II). It should be emphasized that in this section, x ( x i ) i = 1 , 2 , 3 are general spatial coordinates which can be either Cartesian or any other curved coordinates. However, the spatial position vector r = r ( x ) is assumed to be a function of only the spatial coordinates x ( x i ) i = 1 , 2 , 3 and it is independent of time t. The gyrocenter distribution function in the (x, v)-space is denoted by Fa, and the number of particles of species a in the phase-space volume element d 3 x d 3 v d x 1 d x 2 d x 3 dUd μ d ϑ at time t is given by F a ( x , v , t ) d 3 x d 3 v. This paper employs the summation convention that the same symbol used for a pair of upper and lower indices within a term [such as seen in Eq. (40) as well as in the equations shown below] indicates summation over the range { 1 , 2 , 3 } of the symbol index. The contravariant metric tensor components gij in the general spatial coordinates x ( x i ) are related to the covariant components gij by g i k g k j = δ j i, where δ j i represents the Kronecker delta. The determinant of the covariant metric tensor matrix is denoted by g ( x ) det [ g i j ( x ) ]. As the spatial position vector r is a function of only the spatial coordinates x ( x i ), g i j ( x ) , g i j ( x ), and g(x) are all independent of time t.
The gyrocenter Lagrangian LGYa, which enters L G K in Eq. (40), is represented in the Eulerian picture by
L GYa ( e a c A j ( x , t ) + m a U b i ( x , t ) g i j ( x ) ) u a x j ( x , v , t ) + m a c e a μ u a ϑ ( x , v , t ) H GYa ( x , U , μ , t ) ,
(41)
where b i B i / B is the ith contravariant component of the unit vector parallel to the background magnetic field and the background field strength is given by B ( x , t ) g i j ( x ) B i ( x , t ) B j ( x , t ). The contravariant components of the background and perturbed magnetic fields are expressed in terms of the covariant components of the vector potentials as
B i ( x , t ) ε ijk g ( x ) A k ( x , t ) x j , B ̂ i ( x , t ) ε ijk g ( x ) A ̂ k ( x , t ) x j ,
(42)
where the Levi-Civita symbol is denoted by
ε ijk ε ijk { 1 ( ( i , j , k ) = ( 1 , 2 , 3 ) , ( 2 , 3 , 1 ) , ( 3 , 1 , 2 ) ) 1 ( ( i , j , k ) = ( 1 , 3 , 2 ) , ( 2 , 1 , 3 ) , ( 3 , 2 , 1 ) ) 0 ( otherwise ) .
(43)
The gyrocenter Hamiltonian is written here as
H GYa ( x , U , μ , t ) 1 2 m a U 2 + μ B ( x , t ) + e a Ψ a ( x , U , μ , t ) ,
(44)
and the fluctuation part is given by
Ψ a ( x , μ , t ) ϕ ( x , t ) + Ψ E 1 a ( x , μ , t ) + Ψ A ̂ 1 a ( x , U , μ , t ) + Ψ E 2 a ( x , μ , t ) + Ψ E A ̂ a ( x , U , μ , t ) + Ψ A ̂ 2 a ( x , U , μ , t ) ,
(45)
where
Ψ E 1 a ( x , μ , t ) = n = 1 α a j 1 j n ( x , μ , t ) n ! j 1 j n ϕ ( x , t ) = n = 1 α a j 1 j n ( x , μ , t ) n ! j 1 j n 1 ( E L ) j n ( x , t )
(46)
and
Ψ A ̂ 1 a ( x , U , μ , t ) = 1 c n = 0 1 n ! [ α a j 1 j n ( U b i + v B a i ) + Ω a × g ε klm α a j 1 j n k b l g i m ] j 1 j n A ̂ i ( x , t ) .
(47)
Here, α a j 1 j n ( x , μ , t ) is defined by Eqs. (A9)–(A11) in  Appendix A by using h i j g i j b i b j. The second-order parts Ψ E 2 a ( x , μ , t ), Ψ E A ̂ a ( x , U , μ , t ), and Ψ A ̂ 2 a ( x , U , μ , t ) on the right-hand side of Eq. (45) are obtained from Eqs. (A15), (A17), and (A18), respectively, with the partial derivative j replaced by the covariant derivative j.
The temporal change rates of the gyrocenter coordinates in Eq. (19) are denoted here by u a x i ( x , v , t ), u a U ( x , v , t ), u a μ ( x , v , t ), and u a ϑ ( x , v , t ). Then, Eq. (20) for the distribution function F a ( x , v , t ) is written as
F a t + x j ( F a u a x j ) + U ( F a u a U ) + μ ( F a u a μ ) + ϑ ( F a u a ϑ ) = 0.
(48)
The gyrocenter Hamiltonian HGYa given in Eq. (44) takes a functional form
H GYa = H GYa [ v , j A i ( x , t ) , j k A i ( x , t ) , { J ϕ ( x , t ) } , { J A ̂ i ( x , t ) } , { J g i j ( x ) } ] ,
(49)
which depends on the velocity space coordinates (except for ϑ) as well as the general spatial coordinates x ( x i ) i = 1 , 2 , 3 through the field variables [ j A i ( x , t ) , j k A i ( x , t ) , { J ϕ ( x , t ) } , { J A ̂ i ( x , t ) } , { J g i j ( x ) } ]. Here, the notation J ( j 1 , j 2 , , j n ) ( n = 0 , 1 , 2 , ; j 1 , j 2 , , j n = 1 , 2 , 3) is used to write
J F { F ( n = 0 ) j 1 j 2 j n F n F / x j 1 x j 2 x j n ( n 1 ) ,
(50)
where F is an arbitrary function of x = ( x i ) i = 1 , 2 , 3. Then, { J ϕ } { ϕ , j ϕ , j k ϕ , jkl ϕ , }, and the definitions of other compact notations { J A ̂ i }, and { J g i j ( x ) } in Eq. (49) are understood in the same way. One obtains ( E L ) i i ϕ from Eq. (23) which is used to replace { J ϕ } with { ϕ , { J ( E L ) i } } where { J ( E L ) i } { ( E L ) i , j ( E L ) i , j k ( E L ) i , jkl ( E L ) i , }. Note that high-order spatial derivative terms due to finite gyroradii enter the gyrocenter Hamiltonian HGYa as seen in Eqs. (46) and (47) where the covariant derivatives contain the spatial derivatives of gij through the Christoffel symbols [see Eq. (A4) in Ref. 37].
It is found from Eqs. (41) and (49) that the functional form of the gyrocenter Lagrangian LGYa is written as
L GYa = L GYa [ v , u a x i ( x , v , t ) , u a ϑ ( x , v , t ) , A i ( x , t ) , j A i ( x , t ) , j k A i ( x , t ) , { J ϕ ( x , t ) } , { J A ̂ i ( x , t ) } , { J g i j ( x ) } ] ,
(51)
where u a x i ( x , v , t ) , u a ϑ ( x , v , t ) , ϕ ( x , t ), and A ̂ i ( x , t ) are the functions, the governing equations of which are derived from the variation principle in Sec. IV C while the dependence of LGYa on A i ( x , t ) , j A i ( x , t ), j k A i ( x , t ), and J g i j ( x , t ) is also explicitly shown because their variations need to be taken into account to evaluate the variation of LGYa in Sec. IV where the local momentum balance is derived using the general spatial coordinate transformation which causes the variations in the functional forms of both ( u a x i , u a ϑ , ϕ , A ̂ i ) and ( A i , g i j ).
It is found from Eqs. (40), (41), and (44) that the part of the Lagrangian density including Ψ a is given by
L Ψ a L Ψ a a d 3 v F a e a Ψ a = ρ c ( g ) ϕ + L E 1 + L A ̂ 1 + L E 2 + L E A ̂ + L A ̂ 2 ,
(52)
which determines the electromagnetic interaction of particles. Here, the gyrocenter charge density ρ c ( g ) is defined by
ρ c ( g ) a e a N a ( g ) a e a d 3 v F a ,
(53)
where N a ( g ) represents the gyrocenter density. The terms L E 1 , L A ̂ 1, L E 2 , L E A ̂, and L A ̂ 2 on the right-hand side of Eq. (52), are defined by
L E 1 a L E 1 a = k = 1 Q 0 j 1 j 2 k j 1 j 2 k 1 ( E L ) j 2 k ,
(54)
L A ̂ 1 a L A ̂ 1 a = n = 1 R 0 j 1 j n j 1 j n 1 A ̂ j n ,
(55)
L E 2 a L E 2 a = 1 2 m = 1 n = 1 χ E i 1 i m ; j 1 j n i 1 i m 1 ( E L ) i m × j 1 j n 1 ( E L ) j n ,
(56)
L E A ̂ a L E A ̂ a = m = 1 n = 1 χ E A ̂ j 1 j m ; k 1 , k n j 1 j m 1 ( E L ) j m × k 1 k n 1 A ̂ k n ,
(57)
and
L A ̂ 2 a L A ̂ 2 a = 1 2 m = 1 n = 1 χ A ̂ j 1 j m ; k 1 , k n j 1 j m 1 A ̂ j m × k 1 k n 1 A ̂ k n ,
(58)
respectively. Here, Q 0 j 1 j 2 k , R 0 j 1 j n, χ E i 1 i m ; j 1 j n, χ E A ̂ j 1 j m ; k 1 , k n, and L A ̂ 2 included in Eqs. (54)–(58) are given by
[ Q 0 j 1 j 2 k , R 0 j 1 j n ] a [ Q 0 a j 1 j 2 k , R 0 a j 1 j n ] , [ χ E i 1 i m ; j 1 j n , χ E A ̂ j 1 j m ; k 1 , k n , χ A ̂ j 1 j m ; k 1 , , k n ] a [ χ E a i 1 i m ; j 1 j n , χ E A ̂ a j 1 j m ; k 1 , , k n , χ A ̂ a j 1 j m ; k 1 , , k n ] ,
(59)
where Q 0 a j 1 j 2 k , R 0 a j 1 j n, χ E a i 1 i m ; j 1 j n, χ E A ̂ a j 1 j m ; k 1 , , k n, and χ A ̂ a j 1 j m ; k 1 , , k n are defined by Eqs. (B4), (B6), and (B12) in  Appendix B. As seen in  Appendixes C and  D, the charge and current densities in the gyrokinetic Poisson and Ampère equations are derived from L ψ, of which the components L E 1 , L A ̂ 1, L E 2 , L E A ̂, and L A ̂ 2 are associated with the polarization charge and the magnetization current.
Here, general spatial coordinates x = ( x i ) i = 1 , 2 , 3 are used for the Eulerian variational derivation of the gyrokinetic Vlasov–Poisson–Ampère equations for the electromagnetic gyrokinetic system. The virtual variations in the phase-space trajectory are now represented in the Eulerian picture by δ x a E i ( x , v , t ), δ U a E ( x , v , t ), δ μ a E ( x , v , t ), and δ ϑ a E ( x , v , t ). Then, Eqs. (25), (26), and (27) are rewritten as
[ δ u a x i , δ u a U , δ u a μ , δ u a ϑ ] = ( t + u a x j x j + u a U U + u a μ μ + u a ϑ ϑ ) [ δ x a E i , δ U a E , δ μ a E , δ ϑ a E ] ( δ x a E j x j + δ U a E U + δ μ a E μ + δ ϑ a E ϑ ) [ u a x i , u a U , u a μ , u a ϑ ] ,
(60)
δ F a = x j ( F a δ x a E j ) U a ( F a δ U a E ) μ ( F a δ μ a E ) ϑ ( F a δ ϑ a E ) ,
(61)
and
δ I GKF = a t 1 t 2 d t V d 3 x d 3 v F a × [ { ( L GYa x i ) u ( d d t ) a ( L GYa u a x i ) } δ x a E i + ( L GYa U ) u δ U a E + ( L GYa μ ) u δ μ a E + { ( L GYa ϑ ) u ( d d t ) a ( L GYa u a ϑ ) } δ ϑ a E ] + t 1 t 2 d t V d 3 x ( δ L G K δ ϕ δ ϕ + δ L G K δ A ̂ i δ A ̂ i + δ L G K δ λ δ λ ) + B . T . ,
(62)
respectively, where ( L GYa / x i ) u, ( L GYa / U ) u, ( L GYa / μ ) u, and ( L GYa / ϑ ) u denote the derivatives of LGYa in x i, U, μ, and ϑ, respectively, with ( u a x i , u a ϑ ) kept fixed in LGYa, and the time derivative along the phase-space trajectory is represented by
( d d t ) a t + u a x k x k + u a U U + u a μ μ + u a ϑ ϑ .
(63)
Using Eq. (27) and δ I GKF = 0, one first obtains
( d d t ) a p a i = ( L GYa x i ) u ,
(64)
where p a i L GYa / u a x i = ( e a / c ) A i ( x , t ) + m a U b i ( x , t ) ( e a / c ) A a i * ( x , U , t ) represents the covariant vector component of the canonical momentum. Equation (64) can be deformed to obtain
m a u a U b i = e a ( Ψ a x i 1 c A a i * t + 1 c g ε ijk u x j B * k ) μ B x i ,
(65)
where B a * i ( ε ijk / g ) ( A a k * / x j ). Next, combining Eq. (65) with ( L GYa / U ) u = m a ( u a x i b i U ) = 0, ( L GYa / μ ) u = ( m a c / e a ) u a ϑ B e a Ψ a / μ = 0, and ( d / d t ) a ( L GYa / u a ϑ ) = ( m a c / e a ) u a μ = ( L GYa / ϑ ) u = 0, the gyrocenter motion equations are derived as
u a x i = 1 B a * [ ( U + e a m a Ψ a U ) B a * i + c ε ijk g b j ( μ e a B x k + Ψ a x k + 1 c A a k * t ) ] ,
(66)
m a u a U = B a * i B a * [ μ B x i + e a ( Ψ a x i + 1 c A a i * t ) ] ,
(67)
u a μ = 0 ,
(68)
and
u a ϑ = Ω a + e a 2 m a c Ψ a μ ,
(69)
where B a * B a * i b i. Substituting Eqs. (66)–(69) into Eqs. (48) and taking its average with respect to the gyrophase ϑ, the gyrokinetic Vlasov equation is derived as
F ¯ a t + x i [ F ¯ a 1 B a * { ( U + e a m a Ψ a U ) B a * i + c ε ijk g b j × ( μ e a B x k + Ψ a x k + 1 c A a k * t ) } ] + U [ F ¯ a B a * i m a B a * { e a ( Ψ a x i + 1 c A a i * t ) μ B x i } ] = 0 ,
(70)
where F ¯ a F a ϑ F a d ϑ / 2 π is the gyrophase-averaged distribution function.
The remaining conditions for δ I GKF = 0 are given by δ L G K / δ ϕ = δ L G K / δ A ̂ i = δ L G K / δ λ = 0. The Coulomb gauge condition is obtained as δ L GKF / δ λ = ( 2 / c ) i A ̂ i = 0. The gyrokinetic Poisson equation is given by
δ L GKF δ ϕ = ρ c + 1 4 π x i [ g ( E L ) i ] = 0 ,
(71)
where the charge density ρc is written as
ρ c = δ L G K δ ϕ = ρ ( g c ) i P G i ,
(72)
with the generalized polarization vector density P G i defined by
P G i n = 0 ( 1 ) n i 1 i n Q i i 1 i n .
(73)
Here, the multipole moments Q i 1 i m are given by Eq. (C7) in  Appendix C. The gyrokinetic Ampère's law is derived from the condition δ L G K / δ A ̂ i = 0 which is written as
δ L GKF δ A ̂ i = 1 c j i 1 4 π ε ijk ( B k + B ̂ k ) x j g 4 π c g i j λ x j = 0.
(74)
Here, the current density is written as
j l = ( j 0 ) l + c k N k l ,
(75)
where ( j 0 ) l is defined by Eq. (D9) in  Appendix D and Nkl is given by
N k l n = 0 ( 1 ) n + 1 j 1 j n R j 1 j n k l ,
(76)
with R j 1 j n k defined in Eq. (D5).

In this section, the invariance of the Lagrangian under arbitrary infinitesimal transformations of spatial coordinates is used to derive the local momentum balance equations for the single-particle-species system and for the whole system including particles of multiple species and electromagnetic fields.

An arbitrary infinitesimal transformation of spatial coordinates from x = ( x i ) i = 1 , 2 , 3 to x = ( x i ) i = 1 , 2 , 3, is written as
x i = x i + ξ i ( x ) ,
(77)
where the infinitesimal variation in the spatial coordinate xi is denoted by ξ i ( x ) which is an arbitrary function of only the spatial coordinates x = ( x i ) i = 1 , 2 , 3 and independent of time t.
The gyrokinetic Lagrangian LGKa is written as
L GKa V d 3 x L GKa V d 3 x d 3 v F a L GYa ,
(78)
where Fa and LGYa defined in Eq. (1) behave as a scalar density field and a scalar field, respectively, under the transformation of the spatial coordinates. The variation in LGKa under the infinitesimal spatial coordinate transformation in Eq. (77) is written as
δ ¯ L GKa V d 3 x ( ( ξ i L GKa ) x i + δ ¯ L GKa ) .
(79)
Here and hereafter, the notation δ ¯ represents the variation caused by the infinitesimal coordinate transformation in Eq. (77) and it should be distinguished from the variation δ due to the virtual displacement used in Secs. III and IV C. The expression of the integral in Eq. (79) takes the form often found in conventional textbooks (see, e.g., Ref. 55) to give the change in the integral caused by the infinitesimal transformation. In the integrand in Eq. (79), the divergence term ( ξ i L GKa ) / x i is obtained using Gauss's theorem for the difference between the domains of integrations in x = ( x i ) i = 1 , 2 , 3 and x = ( x i ) i = 1 , 2 , 3 while δ ¯ L GKa is written using the Leibniz rule for the derivative operation by δ ¯ as
δ ¯ L GKa = d 3 v δ ¯ ( F a L GYa ) = d 3 v ( δ ¯ F a · L GYa + F a · δ ¯ L GYa ) ,
(80)
where δ ¯ F a and δ ¯ L GYa represent the variations in the spatial functional forms of Fa and LGYa under the infinitesimal spatial coordinate transformation.
Then, applying the chain rule formula for the derivative operation δ ¯ L GYa [ u a x i , u a ϑ , { J A i } , { J A ̂ i } , { J ϕ } , { J g i j } ] yields
δ ¯ L GYa = L GYa u a x i δ ¯ u a x i + L GYa u a ϑ δ ¯ u a ϑ + J L GYa ( J A i ) δ ¯ ( J A i ) + J L GYa ( J A ̂ i ) δ ¯ ( J A ̂ i ) + J L GYa ( J ϕ ) δ ¯ ( J ϕ ) + J L GYa ( J g i j ) δ ¯ ( J g i j ) ,
(81)
where L GYa / ( J A i ) = 0 when the order n of J = ( j 1 , , j n ) is greater than or equal to three [see Eq. (51)].
As shown in Ref. 37, the variation in the functional form under the infinitesimal spatial coordinate transformation can be represented by δ ¯ = L ξ, where L ξ is the Lie derivative56 associated with the vector field given by ( ξ i ) i = 1 , 2 , 3 and it acts on an arbitrary tensor field (as well as an arbitrary tensor field density). Using the fact that L ξ L GKa = ( ξ i L GKa ) / x i holds from the definition of the Lie derivative acting on a scalar density field, one finds that the integrand in Eq. (79) is written as L ξ L GKa + δ ¯ L GKa which is found to vanish from δ ¯ = L ξ. Then, the integral in Eq. (79) also vanishes and accordingly δ ¯ L GKa = 0, which means that LGKa is a scalar constant which is invariant under the coordinate transformation. Using Eqs. (79)–(81) and δ ¯ = L ξ, δ ¯ L GKa = 0 can also be written as
δ ¯ L GKa = V d 3 x d 3 v F a ( ξ i L GYa x i + δ ¯ L GYa ) = V d 3 x d 3 v F a ( ξ i L GYa x i + L GYa u a x i δ ¯ u a x i + L GYa u a ϑ δ ¯ u a ϑ + J L GYa ( J A i ) δ ¯ ( J A i ) + J L GYa ( J A ̂ i ) δ ¯ ( J A ̂ i ) + J L GYa ( J ϕ ) δ ¯ ( J ϕ ) + J L GYa ( J g i j ) δ ¯ ( J g i j ) ) = 0.
(82)

Recall that L ξ annihilate any scalar constant. Therefore, when LGKa is a scalar constant, one can naturally write δ ¯ L GKa = L ξ L GKa ( = 0 ). Thus, one can confirm that the relation δ ¯ = L ξ is consistent and useful when treating all tensor variables (including scalar fields and scalar constants) and deriving the invariance properties associated with the coordinate transformation. This relation δ ¯ = L ξ can be applied under the condition that all quantities, in which the variations due to the coordinate transformation are considered, can be written in terms of tensor fields (or tensor field densities) on which the operation of the Lie derivative can be defined. This condition means that integrals using such tensor fields yield scalar constants which represent geometric quantities and take invariant values independent of the choice of the spatial coordinates. It should be stressed that δ ¯ L GKa = L ξ L GKa = 0 and Eq. (82), which are derived only from the above-mentioned invariance property under the spatial coordinate transformation, are valid whether the gyrokinetic equation derived from the variational principle associated with the virtual variation in the phase-space trajectory in Secs. II and III holds or not.

In the same way as in Eq. (79), the invariance of the Lagrangian LGKF of the whole system under the infinitesimal spatial coordinate transformation can be written as
δ ¯ L GKF V d 3 x ( ( ξ i L GKF ) x i + δ ¯ L GKF ) = a δ ¯ L GKa + V d 3 x ( ( ξ i L F ) x i + δ ¯ L F ) = 0 ,
(83)
where LGKF is defined by Eq. (39) with Eq. (40). The variation δ ¯ L F of the field Lagrangian density L F defined in Eq. (40) can be written as
δ ¯ L F = ( ξ i L F ) x i = L F ( j A i ) δ ¯ ( j A i ) + J L F ( J A ̂ i ) δ ¯ ( J A ̂ i ) + L F ( i ϕ ) δ ¯ ( i ϕ ) + J L F ( J g i j ) δ ¯ ( J g i j ) + L F λ δ ¯ λ ,
(84)
where δ ¯ = L ξ, L ξ L F = i ( ξ i L F ), and the chain rule for δ ¯ L F [ { j A i } , { J A ̂ i } , { i ϕ } , { J g i j } , λ ] are used. We now use Eqs. (82), and (84) to rewrite Eq. (83) as
δ ¯ L GKF = V d 3 x [ a d 3 v F a ( ξ i L GYa x i + L GYa u a x i δ ¯ u a x i + L GYa u a ϑ δ ¯ u a ϑ + J L GYa ( J A i ) δ ¯ ( J A i ) + J L GYa ( J A ̂ i ) δ ¯ ( J A ̂ i ) + J L GYa ( J ϕ ) δ ¯ ( J ϕ ) + J L GYa ( J g i j ) δ ¯ ( J g i j ) ) + ( ξ i L F ) x i + L F ( j A i ) δ ¯ ( j A i ) + J L F ( J A ̂ i ) δ ¯ ( J A ̂ i ) + L F ( i ϕ ) δ ¯ ( i ϕ ) + J L F ( J g i j ) δ ¯ ( J g i j ) + L F λ δ ¯ λ ] = 0.
(85)

As seen from Eqs. (81), (82), (84), and (85), the invariance of the scalar constants LGKa and LGKF under the infinitesimal spatial transformation can be confirmed using the chain rule formulas for the derivative operation δ ¯ = L ξ acting on the scalar field L GYa [ u a x i , u a ϑ , { J A i } , { J A ̂ i } , { J ϕ } , { J g i j } ] and the scalar field density L F [ { j A i } , { J A ̂ i } , { i ϕ } , { J g i j } , λ ] which are given as composite functions. In Secs. V B and V C, Eqs. (82) and (85) are used to derive the local momentum balance equation for the single-particle-species system and that for the whole system consisting of particles of all species and electromagnetic fields.

We now use δ ¯ = L ξ, L ξ u a x j = ξ i i u a x j u a x i i ξ j, L ξ u a ϑ = ξ i i u a ϑ, and the Euler–Lagrange equations for gyrocenter motion [Eq. (64), ( L GYa / U ) u = 0, ( L GYa / μ ) u = 0, and ( d / d t ) a ( L GYa / u a ϑ ) = ( L GYa / ϑ ) u = 0] to write the first three terms in the integrand on the right-hand side of Eq. (82) as
F ( ξ i L GYa x i + L GYa u a x j δ ¯ u a x j + L GYa u a ϑ δ ¯ u a ϑ ) = ξ i [ t ( F a L GYa u a x i ) + U ( F a u a U L GYa u a x i ) L GYa u a x i { F a t + x j ( F a u a x j ) + U ( F a u a U ) + ϑ ( F a u a ϑ ) } ] + x j ( ξ i F a u a x j L GYa u a x i ) .
(86)
Then, substituting Eq. (86) into Eq. (82), using δ ¯ J = J δ ¯, and performing partial integrals yield
δ ¯ L GKa = V d 3 x [ ξ i d 3 v { t ( F a L GYa u a x i ) D t F a L GYa u a x i } + δ L GKa δ A i δ ¯ A i + δ L GKa δ A ̂ i δ ¯ A ̂ i + δ L GKa δ ϕ δ ¯ ϕ + δ L GKa δ g i j δ ¯ g i j ] + B . T . = 0 ,
(87)
where
D t F a F a t + x j ( F a u a x j ) + U ( F a u a U ) + μ ( F a u a μ ) + ϑ ( F a u a ϑ ) .
(88)
Furthermore, substituting δ ¯ A i = ξ j ( j A i ) ( i ξ j ) A j = ξ j ( j A i ) ( i ξ j ) A j, δ ¯ A ̂ i = ξ j ( j A ̂ i ) ( i ξ j ) A ̂ j = ξ j ( j A ̂ i ) ( i ξ j ) A ̂ j, δ ¯ ϕ = ξ j j ϕ = ξ j j ϕ, and δ ¯ g i j = i ξ j j ξ i into Eq. (87) and performing partial integrals, one obtains
δ ¯ L GKa = V d 3 x ξ j ( J GKa ) j + B . T . = 0 ,
(89)
where
( J GKa ) j t ( d 3 v F a L GYa u a x j ) d 3 v D t F a L GYa u a x j + 2 i ( g j k δ L GKa δ g i k ) δ L GKa δ ϕ j ϕ δ L GKa δ A k j A k δ L GKa δ A ̂ k j A ̂ k + k ( δ L GKa δ A k A j + δ L GKa δ A ̂ k A ̂ j ) .
(90)
Here, it should be noted that
( J GKa ) j = 0
(91)
holds because Eq. (89) is valid for an arbitrary infinitesimal vector field represented by ξj which vanishes on the boundary of V.
Recalling that the canonical momentum of a single particle is given by p a j L GYa / u a x j, it is found that the first term on the right-hand side of Eq. (90) represents the rate of change in the momentum of particles per volume. The pressure tensor of the particle species a is given in terms of the variational derivative δ L GKa / δ g i j as
P a i j 2 δ L GKa δ g i j 2 J ( 1 ) # J J ( d 3 v F a L GYa ( J g i j ) ) = P CGL a i j + π a i j + π Ψ a i j + P Ψ a i j ,
(92)
where # J = n represents the order of J ( j 1 , j 2 , , j n ),
P CGL a i j = d 3 v F a [ m a U 2 b i b j + μ B ( g i j b i b j ) ] ,
(93)
π a i j d 3 v F a m a U [ b i ( u a x ) j + ( u a x ) i b j ] ,
(94)
π Ψ a i j e a d 3 v F a U Ψ a U b i b j ,
(95)
and
P Ψ a i j 2 e a δ δ g i j ( d 3 x d 3 v F a Ψ a ) .
(96)
The pressure tensor P CGL a i j defined in Eq. (93) takes the Chew–Goldberger–Low (CGL) form49 and it plays an essential role in the neoclassical transport. Effects of turbulent fluctuations on the momentum transport are included in π a i j, π Ψ a i j, and P Ψ a i j, which are investigated in detail in Sec. VIII. It is found from Eqs. (3), (4), (6), (7), and (95) that π Ψ a i j vanishes when A ̂ = 0.
The particle density N a ( p ) and the particle flux Γ a i of species a are given from the functional derivatives δ L GKa / δ ϕ and δ L GKa / δ A ̂ i by
e a N a ( p ) δ L GKa δ ϕ J ( 1 ) # J J ( d 3 v F a L GYa ( J ϕ ) ) d 3 v F a L GYa ϕ + n = 1 ( 1 ) n × j 1 j n 1 ( d 3 v F a L GYa ( j 1 j n 1 ( E L ) j n ) ) ,
(97)
and
e a c Γ a i δ L GKa δ A ̂ i J ( 1 ) # J J ( L GKa ( J A ̂ i ) ) ,
(98)
respectively. The electric current density defined by j i a e a Γ a i with the particle flux Γ a i in Eq. (98) enters the gyrokinetic Ampère's law, which is derived from the variational principle δ L GKF / δ A ̂ i = 0 in Eq. (74). On the other hand, the variational derivative δ L GKa / δ A i gives
e a c Γ # a i δ L GKa δ A i J ( 1 ) # J J ( L GKa ( J A i ) ) = L GKa A i x j ( L GKa ( j A i ) ) + 2 x j x k ( L GKa ( j k A i ) ) ,
(99)
from which another type of the particle flux Γ # a i of species a is derived as
Γ # a i d 3 v F a u a x i + c e a ε ijk x j ( d 3 v F a g × [ μ b k + m a U B { ( u a x ) k ( u a x ) l b l b k } e a { Ψ a B k 1 F a x l ( F a Ψ a ( l B k ) ) } ] ) .
(100)
Here, one note that neither of the last two terms in the last line of Eq. (100) is seen to take the form of a vector component in the general spatial coordinate system. However, the sum of these two terms is shown to be rewritten by
Ψ a B k 1 F a x l ( F a Ψ a ( l B k ) ) = ( Ψ a B k ) B 1 F a l ( F a Ψ a ( l B k ) ) .
(101)
On the left-hand side of Eq. (101), derivatives are performed with regarding Ψ a as taking a functional form of Ψ a ( x i , B k ( x i ) , j B k ( x i ) ) while, on the right-hand side, a functional form of Ψ a ( x i , B k ( x i ) , j B k ( x i ) ) is considered. Then, each term of the right-hand side of Eq. (101) is found to have the form of a vector component, and the sum of the two terms on the left-hand side also becomes a vector component. Therefore, it is understood that Γ # a i give by Eq. (100) is transformed as a vector density under the general spatial coordinate transformation by noting that Fa and εijk are a scalar density and a tensor density.

The two types of particle fluxes Γ a i = ( c / e a ) δ L GKa / δ A ̂ i and Γ # a i = ( c / e a ) δ L GKa / δ A i in Eqs. (98) and (99) arise because the separation of the magnetic field into the average and fluctuating parts are done in the case of electromagnetic turbulence. As described in Ref. 31, it is Γ a i that accurately represents both the average and fluctuating parts of the particle flux and is used to evaluate the current density in Ampère's law as shown in Eqs. (36), (74), and (D1). It is found from comparing Γ # a i with Γ a i that the average part of Γ # a i equals that of Γ a i to the lowest order in δ = ρ / L while their fluctuating parts differ from each other.

It is emphasized here that ( J GKa ) j = 0 in Eq. (91) is valid for ( J GKa ) j defined in Eq. (90) where the gyrocenter distribution function Fa can be arbitrarily chosen and it does not need to be a solution of the gyrokinetic Vlasov equation given by Eq. (48) or Eq. (70). It is recalled that the variation associated with the spatial coordinate transformation should be clearly distinguished from the variation (or virtual displacement) used for deriving the gyrokinetic Vlasov equation; the fact that the former variation of the Lagrangian vanishes can be used to derive the momentum balance equation even in a more general case where the governing kinetic equation differ from the gyrokinetic Vlasov equation derived using the latter variational principle. We now assume Fa to satisfy not the gyrokinetic Vlasov equation, Eq. (48), but a more general one, that is, the gyrokinetic Boltzmann equation given by
D t F a F a t + x j ( F a u a x j ) + U ( F a u a U ) + μ ( F a u a μ ) + ϑ ( F a u a ϑ ) = K a ,
(102)
where K a represents the rate of temporal change in Fa due to collisions and/or external sources for the species a. It is assumed in the present work that
a e a d 3 v K a = 0
(103)
is satisfied by K a. Therefore, the charge density is not changed by K a. Using Eqs. (90), (92), (97)–(99), and (102), Eq. (91) implies that the solution of the gyrokinetic Boltzmann equation, Eq. (102), satisfies the canonical momentum balance equation
t ( d 3 v F a p a j ) d 3 v K a p a j + i ( P a ) j i = e a N a ( p ) j ϕ + e a c [ Γ # a k j A k + Γ a k j A ̂ k k ( Γ # a k A j + Γ a k A ̂ j ) ] ,
(104)
which can be written in the conventional dyadic notation representing vectors and tensors in terms of boldface letters as
t ( d 3 v F a p a ) d 3 v K a p a + · P a = e a N a ( p ) ϕ + e a c [ ( A ) · Γ # a + ( A ̂ ) · Γ a · ( Γ # a A + Γ a A ̂ ) ] .
(105)
Here, as expected from Noether's theorem, one can confirm that Eq. (105) takes the conservation form of the canonical momentum in the direction of the constant vector e if the term including K a vanishes and the electric and magnetic fields satisfy the symmetry conditions e · ϕ = 0 and e · A = e · A ̂ = 0. In a case where the electric and magnetic fields are axisymmetric, conservation of the toroidal angular momentum can also be derived from Eq. (105) in the same manner as shown in Sec. VI.
The canonical momentum balance equation, Eq. (105), is also written as
t ( d 3 v F a p a ) d 3 v K a p a + · P a = e a N a ( p ) ϕ + e a c [ ( Γ # a × B + Γ a × B ̂ ) A ( · Γ # a ) A ̂ ( · Γ a ) ] .
(106)
Furthermore, Eq. (106) is deformed to
t ( m a N a ( g ) V a g b ) d 3 v K a m a U b + · P a = e a ( N a ( p ) E L + N a ( g ) E T ) + e a c [ ( Γ # a × B + Γ a × B ̂ ) A ̂ ( · Γ a ) ] ,
(107)
where E L = ϕ and E T = c 1 A / t are used. Here,
N a ( g ) d 3 v F a and N a ( g ) V a g d 3 v F a U
(108)
represent the density and the parallel flux of the gyrocenters, respectively. From Eqs. (100), (102), and (108), one can obtain
N a ( g ) t + · Γ # a = d 3 v K a ,
(109)
which is used to derive Eq. (107) from Eq. (106). The first term on the left-hand side of Eq. (107) is the change rate of the density of the kinetic momentum which is obtained by extracting the vector potential term p a ( e a / c ) A + m a U b. The second and third terms on the left-hand side of Eq. (107) represent the effects of collisions (or external sources) and the pressure tensor, respectively, while the right-hand side contains Lorentz forces due to the electric and magnetic fields. The last term on the right-hand side appears due to the perturbed vector potential and it is carried over from Eq. (106).

It is noted that the time derivative terms in Eqs. (105)–(107) are missing the perpendicular kinetic momentum part. This originates from the fact that the canonical momentum p a associated with the gyrocenter Lagrangian does not include the perpendicular kinetic moment due to the perpendicular velocity v which depends on the gyrophase angle. The perpendicular part of the kinetic momentum density is given by m a Γ a and its time derivative m a Γ a / t is considered as neglected in the gyrokinetic momentum balance [Eqs. (105)–(107)], which the gyrokinetic Boltzmann equation, Eq. (102), satisfies. The leading order of the magnitude of terms in the perpendicular part of Eqs. (105)–(107) is given from the order of the Lorentz force terms and it is estimated to be O ( e a | Γ a | B / c ) = O ( m a Ω a | Γ a | ) where | Γ # a | | Γ a | N a ( g ) c | E L | / B ( ρ a / L ) N a ( g ) v T a is regarded as valid for both the average and fluctuating parts. Thus, the neglected term m a Γ a / t in the perpendicular momentum balance is smaller than the leading-order terms by a factor O ( Ω a 1 / t ). Here, the transport timescale ordering gives Ω a 1 / t ( ρ a / L ) 3 for the ensemble-averaged part while Ω a 1 / t ρ a / L is obtained for the fluctuating part from the gyrokinetic ordering. Thus, neglecting m a Γ a / t is not considered to give a significant influence on the perpendicular part of the local momentum balance although this higher-order term should be correctly included for accurately describing the flux-surface-averaged momentum balance along the symmetry direction in up-down symmetric tokamaks and stellarator-symmetric quasisymmetric stellarators.51–53 

One can follow the same procedures as used in deriving Eq. (87) to deform Eq. (85) to
δ ¯ L GKF = V d 3 x [ ξ i a d 3 v { t ( F a L GYa u a x i ) D t F a L GYa u a x i } + δ L GKF δ A i δ ¯ A i + δ L GKF δ A ̂ i δ ¯ A ̂ i + δ L GKF δ ϕ δ ¯ ϕ + δ L GKF δ g i j δ ¯ g i j ] + B . T . = 0.
(110)
Substituting δ ¯ A i = ξ j ( j A i ) ( i ξ j ) A j, δ ¯ A ̂ i = ξ j ( j A ̂ i ) ( i ξ j ) A ̂ j, δ ¯ ϕ = ξ j j ϕ into Eq. (110), and performing a partial integral, one obtains
δ ¯ L GKF = V d 3 x ξ j ( J GKF ) j + B . T . = 0 ,
(111)
where
( J GKF ) j t ( a d 3 v F a L GYa u a x i ) a d 3 v D t F a L GYa u a x i δ L GKF δ ϕ j ϕ δ L GKF δ A i j A i δ L GKF δ A ̂ i j A ̂ i + i ( δ L GKF δ A i A j + δ L GKF δ A ̂ i A ̂ j ) + 2 i ( g j k δ L GKF δ g i k ) .
(112)
Since Eq. (111) is valid for any ξj which vanishes on the boundary of V,
( J GKF ) j = 0
(113)
holds for ( J GKF ) j defined in Eq. (112) where Fa, ϕ, Ai, and A ̂ i can be arbitrarily chosen and they do not need to be determined by any governing equations.
When Fa, ϕ, and A ̂ i satisfy the gyrokinetic Boltzmann equation shown in Eq. (102) and the gyrokinetic Poisson-Ampère equations given by δ L GKF / δ ϕ = 0 and δ L GKF / δ A ̂ i = 0, Eqs. (112) and (113) lead to the total canonical momentum balance equation,
t ( a d 3 v F a p a j ) a d 3 v K a p a j + i Θ j i = δ L GKF δ A i j A i i ( δ L GKF δ A i A j ) ,
(114)
where p a j L GYa / u a x j ( e a / c ) A j ( x , t ) + m a U b j ( x , t ) and the total pressure tensor density Θ i j Θ k i g j k is defined by
Θ i j 2 δ L GKF δ g i j Θ G K i j + Θ F i j .
(115)
The gyrokinetic part Θ G K i j of Θ i j is written as
Θ G K i j 2 δ L G K δ g i j 2 a δ L GKa δ g i j = P CGL i j + π i j + π Ψ i j + P Ψ i j ,
(116)
where P CGL i j , π i j, π Ψ i j, and P Ψ i j are defined using Eqs. (93)–(96) as
[ P CGL i j , π i j , π Ψ i j , P Ψ i j ] a [ P CGLa i j , π a i j , π Ψ a i j , P Ψ a i j ] .
(117)
The field part Θ F i j of Θ i j is given by
Θ F i j 2 δ L F δ g i j 2 J ( 1 ) # J L F ( J g i j ) = g 4 π [ g i j 2 { ( E L ) k ( E L ) k + ( B k + B ̂ k ) ( B k + B ̂ k ) } { ( E L ) i ( E L ) j + ( B i + B ̂ i ) ( B j + B ̂ j ) } + 1 c { g i j ( k λ ) A ̂ k + ( i λ ) A ̂ j + ( j λ ) A ̂ i } ] ,
(118)
which contains the well-known Maxwell stress tensor and the additional terms in the same form as found in the Vlasov–Darwin model.45 
Now, using Cartesian spatial coordinates and the conventional dyadic notation representing vectors and tensors in terms of boldface letters, Eq. (114) is rewritten as
t ( a d 3 v F a p a ) a d 3 v K a p a + · Θ = ( A ) · δ L GKF δ A · ( δ L GKF δ A A ) .
(119)
Here, δ L GKF / δ A is given by
δ L GKF δ A = j # c 1 4 π × ( B + B ̂ ) ,
(120)
where j # is defined from Γ # a in Eq. (100) by
j # c δ L G K δ A a e a Γ # a j ( g c ) + c × M # .
(121)
Here, j ( g c ) is the gyrocenter current defined by Eq. (C12) in  Appendix C and M # is given by
M # a e a d 3 v F a [ μ b + m a U B ( u a x ) e a { Ψ a B 1 F a x l ( F a Ψ a ( l B ) ) } ] .
(122)
Now, one can confirm the validity of Noether's theorem again from Eq. (119), which takes the conservation form of the total canonical momentum in the direction specified by the constant vector e when the background magnetic field satisfies the symmetry condition e · A = 0 and a d 3 v K a p a can be ignored. In Sec. VI, toroidal angular momentum conservation is derived in the case of the axisymmetric background field. It should be noted that no specific conditions to determine A are imposed from the variational principle in contrast to ϕ and A ̂ which are variationally determined. Thus, δ L GKF / δ A given in Eq. (120) does not vanish generally.
The total canonical momentum balance equation in Eq. (119) can be rewritten as
t ( a d 3 v F a p a ) a d 3 v K a p a + · Θ = δ L GKF δ A × B A · ( δ L GKF δ A ) ,
(123)
which is also deformed to
t ( a d 3 v F a m a U b ) a d 3 v K a m a U b + · Θ = ρ c ( g c ) E T + δ L GKF δ A × B ,
(124)
where ρ c ( g c ) a e a N a ( g ) and E T c 1 A / t. Equation (124) represents the total balance equation of the kinetic momentum instead of the canonical one. The effects of collisions (or external sources) and the total pressure tensor are shown on the left-hand side of Eq. (124) while the Lorentz forces due to the back ground inductive field and the background magnetic field appear on the right-hand side.
Finally, Eq. (124) can be deformed through vector calculus to
t ( a d 3 v F a m a U b + 1 4 π c ( D L × B ) ) + · ( Θ + D L E T + E T D L 4 π ) + ( E T · D L 4 π ) = a d 3 v K a m a U b + ( δ L GKF δ A ) T × B ,
(125)
which is written in more detail as
t ( a d 3 v F a m a U b + 1 4 π c ( D L × B ) ) + · ( P CGL + π + π Ψ + P Ψ ) + ( | E L | 2 8 π + E T · D L 4 π ) · ( E L E L + D L E T + E T D L 4 π ) + ( | B + B ̂ | 2 8 π ) · ( ( B + B ̂ ) ( B + B ̂ ) 4 π ) ( λ · A ̂ 4 π c ) + · ( ( λ ) A ̂ + A ̂ ( λ ) 4 π c ) = a d 3 v K a m a U b + ( ( j # ) T c × ( B + B ̂ ) 4 π ) × B .
(126)
In Eq. (126), D L is the longitudinal part of the displacement vector defined by Eq. (C4) and j # is defined in Eq. (121). The change rate of the kinetic momentum density plus the electromagnetic momentum density ( D L × B ) / ( 4 π c ) is described by Eq. (126). The left-hand side of Eq. (126) shows all terms of momentum transport written as the divergence of pressure tensors due to particles' motion and Maxwell stresses including both average and fluctuating parts of the electromagnetic field. Except for the terms on the right-hand side, Eq. (126) takes the conservation form similar to that of the total momentum conservation equation of the Vlasov-Darwin model derived in Ref. 45. Since j # = a e a Γ # a does not accurately represents the fluctuating part of the current density, ( j # ) T / c ( 4 π ) 1 × ( B + B ̂ ) cannot be neglected as far as turbulent electromagnetic fields exist. In Sec. VII, the self-consistency condition to determine the average field B is considered to make the ensemble average of Eq. (126) take the conservation form.
In this section, we investigate the momentum balance in toroidal systems based on the results obtained in Sec. V. The background magnetic field B with and without symmetry are considered and the momentum balance equations averaged over the ensemble and the flux surface are examined. Here, the background magnetic field B is assumed to satisfy the toroidal MHD equilibrium equation
1 4 π ( × B ) × B = P 0 ,
(127)
where P0 is a magnetic flux surface function representing equilibrium pressure.
The axisymmetric toroidal background magnetic field is represented by
B = I ζ + ζ × χ ,
(128)
where ζ and χ represents the toroidal angle and the poloidal flux (divided by 2 π), respectively, and the covariant toroidal component I is a flux surface function which is independent of the toroidal and poloidal angles. Denoting the major radius by R and writing the contravariant basis vector e ζ in the toroidal direction as
e ζ R 2 ζ x ζ ,
(129)
one obtains the following relation:
e ζ = R 1 [ ( R ) e ζ e ζ ( R ) ] = n × I = I × n ,
(130)
where n R 1 ( e ζ × R ) is the unit vector parallel to the direction of the major axis and I is the unit tensor. It is shown from Eq. (130) that an arbitrary symmetric tensor S ( S i j = S j i ) satisfies
( · S ) · e ζ = · ( S · e ζ ) .
(131)
It is also noted that · e ζ = 0 and
e ζ · S = · ( S e ζ ) ,
(132)
where S is an arbitrary scalar function.
In the axisymmetric background field B, A can also be given by the axisymmetric field which satisfies
e ζ · A = 1 R A × n .
(133)
Then, the inner product of e ζ and Eq. (114) is taken and Eqs. (130), (131), and (133) are used to derive
t ( a d 3 v F a p a · e ζ ) + · [ ( Θ + δ L GKF δ A A ) · e ζ ] = a d 3 v K a p a · e ζ .
(134)
Except for the right-hand side, Eq. (134) takes the conservation form of the canonical momentum conjugate to the toroidal angle as expected from Noether's theorem. It is found from the assumption given in Eq. (103) that a d 3 v K a p a · e ζ = a d 3 v K a m a U b ζ where b ζ = b · e ζ. In the zero-gyroradius limit, when using a particle collision operator for K a, one finds that a d 3 v K a m a U b ζ to vanish because the momentum of particles is conserved in collisions. Furthermore, it can be shown that when K a is given by the collision operator which appropriately includes the finite gyroradius effect,57–59, a d 3 v K a m a U b ζ can be written as a divergence of the sum of classical momentum transport fluxes. Therefore, without external momentum sources, Eq. (134) keeps the conservation form even though the collision term is present. In addition to the case of axisymmetry, the canonical momentum conservation is confirmed in other cases of symmetry under continuous isometric transformations such as a translational symmetry and a helical (or screw) symmetry.37 
Next, the toroidal component of the momentum balance equation in Eq. (126) is considered. The transverse part of j # on the right-hand side of Eq. (126) can be written in terms of a certain field B # as
( j # ) T = c 4 π × B # ,
(135)
which is combined with B × e ζ = χ to derive
( ( j # ) T × B ) · e ζ = ( j # ) T · χ = · ( c 4 π B # × χ ) .
(136)
One also obtains
( ( × ( B + B ̂ ) ) × B ) · e ζ = ( × ( B + B ̂ ) ) · χ = · ( ( B + B ̂ ) × χ ) .
(137)
Now, taking the inner product of Eq. (126) and e ζ, it is found that the toroidal angular momentum balance equation can be written in the following form:
t ( a d 3 v F a m a U b ζ + 1 4 π c D L · χ ) + · [ ( Θ + D L E T + E T D L 4 π ) · e ζ + E T · D L 4 π e ζ ] + · [ 1 4 π ( B # B B ̂ ) × χ ] = a d 3 v K a m a U b ζ ,
(138)
where ( D L × B ) · e ζ = D L · χ is used. The expression of the divergence term on the left-hand side of Eq. (138) is straightforwardly derived from Eq. (126) using Eqs. (131), (132), (136), and (137). Without external momentum sources, Eq. (138) keeps the conservation form in the same way as Eq. (134).

It is shown above that the symmetry of the pressure tensor is essential in deriving the equation of the toroidal angular momentum conservation in an axisymmetric system. From this point of view, the derivation of the symmetric pressure tensor from the variational derivative of the Lagrangian with respect to the metric tensor is useful.

In non-axisymmetric toroidal systems, the background field B is expressed by
B = ψ ( s ) × θ + ζ × χ ( s ) ,
(139)
where θ and ζ are the poloidal and toroidal angles, respectively, and s is an arbitrarily chosen flux-surface label. Here, we assume the background B to be stationary, B / t = 0, for simplicity, and use the Hamada coordinates60, ( s , θ , ζ ), in which the Jacobian g [ s · ( θ × ζ ) ] 1, the poloidal field B θ B · θ, and the toroidal field B ζ B · ζ are flux-surface functions. Then, the contravariant basis vector e ζ in the toroidal direction is written as
e ζ x ζ g s × θ ,
(140)
which satisfies · e ζ = 0 and B × e ζ = χ. Therefore, equations in the same form as Eqs. (132), (136), and (137) hold while Eq. (131) does not because Eq. (130) is not valid in the non-axisymmetric case. Now, taking the inner product of Eq. (126) and e ζ gives the toroidal angular momentum balance equation in the following form:
t ( a d 3 v F a m a U b ζ + 1 4 π c D L · χ ) + · ( ) + e ζ · ( · ( P CGL + ) ) = 0 ,
(141)
where the effects of external momentum sources are ignored. The conservation form is broken in Eq. (141) because of e ζ · ( · ( P CGL + ) ) on the left-hand side of Eq. (141). Here, the CGL and other pressure tensors and anisotropic Maxwell stress terms are included in e ζ · ( · ( P CGL + ) ) while, using Eq. (132), the tensors proportional to the unit tensor can be transferred to the inside of the divergence term · ( ) in Eq. (141). Using the flux-surface average defined by d θ d ζ g / V ( s ) with V ( s ) d θ d ζ g, it is found that the divergence term · ( ) in Eq. (141) is annihilated by the flux-surface average because, as seen from Eqs. (132), (136), (137), and (140), the inner products of the vectors in ( ) and s vanish and
· T = 1 V ( s ) d d s ( V ( s ) T · s )
(142)
holds for any vector T. Then, taking the flux-surface average, Eq. (141) is reduced to the more compact form,
t a d 3 v F a m a U b ζ + 1 4 π c D L · χ + e ζ · ( · ( P CGL + ) ) = 0.
(143)
It is recalled that, in neoclassical theory for non-axisymmetric systems,63 the lowest-order toroidal viscosity is given by e ζ · ( · P CGL ) . It is shown in Sec. VII that, when using the ensemble average and the gyroradius expansion of Eq. (143) in general, non-axisymmetric toroidal systems, this neoclassical toroidal viscosity becomes a dominant term.
In this subsection, quasi-axisymmetric toroidal systems35 are considered using Eq. (143) with the Hamada coordinates ( s , θ , ζ ) to represent the equilibrium magnetic field B. The quasi-axisymmetry is characterized by B = | B | being independent of the toroidal angle, B / ζ = 0, which is equivalent to B / ζ B = 0 in the Boozer coordinates64, ( s , θ B , ζ B ) as proved in Ref. 65. In the quasi-axisymmetric equilibrium field B, the CGL-type pressure tensor P CGL P b b + P ( I bb ) is shown to satisfy
e ζ · ( · P CGL ) = ( P P ) ln B ζ = 0 ,
(144)
which implies that the dominant neoclassical toroidal viscosity which exists in general non-axisymmetric systems vanishes as in the axisymmetric systems. Thus, one of the factors preventing conservation of the toroidal angular momentum in Eq. (143) disappears even though perfect conservation is not allowed. Because of Eq. (144), the magnitude of the dominant terms in Eq. (143) becomes of higher order. Then, as mentioned in Sec. II, basic gyrokinetic equations including higher-order terms, which are not considered in the present study, are required to accurately describe the flux-surface-averaged momentum balance along the symmetry direction in stellarator-symmetric quasisymmetric stellarators as well as up-down symmetric tokamaks.51–53 

In this section, the momentum balance equation in Eq. (126) is ensemble-averaged, by which all terms in the equation are smoothed to make their space-time scales of variations much larger than those of fluctuations. The ensemble average is used as the basic method of statistical mechanics to obtain the macroscopic mean values of physical valuables and it can also be considered to equal the local space–time average, the definition of which is described in detail in Ref. 66. The fluctuations are assumed to have wavelengths of the order of gyroradii in the directions perpendicular to the background magnetic field, and they are treated by the WKB representation in Sec. VIII.

Since the background magnetic field B is considered to include no fluctuations, one can write B = B ens, where ens represents the ensemble average. It should also be noted here that no equation to determine the background magnetic field B is given from the variational principle while B is allowed to change with time in the present model. If the variational condition δ L GKF / δ A = 0 was employed, B would include fluctuation components as seen from Eq. (120). Then, ( δ L GKF / δ A ) T ens = 0 is assumed here instead of δ L GKF / δ A = 0 as the condition for determining B. Using Eq. (120), ( δ L GKF / δ A ) T ens = 0 is written as
× B = 4 π c ( j # ) T ens ,
(145)
which seems to be appropriate because the equilibrium part of j # c δ L G K / δ A [see Eq. (121)] is equal to that of j c δ L G K / δ A ̂ [see Eq. (D1)] to the lowest order in the gyroradius expansion and it represents the current density consistent with the MHD equilibrium. In the neoclassical transport theory,49 the self-consistency condition for the background field B evolving in the transport timescale is given from the MHD equilibrium equation. In passing, the condition for B can be derived as the variational equation in the drift kinetic model44,61 and the variational derivation of the time-evolving axisymmetric background field is considered in the gyrokinetic model.21,62
With the help of Eq. (121), Eq. (145) can be rewritten as
× H # ens = 4 π c ( j ( g c ) ) T ens .
(146)
Here, j ( g c ) is the gyrocenter current given by Eq. (C12) and H # is defined using M # in Eq. (122) as
H # B + B ̂ 4 π M # ,
(147)
from which one has
H # ens = B 4 π M # ens .
(148)
Using Eq. (145), the ensemble average of Eq. (126) is written as
t a d 3 v F a m a U b + 1 4 π c ( D L × B ) ens + · P CGL + π + π Ψ + P Ψ ens + | E L | 2 8 π + E T · D L 4 π ens · E L E L + D L E T + E T D L 4 π ens + | B + B ̂ | 2 8 π ens · ( B + B ̂ ) ( B + B ̂ ) 4 π ens λ · A ̂ 4 π c ens + · ( λ ) A ̂ + A ̂ ( λ ) 4 π c ens = a d 3 v K a m a U b ens .
(149)
The ensemble-averaged momentum balance equation, Eq. (149), takes the conservation form when no external sources of momentum exist and the right-hand side is written as a divergence of the tensor representing classical momentum transport. It is emphasized here that the ensemble-averaged momentum conservation described above is satisfied even in non-axisymmetric systems when the background field B is determined by the equilibrium condition given in Eq. (145). It is interesting to compare Eq. (149) with the momentum conservation law shown by Eqs. (31)–(33) in Ref. 45 for the Vlasov–Poisson–Ampère (or Vlasov–Darwin) system in which collisional effects are ignored and the magnetic field is not divided into background and turbulent parts. One can see that kinetic and electromagnetic momenta, kinetic and electromagnetic pressure tensors, and longitudinal and transverse electric fields in the momentum conservation equation of the Vlasov–Darwin system appear in Eq. (149) in a similar manner and that Eq. (149) additionally includes polarization, magnetization, and other higher-order terms due to finite-gyroradius effects and electromagnetic microturbulence. The similarities and differences described above are regarded as natural results because the electromagnetic gyrokinetic systems are derived from the Vlasov–Darwin system through ordering assumptions regarding gyroradius scales and fluctuation amplitudes.
Hereafter in this section, Eq. (149) is expanded using the ordering parameter given by the normalized gyroradius δ = ρ / L which is the ratio of the gyroradius ρ to the equilibrium scale length L. The zeroth-order part F a 0 of the distribution function is assumed to be given by the local Maxwellian as
F a 0 = F a 0 ens = F a M D a 0 f a M N a 0 D a 0 ( m a 2 π T a 0 ) 3 / 2 exp [ 1 T a 0 ( 1 2 m a U 2 + μ B ) ] ,
(150)
where N a 0 and T a 0 are the background density and temperature of the particle species a, respectively, and D a 0 B / m a is the zeroth-order part of D a B a * / m a. Here, the transport timescale ordering is used for the ensemble-averaged variables in Eq. (149) which means that ens / t δ 2 ( v T / L ) ens where vT is thermal velocity. The CGL pressure tensor defined in Eq. (93) is expanded in δ = ρ / L as
P CGL a = P a 0 I + ( P CGL a ) 1 + O ( δ 2 ) ,
(151)
where the zeroth-order part is isotropic and expressed by the scalar pressure, P a 0 N a 0 T a 0. Then, it is found that the zeroth-order part of the ensemble-averaged momentum conservation equation, Eq. (149), is given by
( P 0 + B 2 8 π ) 1 4 π · ( BB ) = 0 ,
(152)
where the equilibrium pressure is defined by P 0 a P a 0 a N a 0 T a 0. Equation (152) is easily confirmed to be equivalent to Eq. (127) representing the MHD equilibrium.
The first-order part of Eq. (149) comes only from the CGL pressure tensor because the other pressure tensors and the turbulent Maxwell stress tensors are of the second order. Thus, one obtains
· ( P CGL ) 1 ens = 0.
(153)
The turbulent part F ̂ a of the distribution function Fa has no contribution to ( P CGL ) 1 and to Eq. (153) because F ̂ a ens = 0. In neoclassical transport theory,48–50 the parallel component of · ( P CGL ) 1 ens in Eq. (153) automatically vanishes because of a quasineutrality condition and the momentum conservation property of the collision term. Also, the flux-surface average of the toroidal component of Eq. (153),
e ζ · ( · ( P CGL ) 1 ) = 0 ,
(154)
automatically holds in axisymmetric and quasi-axisymmetric toroidal systems as described in Sec. VI. Here, denotes the average over the flux surface and the ensemble. In general, non-axisymmetric systems such as stellarator and heliotron plasmas, Eq. (154) is not automatically satisfied but it imposes an ambipolarity condition on neoclassical particle fluxes, from which the background radial electric field can be determined.63 

The effects of electromagnetic microturbulence with perpendicular wavelengths on the gyroradius scale appear on the second order in Eq. (149). In Sec. VIII, the turbulence contributions to the momentum transport are investigated in detail using the WKB representation for turbulent fluctuations.

The WKB (or ballooning) representation10 is useful for treating turbulent fluctuations which have small wavelengths of the order of the gyroradius ρ in the directions perpendicular to the background magnetic field. The WKB representation for the fluctuation part Q ̂ of an arbitrary function Q ( x , t ) takes the form
Q ̂ ( x , t ) = k Q ̂ k ( x , t ) exp [ i S k ( x , t ) ] ,
(155)
where Q ̂ k ( x , t ) has the equilibrium gradient scale length L while the eikonal S k ( x , t ) represents a rapid variation with the wave number vector k S k ( ρ 1 ) which satisfies k · b = 0.
In the gyroradius expansion using δ = ρ / L 1, the zeroth-order gyrocenter distribution function is assumed to be given by the local Maxwellian, F a 0 = F a M = D a 0 f a M, as shown in Eq. (150). The fluctuation part of Fa appears in the first order and it is given by the WKB representation as F ̂ a 1 = k F ̂ a 1 k exp [ i S k ( X , t ) ], where X is the gyrocenter position vector. As shown in Ref. 31, the k -component F ̂ a 1 k is written as F ̂ a 1 k = D a 0 f ̂ a 1 k , where f ̂ a 1 k is given by
f ̂ a 1 k = f a M e a T a ψ ̂ a ξ k + h ̂ a k .
(156)
Here, h ̂ a k is the nonadiabatic part of the turbulent distribution function and the gyrophase-averaged potential ψ ̂ a ξ k is defined by
ψ ̂ a ξ k = J 0 ( k v Ω a ) ( ϕ ̂ k U c A ̂ k ) + J 1 ( k v Ω a ) v c B ̂ k k ,
(157)
where J0 and J1 are the Bessel functions. In addition, another kind of gyrophase-averaged potential is defined by
χ ̂ a ξ k = k v Ω a J 1 ( k v Ω a ) ( ϕ ̂ k U c A ̂ k ) + [ k v Ω a J 0 ( k v Ω a ) J 1 ( k v Ω a ) ] v c B ̂ k k ,
(158)
which is used later to express the turbulent pressure tensor in Eq. (166).
It is now recalled that the ensemble-averaged quantities are smooth spatial functions with the gradient scale length L. For arbitrary real-valued turbulent fluctuations Q ̂ and Q ̂ , Q ̂ k * Q ̂ k ens = 0 for k k and Q ̂ Q ̂ ens = k Q ̂ k * Q ̂ k ens hold. In the ensemble-averaged momentum balance equation given by Eq. (149), the effects of the electromagnetic turbulence on the momentum transport enter π Ψ ens π ens, and P Ψ ens through the correlation between the turbulent distribution function and the turbulent potential. Using Eqs. (94), (95), and (117), and neglecting terms of higher orders in δ = ρ / L, one finds that
π Ψ ens = bb a ( n a 0 e a 2 m a c 2 ( A ̂ ) 2 ens + e a c d 3 v U h ̂ a A ̂ ϑ ens ) = bb ( ω p 2 4 π c 2 ( A ̂ ) 2 ens + 1 c j ̂ A ̂ ens )
(159)
and
π ens = a d 3 v F a ens m a U [ b i ( u a x ) j ens + ( u a x ) i ens b j ] + π turb ens ,
(160)
where the turbulent part π turb ens is given by
π turb ens = a d 3 v F ̂ a m a U [ b i ( u ̂ a x ) j + ( u ̂ a x ) i b j ] = c B k [ b ( k × b ) + ( k × b ) b ] × a d 3 v m a U Im [ h ̂ a k * ψ ̂ a ϑ k ens ] .
(161)
Using Eqs. (96) and (117), P Ψ ens is written as
P Ψ ens = P ϕ ens ens + P Ψ turb ens ,
(162)
where the effects of the ensemble-averaged (or background) electric field E L ens ϕ ens and the turbulent electromagnetic field are included in
P ϕ ens ens = n a 0 m a c 2 B 2 ( b × ϕ ens ) ( b × ϕ ens ) m a c 2 2 e a B 2 × [ ( n a 0 T a 0 ) ϕ ens + ϕ ens ( n a 0 T a 0 ) ( I bb ) ( n a 0 T a 0 ) · ϕ ens ] + n a 0 m a 0 c 2 T a 0 2 e a B 2 × [ ( I 3 b b ) bb : ( ϕ ens ) + b ( b · ϕ ens ) + ( b · ϕ ens ) b + · ( bb ) ϕ ens + ϕ ens · ( bb ) ϕ ens · ( bb ) + 2 ln B ϕ ens + 2 ϕ ens ln B 2 ( ϕ ens · ln B ) ( I bb ) ]
(163)
and
P Ψ turb ens = P Ψ ad ens + P Ψ nad ens ,
(164)
respectively, where
P Ψ ad ens = ( I bb ) a k n a 0 e a 2 2 T a [ | ϕ ̂ k | 2 + T a m a c 2 | A ̂ k | 2 ens × { 1 Γ 0 ( b a k ) } + T a m a c 2 B ̂ k | 2 k 2 ens [ 1 2 b a k × { Γ 0 ( b a k ) Γ 1 ( b a k ) } ] v T a c Re ϕ ̂ k * B ̂ k k ens × ( 2 b a k ) 1 / 2 { Γ 0 ( b a k ) Γ 1 ( b a k ) } ]
(165)
and
P Ψ nad ens = a k [ 1 2 ( I bb ) k k k 2 ] × e a d 3 v Re [ h ̂ a k * χ ̂ a ϑ k ens ] 1 c j ̂ A ̂ ens 1 2 c j ̂ · A ̂ ens ( I bb )
(166)
are derived from the adiabatic and nonadiabatic parts of the distribution function in Eq. (156), respectively. On the right-hand side of Eq. (165), b a k k 2 T a / ( m a Ω a 2 ) is used, and the functions Γ 0 and Γ 1 are defined by Γ 0 ( b ) I 0 ( b ) exp ( b ) and Γ 1 ( b ) I 1 ( b ) exp ( b ), respectively, where I0 and I1 are the modified Bessel functions.
One finds that the pressure tensor terms consisting of only bb and I parts cannot produce transport of toroidal and poloidal momenta across flux surfaces in toroidal magnetically confined systems. In the axisymmetric toroidal system, in which the magnetic field is given by Eq. (128), one can use ( k × b ) · χ = B R 2 ζ · k and Eqs. (161), (164), (165), and (166) to obtain the radial transport of the toroidal angular momentum due to the interaction of the nonadiabatic distribution function and the turbulent electromagnetic potential as
χ · π turb + P Ψ turb ens · R 2 ζ = a k d 3 v [ c I B m a U Im [ h ̂ a k * ψ ̂ a ϑ k ens ] ( k · R 2 ζ ) + e a Re [ h ̂ a k * χ ̂ a ϑ k ens ] ( k · χ ) ( k · R 2 ζ ) k 2 ] .
(167)
The flux-surface average of Eq. (167) agrees with the low-flow ordering limit of the result given in Eq. (53) of Ref. 67 where the turbulent radial transport of the toroidal angular momentum double-averaged over the ensemble and the flux surface is presented for the case of high-flow ordering.21,66,67 The radial turbulent transport of the toroidal angular momentum in Eq. (167) is not the flux-surface average but a spatially-local expression. It is shown in the case of the low-flow ordering51,68 that, in the axisymmetric configuration with up-down symmetry, the flux-surface average of Eq. (167) vanishes even though the local value of Eq. (167) does not. It can also be shown from Eq. (163) that the flux-surface average of χ · P ϕ ens ens · R 2 ζ vanishes as well.
Next, let us consider the Maxwell stress terms in the ensemble-averaged momentum balance equation in Eq. (149). The Maxwell stress due to the electric field is dominantly given by E L because the magnitude of E T c 1 A / t is smaller than that of E L by a factor of δ. On the left-hand side of Eq. (149), the turbulent magnetic pressure tensor also appears as
1 8 π | B ̂ | 2 ens I 1 4 π B ̂ B ̂ ens ,
(168)
which has the opposite sign to that of the Maxwell stress tensor due to turbulent magnetic fields.

It is also noted that the terms proportional to λ and A ̂ in Eq. (149) are negligible compared with the other magnetic Maxwell stress terms because c 1 | λ | | A ̂ | / | B ̂ | 2 | j L | / ( c k | B ̂ | ) ( E L / t ) / ( c k | B ̂ | ) ( v T / L ) ( c T / e B ) ( e ϕ ̂ / T ) / ( c 2 | B ̂ | / B ) ( ρ / L ) ( v T 2 / c 2 ) 1, where Eq. (37), · j L = ρ c / t = ( 4 π ) 1 ( · E L ) / t, / t v T / L, and e ϕ ̂ / T | B ̂ | / B ρ / L are used.

One can also find that the nonadiabatic distribution function produces the turbulent current j ̂ which correlates with the turbulent vector potential A ̂ and produces the pressure tensor given by
1 c j ̂ A ̂ ens b b 1 c j ̂ A ̂ ens 1 2 c j ̂ · A ̂ ens ( I bb ) ,
(169)
which are included in Eqs. (159) and (166). In the wavenumber representation, the turbulent magnetic field is given by
B ̂ k = i k × A ̂ k = B ̂ k b + B ̂ k ,
(170)
where B ̂ k = i b · ( k × A ̂ k ) and B ̂ k = i ( k × b ) A ̂ k . The turbulent vector potential is written by
A ̂ k = A ̂ k b + A ̂ k ,
(171)
where A ̂ k = b · A ̂ k and A ̂ k = k 2 ( b × k ) [ ( b × k ) · A ̂ k ]. Here, the Coulomb gauge condition k · A ̂ k = 0 is used. Then, one has
| B ̂ k | 2 = | B ̂ k | 2 + | B ̂ k | 2 = k 2 ( | A ̂ k | 2 + | A ̂ k | 2 ) ,
(172)
and Ampère's law is given by
j ̂ k = c 4 π k 2 A ̂ k .
(173)
Using Eqs. (170)–(173), one obtains
1 4 π B ̂ B ̂ ens 1 c j ̂ A ̂ ens = 1 4 π k | A ̂ k | 2 ens ( k k k 2 I ) = 1 4 π k | B ̂ k | 2 ens ( k k k 2 I ) ,
(174)
and the summation of Eqs. (168) and (169) are written as
Eq . (168) + Eq . (169) = 1 8 π k [ | B ̂ k | 2 + 2 | B ̂ k | 2 ens b b | B ̂ k | 2 ens k k k 2 | B ̂ k | 2 + | B ̂ k | 2 ens ( b × k ) ( b × k ) k 2 ] ,
(175)
which is given in terms of only the turbulent magnetic field.
Again, in considering the axisymmetric toroidal system, the radial transport of the toroidal angular momentum due to the turbulent magnetic field is represented by a component of Eq. (175) obtained from a double-dot product with a dyad ( χ ) ( R 2 ζ ), which is equivalent to that of Eq. (174) and written as
χ · [ 1 4 π B ̂ B ̂ ens 1 c j ̂ A ̂ ens ] · R 2 ζ = 1 4 π k | A ̂ k | 2 ens ( k · χ ) ( k · R 2 ζ ) = 1 4 π k | B ̂ k | 2 ens ( k · χ ) ( k · R 2 ζ ) k 2 .
(176)
Taking the flux-surface average of Eq. (176) yields the same expression of the toroidal angular momentum transport across the flux surface caused by the turbulent magnetic field as derived in Refs. 67 and 66. However, in the same manner as in the case of Eq. (167), the flux-surface average of Eq. (176) is shown to vanish in the axisymmetric configuration with up-down symmetry under the low-flow ordering.51,68

It is also found from Eqs. (159) and (165) that the adiabatic part of the perturbed distribution function in Eq. (156) produces pressure tensor terms in π Ψ ens and P Ψ ad ens, which are given in terms of turbulent electrostatic and vector potentials although these terms consist of only the bb and I parts so that they cannot produce toroidal and poloidal momentum transport across flux surfaces in toroidal plasmas.

In this paper, the Eulerian (or Euler-Poincaré) variational formulation is presented to obtain the governing equations of the electromagnetic turbulent gyrokinetic system, for which the local momentum balance equation is derived from the invariance of the Lagrangian of the system under an arbitrary spatial coordinate transformation. In addition, the effects of collisions and external sources are taken into account in the momentum balance equation.

Using the gyrokinetic Lagrangian, which retains proper electromagnetic potential terms and taking the variational derivatives of the Lagrangian with respect to the electrostatic and vector potentials of the perturbed magnetic field, one can obtain the gyrokinetic Poisson equation and Ampère's law where the effects of the polarization and magnetization due to finite gyroradii and electromagnetic microturbulence are correctly included. Especially, the derived gyrokinetic Ampère's law can accurately express the current density from the microscopic gyroradius scale to the macroscopic equilibrium scale so that it is useful for long-time and global gyrokinetic turbulence simulations of high beta plasmas.

The local momentum balance equation obtained in the present work contains the symmetric pressure tensor which is derived from the variational derivative of the Lagrangian with respect to the metric tensor. It is shown that the pressure tensor obtained for the whole system consisting of all particles and fields involves the gyrokinetic and field parts; the neoclassical and turbulent momentum transport processes are described by the former part while the Maxwell stress is by the latter.

One can confirm from the momentum balance equation that, when the background magnetic field has a symmetry such as a translational one and an axisymmetry, the canonical momentum conjugate to the coordinate in the symmetry direction is conserved as predicted by Noether's theorem. The symmetry of the pressure tensor is found to be an important property for derivation of the momentum conservation in the symmetric background field. When the background field is assumed to satisfy the appropriate condition representing the macroscopic Ampère's law, the ensemble-averaged total momentum balance equation is found to take the conservation form even in the asymmetric background field. Thus, this condition can be conveniently applied to long-time gyrokinetic simulations in which the change in the background field occurs with the relaxation of high-beta plasmas. It is also shown that, in the toroidal systems with the quasi-axisymmetric background field, the toroidal angular momentum is not rigorously conserved although the flux-surface-averaged neoclassical toroidal viscosity, which is a dominant component for breaking the toroidal momentum conservation in general, non-axisymmetric systems, vanishes.

The WKB representation is employed to derive detailed expressions of the ensemble-averaged pressure tensor due to the electromagnetic microturbulence, which provide a means for evaluating the local turbulent momentum transport by the local flux-tube gyrokinetic simulation. The radial transport fluxes of the toroidal angular momentum caused by the nonadiabatic distribution function and the turbulent electromagnetic fields in the axisymmetric system are represented as a non-diagonal component of the pressure tensor, which are shown to agree with the results from the previous works based on the classical gyrokinetic formulation. The local pressure tensor represented by a symmetric 3 × 3 matrix contains further information on momentum transport which is useful for more detailed analyses of transport processes by gyrokinetic simulations.

This work is supported in part by the JSPS Grants-in-Aid for Scientific Research (Grant No. 19H01879) and in part by the NIFS Collaborative Research Program (No. NIFS23KIPT009).

The author has no conflicts of interest to disclose.

Hideo Sugama: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Writing – original draft (lead).

The authors have no conflicts to disclose.

Hideo Sugama: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Writing – original draft (lead).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

The potential field defined in Eq. (3) is decomposed here as
Ψ a ( Z , t ) ϕ ( X , t ) + Ψ E 1 a ( X , μ , t ) + Ψ A ̂ 1 a ( X , U , μ , t ) + Ψ E 2 a ( X , μ , t ) + Ψ E A ̂ a ( X , U , μ , t ) + Ψ A ̂ 2 a ( X , U , μ , t ) ,
(A1)
where
Ψ E 1 a ( X , μ , t ) ϕ ( X + ρ a , t ) ϑ ϕ ( X , t ) ,
(A2)
Ψ A ̂ 1 a ( X , U , μ , t ) 1 c [ v · A ̂ ϑ + v B a · A ̂ ϑ ]
(A3)
Ψ E 2 a ( X , μ , t ) e a 2 B μ ( ϕ ̃ ) 2 ϑ ,
(A4)
Ψ E A ̂ a ( X , U , μ , t ) e a c B μ ϕ ̃ ( v · A ̂ ) ϑ ,
(A5)
and
Ψ A ̂ 2 a ( X , U , μ , t ) e a 2 m a c 2 | A ̂ | 2 ϑ e a 2 c 2 B μ [ ( v · A ̂ ) v · A ̂ ϑ ] 2 ϑ .
(A6)
It should be noted that, in this appendix, the Cartesian spatial coordinates and the conventional dyadic notation representing vectors and tensors in terms of boldface letters are used. Then, the electrostatic potential ϕ ( X + ρ a ) and the perturbed vector potential A ̂ ( X + ρ a ) are Taylor expanded about the gyrocenter position X as
[ ϕ ( X + ρ a ) A ̂ ( X + ρ a ) ] = n = 0 1 n ! ρ a j 1 ρ a j n j 1 j n [ ϕ ( X ) A ̂ ( X ) ] ,
(A7)
where Xj and ρj (j = 1, 2, 3) and the Cartesian spatial coordinates of the gyrocenter position vector X and the gyroradius vector ρ, respectively, and the partial derivatives are represented using the simplified notation, j 1 j n n / X j 1 X j n. Here, for simplicity, we omit the t-dependence of ϕ and employ the summation convention that the same symbol used for a pair of upper and lower indices indicates a summation over the range { 1 , 2 , 3 }. Therefore, summation notations j 1 = 1 3 j n = 1 3 are dropped in Eq. (A7).
Using Eqs. (A2) and (A7), one can write the part of the potential function which linearly depends on the electric field and its derivatives as
Ψ E 1 a = n = 1 α a j 1 j n n ! j 1 j n ϕ ( X ) = n = 1 α a j 1 j n n ! j 1 j n 1 ( E L ) j n ( X ) ,
(A8)
where the gyrophase average of a product of n gyroradius vector components is denoted by
α a j 1 j n ρ a j 1 ρ a j n ϑ .
(A9)
We see that α a j 1 j n is symmetric with respect to arbitrary permutations of the indices j 1 , , j n. We also find that α a j 1 j n = 0 for odd n and
α a j 1 j 2 l = 1 ( 2 l ) ! σ S 2 l η a j σ ( 1 ) j σ ( 2 l ) ,
(A10)
where S 2 l is the symmetric group of permutations of the set { 1 , 2 , , 2 l } and η a j 1 j 2 l is defined by
η a j 1 j 2 l = ( 2 l ) ! ( l ! ) 2 ( ρ a 2 ) 2 l h j 1 j 2 h j 3 j 4 h j 2 l 1 j 2 l ,
(A11)
with ρ a ( c / e a ) 2 m a μ / B and h i j δ i j b i b j. Here, bi is the ith component of b B / B and δij represents the Kronecker delta; δ i j = 1 (for i = j), 0 (for i j).
Next, using Eqs. (A3), the part which linearly depends on the vector potential and its derivatives is written as
Ψ A ̂ 1 a = 1 c [ ( U b + v B a ) · A ̂ ( X + ρ a ) ϑ + v · A ̂ ( X + ρ a ) ϑ ] = 1 c [ ( U b + v B a ) · A ̂ ( X + ρ a ) ϑ + Ω a ( ρ a × b ) · A ̂ ( X + ρ a ) ϑ ] = 1 c n = 0 1 n ! [ α a j 1 j n ( U b i + v B a i ) + Ω a ε klm α a j 1 j n k b l δ i m ] j 1 j n A ̂ i ( X ) .
(A12)
Now, using Eqs. (A7), one has
( ϕ ̃ ) 2 ϑ = m = 1 n = 1 β i 1 i m ; j 1 j n m ! n ! ( i 1 i m ϕ ( X ) ) ( j 1 j n ϕ ( X ) ) .
(A13)
Here, β i 1 i m ; j 1 j n is defined by
β i 1 i m ; j 1 j n α a i 1 i m j 1 j n α a i 1 i m α a j 1 j n ,
(A14)
which satisfies β i 1 i m ; j 1 j n = 0 for odd ( m + n ). Substituting Eq. (A13) into Eq. (A4), the part which is quadratically dependent on { i 1 i m 1 ( E L ) i m } is derived as
Ψ E 2 a ( X , μ , t ) = 1 2 m = 1 n = 1 C E 2 a i 1 i m ; j 1 j n i 1 i m 1 ( E L ) i m j 1 j n 1 ( E L ) j n ,
(A15)
where C E 2 a i 1 i m ; j 1 j n is given by
C E 2 a i 1 i m ; j 1 j n e a 2 μ B ( m + n ) β i 1 i m ; j 1 j n m ! n ! .
(A16)
It is also found from Eqs. (A5) and (A6) that the remaining potential functions Ψ E A ̂ a and Ψ A ̂ 2 a take bilinear and quadratic forms which are given by
Ψ E A ̂ a ( X , μ , t ) = m = 1 n = 1 C E A ̂ a i 1 i m ; j 1 j n i 1 i m 1 ( E L ) i m j 1 j n 1 A ̂ j n ,
(A17)
and
Ψ A ̂ 2 a ( X , μ , t ) = 1 2 m = 1 n = 1 C A ̂ 2 a i 1 i m ; j 1 j n i 1 i m 1 A ̂ i m j 1 j n 1 A ̂ j n ,
(A18)
respectively. One can use Eq. (A9) to derive the coefficients C E A ̂ a i 1 i m ; j 1 j n and C A ̂ 2 a i 1 i m ; j 1 j n in Eqs. (A17) and (A18) as functions of α a j 1 j n. The expressions given in Eqs. (A8), (A12), (A15), (A17), and (A18) are valid in the Cartesian spatial coordinates although they can be easily transformed into those in general, spatial coordinates as shown in Sec. III.
The gyrokinetic Lagrangian given by Eq. (18) in Sec. III is written as L G K d 3 X L G K where the gyrokinetic Lagrangian density L G K is defined as a function of ( X , t ) by L G K a d 3 v F a L GYa. The part of L G K including the potential field Ψ a is represented by
L Ψ a L Ψ a = a d 3 v F a e a Ψ a = ρ c ( g ) ( X , t ) ϕ ( X , t ) + L E 1 + L A ̂ 1 + L E 2 + L E A ̂ + L A ̂ 2 ,
(B1)
where Eq. (A1) is used. Equation (B1) describes the electromagnetic interaction of charged particles and is used to derive gyrokinetic expressions for polarization and magnetization as shown in Appendixes C and D. In Eq. (B1), the gyrocenter charge density ρ c ( g ) ( X , t ) is given by ρ c ( g ) ( X , t ) a e a N a ( g ) ( X , t ) a e a d 3 v F a ( X , U , μ , t ) and the other components of the Lagrangian density are defined by
[ L E 1 , L A ̂ 1 , L E 2 , L E A ̂ , L A ̂ 2 ] a [ L E 1 a , L A ̂ 1 a , L E 2 a , L E A ̂ a , L A ̂ 2 a ] a d 3 v F a e a [ Ψ E 1 a , Ψ A ̂ 1 a , Ψ E 2 a Ψ E A ̂ a , Ψ A ̂ 2 a ] .
(B2)
Here, using Eqs. (A8) and (B2), L E 1 a can be represented in the linear form of E L and its spatial derivatives,
L E 1 a = k = 1 Q 0 a j 1 j 2 k j 1 j 2 k 1 ( E L ) j 2 k ,
(B3)
where Q 0 a j 1 j 2 k represents the multipole moment of the electric charge distribution54 of species a induced by finite gyroradius,
Q 0 a j 1 j 2 k e a d 3 v F a α a j 1 j 2 k ( 2 k ) ! .
(B4)
Substituting Eq. (A12) into Eq. (B2) yields the linear form of A ̂ and its spatial derivatives,
L A ̂ 1 a = n = 1 R 0 a j 1 j n j 1 j n 1 A ̂ j n ,
(B5)
where
R 0 a j 1 j n i e a c d 3 v F a 1 n ! [ α a j 1 j n ( U b i + v B a i ) + Ω a ε klm α a j 1 j n k b l δ i m ] .
(B6)
Especially, in the cases of n = 0 and 1, Eq. (B6) is written as
R 0 a i e a c d 3 v F a ( U b i + v B a i ) = [ e a c N a ( g ) V a g b + b B × ( P a b · b + P a ln B ) ] i ,
(B7)
and
R 0 a j i e a c d 3 v F a Ω a ε klm α a j k b l δ i m = P a B ( I × b ) j i ,
(B8)
respectively, where [ P a , P a ] d 3 v F a [ m a U 2 , μ B ]. From Eqs. (A15), (A17), (A18), and (B2), one obtains the quadratic forms,
L E 2 a = 1 2 m = 1 n = 1 χ E a i 1 i m ; j 1 j n i 1 i m 1 ( E L ) i m j 1 j n 1 ( E L ) j n ,
(B9)
L E A ̂ a = n = 1 Q A ̂ a j 1 j n j 1 j n 1 ( E L ) j n = n = 1 R E a k 1 k n k 1 k n 1 A ̂ k n = m = 1 n = 1 χ E A ̂ a j 1 j m ; k 1 , k n j 1 j m 1 ( E L ) j m k 1 k n 1 A ̂ k n ,
(B10)
and
L A ̂ 2 a = 1 2 n = 1 R A ̂ a j 1 j n j 1 j n 1 A ̂ j n = 1 2 m = 1 n = 1 χ A ̂ a j 1 j m ; k 1 , k n ( j 1 j m 1 A ̂ j m ) ( k 1 k n 1 A ̂ k n ) .
(B11)
Here, χ E a i 1 i m ; j 1 j n, χ E A ̂ a i 1 i m ; j 1 j n, and χ A ̂ a i 1 i m ; j 1 j n are regarded as the generalized electromagnetic susceptibilities which are defined by
[ χ E a i 1 i m ; j 1 j n , χ E A ̂ a i 1 i m ; j 1 j n , χ A ̂ a i 1 i m ; j 1 j n ] e a d 3 v F a [ C E 2 a i 1 i m ; j 1 j n , C E A ̂ a i 1 i m ; j 1 j n , C A ̂ 2 a i 1 i m ; j 1 j n ] .
(B12)
In this appendix, it is shown in detail that the charge density of particles consists of the charge density of gyrocenters and other components including multipole moments which appear due to finite gyroradii of charged particles and turbulent electromagnetic fields. The charge density ρc which appears in Poisson's equation, Eq. (33), is given by the variational derivative of the gyrokinetic Lagrangian LGK with respect to the electrostatic potential ϕ as shown in Eq. (34) where the delta function δ 3 ( X + ρ a x ) is used. Here, to obtain another expression of ρc, the variational derivative is written as
ρ c = δ L G K δ ϕ = n = 0 ( 1 ) n j 1 j n ( L Ψ ( j 1 j n ϕ ) ) = L Ψ ϕ + n = 1 ( 1 ) n j 1 j n L Ψ ( j 1 j n 1 ( E L ) j n ) = ρ c ( g c ) · P G ,
(C1)
where ρ c ( g c ) and P G are the gyrocenter charge density and the polarization density vector, defined by
ρ c ( g c ) L Ψ ϕ a e a d 3 v F a a e a N a ( g ) ,
(C2)
and
P G δ L G K δ E L n = 0 ( 1 ) n j 1 j n L Ψ ( j 1 j n E L ) ,
(C3)
respectively. Then, the electric displacement field D is given by
D E + 4 π P G ,
(C4)
in terms of which the gyrokinetic Poisson equation is written as
· D = 4 π ρ c ( g c ) .
(C5)
From Eq. (C3), the ith component of the gyrokinetic polarization density vector is written as
P G i = n = 0 ( 1 ) n i 1 i n Q i i 1 i n ,
(C6)
where the multipole moments Q i i 1 i n ( n = 0 , 1 , 2 , ) are given using Eqs. (B1), (B2), (B3), (B9), and (B10) as
Q i i 1 i n = L Ψ ( i 1 i n ( E L ) i ) = a ( Q 0 a i i 1 i n + Q E a i i 1 i n + Q A ̂ a i i 1 i n ) .
(C7)
Here, Q 0 a i i 1 i n is defined in Eq. (B4). The other multipole moments Q E a i i 1 i m and Q A ̂ a i i 1 i n of the electric charge distribution of species a are written in terms of in the linear forms with respect to ( E L ) i , A ̂ i, and their spatial derivatives as
Q E a i i 1 i m n = 1 χ E a i i 1 i m ; j 1 j n j 1 j n 1 ( E L ) j n ,
(C8)
and
Q A ̂ a i i 1 i n = m = 1 χ E A ̂ a i i 1 i n ; k 1 k m k 1 k m 1 A ̂ k m ,
(C9)
respectively, where χ E a i i 1 i n ; j j 1 j m and χ E A ̂ a i i 1 i n ; k 1 k m are defined in Eq. (B12). It is found from Eqs. (A14), (A16), (B4), and (B12) that Q 0 a i 1 i m and Q E a i 1 i m are both symmetric with respect to arbitrary permutations of the indices i 1 , , i m because α a i 1 i m defined in Eq. (A9) has the same symmetry.
When retaining only the n = 0 term in Eq. (C9) and using the lowest order distribution function F a 0 given by the local Maxwellian, the polarization density vector is approximated by
P G i a Q E a i n n a 0 m a c 2 B 2 E L = c 2 4 π v A 2 E L ,
(C10)
where n a 0 d 3 v F a 0 and v A B 2 / ( 4 π a n a 0 m a ) represent the equilibrium density and the Alfvén velocity, respectively. Equation (C10) presents a well-known expression of polarization. It should be noted that Eq. (C1) with Eq. (C3) including all multipole moments gives the charge density which is equivalent to that presented in Eq. (34). Then, as shown in Ref. 31, the gyrokinetic Poisson equation given by the classical gyrokinetic theory1–3 based on the WKB formalism can be derived as well from the turbulent part of Eq. (33) using the charge density given by Eq. (34) [or Eq. (C1) with Eqs. (C2) and (C3)].
It is remarked here that the gyrocenter charge density ρ c ( g c ) defined in Eq. (C2) satisfies
ρ c ( g c ) t + · j ( g c ) = 0 ,
(C11)
where the gyrocenter current density j ( g c ) is given by
j ( g c ) a e a Γ a ( g c ) a e a d 3 v F a u a x .
(C12)
Equation (C11) is derived by taking the velocity-space integral and the species summation of Eq. (102) with the help of Eq. (103). Combining Eqs. (C5) and (C11), one obtains
· ( D t + 4 π j ( g c ) ) = 0 ,
(C13)
from which one can write
4 π j L ( g c ) = D L t = t ( ϕ D ) .
(C14)
In Eq. (C14), the subscript L denotes the longitudinal (or irrotational) vector part, and the potential ϕ D for D L is defined such that D L = ϕ D.
In a similar way to  Appendix C, this appendix shows how the current density of particles is expressed by the sum of the current density of gyrocenters and other components induced by finite gyroradii of charged particles and turbulent electromagnetic fields. The gyrokinetic Ampère's law is presented in Eq. (38) which contains the transverse part j T of the current density j given by the variational derivative of the gyrokinetic Lagrangian LGK with respect to the perturbed vector potential A ̂ as shown in Eq. (36). Here, another expression of j is obtained by writing the variational derivative as
1 c j = δ L G K δ A ̂ = n = 0 ( 1 ) n j 1 j n ( L Ψ ( j 1 j n A ̂ ) ) = a e a d 3 v F a n = 0 ( 1 ) n + 1 j 1 j n ( Ψ a ( j 1 j n A ̂ ) ) = a e a c Γ a .
(D1)
The particle flux Γ a of species a in Eq. (D1) is written as
Γ a = n = 0 Γ a ( n ) ,
(D2)
where Γ a ( n ) is defined by
Γ a ( n ) = c d 3 v F a ( 1 ) n + 1 j 1 j n ( Ψ a ( j 1 j n A ̂ ) ) .
(D3)
The zeroth-order flux Γ a ( 0 ) is written as
Γ a ( 0 ) c d 3 v F a ( Ψ a A ̂ ) = d 6 Z δ 3 ( X x ) [ F a ( Z , t ) ( v e a m a c A ̂ + v B a ) + e a ψ ̃ a B F a μ v ] = d 3 v F a [ ( U e a m a c A ̂ ) b + v B a + c B b × ψ a ϑ ] + O ( δ 2 ) ,
(D4)
which is equivalent to Γ a ( g c ) d 3 v F a u a x to the lowest order in δ.
Using Eqs. (B1), (B5), (B10), and (B11), one can represent the derivatives L Ψ / ( j 1 j n A ̂ k ) by
R j 1 j n k L Ψ ( j 1 j n A ̂ k ) = a e a d 3 v F a Ψ a ( j 1 j n A ̂ k ) = a ( R 0 a j 1 j n k + R E a j 1 j n k + R A ̂ a j 1 j n k ) ,
(D5)
where
R E a j 1 j n k = m = 0 χ E A ̂ a i 1 i m i ; j 1 j n k i 1 i m ( E L ) i ,
(D6)
R A ̂ a j 1 j n k = m = 0 χ A ̂ a j 1 j n k ; l 1 l m l l 1 l m A ̂ l ,
(D7)
and R 0 a j 1 j n k is defined by Eq. (B6). The coefficients χ E A ̂ a i 1 i m i ; j 1 j n k and χ A ̂ a j 1 j n k ; l 1 l m l are given by Eq. (B12).
Now, it is found from Eqs. (D1)–(D4) that the lth component of j can be expressed as
j l = ( j ( 0 ) ) l + c k N k l ,
(D8)
where
j ( 0 ) c L Ψ A ̂ = c a e a d 3 v F a Ψ a A ̂ = a e a Γ a ( 0 )
(D9)
is regarded as the current of gyrocenters and
N k l n = 0 ( 1 ) n + 1 j 1 j n R j 1 j n k l .
(D10)
Here, it should be recalled that j ( 0 ) a e a Γ a ( 0 ) defined above equals j ( g c ) a e a Γ a ( g c ) a e a d 3 v F a u a x to the lowest order in δ although the equality does not rigorously holds. When assuming | ρ a · | < 1 and retaining only the lowest order of Nkl in the expansion with respect to | ρ a · |, one just has the nonturbulent contribution to N k l as
N k l R k l a R 0 a k l .
(D11)
Then, using Eqs. (B8) and (D11) leads to
c k N k l a [ × ( c P a B b ) ] l ,
(D12)
where P a is defined after Eq. (B8). Equation (D12) is a well-known expression of a magnetization current [see Ref. 69]. Then, from using Eqs. (D4), (D8), (D9), and (D12) with the lowest order distribution function F a 0 given by the local Maxwellian, the perpendicular component of the equilibrium current can be derived as
j = a e a d 3 v F a 0 ( v B a + c B b × ϕ ens ) [ × ( c P 0 B b ) ] = c B b × P 0 ,
(D13)
where P 0 a P a 0 a d 3 v F a 0 μ B denotes the equilibrium pressure and a e a d 3 v F a 0 = 0 is used. Equation (D13) presents the magnetization law69 and one can see that, as pointed out in Ref. 31, the diamagnetic current consistent with the MHD equilibrium c 1 j × B = P 0 is correctly derived from the variational formulation with the v B a term retained in the potential part of the Hamiltonian given by Eq. (2) with Eq. (3). It is also shown in Ref. 31 that the gyrokinetic Ampère's law given by the classical gyrokinetic theory1–3 based on the WKB formalism can be derived from the turbulent part of Eq. (38) using Eq. (36), which is equivalent to Eq. (D1) with all the gyroradius expansion terms retained.
Here, j ( 0 ) defined in Eq. (D9) is divided into the longitudinal and transverse parts as j ( 0 ) = j L ( 0 ) + j T ( 0 ), where the longitudinal part can be written in terms of the scalar function λ ( 0 ) as j L ( 0 ) = 1 4 π λ ( 0 ). Similarly, N k l in Eq. (D10) is represented by the sum of two parts,
N k l = N L k l + N T k l ,
(D14)
where N L k l and N T k l are defined such that ε mnl n N L k l = 0 and and l N T k l = 0 are satisfied. Then, there exist Vk and W n k in terms of which N L k l and N T k l are given by
N L k l = l V k , N T k l = ε lmn m W n k .
(D15)
Using Eq. (37) and j L ( 0 ) 1 4 π λ ( 0 ), one obtains
λ 4 π = λ ( 0 ) 4 π + c · V
(D16)
and
j T = j T ( 0 ) + c × M W ,
(D17)
where V is the vector with the components Vk (k = 1, 2, 3), and the jth component of the vector M W is defined by
( M W ) j = k W j k .
(D18)
It is found from Eq. (D17) that the magnetization field can be represented by M W up to a gradient of an arbitrary scalar function. Using M W, the magnetic intensity field H is defined by
H = B + B ̂ 4 π M W ,
(D19)
which is distinguished from H # given in Eq. (147). Then, using Eq. (D19), the gyrokinetic Ampère's law, Eq. (38) can be rewritten as
× H = 4 π c j T ( 0 ) .
(D20)
This appendix presents energy balance equations in electromagnetic gyrokinetic turbulence. Here, the Cartesian coordinate system is used and three-dimensional vectors are written in terms of boldface letters. The energy of a single charged particle of species a is denoted by E a which is equal to the gyrocenter Hamiltonian HGYa in Eq. (2) and written as
E a 1 2 m a U 2 + μ B + e a Ψ a H GYa = L GYa u a Z · u a Z L GYa .
(E1)
It can be shown from Eqs. (28), (32), and (E1) that the total derivative of HGYa is written as
E ̇ a ( d d t ) a H GYa ( t + u a Z · Z ) H GYa = ( L GYa t ) u = e Ψ a t + μ B t e a c u a x · A a * t ,
(E2)
where ( L GYa / t ) u denotes the time derivatives of LGYa with u a Z kept fixed in LGYa. Multiplying Eq. (E2) with Fa and taking its velocity-space integral, the local energy balance equation for the system of the single particle species is obtained as
t ( d 3 v F a E a ) + · ( d 3 v F a E a u a x ) = d 3 v ( F a E ̇ a + K a E a ) ,
(E3)
where the gyrokinetic Boltzmann equation shown in Eq. (102) is used and u a x = ( d / d t ) a X represents the gyrocenter velocity defined at the right-hand side of Eq. (11).
Next, the energy balance in the whole system including particles of all species and the turbulent electromagnetic fields is considered. From Eq. (40), one finds
a d 3 v F a L GYa t + L F t = J { L GKF ( J ϕ ) ( J ϕ ) t + L GKF ( J A ) · ( J A ) t + L GKF ( J A ̂ ) · ( J A ̂ ) t } + L GKF λ λ t .
(E4)
Then, taking the species summation of Eq. (E3) and using Eqs. (E1), (E2), and (E4), they yield
t ( a d 3 v F a H GYa L F ) + · ( a d 3 v F a H GYa u a x ) + J { L GKF ( J ϕ ) ( J ϕ ) t + L GKF ( J A ) · ( J A ) t + L GKF ( J A ̂ ) · ( J A ̂ ) t } + L GKF λ λ t = a d 3 v K a H GYa .
(E5)
Now, the variational equations, δ L G K / δ ϕ = 0, δ L G K / δ A ̂ = 0, and δ L G K / δ λ = 0, which are equivalent to the gyrokinetic Poisson equation, Ampère's law, and the Coulomb gauge condition, respectively, are used to rewrite Eq. (E5) as
t ( a d 3 v F a H GYa L F ) + · ( a d 3 v F a H GYa u a x ) + n = 1 k = 1 n ( 1 ) k 1 j k { j 1 j k 1 ( L GKF ( j 1 j n ϕ ) ) ( j k + 1 j n ϕ ) t + j 1 j k 1 ( L GKF ( j 1 j n A ̂ l ) ) ( j k + 1 j n A ̂ l ) t } + · [ c 4 π E T × ( B + B ̂ 4 π M # ) ] + i [ Ξ j i B j t ] = [ j # c 4 π × ( B + B ̂ ) ] · E T + a d 3 v K a H GYa ,
(E6)
where Ξ j i on the left-hand side is defined by
Ξ j i = L G K ( i B j ) = a d 3 v F a L GYa ( i B j ) .
(E7)
Equation (E6) can be further deformed to obtain the local energy balance equation of the whole system as
t [ a d 3 v F a ( 1 2 m a | v e a m a c A ̂ | 2 e a c v B a · A ̂ e a 2 2 c 2 B μ [ ( v · A ̂ ̃ ) 2 ] ) + | E L | 2 8 π + E L · P D 2 + 1 2 n = 1 ( j 1 j n ( E L ) i ) Q E i j 1 j n + E T · D L 4 π + | B + B ̂ | 2 8 π ] + · [ a d 3 v F a { 1 2 m a U 2 + μ B + e a ( Ψ a ϕ ( x ) ) } u a x + c 4 π E × H # + 1 4 π ϕ D E T t + ϕ ( P G ) T t ] + i [ Ξ j i B j t ] + i [ n = 1 k = 0 n 1 ( 1 ) n k { ( j n k j n 1 ( E L ) j n ) j 1 j n k 1 ( Q i j 1 j n t ) + c ( j n k j n ( E ̂ T ) l ) j 1 j n k 1 R i j 1 j n l } ] + · [ c 4 π E ̂ T × ( B + B ̂ ) λ 4 π E ̂ T ] + i [ c N i j ( E ̂ T ) j ] = [ ( j # ) T c 4 π × ( B + B ̂ ) ] · ( E T + E L ) + a d 3 v K a { 1 2 m a U 2 + μ B + e a ( Ψ a ϕ ( x ) ) } ,
(E8)
where H # is defined in Eq. (147). The rate of change in the sum of the kinetic and electromagnetic energy densities is described by Eq. (E8). There appear the effects of polarization and magnetization including all multipole moments. The energy flux on the left-hand side of Eq. (E8) contains the kinetic energy flow due to the gyrocenter motion, the Poynting vector, and the extra contributions due to the electromagnetic microturbulence. The last terms on the left-hand side of Eq. (E8) can be deformed into
· [ c 4 π E ̂ T × ( B + B ̂ ) λ 4 π E ̂ T ] + i [ c N i j ( E ̂ T ) j ] = · [ c 4 π E ̂ T × H λ ( 0 ) 4 π E ̂ T + c V · E ̂ T ] + i [ c ε ijk W j l l ( E ̂ T ) k + W j i B ̂ j t ] ,
(E9)
where H is given in Eq. (D19). It is seen from Eqs. (E8) and (E9) that the energy balance equation two types of the Poynting vector, ( c / 4 π ) E × H # and ( c / 4 π ) E ̂ T × H, where E E L + E T, E L = ϕ, E T = c 1 A / t and E ̂ T = c 1 A ̂ / t. On the right-hand side of Eqs. (E6) and (E8), the effects of collisions and/or external sources are represented by terms including K a. It is noted that these terms can be written as the divergence of classical energy transport flux when K a is given by the collision operator including the finite gyroradius effect. In addition, terms including E L ϕ on the right-hand side of Eq. (E8) are given in the divergence form as
[ ( j # ) T c 4 π × ( B + B ̂ ) ] · E L = · [ c 4 π ϕ × ( B # B B ̂ ) ] .
(E10)
Therefore, in the stationary background magnetic field where E T c 1 A / t = 0 and with no external energy sources, Eqs. (E6) and (E8) take the conservation form. Furthermore, even when E T c 1 A / t 0, the ensemble average of Eq. (E6) and (E8) can take the form of total energy conservation on macroscopic spatiotemporal scales in the background field determined by the condition in Eq. (145). It can also be confirmed from comparison with Eq. (22) in Ref. 45 that the kinetic and electromagnetic energies, the kinetic energy flux, the Poynting vector, and the longitudinal and transverse electric fields in the energy conservation equation of the Vlasov–Darwin system are retained in Eq. (E8). There, additional terms due to finite-gyroradius effects and electromagnetic microturbulence are included as well.
1.
T. M.
Antonsen
, Jr.
and
B.
Lane
,
Phys. Fluids
23
,
1205
(
1980
).
2.
P. J.
Catto
,
W. M.
Tang
, and
D. E.
Baldwin
,
Plasma Phys.
23
,
639
(
1981
).
3.
E. A.
Frieman
and
L.
Chen
,
Phys. Fluids
25
,
502
(
1982
).
4.
A. J.
Brizard
and
T. S.
Hahm
,
Rev. Mod. Phys.
79
,
421
(
2007
).
5.
H.
Sugama
,
Phys. Plasmas
7
,
466
(
2000
).
6.
A. A.
Schekochihin
,
S. C.
Cowley
,
W.
Dorland
,
G. W.
Hammett
,
G. G.
Howes
,
E.
Quataert
, and
T.
Tatsuno
,
Astrophys. J., Suppl.
182
,
310
(
2009
).
7.
J. A.
Krommes
,
Ann. Rev. Fluid Mech.
44
,
175
(
2012
).
8.
W.
Horton
,
Turbulent Transport in Magnetized Plasmas
, 2nd ed. (
World Scientific
,
Singapore
,
2018
), Chap. 12.
9.
H.
Sugama
,
M.
Okamoto
,
W.
Horton
, and
M.
Wakatani
,
Phys. Plasmas
3
,
2379
(
1996
).
10.
R. D.
Hazeltine
and
J. D.
Meiss
,
Plasma Confinement
(
Addison-Wesley
,
Redwood City, CA
,
1992
), Chap. 7.10.
11.
A. M.
Dimits
,
G.
Bateman
,
M. A.
Beer
,
B. I.
Cohen
,
W.
Dorland
,
G. W.
Hammett
,
C.
Kim
,
J. E.
Kinsey
,
M.
Kotschenreuther
,
A. H.
Kritz
,
L. L.
Lao
,
J.
Mandrekas
,
W. M.
Nevins
,
S. E.
Parker
,
A. J.
Redd
,
D. E.
Shumaker
,
R.
Sydora
, and
J.
Weiland
,
Phys. Plasmas
7
,
969
(
2000
).
12.
F.
Jenko
,
W.
Dorland
,
M.
Kotschenreuther
, and
B. N.
Rogers
,
Phys. Plasmas
7
,
1904
(
2000
).
13.
J.
Candy
and
R. E.
Waltz
,
J. Comp. Phys.
186
,
545
(
2003
).
14.
T.-H.
Watanabe
and
H.
Sugama
,
Nucl. Fusion
46
,
24
(
2006
).
15.
A. G.
Peeters
,
Y.
Camenen
,
F. J.
Casson
,
W. A.
Hornsby
,
A. P.
Snodin
,
D.
Strintzi
, and
G.
Szepesi
,
Comp. Phys. Comm.
180
,
2650
(
2009
).
16.
R. G.
Littlejohn
,
J. Math. Phys.
23
,
742
(
1982
).
17.
B.
Scott
and
J.
Smirnov
,
Phys. Plasmas
17
,
112302
(
2010
).
18.
A. J.
Brizard
and
N.
Tronko
,
Phys. Plasmas
18
,
082307
(
2011
).
19.
F. I.
Parra
and
I.
Calvo
,
Plasma Phys. Control. Fusion
53
,
045001
(
2011
).
20.
H.
Qin
,
J. W.
Burby
, and
R. C.
Davidson
,
Phys. Rev. E
90
,
043102
(
2014
).
21.
H.
Sugama
,
M.
Nunami
,
M.
Nakata
, and
T.-H.
Watanabe
,
Phys. Plasmas
24
,
020701
(
2017
);
H.
Sugama
,
Rev. Mod. Plasma Phys.
1
,
9
(
2017
).
22.
P.
Fan
,
H.
Qin
, and
J.
Xiao
,
Plasma Sci. Technol.
23
,
105103
(
2021
).
23.
Z.
Lin
,
T. S.
Hahm
,
W. W.
Lee
,
W. M.
Tang
, and
R. B.
White
,
Science
281
,
1835
(
1998
).
24.
Y.
Idomura
,
Phys. Plasmas
24
,
080701
(
2017
).
25.
S.
Ku
,
C. S.
Chang
, and
P. H.
Diamond
,
Nucl. Fusion
49
,
115021
(
2009
).
26.
W. X.
Wang
,
T. S.
Hahm
,
S.
Ethier
,
G.
Rewoldt
,
W. W.
Lee
,
W. M.
Tang
,
S. M.
Kaye
, and
P. H.
Diamond
,
Phys. Rev. Lett.
102
,
035005
(
2009
).
27.
V.
Grandgirard
,
J.
Abiteboul
,
J.
Bigot
,
T.
Cartier-Michaud
,
N.
Crouseilles
,
G.
Dif-Pradalier
,
C.
Ehrlacher
,
D.
Esteve
,
X.
Garbet
,
P.
Ghendrih
,
G.
Latu
,
M.
Mehrenberger
,
C.
Norscini
,
C.
Passeron
,
F.
Rozar
,
Y.
Sarazin
,
E.
Sonnendrücker
,
A.
Strugarek
, and
D.
Zarzoso
,
Comput. Phys. Commun.
207
,
35
(
2016
).
28.
A.
Bottino
and
E.
Sonnendrücker
,
J. Plasma Phys.
81
,
435810501
(
2015
).
29.
J. A.
Heikkinen
,
S. J.
Janhunen
,
T. P.
Kiviniemi
, and
F.
Ogando
,
J. Comput. Phys.
227
,
5582
(
2008
).
30.
N. R.
Mandell
,
A.
Hakim
,
G. W.
Hammett
, and
M.
Francisquez
,
J. Plasma Phys.
86
,
905860109
(
2020
).
31.
H.
Sugama
,
S.
Matsuoka
, and
M.
Nunami
,
Phys. Plasmas
29
,
052509
(
2022
).
32.
K. H.
Burrell
,
Phys. Plasmas
27
,
060501
(
2020
).
33.
A. H.
Boozer
,
Phys. Fluids
26
,
496
(
1983
).
34.
J.
Nührenberg
and
R.
Zille
,
Phys. Lett.
129
,
113
(
1988
).
35.
D. A.
Spong
,
Phys. Plasmas
12
,
056114
(
2005
).
36.
37.
H.
Sugama
,
S.
Matsuoka
,
M.
Nunami
, and
S.
Satake
,
Phys. Plasmas
28
,
022312
(
2021
).
38.
W. A.
Newcomb
,
Nucl. Fusion Suppl. Pt
2
,
451
(
1962
).
39.
A. J.
Brizard
and
C.
Tronci
,
Phys. Plasmas
23
,
062107
(
2016
).
40.
E.
Hirvijoki
,
J. W.
Burby
,
D.
Pfefferlé
, and
A. J.
Brizard
,
J. Phys. A: Math. Theor.
53
,
235204
(
2020
).
41.
H.
Cendra
,
D. D.
Holm
,
M. J. W.
Hoyle
, and
J. E.
Marsden
,
J. Math. Phys.
39
,
3138
(
1998
).
42.
J. E.
Marsden
and
T. S.
Ratiu
,
Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems
, Texts in Applied Mathematics (
Springer
,
1999
), Chap. 13.5.
43.
J.
Squire
,
H.
Qin
,
W. M.
Tang
, and
C.
Chandre
,
Phys. Plasmas
20
,
022501
(
2013
).
44.
H.
Sugama
,
M.
Nunami
,
S.
Satake
, and
T.-H.
Watanabe
,
Phys. Plasmas
25
,
102506
(
2018
).
45.
H.
Sugama
,
T.-H.
Watanabe
, and
M.
Nunami
,
Phys. Plasmas
20
,
024503
(
2013
).
46.
A. J.
Brizard
,
J. Plasma Phys.
87
,
905870307
(
2021
).
47.
L. D.
Landau
and
E. M.
Lifshitz
,
The Classical Theory of Fields
(
Pergamon
,
New York
,
1975
), Chap. 11.
48.
F. L.
Hinton
and
R. D.
Hazeltine
,
Rev. Mod. Phys.
48
,
239
(
1976
).
49.
S. P.
Hirshman
and
D. J.
Sigmar
,
Nucl. Fusion
21
,
1079
(
1981
).
50.
P.
Helander
and
D. J.
Sigmar
,
Collisional Transport in Magnetized Plasmas
(
Cambridge University Press
,
Cambridge
,
2002
), Chap. 8.
51.
H.
Sugama
,
T.-H.
Watanabe
,
M.
Nunami
, and
S.
Nishimura
,
Plasma Phys. Control. Fusion
53
,
024004
(
2011
).
52.
F. I.
Parra
,
M.
Barnes
,
I.
Calvo
, and
P. J.
Catto
,
Phys. Plasmas
19
,
056116
(
2012
).
53.
I.
Calvo
and
F. I.
Parra
,
Plasma Phys. Control. Fusion
57
,
075006
(
2015
).
54.
J. D.
Jackson
,
Classical Electrodynamics
, 3rd ed. (
Wiley
,
New York
,
1998
), Chap. 6.
55.
H.
Goldstein
,
C.
Poole
, and
J.
Safko
,
Classical Mechanics
, 3rd ed. (
Addison-Wesley
,
San Francisco
, CA,
2002
), p.
593
.
56.
J. E.
Marsden
and
T. S.
Ratiu
,
Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems
, Texts in Applied Mathematics (
Springer
,
1999
), Chap. 4.3.
57.
H.
Sugama
,
T.-H.
Watanabe
, and
M.
Nunami
,
Phys. Plasmas
22
,
082306
(
2015
).
58.
A. J.
Brizard
,
Phys. Plasmas
11
,
4429
(
2004
).
59.
J. W.
Burby
,
A. J.
Brizard
, and
H.
Qin
,
Phys. Plasmas
22
,
100707
(
2015
).
61.
A. J.
Brizard
,
Phys. Plasmas
30
,
102106
(
2023
).
62.
H.
Sugama
,
T.-H.
Watanabe
, and
M.
Nunami
,
Phys. Plasmas
21
,
012515
(
2014
).
63.
M.
Wakatani
,
Stellarator and Heliotron Devices
(
Oxford University Press
,
New York
,
1998
), Chap. 7.
64.
A. H.
Boozer
,
Phys. Fluids
26
,
1288
(
1983
).
65.
H.
Sugama
and
S.
Nishimura
,
Phys. Plasmas
9
,
4637
(
2002
).
66.
I. G.
Abel
,
G. G.
Plunk
,
E.
Wang
,
M.
Barnes
,
S. C.
Cowley
,
W.
Dorland
, and
A. A.
Schekochihin
,
Rep. Prog. Phys.
76
,
116201
(
2013
).
67.
H.
Sugama
and
W.
Horton
,
Phys. Plasmas
5
,
2560
(
1998
).
68.
F. I.
Parra
,
M.
Barnes
, and
A. G.
Peeters
,
Phys. Plasmas
18
,
062501
(
2011
).
69.
R. D.
Hazeltine
and
J. D.
Meiss
,
Plasma Confinement
(
Addison-Wesley
,
Redwood City, CA
,
1992
), Chap. 4.5.