Effects of collisions on conservation laws in gyrokinetic field theory

Effects of collisions on conservation laws for toroidal plasmas are investigated based on the gyrokinetic field theory. Associating the collisional system with a corresponding collisionless system at a given time such that the two systems have the same distribution functions and electromagnetic fields instantaneously, it is shown how the collisionless conservation laws derived from Noether's theorem are modified by the collision term. Effects of the external source term added into the gyrokinetic equation can be formulated similarly with the collisional effects. Particle, energy, and toroidal momentum balance equations including collisional and turbulent transport fluxes are systematically derived using a novel gyrokinetic collision operator, by which the collisional change rates of energy and canonical toroidal angular momentum per unit volume in the gyrocenter space can be given in the conservative forms. The ensemble-averaged transport equations of particles, energy, and toroidal momentum given in the present work are shown to include classical, neoclassical, and turbulent transport fluxes which agree with those derived from conventional recursive formulations.


I. INTRODUCTION
Gyrokinetic theories and simulations are powerful means to investigate microinstabilities and turbulent transport processes in magnetically confined plasmas. [1][2][3][4] Originally, gyrokinetic equations are derived by recursive techniques combined with the WKB or ballooning representation. [5][6][7][8][9][10] On the other hand, modern derivations of the gyrokinetic equations are based on the Lagrangian and/or Hamiltonian formulations, 11 in which conservation laws for the phase-space volume and the magnetic moment are automatically ensured by Liouville's theorem and Noether's theorem, respectively. 12 Besides, conservation of the total energy and momentum is naturally obtained in the gyrokinetic field theory, where all governing equations for the distribution functions and the electromagnetic fields are derived from the Lagrangian which describes the whole system consisting of particles and fields. [13][14][15][16][17] A subtle point regarding the Lagrangian/Hamiltonian gyrokinetic formulations is that they basically treat collisionless systems so that Noether's theorem and conservation laws do not hold directly for collisional systems. In this paper, we examine how the collision and external source terms added into the gyrokinetic equations influence the conservation laws derived from Noether's theorem in the gyrokinetic field theory for collisionless systems.
For a given collisional kinetic system, we can imagine a corresponding collisionless kinetic system such that the two systems have the same distribution functions and electromagnetic fields instantaneously. As an example of two such systems, the Boltzmann-Poisson-Ampère system and the Vlasov-Poisson-Ampère system are considered in Sec. II, where we express the variation of the action integral for the latter collisionless system in terms of the distribution functions and the electromag-netic fields for the former collisional system to show how the conservation laws derived from Noether's theorem in the collisionless system are modified in the collisional system with external sources of particles, energy, and momentum. There, we confirm the natural result that, when adding no external sources but only the collision term that conserves the energy and momentum, the energy and momentum conservation laws for the Boltzmann-Poisson-Ampère system take the same forms as those for the Vlasov-Poisson-Ampère system. The above-mentioned procedures are repeated in Sec. III to treat the collisional and collisionless gyrokinetic systems. In our previous work, 18 using the gyrokinetic Vlasov-Poisson-Ampère system of equations, conservation laws of particles, energy, and toroidal angular momentum are obtained for collisionless toroidal plasmas, in which the slow temporal variation of the background magnetic field is taken into account in order to enable self-consistent treatment of physical processes on transport time scales. Based on these results, the particle, energy, and toroidal angular momentum balance equations for the collisional plasma are derived from the gyrokinetic Boltzmann-Poisson-Ampère system of equations in Secs. IV and V. In Sec. VI, it is shown by taking the ensemble average of these balance equations that the particle, energy, and toroidal angular momentum transport fluxes are given by the sum of the conventional expressions of the classical, neoclassical, and turbulent transport fluxes to the lowest order in the normalized gyroradius parameter. Conclusions are given in Sec. VII and formulas for transformation from particle to gyrocenter coordinates are presented in Appendix A.
Regarding the collision operator for the gyrokinetic equation, several works have been done, which take account of finite-gyroradius effects to modify the Landau collision operator defined in the particle coordinates. [19][20][21][22][23][24][25][26][27] The relation of the collision operator in the gyrocenter coordinates to that in the particle coordinates is explained in Appendix B. The Landau operator for Coulomb collisions conserves particles' number, kinetic energy, and momentum locally at the particle position although, in the gyrocenter position space, collisions induce transport fluxes of particles, energy, and momentum. Besides, it is emphasized in this work that the collisional change rates of the gyrocenter Hamiltonian (which includes not only the kinetic energy but also the potential energy) and of the canonical momentum (instead of the kinetic momentum) per unit volume in the gyrocenter space should take the conservative (or divergence) forms in order to properly derive the energy and momentum conservation laws for the collisional gyrokinetic system. The approximate collision operator which keeps these conservation properties of the gyrocenter energy and canonical toroidal angular momentum is shown in Appendix C. [It is noted that another form of the gyrokinetic collision operator, which satisfies the energy and momentum conservation laws, has recently been presented by Burby et al. 26 ] Appendix D is given to describe how to derive the formula for the toroidal angular momentum transport flux due to the collision term.

