The Eulerian variational formulation of the gyrokinetic system in general spatial coordinates

The Eulerian variational formulation of the gyrokinetic system with electrostatic turbulence is presented in general spatial coordinates by extending our previous work [H. Sugama, {\it et al}., Phys.\ Plasmas {\bf 25}, 102506 (2018)]. The invariance of the Lagrangian of the system under an arbitrary spatial coordinate transformation is used to derive the local momentum balance equation satisfied by the gyrocenter distribution functions and the turbulent potential which are given as solutions of the governing equations. This derivation is in contrast with the conventional method using the spatial translation in which the asymmetric canonical pressure tensor generally enters the momentum balance equation. In the present study, the variation of the Lagrangian density with respect to the metric tensor is taken to directly obtain the symmetric pressure tensor which includes the effect of turbulence on the momentum transport. In addition, it is shown in this work how the momentum balance is modified when the collision and/or external source terms are added to the gyrokinetic equation. The results obtained here are considered useful for global gyrokinetic simulations investigating both neoclassical and turbulent transport processes even in general non-axisymmetric toroidal systems.


I. INTRODUCTION
Gyrokinetics [1][2][3][4][5][6][7] has been used for several decades as a basic theoretical framework to study microinstabilities and turbulent processes in magnetized plasmas. 8 significant number of large-scale simulations are being performed based on gyrokinetic equations to analyze and predict turbulent transport fluxes of particles, heat, and momentum. 9,10The gyrokinetic equations derived from the Lagrangian and/or Hamiltonian describing the gyrocenter motion possess conservation properties 4,5 which are suitable for long-time and global transport simulations.2][13] The variational formulations are also applied to derive the useful conservative numerical schemes [14][15][16][17][18] in plasma physics for solving the guiding center equations and the ideal magnetohydrodynamics equations as well as the Vlasov-Poisson and Vlasov-Maxwell equations.
5][26][27][28][29] In our previous work, 30 the Eulerian formulation, 31 which is also called the Euler-Poincaré reduc-tion procedure, 11,13,[32][33][34] is applied to derive the governing equations of the Vlasov-Poisson-Ampère (or Vlasov-Darwin) system and those of the drift kinetic system in the general spatial coordinates, and the momentum balance equations for these systems are obtained by using the invariance of the action integrals for the systems under arbitrary spatial coordinate transformations.In this paper, the previous work is extended to present the Eulerian variational formulation of the gyrokinetic system with electrostatic turbulence in general spatial coordinates and to derive the local momentum balance equation with the symmetric pressure tensor including the effects of electrostatic turbulent fluctuations on the momentum transport.
In the present work, the governing equations of the gyrokinetic system are represented in general spatial coordinates so that they are useful for application to systems with complex geometries such as stellarator and heliotron plasmas 35 in which it is convenient to employ flux surface coordinates (e.g., Hamada coordinates 36 and Boozer coordinates 37 ) to express these equations for analytic and numerical studies.In our formulation, the momentum balance equation is derived using the invariance of the Lagrangian under arbitrary spatial coordinate transformations, which is analogous to the derivation of energymomentum conservation laws from the invariance of the action integral under arbitrary transformations of spatiotemporal coordinates in the theory of general relativity. 38In the same way as in the previous work, 30 the symmetric pressure tensor entering our momentum balance equation is directly obtained by taking the variational derivative of the Lagrangian density with respect to the metric tensor, which is in contrast to the conventional technique using the spatial translation transformation for the derivation of the canonical momentum balance including the asymmetric pressure tensor in the presence of the magnetic field.The canonical pressure tensor is asymmetric and gauge-dependent because of the vector potential included in the canonical momentum.In order to obtain the symmetric pressure tensor from the asymmetric canonical pressure tensor, additional complicated procedures of the Belinfante-Rosenfeld type using the angular momentum conservation law derived from the rotational symmetry are required. 39,40On the other hand, our method can more directly derive the symmetric pressure tensor which is clearly shown in this work to describe both neoclassical 41,42 and turbulent transport of the momentum in the gyrokinetic system.
We should note here that while employing the Hamiltonian gyrocenter motion equations, an irreversible collision term must be included into the gyrokinetic equation in order to treat neoclassical and turbulent transport simultaneously. 43,44It is shown in our formulation how the momentum balance equations are modified when the collision and/or external source terms are added into the gyrokinetic equation.This is possible because the momentum balance equations are derived from the invariance of the Lagrangian under the general spatial coordinate transformations, with the help of the gyrocenter motion equations and the gyrokinetic Poisson's equation while we have freedom in choosing the governing equation for the gyrocenter distribution function.The momentum balance equations including the effects of collisions and external sources are useful for checking and analyzing results of global gyrokinetic simulations using the Lagrangian and Hamiltonian equations to investigate neoclassical and turbulent transport in plasmas with particle, momentum, and/or heat sources.
It is valuable to make comparisons of the present work with recent works 45,46 on similar subjects.In Ref. 45 , a constrained Eulerian variational principle presented by Brizard 47 is used to derive the governing equations and local energy-momentum conservation laws of the electromagnetic gyrokinetic Vlasov-Maxwell system, for which the gyrocenter motion equations are expressed only in terms of the perturbed electric and magnetic fields by including the perturbed fields in the Poisson-bracket structure.In Ref. 46 , the field theory presented by Qin et al. 48or particle-field systems on heterogenous manifolds is applied to derive local energy and momentum conservation laws in the electromagnetic gyrokinetic system.Except for dropping magnetic fluctuations and some of the second or higher order terms with respect to the small perturbation amplitude, the gyrokinetic system treated here is basically the same as seen in other earlier works 4,5 , and it contains full finite gyroradius effects due to electrostatic fluctuations with high wavenumbers which are not included in Ref. 46 .Also, in the same way as in our previous work 30 , the governing equations for the gyrokinetic system are derived in this paper based on the Eule-rian (or the Euler-Poincaré) formulation which is historically older than the methods employed in Refs. 45,46.3][34] .Here, the Eulerian formulation implies that, for the present gyrokinetic case, the phase-space velocity (or the temporal change rate of the gyrocenter coordinates in the phase space) of the gyrocenter is regarded as a field function of time and phase-space coordinates of the gyrocenter at the time when the gyrocenter passes through the considered point.In the formulation presented in Refs. 46,48, the gyrocenter (or particle) phase-space velocity is described by not the Eulerian but Lagrangian view point and is coupled with the Eulerian description of the electromagnetic fields.In Refs. 45,46, not general but isometric transformations such as spatial translation and rotation are considered to derive local conservation laws of canonical linear and angular momentum in collisionless systems.In the present work, recognizing that the Lagrangian and the governing equations of the gyrokinetic system can be expressed in the invariant form under arbitrary spatial coordinate transformations in the Eulerian framework, this invariance property is used to derive the local momentum balance equation even in the presence of collisions and external sources, and the derived balance equation is shown to give the local momentum conservation law in the direction of symmetry for collisionless systems without external sources.
The rest of this paper is organized as follows.In Sec.II, the Lagrangian for describing the single-particle gyrocenter motion is given, in which the electrostatic potential fluctuation with the wavelength of the order of the gyroradius is included.Next, in Sec.III, we use the general spatial coordinates and define the Lagrangian of the whole system including particles of all species and turbulent electrostatic fields to present the Eulerian variational principle, from which the collisionless gyrokinetic equations for the gyrocenter distribution functions and the gyrokinetic Poisson's equation for the electrostatic potential are derived.In Sec.IV, we make use of the invariance of the Lagrangians for the single-species and whole systems under arbitrary spatial coordinate transformations to derive the momentum balance equations for both systems while allowing the collision and/or external source terms to be included in the gyrokinetic equations.Then, the symmetric pressure tensors, which enter these momentum balance equations, are obtained from the variational derivatives of the Lagrangians with respect to the metric tensor components and they are verified to describe both neoclassical and turbulent transport of the momentum.Finally, conclusions are given in Sec.V.In addition, Appendix A is given to briefly explain covariant derivatives and Christoffel symbols, which are used in the general spatial coordinates.Variations in the functional forms of vector and tensor fields under the infinitesimal transformation of the spatial coordinates and variational derivative with respect to the metric tensor components are described in Appendices B and C, respectively.In Appendix D, the WKB representation 49 is used to express the turbulent part of the pressure tensor in the form comparable with the one obtained in the past work on the momentum transport.Another derivation of the momentum balance equations, in which the asymmetric canonical pressure tensors appear, is shown in Appendix E and the energy balance equations are derived in Appendix F. The case of the symmetric background magnetic field is considered in Appendix G where it is shown how the local conservation law of the canonical momentum in the direction of symmetry is derived in the present formulation.

II. THE GYROCENTER LAGRANGIAN
The Lagrangian for describing the gyrocenter motion of the particle with mass m and charge e is given by 4 where the gyrocenter phase-space coordinates X, v , µ ≡ mv 2 ⊥ /(2B), and ϑ denote the gyrocenter position, the velocity component parallel to the magnetic field, the magnetic moment, and the gyrophase angle, respectively, and ˙≡ d/dt represents the time derivative along the trajectory in the phase space.The vector potential and the unit vector parallel to the magnetic field B are written by A and b ≡ B/B, respectively.Here, we suppose 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 −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 defined by and its electrostatic fluctuation part is written as where the gyroradius vector is denoted by ρ ≡ b × v/Ω, v is the particle's velocity, Ω ≡ eB/(mc) is the gyrofrequency, and the average of the electrostatic potential φ at the particle position X + ρ over the gyrophase ϑ is given by and the gyrophase-dependent part is denoted by The last term on the right-hand side of Eq. ( 3) is of the second order in the gyrokinetic ordering parameter ǫ where ǫ ∼ e φ/(m|v| 2 ) ∼ ρ/L with L representing the equilibrium gradient scale length is assumed.This second-order term is retained because it is necessary for deriving the gyrokinetic Poisson's equation correctly including the polarization effect as shown in Sec.III.C.However, other second-order terms shown in Ref. 50 are neglected in the gyrocenter Hamiltonian given by Eqs. ( 2) and (3).Therefore, rigorously speaking, the accuracy of the present model is up to the first order.The turbulent fluctuations are assumed to have the characteristic wavelength ∼ (k ⊥ ) −1 ∼ ρ.Then, the fluctuation terms in Eq.( 3) are considered to contain all-order terms in k ⊥ ρ(∼ 1) even though small amplitude terms of higher order in ǫ ∼ ρ/L(≪ 1) are neglected.In Eqs. ( 4) and ( 5) as well as in the equations shown below, the time variable t on which φ depends is omitted for simplicity.
In this section, the Cartesian spatial coordinates are used and three-dimensional vectors are represented in terms of boldface letters.Then, the electrostatic potential φ(X + ρ) is Taylor expanded about the gyrocenter position X as where X j and ρ j (j = 1, 2, 3) and the Cartesian spatial coordinates of the gyrocenter position vector X and the gyroradius vector ρ, respectively.Substituting Eq. ( 6) into Eq.( 4), we obtain where the gyrophase average of a product of n gyroradius vector components is denoted by Obviously, α and where S 2l is the symmetric group of permutations of the set {1, 2, • • • , 2l} and η j1•••j 2l is defined by with and Here, b i is the ith contravariant component of b ≡ B/B and δ ij represents the Kronecker delta; δ ij = 1 (for i = j), 0 (for i = j).In Eq. ( 3), we also find that where is defined and it satisfies The expressions given in Eqs. ( 7), (13), and ( 14) are valid in the Cartesian spatial coordinates although they can be easily transformed into those in general spatial coordinates as shown in Sec.III.

III. EULERIAN VARIATIONAL PRINCIPLE FOR DERIVATION OF THE COLLISIONLESS GYROKINETIC SYSTEM OF EQUATIONS IN GENERAL SPATIAL COORDINATES
In this section, the governing equations for the distribution functions and the turbulent electrostatic potential in the collisionless gyrokinetic system are derived from the Eulerian variational principle using the Lagrangian density represented in general spatial coordinates.

A. The Lagrangian density represented in general spatial coordinates
The action integral I GKF for the Eulerian variational principle to derive all of the governing equations of the collisionless electrostatic gyrokinetic turbulent system is written as where the Lagrangian density L GKF is given by Here, d 3 x ≡ dx 1 dx 2 dx 3 and d 3 v ≡ dv dµdϑ are used and the subscript a represents the particle species with mass m a and charge e a .We use x ≡ (x i ) i=1,2,3 and v ≡ (v , µ, ϑ) as the gyrocenter phase-space coordinates.It should be emphasized that, in this section, x ≡ (x i ) i=1,2,3 represent general spatial coordinates of the gyrocenter position which can be either Cartesian or any other curved coordinates.However, here we assume that the spatial position vector r = r(x) is 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 on the (x, v)-space is denoted by F a , and the number of particles of species a in the phase-space volume element d 3 xd 3 v ≡ dx 1 dx 2 dx 3 dv dµdϑ at time t is given by F a (x, v, t)d 3 xd 3 v.The field part L F of the Lagrangian density L GKF in Eq. ( 18) is written in terms of the covariant components of the electric field due to the electrostatic potential, In this paper, we employ the summation convention that the same symbol used for a pair of upper and lower indices within a term [such as seen in Eq. ( 18) 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 g ij in the general spatial coordinates x ≡ (x i ) are related to the covariant components g ij by g ik g kj = δ i j , where δ i j represents the Kronecker delta.The determinant of the covariant metric tensor matrix is denoted by Note that g ij (x), g ij (x), and g(x) are all independent of time t because the spatial position vector r is assumed to be given by a function of only the spatial coordinates x ≡ (x i ) i=1,2,3 , independently of time t as mentioned earlier.
The gyrocenter Lagrangian for species a, which is multiplied by F a to define L GK in Eq. (18), is represented in the Eulerian picture by where A j is the jth covariant component of the vector potential, b i ≡ B i /B is the ith contravariant component of the unit vector parallel to the background magnetic field, and the magnetic field strength is given by The contravariant components (B i ) i=1,2,3 of the magnetic field are expressed in terms of the covariant components (A i ) i=1,2,3 of the vector potential as where the Levi-Civita symbol is denoted by The gyrocenter Hamiltonian is written here as and the electrostatic fluctuation part is given by Ψ a (x, µ, t) ≡ φ(x, t) + Ψ E1a (x, µ, t) + Ψ E2a (x, µ, t), (26)   where µ, t) are defined by Eqs. ( 8)-( 12) and (15) in which h ij is represented in the general spatial coordinates x ≡ (x i ) i=1,2,3 by and the covariant derivative ∇ j is defined in Appendix A.
The Eulerian representation of the gyrocenter Lagrangian L GY a shown in Eq. ( 21) contains the temporal change rates of the gyrocenter position coordinates, parallel velocity, magnetic moment, and gyrophase at time t which are represented by the functions u i ax (x, v, t), u av (x, v, t), u aµ (x, v, t), and u aϑ (x, v, t), respectively, in the same way as in Ref. 30 .Then, the distribution function F a (x, v, t) satisfies We find here that the gyrocenter Hamiltonian H GY a given in Eq. ( 25) takes a functional form, 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 φ(x, t)}, {∂ J g ij (x)}].Here, we use the notation where F is an arbitrary function of x = (x i ) i=1,2,3 .
Then, the compact notations {∂ J φ(x, t)} and {∂ J g ij (x)} in Eq. (31) imply and respectively.Note that the high-order spatial derivatives included in Eqs.(33) and ( 34) enter the gyrocenter Hamiltonian H GY a due to finite gyroradius as seen in Eqs.(27) and Eqs.(28) where the covariant derivatives contain the spatial derivatives of g ij through the Christoffel symbols defined by Eq. (A4) in Appendix A.
In the same way as in Eq. ( 31), the functional form of the gyrocenter Lagrangian L GY a is written as where the Eulerian representations of the temporal change rates of the gyrocenter position and the gyrophase, u i ax (x, v, t) and u aϑ (x, v, t), are additionally included.In Eq. ( 35), u ax (x, v, t), u aϑ (x, v, t), and φ(x, t) are the functions, the governing equations of which are derived from the variation principle in Sec.III.C while the dependence of L GY a on A i (x, t), ∂ j A i (x, t), and ∂ J g ij (x, t) is also explicitly shown because their variations need to be taken into account to evaluate the variation of L GY a in Sec.IV where, in order to derive the local momentum balance, we consider the general spatial coordinate transformation which causes the variations in the functional forms of both (u ax , u aϑ , φ) and (A i , g ij ).