II. BOLTZMANN-POISSON-AMPÈRE SYSTEM
In this section, conservation laws are investigated for the Boltzmann-Poisson-Ampère system of equations which provide the basis of approximate description by the collisional electromagnetic gyrokinetic system of equations for strongly magnetized plasmas considered in Secs. III-VI. Time evolution of the distribution function f a (x, v, t) for particle species a is described by the Boltzmann kinetic equation, where K a (x, v, t) denotes the rate of change in the distribution function f a due to Coulomb collisions and it may also include other parts representing external particle, momentum, and/or energy sources if any. The electromagnetic fields E(x, t) and B(x, t) are written as E = −∇φ − c −1 ∂A/∂t and B = ∇ × A, where the electrostatic potential φ and the vector potential A are determined by Poisson's equation, and Ampère's law, respectively. Here, the Coulomb (or transverse) gauge condition ∇ · A = 0 is used and the current density j ≡ a e a n a u a ≡ a e a f a (x, v, t)vd 3 v (or any vector field) is written as represent the longitudinal (or irrotational) and transverse (or solenoidal) parts, respectively. 28 Equations (1), (2), and (3) are the governing equations for the Boltzmann-Poisson-Ampère system. Suppose that f a , φ, and A which satisfy Eqs. (1)-(3) are given. Then, for the electromagnetic fields E = −∇φ − c −1 ∂A/∂t and B = ∇ × A given from φ and A, we consider the distribution function f V a which is the solution of the Vlasov equation, We also assume f V a to coincide instantaneously with f a at a given time . Therefore, equations obtained from Eqs. (2) and (3) with f a replaced by f V a also hold at t 0 . In other words, f V a , φ, and A satisfy the Vlasov-Poisson-Ampère system of equations at t 0 . Note that the Vlasov-Poisson-Ampère system of equations can be derived from the variational principle using the action I defined by Eq. (1) in Ref. 29 where its variation δI associated with infinitesimal transformations of independent and dependent variables [see Eqs. (10) and (15) in Ref. 29] are explicitly shown in order to apply Noether's theorem for obtaining conservation laws of energy and momentum. Now, let us use f V a , φ, and A to define the action integral I over a small time interval, t 0 − h/2 ≤ t ≤ t 0 + h/2, during which the Vlasov-Poisson-Ampère system of equations are approximately satisfied by f V a , φ, and A within the errors of order h. Then, neglecting the errors of higher order in h, the variation δI can be written in the same form as in Eq. (15) of Ref. 29, where δG V 0 and δG V are written as respectively, where the superscript T represents the transpose of the tensor and I denotes the unit tensor. The field variable λ which appears in Eq. (7) Suppose that the variations t E , δx E , δφ, and δA are such that δI = 0 holds for an arbitrary x-integral domain in Eq. (5). Then, taking the small time interval limit h → +0 in Eq. (6), we find that the conservation law, should be satisfied. This is the so-called "Noether's the- We now define E c , P c , Q c , and Π c from E V c , P V c , Q V c , and Π V c , respectively, by replacing f V a with f a in Eq. (7). Correspondingly, δG 0 and δG are defined from δG V 0 and δG V by replacing E V c , P V c , Q V c , and Π V c with E c , P c , Q c , and Π c , respectively, in Eq. (7). These definitions immediately yield δG V (x, t 0 ) = δG(x, t 0 ) and where Eq. (10) is used and δK G0 is defined by Substituting Eq. (11) into Eq. (9), we find that the conservation law is modified for the Boltzmann-Poisson-Ampère system as where t 0 is rewritten as t because t 0 is an arbitrarily chosen time. Equation (13) shows that δK G0 represents effects of K a in Eq. (1) on the conservation law. If K a is given by the Coulomb collision term only, δK G0 defined by Eq. (12) vanishes because the collision term conserves particles' number, momentum, and energy. Energy and momentum balance equations can be derived from Eq. (13) using symmetries of the system under infinitesimal time and space translations as shown later. Before deriving them, we first consider the equation for the particle number density n a ≡ f a d 3 v which is obtained by taking the velocity-space integral of Eq. (1) as We hereafter assume that a e a K a d 3 v = 0, which means that the source terms K a conserve electric charge even if K a d 3 v = 0 for each species a. This seems a reasonable assumption in consistency with Eqs. (2) and (3) in which no external source terms are included. Then, multiplying Eq. (14) with the electric charge e a and performing the summation over species result in the charge conservation law, We also find from a e a K a d 3 v = 0 that, in Eq. (12), the terms including φ and A vanish and make no contribution to K Ec and K P c . As seen from Eqs. (2), (8), and (15), we can put λ = ∂φ/∂t which is also used in Ref. 29 to derive energy and momentum conservation laws for the Vlasov-Poisson-Ampère system. We now note that the action integral is invariant, namely, δI = 0 under the infinitesimal translations in space and time represented by δt E = ǫ 0 , δx a = δx E = ǫ, δv a = 0, δφ = 0, and δA = 0, where ǫ 0 and ǫ are constant in time and space. These invariance properties hold because the integrands in the action integral I where the canonical energy density and flux (E c , Q c ) are given by replacing f V a with f a in the definitions of (E V c , Q V c ) in Eq. (7). In the same way as in Ref. 29, we use the kinetic energy density and flux, (E p , Q p ), defined by to modify Eq. (16) into more familiar forms. Then, the energy balance equation is finally written as where E L ≡ −∇φ and E T ≡ −c −1 ∂A/∂t are the longitudinal and transverse parts of the electric field, respectively.
Next, from the space translational symmetry and Eq. (13), we obtain the canonical momentum balance equation, where the canonical momentum density and tensor (P c , Π c ) are given by replacing f V a with f a in the definitions of (P V c , Π V c ) in Eq. (7). Furthermore, in the same way as in Ref. 29, the invariance of I under the infinitesimal rotation is shown to give the equation for the angular momentum, which is used to modify Eq. (19) into the momentum balance equation, (20) Here, the particle parts (P p , Π p ) of the momentum density and the pressure tensor are defined by and the field parts (P f , Π f ) are given by Equations (18) and (20) take physically familiar forms of energy and momentum balance equations including external source terms. As mentioned earlier, if K a is given by the Coulomb collision term only, K Ec and K P c vanish so that the energy and momentum balance equations for the Boltzmann-Poisson-Ampère system take the same forms as those for the Vlasov-Poisson-Ampère system. 29