B. The Lagrangian density associated with polarization
We find from Eqs. ( 21) and ( 25)-( 28) that the part of the Lagrangian density L GY a which includes the perturbed potential is given by where the gyrocenter density N (g) a is defined by L E1a is given in the linear form with respect to the longitudinal electric field (E L ) i ≡ −∂φ/∂x i and its covariant derivatives, and L E2a is given in the quadratic form, In this paper, the longitudinal (irrotational) and transverse (solenoidal) parts of any vector field are represented by the subscripts L and T , respectively. 51In Eq. ( 40), there appear the multipole moments of the electric charge distribution 51 of species a induced by finite gyroradius, which can exist even without the electric field (E L ) i .The multipole moments

Ea
of the electric charge distribution of species a induced by the longitudinal electric field (E L ) i is shown in Eq. ( 41), and they are given in the linear form with respect to (E L ) i and its covariant derivatives as Here, the coefficients which are regarded as the generalized electric susceptibility produced by finite gyroradius.It is remarked that has the same symmetry.Taking the summation of Eq. ( 38) over species a, we obtain where the gyrocenter charge density ρ (g) is given by and the parts of the Lagrangian density associated with polarization are represented by L E1 and L E2 .Here, L E1 is defined by with the multiple moments induced by finite gyroradius, and L E2 is defined by with the multiple moments induced by the longitudinal electric field (E L ) i and its covariant derivatives, and the generalized electric susceptibility due to finite gyroradius, C. Derivation of the governing equations of the collisionless electrostatic gyrokinetic turbulent system Here we virtually allow the phase-space trajectories for all species and the electrostatic potential to vary infinitesimally.Following the same procedure as in Ref. 30 , the variations in the gyrocenter position, parallel velocity, magnetic moment, and gyrophase of the phasespace trajectory are represented in the Eulerian picture by δx i aE (x, v, t), δv a E (x, v, t), δµ aE (x, v, t), and δϑ aE (x, v, t), respectively.We also denote the variation in the electrostatic potential by δφ.Then, in terms of δx i E , δv E , δµ E , and δϑ E , the variations in the functional forms of u i ax , u av , u aµ , and u aϑ due to the virtual displacement are written as and the variation in the distribution function F a is given by Using Eqs. ( 17), ( 18), ( 45), ( 52) and ( 53), the variation in the action integral I GKF is expressed as where D i represents the electric displacement vector defined later in Eq. ( 73) and B.T. represents boundary terms that appear due to partial integrals.
Here, (∂L GY a /∂x i ) u , (∂L GY a /∂v ) u , (∂L GY a /∂µ) u , and (∂L GY a /∂ϑ) u denote the derivatives of L GY a in x i , v , µ, and ϑ, respectively, with (u i ax , u aϑ ) kept fixed in L GY a , and the time derivative along the phase-space trajectory is represented by We now employ the Eulerian variation principle which implies that the collisionless gyrokinetic equations for the distribution functions of all species and the gyrokinetic Poisson's equation for the electrostatic potential can be derived from the condition that δI GKF = 0 for arbitrary variations δx i aE , δv a E , δµ aE , δϑ aE , and δφ which vanish on the boundaries of the integral region to make B.T. disappear in Eq. (54).
We first use δI GKF /δx i aE = 0 to obtain where p ai represents the covariant vector component of the canonical momentum defined by We should note that the distribution function F a is included as a factor in δI GKF /δx i aE = 0 although it is omitted from Eq. ( 56) for simplicity.This omission of F a is also performed in the other equations obtained below from δI GKF /δv a E = 0, δI GKF /δµ aE = 0, and δI GKF /δϑ aE = 0 although it does not make a difference in deriving the resultant collisionless gyrokinetic equation in Eq. ( 69).We can rewrite Eq. ( 56) as where the modified magnetic field is defined by respectively.
Next, δI GKF /δv a E = 0 is used to obtain from which we have Furthermore, δI GKF /δµ aE = 0 and δI GKF /δϑ aE = 0 yield and Equations ( 58), ( 61), (62), and (63) are rewritten as and where Ω a ≡ e a B/(m a c) and Equations ( 64) and ( 65) are obtained by taking the vector and scalar products between the magnetic field and Eq. ( 58), respectively.We can verify that the right-hand sides of Eqs. ( 64)-( 67) are all independent of ϑ and that the magnetic moment µ is an invariant of motion as seen from Eq. (66).Substituting Eqs. ( 64)-(67) into Eqs.( 30) and taking its average with respect to the gyrophase ϑ, the collisionless gyrokinetic equation is derived as where F a denotes the gyrophase-averaged distribution function, The remaining governing equation of the system, namely, the gyrokinetic Poisson's equation is derived from the condition that the variational derivative of the action integral I GKF with respect to the electrostatic potential φ vanishes, δI GKF /δφ = 0. Since the time derivative of φ never appears in the Lagrangian density L GKF , the above-mentioned condition can be replaced using the Lagrangian L GKF instead of I GKF by where δ 3 x with the subscript x = (x i ) i=1,2,3 represents the function that takes a value δ 3 x (y) = δ(y 1 − x 1 )δ(y 2 − x 2 )δ(y 3 − x 3 ) at y = (y i ) i=1,2,3 .Equation (71) gives the gyrokinetic Poisson's equation, where the electric displacement vector D i is written as Here, the generalized polarization vector density P i G is defined by in which not only the dipole moment but also other multipole moments 51 occurring due to finite gyroradius are included.The multipole moments Q i1•••im in Eq. ( 74) are generally given by the sum of the two parts, where are defined by Eqs. ( 48) and (50), respectively.It should be noted that Q i1•••im 0 vanishes unless m is an even integer.Taking only the contribution of m = 1 to the summation over m in Eq. (75), we obtain the dipole moment density which is denoted by The presence of the infinite series due to the finite gyroradius in Eq. ( 74) and other places has an analogy with that for the case of macroscopic electromagnetism described in detail in Ref. 51 where the macroscopic charge density is evaluated for the system consisting of molecules.The contribution of each molecule to the macroscopic charge density is calculated by spatially averaging the microscopic charge density (given by the point charges constituting the molecule) around the center of mass of the molecule.Then, the resultant expression of the macroscopic charge density is given by the series expansion associated with the multipole moments due to the finite distance of each point charge from the center of mass of the molecule.The local spatial average of the microscopic charge density in the system of molecules is replaced by the phase-space integration in the present case of the gyrokinetic system to represent the macroscopic charge density as which contains the polarization effect due to the finite gyroradius and the microscopic electrostatic fluctuations.In Eq. ( 77), the Cartesian spatial coordinates are used and three-dimensional vectors are represented in terms of boldface letters.In the square brackets on the righthand side of Eq. ( 77), the first and second terms give the point charge density of the particle at x = x ′ + ρ a and its correction due to the electrostatic fluctuation, respectively.We can see that the effect of the finite gyroradius (or the finite distance between the particle and gyrocenter positions) is included in the delta functions, which cause the infinite series expansions to appear in Eq. ( 74) and other equations related to the polarization (or dipole and multipole moments).

IV. DERIVATION OF THE MOMENTUM BALANCE
In this section, we use the invariance of the Lagrangian under arbitrary infinitesimal transformations of spatial coordinates to derive the momentum balance equation for the single-particle-species system and that for the whole system including all species and the field.

A. Invariance of the Lagrangian under arbitrary spatial coordinate transformations
We now consider the infinitesimal transformation of the spatial coordinates from x = (x i ) i=1,2,3 to x ′ = (x ′i ) i=1,2,3 , which is written as Here, the infinitesimal variation in the spatial coordinate x i 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.From Eq. ( 18), the gyrokinetic part of the Lagrangian for species a is found to be given by where the gyrocenter Lagrangian L GY a is defined in Eq. (21).Since L GKa shown in Eq. ( 79) is a scalar con-stant which is invariant under arbitrary spatial coordinate transformations including the infinitesimal transformation given by Eq. (78), we have Note that here and hereafter we use δ • • • to represent the variation associated with the infinitesimal spatial coordinate transformation which should be distinguished from the variation δ • • • due to the virtual displacement treated in Sec.III.The divergence term in the integrand of Eq. ( 80) appears due to 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 in the integrand represents the variation in the spatial functional form of L GKa due to the infinitesimal spatial coordinate transformation.
As shown in Appendix B, the variation in the spatial functional form of an arbitrary tensor field (as well as an arbitrary tensor field density) due to the infinitesimal spatial coordinate transformation given in Eq. ( 78) is written as the opposite sign of its Lie derivative 52 with respect to the generating vector field ξ i , and it is represented by δ = −L ξ .It is also noted that since the Lagrangian density behaves as a scalar field density under arbitrary spatial coordinate transformations, the Lie derivative of L GKa is given by Therefore, the integrand in Eq. ( 80) is simply rewritten as L ξ L GKa − L ξ L GKa (= 0), from which its spatial integral δL GKa is naturally found to vanish as shown in Eq. (80).We now use Eq. ( 79) and the Leibniz rule for the derivative operation by δ = −L ξ to write the variation in the spatial functional form of the Lagrangian density L GKa as Then, using Eqs.( 81)-( 82) and δ = −L ξ , Eq. ( 80) is rewritten as Equation ( 83) is found to hold by noting that L ξ L GY a = ξ i ∂L GY a /∂x i (because L GY a behaves as a scalar field under arbitrary spatial coordinate transformations) and that the chain rule is applied to the derivative operation where ∂L GY a /∂(∂ J A i ) = 0 when the order of J is greater than or equal to two [see Eq. ( 35)].
We next consider the invariance of the Lagrangian L GKF of the whole system under the infinitesimal spatial coordinate transformation, which is written in the same way as in Eq. ( 80) by where L GKF consists of the gyrokinetic and field parts as shown in Eqs. ( 17) and (18).Procedures similar to those leading to Eq. ( 83) can be made to obtain The invariance shown in Eq. (86) can also be confirmed using Eq. ( 84) and the following formula for the derivative of It is summarized from Eqs. (83), (84), (86), and (87) that the invariance of the scalar constants of L GKa and L GKF under the infinitesimal spatial transformation can be verified using the chain rule formulas for the derivative operation δ = −L ξ on the scalar field L GY a and the scalar field density L F .The invariance formulas shown in Eqs. ( 83) and ( 86) are used to derive the momentum balance equation for the single-species system in Sec.IV.B and that for the whole system including all species and the field in Sec.IV.C, respectively.
B. Derivation of the momentum balance for a single particle species We now use the Euler-Lagrange equations for gyrocenter motion [Eqs.(56), ( 60), (62), and (63)] and perform partial integrals to rewrite (83) as where, instead of using Eq. ( 30), the gyrocenter distribution function F a is assumed to satisfy Here, K a represents the rate of temporal change in F a due to collisions and/or external sources for the species a.In the present work, we assume that is satisfied by K a .In fact, it is shown in Ref. 43 that Eq. (90) corresponds to the intrinsic ambipolarity of the classical particle fluxes when K a represents the collision operator in the gyrocenter coordinates.The variational derivative of L GKa with respect to A i is given by where the particle flux of species a is represented by The variational derivatives δL GKa /δφ and δL GKa /δg ij are written as and respectively, where the particle density a and the symmetric pressure tensor density P ij a of species a are defined by e a N (p) a ≡ e a N (g) a − ∇ i P i Ga (95) and respectively, and #J = n represents the order of J ≡ (j 1 , j 2 , • • • , j n ).
On the right-hand side of Eq. ( 95), P i Ga represents the contribution of species a to the generalized polarization vector density P i G defined in Eq. ( 74) and it is written as On the right-hand side of Eq. ( 96), P ij CGLa is given in the Chew-Goldberger-Low (CGL) form, 41 and π ij ∧a is defined by Here, the perpendicular component of the gyrocenter velocity is represented by (u ax ) i ⊥ ≡ u i ax − u k ax b k b i .We find that Eqs.(98) and (99) agree with those given by Eqs. ( 124) and (125) in Ref. 30 except for the effects of the electrostatic fluctuations included in F a and (u ax ) i ⊥ on the right-hand side of Eqs.(98) and (99).We also note that Eq. (99) contains the product of the fluctuation parts of F a and (u ax ) i ⊥ , the ensemble average of which does not disappear but contributes to the turbulent momentum transport.In the neoclassical transport theory, 41,42 it is considered that the CGL pressure tensor shown in Eq. (98) contains the scalar (or isotropic) part, which represents background pressure, and the anisotropic part, the magnitude of which is smaller than the background pressure by the factor ∼ ρ/L where ρ and L represent the gyroradius and the equilibrium gradient scale length, respectively.The anisotropic part of the CGL pressure tensor causes the viscous force which plays an essential role in the neoclassical transport processes when the distribution function deviates from the local Maxwellian under the influence of collisions. 41,42he magnitude of π ij ∧a defined in Eq. ( 99) is regarded as ∼ (ρ/L) 2 .
The last term on the right-hand side of Eq. ( 96) is given by (100) where P ij E1a and P ij E2a are defined in Eqs.(C5) and (C8) of Appendix C, respectively.We note that the effects of the turbulent electrostatic potential are included in the definitions of P ij Ψa explicitly as well as in π ij ∧a through the turbulent drift velocity part of (u ax ) i ⊥ .Then, the nonlinear interaction of the turbulent potential and the fluctuation part of the gyrocenter distribution function included in P ij Ψa and π ij ∧a causes the turbulent momentum transport.In Appendix D, the ensemble-averaged pressure tensor describing the turbulent momentum transport is given by the WKB representation.
Performing further partial integrals in Eq. ( 88) finally gives where and In deriving Eq. ( 101)-(103), we have deformed Eq. (88) using Eqs.(89), (91), ( 93), (94), )] and partial integration.Especially, the term ∇ i P ij a in Eq. ( 102) is derived from the term (δL GKa /δg ij )δg ij which is rewritten as , where P ij a = P ji a is also used and the spatial integral of the last term ∇ i (ξ j P ij a ) becomes one of the boundary terms.We now recall that an arbitrary infinitesimal vector field can be employed as ξ i .Then, in order for Eq.(101) to hold for any ξ i , we need to have J j GKa = 0, which gives the momentum balance equation as As seen in Eqs. ( 126) and (152) of Ref. 30 , the relation between the symmetric and canonical pressure tensors is obtained from the other condition for the sum of the boundary terms (B.T.) in Eq. ( 101) to vanish although its complicated expression is not shown here.We see that the inertia term in the momentum balance equation, Eq. ( 104), contains only the parallel momentum component while the electric current e a Γ k a in the Lorentz force term consists of the guiding-center current and the magnetization current as shown in Eq. (92) [see also Eq. ( 114)].For comparison with Eq.( 104), the canonical momentum balance equation derived by the conventional method is shown in Eq. (E5) of Appendix E where the divergence of the asymmetric canonical pressure tensor appears.In addition, the energy balance equation for the single-species gyrokinetic system is given in Eq. (F2) of Appendix F.