III. GYROKINETIC BOLTZMANN-POISSON-AMPÈRE SYSTEM
Let us start from the gyrokinetic Boltzmann equation written as where F a (Z, t) is the gyrocenter distribution function for species a, C g ab [F a , F b ](Z, t) represents the rate of change in F a (Z, t) due to Coulomb collisions between particle species a and b, and S a (Z, t) denotes other parts including external particle, momentum, and/or energy sources if any. The gyrocenter coordinates are written as Z a = (X a , U a , µ a , ξ a ), where X a , U a , µ a , and ξ a represent the gyrocenter position, parallel velocity, magnetic moment, and gyrophase angle, respectively. Appendix A shows the relation of the gyrocenter coordinates to the particle coordinates in detail. The perturbation expansion parameter in the gyrokinetic theory is denoted by δ which represents the ratio of the gyroradius ρ to the macroscopic scale length L of the background field. It is shown in Appendix B how the collision operator C g ab [F a , F b ] for the gyrocenter distribution functions F a and F b is given from the collision operator C p ab [f a , f b ] for the particle distribution functions f a and f b .
The deviation of each distribution function from the local Maxwellian is regarded as of O(δ), and accordingly the collision term C g ab is considered to be of O(δ). We assume that the source term S a is of O(δ 2 ) so that its effect appears only in the transport time scale. We also assume that a e a dU dµ dξ D a S a (Z, t) = 0 in order to prevent the source term from affecting the charge conservation laws [see Eq. (60)]. Here, D a denotes the Jacobian for the gyrocenter coordinates, D a ≡ det[∂(x a , v a )/∂(X a , U a , ξ a , µ a )], where (x a , v a ) represent the particle coordinates consisting of the particle's position and velocity vectors.
We treat toroidal systems, for which the equilibrium magnetic field is given in the axisymmetric form as where I and χ are constant on toroidal flux surfaces labeled by an arbitrary radial coordinate s and ζ is the toroidal angle. We note that I and χ represent the covariant toroidal component of the equilibrium field B 0 and the poloidal magnetic flux divided by 2π, respectively. The equilibrium field B 0 is allowed to be dependent on time. Then, following Ref. 18, the gyrocenter motion equations are written as where the gyrocenter Hamiltonian, which is independent of ξ a , is defined by and A * a is given by A * a = A 0 (X a , t) + (m a c/e a )U a b(X a , t).
Using the nonvanishing components of the Poisson brackets for pairs of the gyrocenter coordinates given by the gyrocenter motion equations in Eq. (25) are rewritten as and Here, Ω a ≡ e a B 0 /(m a c), b = B 0 /B 0 , B * a = B * a · b, and B * a = ∇ × A * a . The field variable Ψ a is defined by where the field variable ψ a is defined in terms of the electrostatic potential φ and the perturbation part of the vector potential A 1 as The gyroradius vector is given by ρ a = b(X a , t) × v a0 (Z a , t)/Ω a (X a , t) and the zeroth-order particle velocity v a0 is written in terms of the gyrocenter coordinates as where the unit vectors (e 1 , e 2 , b) form a right-handed orthogonal system.
The gyrophase-average and gyrophase-dependent parts of an arbitrary periodic function Q(ξ a ) of the gyrophase ξ a are written as respectively. In the gyrocenter motion equations, effects of the time-dependent background magnetic field and those of the fluctuating electromagnetic fields appear through ∂A * a /∂t and Ψ a , respectively. It should be noted that dZ a /dt on the left-hand side of Eq. (23) is regarded as a function of (Z, t) which is given by the right-hand side of Eq. (25).
We find from Eqs. (A3)-(A4) in Appendix A and Eq. (B1) in Appendix B that the gyrophase-dependent part of the right-hand side of Eq. (23) appears from C g ab and it is of O(δ). Using Ω a = O(δ −1 ), the gyrophasedependent part of the left-hand side of Eq. (23) is written as Ω a ∂ F a /∂ξ to the lowest order in δ. Then, it is concluded that F a = O(δ 2 ). Taking the gyrophase average of Eq. (23), we obtain where F a (Z, t) and S a (Z, t) are both regarded as independent of the gyrophase ξ and · · · ξ are omitted from them for simplicity. It is seen from Eq. (B1) that ef- Here and hereafter, we neglect F a = O(δ 2 ) in both sides of the gyrokinetic Boltzmann equation given by Eq. (35). Even so, its moment equations can correctly include the collisional transport fluxes of particles, energy, and toroidal momentum up to the leading order, that is O(δ 2 ), as confirmed later. In Appendix C, Eq. (C1) combined with Eqs. (C2), (C10), (C15), and (C16) presents the approximate gyrokinetic collision operator, which has favorable conservation properties and correctly describes collisional transport of energy and toroidal angular momentum.
The gyrokinetic Poisson equation and the gyrokinetic Ampère's law are written as 18 and respectively, where (j G ) T is the transverse part of the gyrokinetic current density j G defined by It should also be noted that the Coulomb gauge conditions ∇ · A 0 = 0 and ∇ · A 1 = 0 for the equilibrium and perturbation parts of the vector potential are used here. The equilibrium vector potential A 0 is given by 18, additional governing equations are derived in order to self-consistently determine I and χ for the time-dependent axisymmetric background field B 0 = ∇ × A 0 = I∇ζ + ∇ζ × ∇χ. They are given by and where the toroidal-angle average is represented by · · · ≡ (2π) −1 · · · dζ, the poloidal angle is denoted by θ, and Here, the covariant toroidal component of an arbitrary vector V is written as V ζ . The turbulent part of the magnetic field is given by For the gyrocenter coordinates which have Poisson brackets given by Eq. (27), the Jacobian is given by D a = B * a /m a . It is important to note that the Jacobian D a satisfies the gyrocenter phase-space conservation law, Then, using Eq. (41), the gyrokinetic Boltzmann equation in Eq. (35) can be rewritten as where K a is the gyrophase-independent function given by the right-hand side of Eq. (35), We hereafter derive conservation laws for the gyrokinetic Boltzmann-Poisson-Ampère system of equations following the procedures similar to those shown in Sec. II. For that purpose, suppose that F a , φ, A 1 , I, and χ satisfy Eqs. (35), (36), (37), (39), and (40). Then, we consider the gyrocenter distribution function F V a which obeys the gyrokinetic Vlasov equation, where dZ a /dt is evaluated by using the above-mentioned fields (φ, A 1 , I, χ) obtained from the solution of the gyrokinetic Boltzmann-Poisson-Ampère system of equations. Here, it should be noted that, if the distribution functions F a and F V a , which are given as the solutions of Eqs. (35) and (44), respectively, are initially gyrophase-independent, they are gyrophase-independent at any time. Besides, F V a is assumed to coincide instantaneously with F a at a given time t 0 . Therefore, Eqs. (36), (37), (39), and (40) are all satisfied at that moment even if F a is replaced with F V a in these equations. Thus, the gyrokinetic Vlasov-Poisson-Ampère system of equations are instantaneously satisfied by (F V a , φ, A 1 , I, χ) at t = t 0 . In Ref. 18, the action integral I is defined to derive all the governing equations for the gyrokinetic Vlasov-Poisson-Ampère system based on the variational principle, and its variation δI associated with the infinitesimal variable transformations are given to obtain conservation laws from Noether's theorem. Here, the action integral I can be expressed in terms of (F V a , φ, with the functions δG V 0 and δG V defined by where E V c and P V c are defined by and definitions of other variables S χ , and δT V are shown in Eq. (79) of Ref. 18. The superscript V in the variables (E V c , P V c , · · · ) implies that they are defined using the distribution function F V a instead of F a .
As explained in Ref. 18, the integral domain of Eq. (45) is not an arbitrary local one in the X-space but it can be local only in the radial direction in order for Eq. (45) to be valid. Then, if the variations δt E , δx E , · · · in Eq. (46) are such that δI = 0 holds for a spatiotemporal integral domain defined by where [s 1 , s 2 ] represents an arbitrary spatial volume region sandwiched between two flux surfaces labeled by s 1 and s 2 , then the conservation law is derived as This is Noether's theorem for the gyrokinetic Vlasov-Poisson-Ampère system. In Eq. (48), V ′ ≡ ∂V /∂s represents the derivative of V (s, t) with respect to s and V (s, t) denotes the volume enclosed by the flux surface with the label s at the time t.
where K is defined by Eq. (43). Let us also define δG 0 and δG from δG V 0 and δG V by replacing F V a with F a . Then, we have G V (X, t 0 ) = G(X, t 0 ) and Here, p c a denotes the canonical momentum for species a defined by Substituting Eq. (50) into Eq. (48) and rewriting the arbitrarily chosen time t 0 as t, we obtain the conservation law for the gyrokinetic Boltzmann-Poisson-Ampère system, where δK G0 represents effects of the collision and source terms on the conservation law. Under the nonstationary background field B 0 , flux surfaces may change their shapes and the grid of the flux coordinates moves. Then, Eq. (53) is rewritten as

IV. EQUATIONS FOR GYROCENTER DENSITIES AND POLARIZATION
In this section, we take the velocity-space integral of the gyrokinetic Boltzmann equation in Eq. (35) to consider the particle transport before treating the energy and toroidal angular momentum transport in Sec. V. We define the gyrocenter density n and the gyrocenter flux where u (gc) a represents the gyrocenter fluid velocity and the gyrocenter drift velocity v (gc) a = dX a /dt is given by evaluating the right-hand side of Eq. (28) at (X, U, µ). Then, integrating the gyrokinetic Boltzmann equation, Eq. (42), with respect to the gyrocenter velocity-space coordinates (U, µ, ξ) and using Eq. (43), we obtain (57) Using the approximate collision operator given by Eq. (C1) in Appendix C, the particle flux Γ C a due to collisions and finite gyroradii is defined by Eq. (C6) with putting A p a (z) = 1. [If the collision operator given by As shown in Ref. 18, the gyrokinetic Poisson equation in Eq. (36) is rewritten as a e a n (gc) where E L = −∇φ, ∇ = ∂/∂X, and P (pol) represents the polarization density defined by Here, ρ ai denotes the ith Cartesian component of /Ω a (X, t), and F * a = F a + (e a ψ a /B 0 )(∂F a /∂µ).
As mentioned before Eq. (24) in Sec. III, a e a dU dµ dξ D a S a (Z, t) = 0 is assumed. Then, using Eq. (57), we can obtain the charge conservation law, ∂ ∂t a e a n (gc) where the current density due to the gyrocenter drift and that due to the collisional particle transport are given by and j C = a e a Γ C a , respectively. Note that the magnetization current is solenoidal and, accordingly, it does not contribute to the charge conservation law in Eq. (60). Equation (58) is substituted into Eq. (60) to show where the subscript L is used to represent the longitudinal part of the vector variable. Then, using Eqs. (58), (60) and (61), we find that the useful formula, holds for any function A(X, t). The relation in Eq. (62) is used in Sec. V.B to derive Eq. (74).

V. GYROKINETIC ENERGY AND TOROIDAL ANGULAR MOMENTUM BALANCE EQUATIONS
In this section, energy and toroidal angular momentum balance equations for gyrokinetic systems including collisional processes are derived by using the results shown in Secs. III-IV.

A. Energy balance equation
The variation δI of the action given in Sec. III vanishes under the infinitesimal time translation represented by δt E = ǫ where ǫ is an infinitesimally small constant. Here, all other infinitesimal variations δx E , δφ, · · · are regarded as zero. Then, δG 0 and δG are determined by these conditions for the infinitesimal time translation and they satisfy Eq. (54) which, in the same manner as in Ref. 18, leads to the energy balance equation, (63) Here, the energy density E is defined by and the energy fluxes Q * c and Q * R are given by and where the energy flux Q C due to collisions and finite gyroradii is defined by taking the summation of Eq. (C6) over species a with putting A p a (z) = 1 2 m a v 2 a + µ 0a B 0 (x a ) + e a φ(x a ). The derivation of Eq. (67) requires Eq. (C7) which is satisfied by the collision operator in the form of Eq.(C1) with appropriately choosing ∆x (2) a , ∆v a , and ∆µ 0a as described in Appendix C.
Substituting Eq. (67) into Eq. (63), the energy balance equation is rewritten as where the energy flux Q is given by The right-hand side of Eq. (68) represents the external energy source. It is confirmed later in Sec. VI.B that the ensemble average of Q · ∇s coincides with the wellknown expression of the radial energy transport to the lowest order in the δ-expansion.

B. Toroidal angular momentum balance equation
The toroidal angular momentum balance equation is derived from the fact that δI = 0 under the infinitesimal toroidal rotation represented by δx E = ǫe ζ (X). Here, ǫ is again an infinitesimally small constant, and e ζ (X) is defined by e ζ (X) = ∂X(R, z, ζ)/∂ζ = R 2 ∇ζ where the right-handed cylindrical spatial coordinates (R, z, ζ) are used. We also defineẑ byẑ = R∇ζ × ∇R which represents the unit vector in the z-direction. Then, if putting the origin of the position vector X at (R, z) = (0, 0), we have e ζ (X) = X ×ẑ. Under the infinitesimal toroidal rotation, the variations of the vector variables are given as δA 1 = ǫA 1 ×ẑ and δA 0 = ǫA 0 ×ẑ although the other variations δt E , δφ, · · · , are all regarded as zero. Then, using these variations of the variables associated with the infinitesimal toroidal rotation, the canonical momentum balance equation is derived from Eq. (53) as Here, the density of the canonical toroidal angular momentum is defined by with the toroidal component of the canonical momentum denoted by the variation of the canonical toroidal angular momentum due to collisions and external sources is given by We follow the same procedures as shown in Sec. V.B of Ref. 18 and use Eqs. (70)-(73) and Eq. (62) with A = A 0ζ = −χ to write the toroidal angular momentum balance equation as where Using Eqs. (73) and (C5) in Appendix C with putting A g a = (p c a ) ζ , the right-hand side of Eq. (74) is rewritten as where J C pζ is defined by taking the summation of where the right-hand side represents the external source of the toroidal angular momentum and is the radial flux of the toroidal angular momentum due to collisions and finite gyroradii. In Sec. VI.C, we derive the ensemble-averaged toroidal angular momentum balance equation from Eq. (77) in order to confirm that it is consistent with the conventional result up to the second order in δ.

VI. ENSEMBLE-AVERAGED BALANCE EQUATIONS FOR PARTICLES, ENERGY, AND TOROIDAL ANGULAR MOMENTUM
In this section, the particle, energy and toroidal angular momentum balance equations derived in Secs. IV and V are ensemble-averaged for the purpose of verifying their consistency with those obtained by conventional recursive formulations. 30,[33][34][35] In the same way as shown in Sec. VI of Ref. 18, we divide an arbitrary physical variable Q into the average and turbulent parts as where · · · ens represents the ensemble average, and we immediately find Q ens = 0. We identify the zeroth fields A 0 and B 0 with the ensemble-averaged parts to write A 0 = A ens , A 1 =Â, B 0 = B ens , and B 1 =B. Regarding the electrostatic potential φ, it is written as the sum of the average and fluctuation parts, φ(x, t) = φ(x, t) ens +φ(x, t). Here, assuming that φ(x, t) ens = 0, the background E × B flow is retained and its velocity is regarded as O(δv T ), where δ and v T represent the drift ordering parameter and the thermal velocity, respectively. Then, using Eq. (33), we have ψ a = ψ a ens +ψ a , where We assume that the ensemble average Q ens of any variable Q considered here has a slow temporal variation subject to the so-called transport ordering, ∂ ln Q ens /∂t = O(δ 2 v T /L), and that it has a gradient scale length L which is on the same order as gradient scale lengths of the equilibrium field and pressure profiles. We also impose the constraint of axisymmetry on Q ens that is written as ∂ Q ens /∂ζ = 0 even though Q itself is not axisymmetric. On the other hand, the turbulent partQ of Q is assumed to vary with a characteristic frequency ω = O(v T /L) and have gradient scale lengths L and ρ in the directions parallel and perpendicular to the equilibrium magnetic field B 0 , respectively.
The ensemble-averaged part F a ens of the distribution function F a for species a consists of the local Maxwellian part and the deviation from it, The local Maxwellian distribution function is written as where the equilibrium density n a0 and temperature T a0 are regarded as uniform on flux surfaces. The first-order ensemble-averaged distribution function F a1 ens is determined by the drift kinetic equation, which can be derived by substituting Eq. (81) into the ensemble average of Eq. (23). The derived equation agrees, to O(δ), with the well-known linearized drift kinetic equation, on which the neoclassical transport theory is based. 31,32 The fluctuation partF a is written aŝ Substituting Eq. (82) into the fluctuation part of the gyrokinetic equation in Eq. (23) yields where C L a represents the linear collision term defined by It is shown by using Eq. (82) and the WKB representation that, to the lowest order in δ, the turbulent parts of Eqs. (36) and (37)  A. Ensemble-averaged particle balance equation Taking the ensemble average of Eq. (57) and subsequently its flux surface average, we obtain ∂ n (gc) a ens ∂t where n (gc) a ens = n a0 + O(δ), and · · · represents a double average over the flux surface and the ensemble. Here, n a0 is the equilibrium density which is a fluxsurface function and characterizes the Maxwellian distribution function F aM . On the right-hand side, the source term S a is regarded as of O(δ 2 ) as well as all other terms in Eq. (86), and it is assumed to have no turbulent component so that S a = S a ens .
It is shown in Ref. 18 that the radial gyrocenter particle flux is given by where the nonturbulent part (Γ NA a ) s and the turbulencedriven part (Γ A a ) s are written as and respectively. On the right-hand side of Eq. (88), P CGL a1 represents the first-order part of the pressure tensor in the Chew-Goldberger-Low (CGL) form 31 defined by P CGL a1 = dU dµ dξ D a F a1 ens m a U 2 bb + µB 0 (I − bb) , and the ensemble-averaged electric field is given by E ens = −∇ φ ens − c −1 ∂A 0 /∂t. Thus, Eq. (88) expresses the neoclassical radial particle flux and the radial E × B drift which are well-known by the neoclassical transport theory. 31 We also find that Eq. (89) agrees with the turbulent radial particle flux derived from the conventional gyrokinetic theory based on the WKB formalism. 30 The radial classical particle flux is given by where F a1 ≡ d 3 vm a vC p a is the collisional friction force. It is well-known that the classical transport equation relating (Γ C a ) s to the gradient forces is immediately derived from Eq. (90) because the first-order gyrophasedependent part of the particle distribution function in Eq. (90) is expressed in terms of the gradient of the background Maxwellian distribution function as f a1 = −ρ a · ∇F aM with the gradient operator ∇ taken for the fixed energy variable ε = 1 2 m a v 2 + e φ ens . In the same manner as in deriving Eq. (54) from Eq. (53), the ensemble-averaged particle transport equation can be obtained from Eq. (86) as where the total radial particle flux is given by the sum of the classical, neoclassical, and turbulent parts as As shown above, the well-known expressions of the classical, neoclassical and turbulent particle fluxes are included in (Γ C a ) s , (Γ NA a ) s , and (Γ A a ) s , respectively. The latter two fluxes are evaluated by the solutions F a1 ens andĥ a of the first-order drift kinetic and gyrokinetic equations, respectively.

B. Ensemble-averaged energy balance equation
The ensemble average of the energy density defined by Eq. (64) is written as where the energy density of the electric field is neglected as a small quantity of O(δ 2 ). It is shown in Ref. 18 that the radial components of the first two terms on the right-hand side of Eq. (65) are double-averaged over the ensemble and the flux surface to give Here, the radial particle flux (Γ (gc) a ) s is given by Eqs. (87)-(89) and the radial heat flux (q a ) s is written as which consists of the nonturbulent part, and the turbulence-driven part, In Eq. (96), the heat stress tensor Θ CGL a is defined by The expression of Eq. (96) coincides with that of the neoclassical radial heat flux in terms of the heat stress tensor. 31 The turbulent heat flux in Eq. (97) takes the same form as that given by the conventional gyrokinetic theory. 30 The radial component of Q C in Eq. (69) is ensembleaveraged to yield where the radial classical heat flux is given by Here, F a2 ≡ d 3 v(m a v 2 /2T a − 5/2)m a vC p a is the collisional heat friction. The expression of the classical heat flux (q C a ) s in Eq. (99) agrees with the conventional one 31 and it immediately gives the classical heat transport equation relating (q C a ) s to the gradient forces in the same way as mentioned after Eq. (90) for the classical particle flux (Γ C a ) s .

Now, Eq. (68) is rewritten as
where the total radial heat flux is given by the sum of the classical, neoclassical, and turbulent parts as and S (Poynting) ≡ (c/4π) E ens × B 0 represents the nonturbulent part of the Poynting vector. Using the relation Equations (100) and (102) take the well-known forms of the energy balance equations 32 except that the terms associated with the electric field energy and the kinetic energies due to the fluid velocities are neglected here as small quantities of higher order in δ.
C. Ensemble-averaged toroidal angular momentum balance equation The ensemble-averaged toroidal angular momentum balance equation is written as where u a represents the nonturbulent part of the parallel fluid velocity for particle species a defined by n a0 u a ≡ dU dµ dξ F a1 ens U and u E ≡ c E ens × b/B 0 is the nonturbulent part of the E × B drift velocity. Equation (103) is derived from Eq. (77) following the same procedures as shown in Ref. 18 except that the additional transport flux Π C a defined in Eq. (107) and the external momentum source are newly included in the present case.
On the left-hand side of Eq. (103), the terms including (u a b + u E ) and S (Poynting) are of O(δ 3 ) although they are written down to explicitly show the inertia-term part. The nonturbulent and turbulence-driven parts of the radial flux of the toroidal angular momentum are defined by and respectively. It is shown in Appendix D that Eq. (78) is ensemble-averaged to give where the radial transport flux of the toroidal angular momentum for species a due to the collision term and finite gyroradii is defined by The expressions for the toroidal momentum fluxes shown in Eqs. (104)-(107) agree with those given by conventional recursive formulations in Refs. 33-35. [Since the so-called high-flow ordering is used in Refs. 33 and 34, the expressions for the toroidal momentum fluxes in it reduce to those in the present work in the low-flow-speed limit.] As argued in Refs. 18 and 35, when there exists the up-down symmetry of the background magnetic field, all toroidal momentum fluxes vanish to O(δ 2 ) and the nontrivial toroidal momentum balance equation is of O(δ 3 ). In this case, gyrokinetic systems equations of higher-order accuracy in δ are required for the correct derivation of this O(δ 3 ) toroidal momentum balance equation to determine the profile of the radial electric field 36 although we should note, at the same time, that the radial electric field is not necessary to determine the particle and energy transport fluxes to the lowest order in δ. 35

VII. CONCLUSIONS
In this paper, particle, energy, and toroidal momentum balance equations including collisional and turbulent transport fluxes are systematically derived from the gyrokinetic Boltzmann-Poisson-Ampère system of equations. Considering an imaginary collisionless system, for which the distribution functions and electromagnetic fields coincide instantaneously with those for the considered collisional system, and expressing the variation of the action integral for the collisionless system in terms of the solution to the governing equations for the collisional system clarify effects of the collision and external source terms on the collisionless conservation laws derived from Noether's theorem. The gyrokinetic collision operator is newly presented, by which the collisional changes in the velocity-space integrals of the gyrocenter Hamiltonian and the canonical toroidal angular momentum can be written in the conservative (or divergence) forms. It is confirmed that, to the lowest order in the normalized gyroradius, the ensemble-averaged fluxes in the derived particle, energy, and toroidal angular momentum balance equations can be written by the sum of conventional expressions of classical, neoclassical, and turbulent transport fluxes. The extension of the present work to the case of the high-flow ordering remains as a future task. We consider the transformation of the phase-space coordinates in this Appendix, where the subscript representing the particle species is omitted as far as it is unnecessary. In terms of the position x and the velocity v of a given particle, we define the parallel velocity v = v · b(x, t), the perpendicular velocity v ⊥ = v − v b, and the zeroth-order magnetic moment, where the equilibrium field at position x and time t is denoted by B 0 (x, t) = B 0 (x, t)b(x, t). We also define the zeroth-order gyrophase by where (e 1 , e 2 , b) are unit vectors which form a righthanded orthogonal system at (x, t). Then, the gyrocenter coordinates Z = (X, U, µ, ξ) are represented in terms of the particle coordinates z = (x, v , µ 0 , ξ 0 ) as The formulas for ∆v ∆µ 0 , and ∆ξ 0 in Eq. (A3) are obtained by combining the guiding center and gyrocenter coordinate transformations. 11,13,37 Here, effects of the background electric field and turbulent electromagnetic fields are included through ψ [see Eq. (80)] and A 1 . When the background electric field and turbulent electromagnetic fields vanish, ∆v ∆µ 0 , and ∆ξ 0 in Eq. (A3) agree with the results in Ref. 37. Denoting the coordinate transformation by T , Eq. (A2) is rewritten as An arbitrary scalar field A on the phase space can be expressed in terms of either the gyrocenter coordinates Z = (X, U, µ, ξ) or the particle coordinates z = (x, v , µ 0 , ξ 0 ) as Using Eqs. (A4), (A5), and the Taylor series expansion, we obtain where T * A g denotes the pullback transformation of A g by T . Using the inverse transformation T −1 , we also The Jacobians D p and D g for the two coordinate systems z and Z are related to each other by where ∂(Z)/∂(z) denotes the Jacobian matrix. Then, we use the following formula, and partial integrals to derive the relation between the expressions of the scalar density DA in the gyrocenter and particle coordinate systems as where the replacement of z with Z is represented by Eq. (3.22) in Ref. 32]. Then, the collision term C g ab represented in the gyrocenter coordinates is related to C p ab by where the distribution function for species a (b) in the particle coordinates is written as the pullback of that in the gyrocenter coordinates F a (F b ) by the coordinate transformation T a (T b ) described in Appendix A, and T −1 * a transforms the collision term as a function of the particle coordinates into that of the gyrocenter coordinates.
In order to see collisional effects on conservation laws, it is convenient to represent the collision term in the gyrocenter coordinate using the transformation formula for the scalar density D a C ab rather than that for the scalar C ab shown in Eq. (B1). Using Eq. (A9), we can derive where A a is an arbitrary scalar field depending on particle species and f a = T * a F a is rewritten by using Eq. (A6) as Then, the gyrocenter representation of the collision operator C g ab acting on F a and F b is obtained by Eq. (B2) with putting A g = A p = 1 and using Eq. (B3) to express f a and f b in terms of F a and F b , respectively. Integrating Eq. (B2) with respect to (U, µ, ξ) and taking the summation over species b yield where C g a = b C g ab and ∇ = ∂/∂X are used and d 3 v = dv dµ 0 dξ 0 D p a (z) denotes the velocityspace integral using the particle coordinates. Here, the transport flux J C Aa of the quantity A a due to collisions and finite gyroradii of particles is defined by The integral of an arbitrary scalar field A a over the whole phase space is written in either the gyrocenter or particle coordinate system as (B6) For the case of A a = 1, Eqs. (B4) and (B5) reduce to respectively, where d 3 v C p a (z) = 0 is used. Here, Γ C a is regarded as the classical particle flux which occurs due to collisions and finite gyroradii. In fact, using ∆x a ≃ −ρ a , we see that the primary term of Γ C a shown in the last line of Eq. (B8) is identical to the conventional definition of the classical particle flux Γ cl a ≡ (c/e a B 0 )F a1 × b, where F a1 ≡ d 3 vm a vC p a is the collisional friction force. Thus, . Let us take the kinetic energy of the particle as A a and put A p a = 1 2 m a v 2 a = 1 2 m a v 2 a + µ 0a B 0 (x a ). Then, it is written in terms of the gyrocenter coordinates as where the inverse T −1 a of the transformation T a given by Eq. (A3) is used. In this case, taking the summation of Eq. (B4) over species a and using the conservation where Q C represents the transport flux of the total kinetic energy due to collisions and finite gyroradii defined by To the lowest order in δ, the collisional energy flux Q C is approximately written as Q C ≃ a (q cl a + 5 2 T a Γ cl a ). Here, the classical heat flux for species a is defined by q cl a is the collisional heat friction. We note from Eq. (B9) that the expression of the kinetic energy in the gyrocenter coordinates should be generally given by the infinite series expansion in δ in order for the gyrocenter velocity-space integral of the collisional rate of change in the kinetic energy to take the form of the divergence of the energy flux without any local source or sink terms. In fact, this energy conservation property is broken if we keep only the lowest order terms in Eq. (B9) and evaluate the gyrocenter velocity-space integral a dU dµ dξ D g a (Z)C g a (Z) 1 2 m a U 2 + µB 0 (X) . The above-mentioned subtle relation between expressions of the collisional energy conservation properties in the particle and gyrocenter coordinate systems is also found when considering the collisional momentum conservation. It should be recalled that the perturbative expansions in δ are truncated up to finite orders in deriving gyrokinetic equations as shown in Sec. III although the conservative form of equations for the energy and the toroidal angular momentum are obtained even from these approximate equations for the collisionless case since they are constructed based on the variational principle. Thus, from the viewpoint of practical applications, it is desirable for the approximate collision operator in the gyrocenter coordinates to keep the conservation properties. More rigorously speaking, we want the gyrokinetic collisional velocityspace integrals a dU dµ dξ D g a (Z)C g a (Z)H a (Z) and a dU dµ dξ D g a (Z)C g a (Z)(p c ζ ) g a (Z) to take the divergence forms and include no local source or sink terms where H a (Z) and (p c ζ ) g a (Z) ≡ (e a /c)A * aζ (Z) are the gyrocenter Hamiltonian and the canonical toroidal angular momentum defined by Eqs. (26) and (72), respectively. Here, it should be noted that not only kinetic parts of energy and toroidal momentum but also contributions from scalar and vector potentials are included in H a (Z) and (e a /c)A * aζ (Z). In Appendix C, we find how to construct the approximate gyrokinetic collision operator, by which the two integrals mentioned above are written in the divergence forms.
We now consider the entropy per unit volume defined in terms of the gyrocenter distribution functions as S g ≡ − a dU dµ dξD g a (Z) log F a (Z), in which the rate of change is given by dS g /dt = − a dU dµ dξ D g a (Z) [log F a (Z) + 1] (dF a /dt). Then, the rate of change in S g due to collisions is obtained by putting A g a = −[log F a (Z) + 1] in Eq. (B4) and taking the summation over species a as where d 3 v C p a (z) = 0 is used although we should recall that dU dµ dξD g a (Z)C g a (Z) does not vanish generally as seen from Eq. (B7). It is well-known that, when Landau's collision operator is used for C p a , the collisional entropy production rate given by the first term on the right-hand side of Eq. (B12) is nonnegative. This is Boltzmann's H-theorem which proves the second law of thermodynamics. The collisional transport flux J C S of the entropy in Eq. (B12) is defined by It is shown that, to the lowest order in δ, the collisional entropy transport flux is written as J C S = a (S a0 u cl a + q cl a /T a ) where the lowest-order entropy density S a0 for species a is given in terms of the local Maxwellian distribution function F aM as S a0 ≡ − dU dµ dξ F aM log F aM , and u cl a is defined by u cl a ≡ Γ cl a /n a . Here, we note again that the infinite series expansion in δ as given by Eq. (B2) is used in deriving Eq. (B12). When the expansion is truncated to finite order, the collisional entropy production term is represented by − a d 3 v C p a log f a plus residual error terms of higher order in δ, and thus the H-theorem is only approximately satisfied. In this Appendix, we consider an approximate gyrokinetic collision operator instead of the one given by Eq. (B1) [or Eq. (B2) with A g = A p = 1] in order to get the gyrokinetic collisional velocity-space integrals of energy and canonical toroidal momentum to take desirable conservative (or divergence) forms. The approximate collision operator is written in the gyrocenter coordinates as where ∆v a and ∆µ 0a are written as a , ∆µ 0a = ∆µ (1) 0a + ∆µ (2) 0a , and f a is given from F a by f a (z a ) = F a (x a + ∆x a , v 0a + ∆v a , µ 0a + ∆µ 0a ). Here, ∆v (1) a , and ∆µ (1) 0a are the O(δ) parts of ∆v a and ∆µ 0a given in Eq. (A3). In this Appendix, we do not derive expressions for the O(δ 2 ) parts ∆x (2) a , ∆v (2) a , and ∆µ (2) 0a by the Lie perturbation expansion method which is used to define the gyrocenter coordinates with the well-conserved magnetic moment because it would unnecessarily give higher-order accuracy to the coordinate transformation than the accuracy of the gyrocenter motion equations themselves shown in Eqs. (28)- (31). Instead, we determine these O(δ 2 ) terms from the conditions that the collisional change rates of energy and canonical toroidal angular momentum per unit volume in the gyrocenter space can be given in the conservative forms as shown below. Thus, the O(δ 2 ) terms are introduced not for accuracy of higher order in δ but for satisfying the conservation property of the collision operator.
In Eq. (C1), the expansions in (∆x (2) a , ∆v a , ∆µ 0a ) are truncated to the first order while the infinite series expansion in ∆x (1) a ≡ −ρ a is retained because fluctuations' wavelengths in the directions perpendicular to the equilibrium magnetic field can be of order of the gyroradius ρ a . In the WKB (or ballooning) representation, the above-mentioned infinite series expansion can be treated using the Bessel functions of the gyroradius normalized by the perpendicular wavelength. 6,8,30 We should also note that the gyrophase average · · · ξa is taken so that the gyrokinetic equation with the collision term is solved only for the gyrophase-averaged part of the gyrocenter distribution function.
For an arbitrary function A g a (Z a ) which is independent of the gyrophase ξ a , we obtain the following formula, where A p a (z a ) and A p0 a (z a ) are defined by A p a (z a ) = A p0 a (z a ) + ∆x (2) a · ∂ ∂x a + ∆v a ∂ ∂v a + ∆µ 0a ∂ ∂µ 0a A g a (z a ), A p0 a (z a ) = A g a (x a − ρ a , v a , µ 0a ).
We should note that the function A p (z a ) defined from A g a (Z a ) in Eq. (C4) does not exactly coincide with that given in Eq. (A6) in Appendix A by the second-and higher-order terms in the series expansion with respect to ∆z a . Integrating Eq. (C3) over the gyrocenter velocity space, we immediately obtain dU dµ dξ D g a (Z)C g a (Z)A g a (Z) where the transport flux J C Aa due to collisions and finite gyroradii is defined by (2) a ]C p a (z)A p a (z) x=X + ∞ n=1 −1 (n + 1)! i1,··· ,in ∂ n ∂X i1 · · · ∂X in × d 3 v ρ a ρ i1 a · · · ρ in a C p a (z)A p a (z) x=X . (C6) To the lowest order in δ, Eqs. (C3), (C5), and (C6) derived from the approximate collision operator in Eq. (C1) agree with Eqs. (B2), (B4), and (B5) given in Appendix B, respectively. The particle flux Γ C a due to collisions and finite gyroradii is given from Eq. (C6) with putting A p a (z) = 1 in the same way as in Eqs. (B8). Now, let us take A g a (Z) = H a (Z) in Eq. (C5). Here, H a (Z) denotes the gyrocenter Hamiltonian defined by Eq. (26). It is desirable that the gyrocenter velocityspace integral a dU dµ dξ D g a (Z)C g a (Z)H a (Z) takes the conservative form, which implies that the integral is expressed by the divergence term only and where ∆H a (z a ) ≡ ∆x (2) a · ∂ ∂x a + ∆v a ∂ ∂v a + ∆µ 0a ∂ ∂µ 0a H a (z a ) It is easily seen that Eq. (C7) is satisfied if ∆H a (z a ) = 0. Then, substituting Eq. (C2) into Eq. (C9) and using ∆H a (z a ) = 0, we have ∆x (2) a · ∂ ∂x a + ∆v (2) a ∂ ∂v a + ∆µ ξa + e a v a b · ∇b · ρ a + 1 We find that the right-hand side of Eq. (C10) is of O(δ 2 ). Then, as remarked after Eq. (C16), we can choose ∆x (2) a , ∆v (2) a , and ∆µ (2) 0a which satisfy Eq. (C10) and are of O(δ 2 ) so as to be consistent with Eq. (A3).
When we use ∆x (2) = 0, ∆v a = ∆v (1) a , and ∆µ 0a = ∆µ (1) 0a for Eq. (C1) by putting ∆v (2) a = ∆µ (2) 0a = 0, we have ∆H a (z a ) = O(δ 2 ) and a d 3 v C p a (z)H p a (z) = O(δ 3 ) because C p a (z) = O(δ) holds for the distribution function, the zeroth order of which is given by the local Maxwellian. Therefore, even for this case where a dU dµ dξ D g a (Z)C g a (Z)H a (Z) is not completely given in the conservative form, the residual term a d 3 v C p a (z)H p a (z) = O(δ 3 ) is smaller by a factor of δ than other transport terms of O(δ 2 ) in the lowest-order energy balance equation given by Eq. (100) in Sec. VI.B.
We next put A g a (Z) = (p c ζ ) g a (Z) in Eq. (C5). Here, (p c ζ ) g a (Z) denotes the canonical toroidal angular momentum defined by where A 0ζ = −χ and b ζ = I/B 0 . We now see that a dU dµ dξ D g a (Z)C g a (Z)(p c ζ ) g a (Z) takes the conservative form if a d 3 v C p a (z)(p c ζ ) p a (z) = 0.
Thus, the collision operator, which has the desired conservation properties as well as the accuracy required for correct description of collisional transport of the energy and the toroidal angular momentum, is given by Eq. (C1), in which ∆x