C. Derivation of the momentum balance for the whole system
In the same way as in deriving Eq. ( 88), performing partial integrals in Eq. ( 86) and using Eqs.(89), we have Now, substituting Eq. (91) into δL GKF /δA i = a δL GKa /δA i and using the gyrokinetic Poisson's equation, δL GKF /δφ = 0, and Eq.(94), Eq. ( 105) is finally rewritten as where Here, the symmetric pressure tensor Θ ij is defined by where the last group of terms including (E L ) i represents the Maxwell stress tensor due to the electrostatic field with the opposite sign.Taking the summation of Eqs. ( 98), (99), and (100) over species defines P ij CGL , π ij ∧ , and P ij Ψ on the right-hand side of Eq. (108) as and respectively.As mentioned after Eq. (100) as well as in Appendix F, the turbulent momentum transport caused by the nonlinear interaction of the turbulent potential and the fluctuation part of the gyrocenter distribution function is included in π ij ∧ and P ij Ψ .From Eqs. ( 106) and (107), we obtain J j DKF = 0 which represents the momentum balance equation for the whole system, For comparison with Eq.( 112), the canonical momentum balance equation derived by the conventional method is shown in Eq. (E11) of Appendix E where the divergence of the asymmetric pressure tensor appears.The condition for the sum of the boundary terms (B.T.) in Eq. ( 106) to vanish results in a complicated expression which is not shown here while it gives the relation between the symmetric pressure tensor Θ ij and the asymmetric canonical pressure tensor in Eq. (E11).Here, using Eq. ( 92), the electric current density can be written as where the covariant components of the magnetization vector are defined by Using Eqs. ( 89) and (90), we have where ρ (g) is defined in Eq. ( 46).Combining Eq. ( 72) and (115), we obtain where the subscript L represents the longitudinal part.It should be noted that the polarization current 4π∂P i G /∂t is included not in J i but in ∂( √ gD i )/∂t.
It can be shown by performing further vector calculations that Eq. ( 112) is rewritten as where the Cartesian spatial coordinates are used and three-dimensional vectors are represented in terms of boldface letters.The longitudinal part of the electric displacement vector defined in Eq. ( 73) is represented by D L .The longitudinal part of the electric field is written as E L ≡ −∇φ while E T ≡ −c −1 ∂A/∂t gives the transverse part under the Coulomb gauge condition ∇ • A = 0. We see that the time derivative of the momentum density (D L × B)/(4πc) due to the electromagnetic field and the spatial divergences of the pressure tensors produced by both the electric and magnetic fields (the opposite sign of the Maxwell stress tensor) appear in the momentum balance equation in Eq. ( 117) where the effects of the polarization P G [see Eq. ( 74)] are included through D L .
Here, recall that F a (x, v, t) represents the gyrocenter distribution function of the gyrocenter (not particle) position coordinates x ≡ (x i ) i=1,2,3 , v ≡ (v i ) i=1,2,3 ≡ (v , µ, ϑ) and time t [see descriptions after Eq.( 18)].Then, from typical gyrokinetic codes using the gyrocenter position coordinates as independent variables for distribution functions, we can directly evaluate the term 117), for which we do not need to specify the transformation from the gyrocenter position coordinates to the particle coordinates.As explained at the end of Sec.III.C, the difference between the gyrocenter and particle positions is taken into account in Eq. (117) through the terms related to the polarization due to the effects of the finite gyroradius and the electrostatic fluctuation.Such polarization terms also appear in the energy balance equations as seen in Eqs.(F7) and (F8) in Appendix F.
When the polarization P G is eliminated from Eq. ( 117), the terms including the electric field are given by ∂ )] which are verified to be the same as given in the momentum conservation law for the Vlasov-Poisson-Ampère (or Vlasov-Darwin) system [see Eq. (32) of Ref. 39 ].In the case using the quasineutrality condition and the self-consistent magnetic field given by ∇ × B = (4π/c)J (J = J T from ∇ • J = 0 due to the quasineutrality) as well as removing polarization effects and K a , we find that (D L × B)/(4πc) and the part of the pressure tensor caused by the electric field disappear from Eq. ( 117) and that Eq. ( 117) agrees with the momentum conservation law for the drift kinetic system shown in Eq. (151) of Ref. 30 .
The energy balance equation for the whole system is derived in Appendix F [see Eqs.(F7) and (F8)] where, in the same way as seen in Eq. ( 117), we can confirm the consistency of the derived energy balance with the energy conservation laws obtained for the Vlasov-Poisson-Ampère (or Vlasov-Darwin) system 39 and the drift kinetic sytem 30 .
We note that, if the Lagrangian L GKF for the gyrokinetic system defined by Eqs. ( 17) and ( 18) is modified to √ gB 2 /(8π), the variational equation δL GKF * /δA = 0 yields ∇ × B = (4π/c)J T which makes Eq. ( 117) take the form of the total momentum conservation law in the case of a d 3 v K a m a v = 0. Nevertheless, in the gyrokinetic turbulent system, this condition is not generally imposed on the given equilibrium magnetic field B. It is because B is considered not to contain the fluctuation part while J T can have fluctuations.However, when the background magnetic field satisfies spatial translation, rotation, or helical symmetry and the effect of K a is neglected, the local conservation law of the canonical momentum in the direction of symmetry can be derived from δL GKF = 0 with Eq. (105) as shown in Appendix G.

V. CONCLUSIONS
In this paper, the governing equations of the gyrokinetic system with electrostatic turbulence are derived in the general spatial coordinates based on the Eulerian variational principle.The local momentum balance equation for each particle species and that for the whole system, which the gyrocenter distribution functions and the potential field satisfy, are obtained from the invariance of the Lagrangians of these systems under arbitrary spatial coordinate transformations.
It is shown that, when the background magnetic field satisfies the consistency condition that its rotation is given by the solenoidal part of the current density as in the Darwin model, the momentum and energy balance equations for the whole system are rewritten in the complete conservative forms where contributions of the turbulent electric field and the background magnetic field are clearly given in the expressions of the momentum and energy densities, the pressure tensor, and the energy flux.The effects of the collision and/or external source terms added into the gyrokinetic equation on the momentum and energy balance equations are clarified as well.
The symmetric pressure tensor is directly obtained by the variational derivative of the Lagrangian with respect to the metric tensor components and it is shown to contain the CGL part representing the neoclassical viscosity as well as the turbulent momentum transport part, the ensemble average of which is confirmed to agree with the previous result obtained from the gyrokinetic theory using the WKB representation.
The representation in terms of the general spatial coordinates is useful in treating complex toroidal plasmas in suitable coordinates such as the flux coordinates.The momentum and energy balance equations obtained here are applicable as a reference for verification of long-time global gyrokinetic simulations based on the Lagrangian and Hamiltonian formulations to study neoclassical and turbulent transport in plasmas with external sources.It may seem troublesome for global simulation codes to treat the finite gyroradius effect represented by the infinite series expansions appearing in Eq. ( 27) and other places.However, for those simulations in which the expansions are truncated or approximated by other simpler expressions such that the Lagrangian corresponding to the reduced model given by the truncation or approximation is clearly defined, the same technique as shown in this work can be applied to that Lagrangian to derive the local energy and momentum balance equations which can be compared with those simulation results.Such applications of the present work to the global simulations based on the reduced gyrokinetic model are considered as future works.For comparison with local flux tube gyrokinetic simulations [53][54][55][56][57] treating the full finite gyroradius effect, useful informations such as expressions of local turbulent momentum transport can be obtained from this work using the WKB approximation as shown in Appendix D. The extension of the present work to the case with magnetic microturbulence also remains as a future study.
where we see that the Lie derivative L ξ can be used again to represent δV i (x).The components of a covariant vector field W i (x) are transformed as Next, following the procedure similar to those used in deriving Eqs.(B3) and (B6), the variation in the functional form of W i (x) under the infinitesimal spatial coordinate transformation is derived as It is shown in the same way as shown above that the variation in the functional form of a mixed tensor field is given by the opposite sign of its Lie derivative as It is also shown that the variations in the functional forms of g ij and g ij under the infinitesimal spatial coordinate transformation are expressed as We now apply the chain rule to the derivative operation δ = −L ξ on √ g = det(g ij ) and use Eqs.(A8) and (B10) to obtain We here note that an arbitrary scalar density S is transformed in the same manner as √ g under arbitrary spatial coordinate transformations, and that S(x)/ √ g(x) is regarded as a scalar field.Then, the variation δS(x) in the functional form of S(x) under the infinitesimal spatial coordinate transformation is given by δS Since the distribution function F (x, v, t) behaves as a scalar density under arbitrary spatial coordinate transformations, the variation δF in its spatial functional form due to the infinitesimal spatial coordinate transformation is written using Eq.(B13) as Note that u i ax are the contravariant components of the gyrocenter velocity vector field while u av , u aµ , and u aϑ , which represent the temporal change rates of the parallel velocity, magnetic moment, and gyrophase, respectively, behave as scalar fields under arbitrary spatial coordinate transformations.Then, we find from Eqs. (B6) and (B3) that the variations in the functional forms of u i ax , u av , u aµ , and u aϑ under the infinitesimal spatial coordinate transformation are given by It is also useful to know that even though the Christoffel symbols are not regarded as tensor components, their derivatives with respect to δ can be defined by the variations in their functional forms under the infinitesimal spatial coordinate transformation and written as Equation (B16) can be derived from Eqs. (A4), (A7), (B10), and the commutative property,

Appendix C: VARIATIONAL DERIVATIVE WITH RESPECT TO METRIC TENSOR COMPONENTS
We here consider the parts of the Lagrangian including the polarization effects denoted by L E1a and L E2a , where the L E1a and L E2a defined in Eqs. ( 40) and ( 41) appear due to finite gyroradius for species a and they are given, respectively, in the linear and quadratic forms of the longitudinal electrostatic field and its spatial gradients.For the purpose of deriving the pressure tensor P ij Ψa due to electrostatic gyrokinetic turbulence, we need to evaluate the variational derivatives of L E1a and L E2a with respect to the metric tensor components g ij .It should be noted here that partial derivatives with respect to g ij need to be carefully performed because 3 × 3 metric tensor components g ij are not completely independent of each other due to the constraint g ij = g ji .Here, for an arbitrary function f of g ij , the notation ∂f /∂g ij is defined such that the infinitesimal variations δg ij in g ij give rise to the variation δf = (∂f /∂g ij )δg ij in f where both δg ij and ∂f /∂g ij must be symmetric under exchange of the indices i and j. 38 For example, we have ∂g kl /∂g ij = 1 2 (δ i k δ j l + δ j k δ i l ) according to the abovementioned definition.In the same manner, derivatives with respective to ∂g ij /∂x k are defined taking into account the symmetry under exchange of the indices i and j.
The variation in L E1a caused by the variation in the metric tensor with keeping the gyrocenter distribution function F a fixed in L E1a is written as where integration by parts is repeatedly performed to finally derive (δL E1a /δg ij ) F .Here, (• • • ) F implies that the gyrocenter function F a is fixed when taking the variation with respect to the metric tensor.From Eqs. (A7), we have where We note that, in Eq. (C4), δ g Γ p j l jn is expressed in the same form as in Eq. (B16).Now, the pressure tensor E1a is given by two times of (δL E1a /δg ij ) F as where Eqs.(C2)-(C4) are used and (∂Q j1•••j 2k 0a /∂g ij ) F is given by using Eq. ( 42) and Next, the variation in L E2a caused by the variation in the metric tensor with keeping the gyrocenter distribution function F a fixed is written as Then, using Eqs.(C3), (C4), and (C7), we obtain the pressure tensor P ij E2a which is given by two times of (δL E2a /δg ij ) F as where (∂χ i1•••im;j1•••jn a /∂g ij ) F is given by using Eqs.( 15), ( 44), (C6), and Appendix D: WKB REPRESENTATION Here, we use the WKB (or ballooning) representation 49 for turbulent fluctuations which have small wavelengths of the order of the gyroradius ρ in the directions perpendicular to the background magnetic field.Such rapid spatial variations are represented using the perpendicular wavenumber vector k ⊥ .We assume that k ⊥ ρ = O(1) and ρ/L ≪ 1 where L is the gradient scale length of equilibrium variables.
The gyrocenter distribution function F a for species a is given by the sum of the zeroth and first-order parts in ρ/L as F a (x, v) = F a0 +F a1 + Fa1 where the zeroth-order part F a0 is the local Maxwellian equilibrium distribution function, and the first-order part representing the deviation from the local Maxwellian consists of the non-turbulent and turbulent functions denoted by F a1 and Fa1 , respectively.Then, the turbulent gyrocenter distribution function Fa1 is expanded as Fa1 = k ⊥ Fa1k ⊥ exp(ik ⊥ • x), where x is the gyrocenter position vector and the k ⊥component F a1k ⊥ is given by the sum of the adiabatic and nonadiabatic parts as 59 Here, φk ⊥ and ĥak ⊥ are the k ⊥ -components of the turbulent electrostatic potential and the nonadiabatic part of the turbulent distribution function, respectively, J 0 is the zeroth-order Bessel function, and T a is the background temperature of the species a.
We now follow the conventional assumption that Fk ⊥ ens = 0 and are satisfied by arbitrary real-valued turbulent fluctuations F and F ′ where (• • • ) * and • • • ens represents the complex conjugate and the ensemble average, respectively.Note that the ensemble-averaged quantities are smooth spatial functions with the gradient scale length L and that F * When taking the ensemble average of the momentum balance equation in Eq. ( 117), we find that the effects of the turbulent fluctuations on the momentum transport are included in π ∧ ens and P Ψ ens through the correlation between the nonadiabatic part of the turbulent distribution function and the turbulent electrostatic potential which are described by ĥ * ak ⊥ φk ⊥ ens .Neglecting terms of higher orders in ρ/L, we find that the turbulent contribution π ∧Ψ ens to π ∧ ens is given by and P Ψ ens is written as where I denotes the unit tensor.We now consider the axisymmetric toroidal system, in which the magnetic field is given by where ζ, I, and χ represent the toroidal angle, the covariant toroidal magnetic field component, and the poloidal flux function, respectively.Then, using (k The flux-surface average of Eq. (D5) agrees with the electrostatic and low-flow ordering limit of the result given in Eq. ( 53) of Ref. 24 where the turbulent radial transport of the toroidal angular momentum double-averaged over the ensemble and the flux surface is presented for the general case allowing the turbulent magnetic field and the high-flow ordering 28,29,60 .It should be emphasized here that Eq. (D5) is not surface-averaged but it presents the spatially-local expression.In the axisymmetric configuration with up-down symmetry, the flux-surface average of Eq. (D5) is shown to vanish in the case of the low-flow ordering 27,61 although it does not imply that the local value of Eq. (D5) itself vanishes as well.

Appendix E: ANOTHER DERIVATION OF THE MOMENTUM BALANCE
In Sec.IV, the momentum balance is derived using the invariance of the Lagrangian of the system under the infinitesimal spatial coordinate transformation induced by the vector field ξ i (x) which has an arbitrary functional form.For comparison with that derivation, another derivation of the momentum balance is given in this appendix in the way closer to the conventional derivation of the canonical momentum conservation law which generally involves the asymmetric canonical pressure tensor.
Equation (83) for the invariance of the gyrokinetic Lagrangian L GKa of species a can be rewritten without separating boundary terms from the spatial integral as where Here, p ai ≡ ∂L GY a /∂u i ax is the canonical momentum [see Eq. ( 57)] and represents the variation in the functional form of the single gyrocenter Lagrangian L GY a under the infinitesimal spatial coordinate transformation with (u i x , u ϑ ) kept fixed in L GY a .The term (∂L GY a /∂x i ) u in Eq. (E3) is written down as Equations ( 56), ( 89), (B15), and (E3) are used for deriving Eq. (E2).Since Eq. (E1) holds for an arbitrary spatial integral domain V , we find J GKa = 0 for any ξ i .Then, the canonical momentum balance equation for the gyrocenters of species a is derived from Eq. (E2) as Equation (E5) also can be derived directly from multiplying Eq. ( 89) by p ai , taking its v-space integral, and using Eq. ( 56).It should be recalled that the general spatial coordinates (x i ) i=1,2,3 are used here.For example, when (x i ) i=1,2,3 represent the Cartesian coordinates (x, y, z), Eq. (E5) represents the linear canonical momentum balance.As another interesting example, we can treat the canonical angular momentum balance in toroidal plasmas such as tokamaks and stellarator/heliotron devices, for which it is convenient to use the cylindrical coordinates and/or the magnetic flux coordinates.In these coordinate systems, the toroidal angle component of Eq. (E5) represents the toroidal angular momentum balance.It is also noted that the momentum balance equations shown in Eqs. ( 104) and (E5) are equivalent to each other although the transformation from Eq. (E5) to Eq. ( 104) requires complicated procedures involving infinite number of times of partial integration.
We next consider the invariance of the Lagrangian L GKF of the whole system under the infinitesimal spatial coordinate transformation shown in Eq. (86).In contrast to Eq. (105), we now rewrite Eq. (86) without separating boundary terms from the spatial integral as with where Eqs. ( 56), ( 72), (89), and (B15) are used.Here, J GKF A and J GKF g originate from the terms including δA j and δg ij in Eq. ( 86), respectively, and they are written as and (E9) The last term J GKF φ included in Eq. (E7) is given in terms of the turbulent electrostatic field as Here, because of the symmetry with respect to permutations of the indices j 1 , • • • , j n , the variables  We now suppose that ξ i represents a Killing vector field or a vector field which generates an isometric transformation so that δg ij = −L ξ g ij = 0, δ(∂ J g ij ) = ∂ J (δg ij ) = 0 and accordingly J GKF g = 0.The isometry is generated by the infinitesimal linear translation and the infinitesimal rotation, and they are represented by the Killing vector fields, ξ = ǫn and ξ = ǫn × r, respectively, where ǫ is an infinitesimal constant, n a constant unit vector, and r the spatial position vector.Here and hereafter, the Cartesian coordinate system is used for (x i ) i=1,2,3 and three-dimensional vectors are represented in terms of boldface letters.For example, the canonical momentum of the gyrocenter of species a is denoted by p a ≡ ∂L GY a /∂u ax = m a v b + (e a /c)A.
Recall again that J GKF defined in Eq. (E7) vanishes because Eq. (E6) holds for any V .Then, substituting ξ = ǫn into Eq.(E7) and noting that J GKF = 0 is satisfied for any direction vector n, we obtain the linear canonical momentum balance equation for the whole system, ∂ ∂t where M is the magnetization vector defined in Eq. ( 114) and Π Ψ is defined using Eq. ( 45) as . (E12) Like P Ψ given in Eq.( 111), Π Ψ contains the effect of the electrostatic field on the momentum transport although P Ψ and Π Ψ are symmetric and asymmetric tensors, respectively.Using Eqs. ( 90) and (115), Eq. (E11) also can be rewritten as The momentum balance equation shown in Eq. (E13) can be transformed into Eq.( 112) although again it is so complicated involving infinite number of times of partial integration.Substituting ξ = ǫn × r into Eq.(E7) and using J GKF = 0 for any direction vector n, the angular momentum balance equation for the whole system is derived as where T Ψ is the asymmetric tensor representing the transport of the angular momentum due to the electrostatic field and defined by .

Appendix F: ENERGY BALANCE IN THE GYROKINETIC SYSTEM
In this Appendix, the energy balance in the electrostatic gyrokinetic turbulent system is derived.In contrast to the case of in Sec.IV where the momentum balance is derived, we here use only the Cartesian coordinate system and represent three-dimensional vectors in terms of boldface letters.Then, the metric tensor components are represented by the Kronecker delta, and they form the 3 × 3 unit matrix with determinant unity.
The partial time derivative of the gyrokinetic Lagrangian density L GKa for particle species a is written as where (∂L GY a /∂t) u denotes the time derivative of L GY a with (u i ax , u aϑ ) kept fixed in L GY a .Here, we consider the gyrocenter distribution function F a satisfying Eq. ( 89) which includes the term K a representing the rate of temporal change in F a due to collisions and/or external sources for the species a. Substituting Eq. (89) into Eq.(F1) and using δI GK /δx aE = 0 [Eq.( 56)], δI GK /δv a E = 0 [Eq.( 60)], δI GK /δµ aE = 0 [Eq.( 62)], and δI GK /δϑ aE = 0 [Eq.( 63)], we obtain the energy balance equation, where the gyrocenter velocity u ax is given by Eq. ( 64) and E a represents the energy of the single particle (or the gyrocenter Hamiltonian H GY a ) defined by The rate of change in the particle's energy is given by (F4) The energy balance equation shown in Eq. (F2) agrees with Eq. (C7) in Ref. 30 except for the effects of the turbulent electrostatic potential included here.It can be seen in Eq. (F2) how the energy balance is modified when the collision (or source) term K a is added into the gyrokinetic equation.
We now consider the energy balance in the extended system consisting of particles of all species and the selfconsistent electrostatic field.The partial time derivative of the Lagrangian density L GKF of this extended system is written as In the same way as in deriving Eq. (F2), we use the Euler-Lagrange equations [Eqs.( 56), ( 60), (62), and (63)] for (u ax , u av , u aµ , u aϑ ), Eq. ( 89) for F a , and Eq. ( 72) for φ in order to rewrite Eq. (F5) as where E T ≡ −c −1 ∂A/∂t represents the electric field induced by the temporal change in the background magnetic field and the magnetization vector M is defined by Eq. (114).It should be recalled that the background magnetic field is allowed to temporally change in this paper so that the long time evolution of the system over the transport time scale can be treated.It is also noted that Eq. ( 90) is used in deriving Eq. (F6).
It can be shown after several analytical manipulations that Eq. (F6) is further rewritten as or where the longitudinal (irrotational) part of the electric displacement vector [see Eq. ( 73)] is denoted by D L ≡ (E + 4πP G ) L ≡ −∇φ D and the magnetic intensity field H is defined by H ≡ B − 4πM with the magnetic induction field B and the magnetization vector field M [see Eq. ( 114)].Equations (F7) and (F8) take the conservative form including no other terms than the time derivative of the total energy density and the divergence of the total energy flux when the right-hand sides vanish.The effect of the collision (or source) term on the total energy balance is shown by the integral term which appears on the right-hand sides of Eqs.(F7) and (F8).This term vanishes when K a represents the collision operator which conserves the summation of the gyrocenter kinetic and polarization energies a d 3 v{ 1 2 m a v 2 + µB + e a (Ψ E1a + Ψ E2a )}.In addition, the right-hand sides of Eqs.(F7) and (F8) contain the difference between the transverse current density J T and (c/4π)∇ × B which vanishes if the self-consistency condition ∇ × B = (4π/c)J T is imposed as is done in the Darwin model 62 and in our previous work by including the magnetic energy with the minus sign into the Lagrangian for the drift kinetic system with self-consistent fields. 30We also find from Eqs. (F7) and (F8) that the terms including E T appear on their left-hand sides in the way consistent with the case of the energy conservation law for the Vlasov-Poisson-Ampère (or Vlasov-Darwin) system shown in Eq.( 22) of Ref. 39 .
All the kinetic energy, the electric energy, and the magnetic energy with the polarization including the dipole and other multipole moment effects are described in Eqs.(F7) and (F8), where the energy flux contains the kinetic energy flow due to the gyrocenter motion, the Poynting vector, and the extra energy flux due to electrostatic turbulent fluctuations with wavelengths of the order of gyroradius.

Appendix G: THE LOCAL CONSERVATION LAW OF THE CANONICAL MOMENTUM IN THE DIRECTION OF SYMMETRY
It is shown in this Appendix how the local conservation law of the canonical momentum in the direction of symmetry can be derived for the gyrokinetic system considered in the present paper.We start from Eq. (105) for which the gyrokinetic equation including the collision (or source) term K a [see Eq. ( 89)] is used.Equation ( 105) is rewritten as where we have used the canonical momentum of the single gyrocenter, the electric current density, and the symmetric pressure tensor defined by p ai ≡ ∂L GY a /∂u i ax , J i ≡ a e a Γ i a = c a δL GKa /δA i = c δL GKF /δA i , and Θ ij ≡ 2 δL GKF /δg ij , respectively, as shown in Eqs. ( 57), (91), (113), and (108).Next, we use δL GKF /δφ = 0, which is equivalent to the gyrokinetic Poisson's equation [see Eq. ( 72)], and substitute )] into Eq.(G1).Then, Eq. (G1) is rewritten after performing partial in-tegrals as where we can also write the last term of the integrand as ξ j ∇ i Θ ij = ξ j ∇ i Θ i j by using ξ j ≡ g jk ξ k , Θ i j ≡ g jk Θ ik , and Eq.(A8).We now recall the invariance of the Lagrangian L GKF under an arbitrary spatial coordinate transformation which implies that δL GKF = 0 for any ξ j .Consequently, the canonical momentum balance equation is obtained from Eq. (G2) as It is emphasized here that Eq. (G3) is valid for arbitrary spatial coordinates (x i ) i=1,2,3 in which the spatial position vector in the Cartesian coordinate system is represented as r = r(x 1 , x 2 , x 3 ).Hereafter, we assume that the background magnetic field has symmetry with respect to continuous isometric transformations (spatial translations, rotations, or screw motions generated by the combination of translations and rotations) in a certain direction.Owing to the use of general spatial coordinates, translation, rotation, and screw (or helical) symmetries can be treated in a unified manner.Under the symmetry assumption for the background magnetic field, there exists a certain coordinate system (x i ) i=1,2,3 where all components of the magnetic field, the vector potential, and the metric tensor become independent of one coordinate for which x 3 is chosen here without loss of generality, ∂ 3 B i = 0, ∂ 3 A i = 0, and ∂ 3 g ij = 0. (G4) For example, such coordinate systems are given by Cartesian, cylindrical (or spherical), and helical coordinate systems for the cases of translation, rotation, and screw (or helical) symmetries, respectively.Then, from Eq. (G3) with j = 3, we have Here, we define the basis vector e ≡ ∂r/∂x 3 associated with the coordinate x 3 .For translation symmetry, e represents the constant vector parallel to the direction of symmetry while, for rotation symmetry, we have e = n×r where n is the direction vector of the rotation axis.For helical symmetry, e is given by the combination of the two types of the vectors described above.Then, the boldface notation p a is used to represent the vector with the covariant components (p aj ) j=1,2,3 , and the inner product of the vectors p a and e is represented by p a •e = p aj e j = p a3 where the jth contravariant component of e ≡ ∂r/∂x 3 is given by e j = ∂x j /∂x 3 = δ j 3 .Thus, p a3 is hereafter treated as a scalar field produced by the inner product of the vectors.In the same way, using A to represent the vector with the covariant components (A j ) j=1,2,3 , we can write A 3 ≡ A • e which is treated as a scalar field produced by the inner product of A and e.
It is noted that, with respect to the spatial coordinate transformation, the distribution function F a is transformed as a scalar density field [see Eq. (B12)], which means F a / √ g is regarded as a scalar field.Since F a enters the definitions of J i [Eq.( 113)] and Θ ij [Eq.( 108)], J i and Θ ij are the components of the contravariant vector density field and the symmetric contravariant tensor density field, respectively.Correspondingly, J i / √ g and Θ ij / √ g represent the vector field and the tensor field.In the Cartesian coordinates where √ g = 1, the scalar, vector, and tensor density fields have the same components as the corresponding scalar, vector, and tensor fields so that we do not need to distinguish these density fields from the corresponding fields.Using the notation J to represent the vector field with the contravariant components (J i / √ g) i=1,2,3 , we can regard J i A 3 / √ g as the ith contravariant component of the vector field J(A • e).We also use Eqs.(A2), (A8), and the formula Γ i ij = (∂ j √ g)/ √ g to write From Eq. (A4) and ∂ 3 g ij = 0 in Eq. (G4), we find that Γ i,j3 = 1 2 (∂ j g i3 − ∂ i g j3 ) = −Γ j,i3 , which is combined with Θ ij = Θ ji to obtain Γ j i3 Θ i j = Γ j,i3 Θ ij = 0. We now use the notation Θ to represent the tensor field with the contravariant components Θ ij / √ g.Also, the tensor-vector contraction Θ • e is used to represent the vector field, the ith component of which is written as (Θ ij / √ g)e j = (Θ i j / √ g)e j = Θ i 3 / √ g.Then, using Eq.(A7) and (A8), we obtain where Γ i ij = (∂ j √ g)/ √ g and Γ j i3 Θ i j = 0 are used.Using Eq. (G5)-(G7), the local canonical momentum balance equation in the direction of symmetry is written  G8)-(G9) with the result derived from the conventional method explained in Appendix E. In deriving Eqs.(G8) and (G9), after the invariance of the Lagrangian under arbitrary spatial coordinate transformations is used to obtain the local momentum balance equation, Eq. (G1), in the general spatial coordinates, Eq. (G1) is represented in the special coordinates, one of which is the coordinate in the direction of symmetry.On the other hand, in Appendix E, the invariance of the Lagrangian under not arbitrary but only translational and rotational coordinate transformations are used to derive the local linear and angular momentum equations in Eqs.(E11) and (E14).In the case where the background magnetic field has translation symmetry along the constant direction vector e, the inner product of Eq. (E11) and e results in the local canonical momentum conservation law, ∂ ∂t where Eqs.( 57), (113), and e • ∇A = 0 are used, and the term including K a is still retained for comparison with Eq. (G9).When the background field has rotational symmetry around the axis which is parallel to the direction vector n and passes through the origin of the position vector r, A satisfies (n× r)•∇A = n× A. Then, we can use the inner product of Eq. (E14) and n to derive the local canonical angular momentum conservation law which has the same form as shown in Eqs.(G10)-(G11) where e is regarded as given by e = n × r.In the same way, it can be shown that, for the case of helical symmetry, the local conservation law of the canonical momentum in the direction of symmetry is written again by Eqs.(G10) and (G11), where we put e = n + αn × r with a constant α to represent the direction of helical symmetry.Now, it is clear that the difference between Eqs. (G9) and (G10) is expressed by that between the symmetric and asymmetric pressure tensors denoted as Θ and Π c , respectively.It is not easy to find out the well-known CGL tensor P CGL [Eq.( 109)] in Π c while P CGL naturally appears in Θ [Eq.( 108)].Incidentally, the procedure of the Belinfante-Rosenfeld type is known for derivation of the symmetric pressure tensor from the asymmetric canonical pressure tensor 39,40 although the symmetric tensor is more directly derived by the present method using the variational derivative of the Lagrangian with the metric tensor.
the radial transport of the toroidal angular momentum due to the turbulent electric field are obtained from Eqs. (D2) and (D3) as

3 = a d 3 vd 3 v
F a p a3 + ∂ i 1 c J i A 3 + Θ i K a p a3 ,(G8)in the spatial coordinates where x 3 is the coordinate in the direction of symmetry.Equation (G8) is also written using the conventional vector and tensor notations in the Cartesian coordinates (with√ g = 1)as∂ ∂t a d 3 v F a (p a • e) + ∇ • 1 c J (A • e) + Θ • e = a K a (p a • e).(G9) Finally, we clearly see that Eqs.(G8) and (G9) represent the local conservation law of the canonical momentum conjugate to the coordinate in the direction of symmetry when the collision (or source) term K a satisfies a d 3 v K a (p a • e) = 0. Let us compare Eqs. (