Eulerian variational formulations and momentum conservation laws for kinetic plasma systems

The Eulerian variational principle for the Vlasov-Poisson-Amp\`{e}re system of equations in a general coordinate system is presented. The invariance of the action integral under an arbitrary spatial coordinate transformation is used to obtain the momentum conservation law and the symmetric pressure in a more direct way than using the translational and rotational symmetries of the system. Next, the Eulerian variational principle is given for the collisionless drift kinetic equation, where particles' phase-space trajectories in given electromagnetic fields are described by Littlejohn's guiding center equations~[R. G. Littlejohn, J. Plasma Phys.\ {\bf 29}, 111 (1983)]. Then, it is shown that, in comparison with the conventional moment method, the invariance under a general spatial coordinate transformation yields a more convenient way to obtain the momentum balance as a three-dimensional vector equation in which the symmetric pressure tensor, the Lorentz force, and the magnetization current are properly expressed. Furthermore, the Eulerian formulation is presented for the extended drift kinetic system, for which, in addition to the drift kinetic equations for the distribution functions of all particle species, the quasineutrality condition and Amp\`{e}re's law to determine the self-consistent electromagnetic fields are given. Again, the momentum conservation law for the extended system is derived from the invariance under the general spatial coordinate transformation. Besides, the momentum balances are investigated for the cases where the collision and/or external source terms are added into the Vlasov and drift kinetic equations.


I. INTRODUCTION
So far, a large number of numerical simulations have been performed to investigate neoclassical and turbulent transport in toroidal plasmas. [1][2][3] As a modern theoretical technique for deriving basic kinetic model equations of such simulations, the variational principle 4-7 is used because the derived equations possess favorable conservation properties for long-time simulations to pursue evolutions of plasma profiles resulting from transport processes. Also, useful numerical schemes for plasma simulation satisfying the conservation properties have been developed by directly utilizing the variational formulation rather than numerically approximating the basic equations derived from the variational principle. [8][9][10][11] In recent years, background flow profiles are regarded as one of key factors which influence magnetic plasma confinement and large-scale gyrokinetic simulations are actively done to investigate momentum transport processes which determine the flow profiles. [12][13][14][15] Thus, pressure tensors or momentum transport fluxes need to be accurately evaluated because they play a critical role for the momentum balance in both neoclassical and turbulent transport theories. [16][17][18][19][20][21][22][23][24] In Ref. 6 , the Lagrangian variational formulation for the electromagnetic gyrokinetic system is presented from an approximate reduction of the Vlasov-Poisson-Ampère system which is equivalent to the Vlasov-Darwin system 25 in which such rapid phenomena as the electromagnetic waves with the speed of light c can be removed from the system (the terminology 'Vlasov-Poisson-Ampère system' has been customarily used instead of 'Vlasov-Darwin system' in the literature on the gyrokinetic theories 6,26 ). It is shown for the Vlasov-Poisson-Ampère system that, in the presence of the magnetic field, the canonical momentum conservation law derived from the space translational symmetry contains the asymmetric pressure tensor. In Ref. 27 , the angular momentum conservation law derived from the rotational symmetry and additional complicated procedures of the Belinfante-Rosenfeld type 28 were used to obtain the symmetric pressure tensor from the asymmetric canonical pressure tensor and to derive the same momentum conservation law as given in Ref. 25 .
In this work, the variational formulations for the Vlasov-Poisson-Ampère system and the drift kinetic system are presented in the invariant forms under general spatial coordinate transformations in analogy with the theory of general relativity. 29 For this purpose, the variational formulations here are completely based on the Eulerian picture [30][31][32][33][34] in which the spatial-coordinate dependence of the particle and field parts of the Lagrangian density can be more equally treated than in another type of formulation using the Lagrangian picture partially for the particle part. 6,27,35 Detailed descriptions about Lagrangian and Eulerian variational formulations are found in a recent paper by Brizard and Tronci 36 . The Eulerian method, which was pioneered by Newcomb 30 to formulate the magnetohydrodynamics equations and is used in the present paper, is also called the Euler-Poincaré reduc-tion procedure recently 31,32,34,36 . Here, in our Eulerian formulation, all the governing equations for these systems also take the invariant forms and the invariance of the action integrals can be utilized to derive the momentum conservation laws and/or the momentum balances as three-dimensional vector equations. The resultant momentum balance equations contain the symmetric pressure tensors which have 3 × 3 symmetric matrix components. These symmetric pressure tensor components are derived from taking the variation of the Lagrangian density with respect to the metric tensor components which appear due to the use of the general spatial coordinate system. The symmetry of the resultant pressure tensor is a natural result because the metric tensor is symmetric. Thus, the derivation of the symmetric pressure tensors shown in the present paper is more direct than the Belinfante-Rosenfeld-type technique and other previous methods. Furthermore, for all systems considered here, not only the momentum conservation laws but also the Belinfante-Rosenfeld type formulas 27,28 relating the symmetric pressure tensors to the asymmetric canonical pressure tensors are simultaneously derived from the invariance of the action integrals under general spatial coordinate transformations.
It is also found that the formulation presented here for deriving the momentum conservation law is more convenient than the conventional method based on taking moments of the basic kinetic equation especially for the drift kinetic system. 37 Normally, only the component of the momentum balance equation in the direction parallel to the magnetic field is derived from the parallel moment of the drift kinetic equation although it is not trivial what moment should be taken for the gyrophase-averaged distribution function to obtain the perpendicular momentum balance. On the other hand, the method based on the invariance with respect to the general spatial coordinate transformation can be applied to derive the momentum balance equations in both parallel and perpendicular directions simultaneously even for the drift kinetic system.
Normally, based on Noether's theorem, 5 the momentum conservation law in a certain direction is derived when a given system has a translational symmetry in that direction. Here, it should be noted that the invariance under the general spatial coordinate transformation holds more generally than the translational symmetry. Even in the case where the latter property is not satisfied, the former property can be valid and used to derive the momentum balance equation which does not take a conservative form. As shown in Sec. III, the drift kinetic system in given electromagnetic fields corresponds to the above-mentioned case. Thus, the momentum balance equation can be obtained for the drift kinetic system with general magnetic geometry. When self-consistent electromagnetic fields are treated as the solutions of the equations given simultaneously with the drift kinetic equations from the variational principle, the explicit dependence on the spatial coordinates is removed from the action integral, and accordingly the momentum conservation law is derived for the total system consisting of the charged particles and fields [see Sec. IV].
The rest of this paper is organized as follows. In Sec. II, the Eulerian formulation of the variational principle for the Vlasov-Poisson-Ampère system is presented. There, the same results as in Ref. 27 are reproduced although the general coordinates are used to write the equations in the invariant form and derive the momentum conservation law in a more direct way than in Ref. 27 . In Sec. III, the Eulerian variational principle is applied to the drift kinetic system, for which the collisionless drift kinetic equation and the momentum balance equation are obtained. In this system, which is immersed in the strong magnetic field, trajectories of charged particles are described by Littlejohn's guiding center equations. 4 In Sec. IV, the variational principle for the drift kinetic system is extended so that the quasineutrality condition and Ampère's law can be derived simultaneously with the drift kinetic equations for all particle species to determine the electromagnetic fields self-consistently with the distribution functions. The momentum conservation law for this extended drift kinetic system is derived as well. In Sec. V, it is shown how the momentum conservation and balance derived in Sec. II-IV are modified when the collision terms are added into the basic kinetic equations there. Finally, conclusions are given in Sec. VI. In Appendix A, the Eulerian variational principle is presented for the Vlasov-Poisson system and its momentum balance is derived. The energy conservation law in the Vlasov-Poisson system is also obtained in Appendix B. In Appendix C, the energy balance equation and the energy conservation law are shown for the drift kinetic systems described in Secs. III and IV.

II. VLASOV-POISSON-AMPÈRE SYSTEM
Here, the Vlasov-Poisson-Ampère system 27 is considered as an example of kinetic systems, for which the Eulerian variational principle is presented. Also, it is shown for this system how to obtain the momentum conservation law from the invariance of the action integral under general coordinate transformations.

A. Eulerian formulation of the variational principle in general coordinates
The distribution function on the phase space for particle species a is denoted by 2,3 are the position and velocity coordinates of the particle, respectively, 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 1 dv 2 dv 3 is given by Here, (x i ) i=1,2,3 represent a general spatial coordinate system which can be either a Cartesian or any other curved coordinate system. However, in the present paper, we assume that the position vector r is a function of only the spatial coordinates (x i ) i=1,2,3 and it is independent of time t. In the given spatial coordinate system, (v i ) i=1,2,3 are defined as contravariant components of the velocity vector by using (∂r/∂x i ) i=1,2,3 as the basis vectors.
In the Lagrangian picture, the motion of a particle of species a in the phase space is described by representing the position and velocity of the particle at time t as the functions, which satisfy the initial conditions at time t 0 , (2) Using the Lagrangian representations of the particle's motion given in Eq. (1), the distribution function at time t is related to that at time t 0 by We next represent the particle's velocity and acceleration in the Eulerian picture by which are related to those in the Lagrangian picture by Here,ḟ = ∂f (x m 0 , v m 0 , t)/∂t stands for the time derivative of an arbitrary function f (x m 0 , v m 0 , t) with (x m 0 , v m 0 ) kept fixed. Using Eqs. (3) and (5), we can show that the distribution function F a (x i , v i , t) satisfies the continuity equation in the six-dimensional phase space, In the present paper, we use the summation convention that an index repeated in a term [such as seen in Eq. (6)] represents summation over the range {1, 2, 3}. The action integral I to describe the Vlasov-Poisson-Ampère system is written as where the Lagrangian L is defined by the spatial integral of the Lagrangian density L over the volume V and L is given by Here, the single-particle Lagrangian L a for species a is written in the Eulerian picture as Note that the single-particle Lagrangian given by Eq. (2) in Ref. 27 is reproduced from Eq. (9) when replacing x i , v i , and u i ax in Eq. (9) with the corresponding Lagrangian representations x i aL , v i aL , andẋ aL . The field Lagrangian density L f on the right-hand side of Eq. (8) is given by Equation (10) is obtained by using the general phasespace coordinates (x i , v i ) to express the field Lagrangian density given by Eq. (3) in Ref. 27 . 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 and the components ∇ i A j (i, j = 1, 2, 3) of the covariant derivative of the covariant vector A j are defined by where the Levi-Civita symbol is denoted by the determinant of the metric tensor matrix is given by and the Christoffel symbols Γ k ij (i, j, k = 1, 2, 3) are defined by 38 The covariant and contravariant components of the metric tensor components are denoted by g ij and g ij , respectively, and they satisfy where δ i j represents the Kronecker delta defined by We now consider the virtual displacement of the particle's trajectory in the phase space, which is represented by the variations of the Lagrangian representations of the particle's position and velocity in Eq. (5) as The variations in the position and velocity are represented in the Eulerian picture by which are related to those in the Lagrangian picture by Making use of Eq. (5) to consider the variations in the particle's velocity and acceleration which result from the virtual displacement denoted by Eq. (18), we obtain where Eq.
We further consider that the spatial functional forms of the electrostatic potential φ, the covariant components A i of the vector potential, and the field λ associated with the Coulomb gauge condition [see Eq. (36)] are virtually varied by δφ, δA i , and δλ in addition to the virtual displacement of the particle's phase-space trajectory. Consequently, the action integral defined by Eq. (7) with where (∂L a /∂x i ) uax and (∂L a /∂v i ) uax represent the derivatives of L a in x i and v i , respectively, with u i ax kept fixed in L a . The definitions of the operators (d/dt) a and ∆ are shown later in Eqs. (30) and (37), respectively, and is the part which can be determined from the values of the variations δx i aE , δφ, and δA i on the boundaries of the integral region because of the divergence theorem.
We now show that the Vlasov-Poisson-Ampère system obeys the Eulerian variation principle. Namely, F a , φ, and A i are determined from the condition that δI = 0 for arbitrary variations δx i aE , δv i aE , δφ, δA i , and δλ which vanish on the boundaries of the integral region. First, it is found from Eq. (23) that δI/δv i aE = 0 gives which is rewritten as Here, we should note that u i ax = v i is derived from Eq. (26) under the condition that F a = 0. However, since u i ax enters Eq. (6) in the form of the product F a u i ax , it doesn't cause any trouble to simply write from now on instead of Eq. (26) even without assuming F a = 0. This simplification of omitting F a will also be done below in the processes where the equation for u i av [see Eq. (34)] is derived.
We next use δI/δx i aE = 0 to obtain where the covariant vector component p ai of the canonical momentum and the time derivative (d/dt) a along the motion of the particle of species a in the phase space are defined by and respectively. From Eq. (30), we also have Equation (28) can be rewritten as the covariant form of Newton's motion equation in the general coordinate system, The contravariant form of Newton's motion equation is obtained from Eq. (32) as It should be noted that the Levi-Civita symbol ǫ ijk ≡ ǫ ijk [see Eq. (13)] can be regarded as either a contravariant tensor density of weight 1 or a covariant tensor density of weight −1. Then, √ gǫ ijk and ǫ ijk / √ g represent covariant and contravariant tensors, respectively, which are used in the Lorentz force terms on the right-hand side of Eqs. (32) and (34).
Substituting Eqs. (27) and (34) into Eq. (6) yields the Vlasov kinetic equation, As noted after Eq. (27), F a appears as a factor in the equations δI/δx i aE = δI/δv i aE = 0 although it is omitted in writing Eqs. (27), (28), (32) and (34). This omission of F a does not make a difference in deriving the Vlasov equation in Eq. (35) by substituting the motion equations, Eqs. (27) and (34), into Eq. (6) because u i ax and u i av enter Eq. (6) in the forms of the products F a u i ax and F a u i av . We use δI/δλ = 0 to obtain the Coulomb (or transverse) gauge condition, Poisson's equation is derived from δI/δφ = 0 as and δI/δA i = 0 gives where j i represents the ith contravariant component of the current density vector defined by The transverse (or solenoidal) part of Eq. (38) is written as Ampère's law, where j i T represents the ith contravariant component of the transverse part of the the current density vector. Note that an arbitrary vector field a can be written as a = a L + a T , where the longitudinal (or irrotational) part a L and the transverse (or solenoidal) part a T satisfy ∇ × a L = 0 and ∇ · a T = 0, respectively. 42 Equations (35), (37), and (40) are the governing equations for the Vlasov-Poisson-Ampère system. Thus, the same system of equations as shown in Ref. 27 are reproduced in the present work although the equations here are represented using the general spatial coordinates (x i ) i=1,2,3 and the contravariant velocity vector components (v i ) i=1,2,3 . Using Eq. (37), the longitudinal part of Eq. (38), and the charge conservation law obtained from Eq. (35), we obtain where E Li = −∂φ/∂x i represents the longitudinal electric field given by the electrostatic potential. Then, we can put 27 which is used hereafter. Then, we find that Eqs. (37), (38), (40), and (41) give the Darwin model 25 as noted in Refs. 6,27 .

B. Transformation of spatial coordinates
We now consider the transformation of the spatial coordinates written as where the infinitesimal variation in the spatial coordinate x i is denoted by ξ i (x n ) which is regarded as an arbitrary function of only the spatial coordinates. Under the transformation of the spatial coordinates, the velocity components (v i ) i=1,2,3 are transformed as the contravariant vector components. Thus, the velocity component v ′i in the transformed coordinate system is written as where the infinitesimal variation δv i in the velocity component is given by 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. II.A. The electrostatic potential is a scalar which is invariant under the transformation of the spatial coordinates, Here, we define the variation δφ in the functional form of φ due to the spatial coordinate transformation by Note that the spatial arguments of φ ′ and φ are the same as each other on the right-hand side of Eq. (47) while they are different in Eq. (46). Then, substituting (46) and using Eq. (47), we obtain where L ξ denotes the Lie derivative 43 associated with the vector field (ξ i ). In the same way as in Eq. (48), the variation δλ in another scalar variable λ is written as In the transformed spatial coordinates, the covariant vector components of the vector potential are written as In the same way as in Eq. (47), we define the variation δA i in the functional form of A i due to the spatial coordinate transformation by (50) and using Eq. (51), we obtain where we see that the Lie derivative L ξ can be used again to represent δA i . The contravariant vector components E i , the covariant metric tensor components g ij , and the contravariant tensor components g ij are transformed as Then, following the procedures similar to those used in deriving Eqs. (48) and (52), the variations in the functional forms of E i , g ij , and g ij due to the spatial coordinate transformation are derived as The transformation of the spatial coordinates given by Eq. (43) changes the Lagrangian representations of the trajectory of the particle's motion in the phase space as , (55) where the particle's position and velocity at time t 0 are written in the transformed coordinates as Since v i is the contravariant vector component, its variation η i caused by the change ξ i in the spatial coordinate x i can be written as In the transformed coordinate system, the distribution function is given by in the transformed and original coordinate systems are related to each other by (59) The variation δF a in the functional form of the distribution function due to the spatial coordinate transformation is defined by Then, using Eqs. (55), (58), (59), and (60), we obtain The relations between the Eulerian and Lagrangian representations of the particle's velocity and acceleration shown in Eq. (5) are rewritten in the transformed coordinate system as We also write to define δu i ax and δu i av as the variations in the Eulerian functional forms of the particle's velocity and acceleration, respectively. Using Eqs. (55), (62), and (63), we find that δu i ax and δu i av are written as (64) C. Derivation of the momentum conservation law i (x ′n , t), and g ′ ij (x ′n , t) in Eqs. (7)- (10) to define the action integral I ′ in the transformed coordinates (x ′n , v ′n ).
Then, using Eqs. (48), (49), (52), (54), (61) and (64), we find that the variation δI ≡ I ′ − I in the action integral is written as where the canonical momentum vector density P j c and the canonical pressure tensor density Π ij c are defined by and respectively. The symmetric tensor density Θ ij and the third-rank tensor density F ijk are defined by and respectively. In deriving Eq. (65), Eqs. (27), (28), (36), (37), and (38), which are derived from the variational principle in Sec. II.A, are also used. The symmetry condition, is naturally confirmed in Eq. (68) because the symmetric metric tensor density components g ij are used for differentiating the Lagrangian density L in the definition of Θ ij . It should be noted that, in Eq. (68), partial derivatives with respect to g ij need to be carefully done 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, 29 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 shown in Eqs. (68) and (69) are defined taking into account the symmetry under exchange of the indices i and j.
We see from Eq. (69) that the third-rank tensor density components F ijk is anti-symmetric with respect to exchanging the superscripts i and k, Using Eq. (71) and the commutation condition, we obtain Note that the commutation condition in Eq. (72) is valid because the three-dimensional real space considered here is a flat one with no curvature. For a general curved space, the Riemann curvature tensor R l mik is used to write 38 where V m is the mth contravariant component of an arbitrary vector field. Then, we find where Eq. (71) and the properties of the Riemann tensor (R a bcd = −R a bdc , R bc ≡ R a bca = R cb ) are used. Thus, we find the interesting fact that Eq. (73) is valid even in the curved space when Eq. (71) is satisfied.
Since the action integral is invariant under an arbitrary transformation of the spatial coordinates, δI shown in Eq. (65) vanishes for any ξ j so that we obtain the momentum conservation law, and the relation of the symmetric pressure tensor density Θ ij to the asymmetric canonical tensor density Π ij c , Equations (76) and (77) which can be used to rewrite the momentum conservation law in Eq. (76) as In Ref. 27 , the momentum conservation law, Eq. (79), including the asymmetric canonical momentum tensor density, Π ij c , is derived from the space translational symmetry of the action integral I before the relation of the Belinfante-Rosenfeld type symmetric pressure tensor Θ ij to Π ij c in Eq. (77) is obtained from the rotational symmetry of I to derive the other momentum conservation law, Eq. (76). On the other hand, in the present work, both the momentum conservation law, Eq. (76), and the relation of Θ ij to Π ij c , Eq. (77), are derived at once from the invariance of I under general spatial coordinate transformations including the space translation and rotation. We should also note that Eq. (76) can be further modified into a more physically familiar form of the momentum conservation law as shown in Eq. (33) of Ref. 27 .
It is shown in Appendix A that, reducing the field Lagrangian density given by Eq. (10) to the more simplified one defined in Eq. (A3) and regarding the vector potential in Eq. (9) as a fixed time-independent field, the governing equations for the Vlasov-Poisson system can be obtained from the Eulerian variational principle in the same manner as shown for the Vlasov-Poisson-Ampère system. As pointed out by Qin et al., 39 when governing equations for a simplified system are obtained by applying a certain approximation to a Lagrangian for another system, the exact energy and momentum conservation laws in the simplified system should be derived from the symmetry properties of the approximate Lagrangian and they generally disagree with those obtained by just making a similar approximation to the conservation laws in the original system. The momentum balance and the energy conservation law in the Vlasov-Poisson system are derived in Appendices A and B, respectively, where they are found to agree with those given by Qin et al. 39

III. DRIFT KINETIC SYSTEM
In this section, the Eulerian variational principle is presented for the collisionless drift kinetic equation which governs the time evolution of the phase-space distribution function of guiding centers of charged particles in the strong magnetic field. The invariance of the drift kinetic system under an arbitrary spatial coordinate transformation is used to obtain the momentum balance as a three-dimensional vector equation in which the symmetric pressure tensor, the Lorentz force, and the magnetization current are properly included.

A. Eulerian variational principle for derivation of the collisionless drift kinetic equation
We here start with defining the action integral for the drift kinetic system by where the Lagrangian density is written as The guiding center position is represented in terms of the general spatial coordinates (x i ) i=1,2,3 , for which the metric tensor is given by g ij . The velocity component of the guiding center along the magnetic field line, the magnetic moment, and the gyrophase angle are denoted by v , µ, and ϑ, respectively. The integral with respect to the velocity space variables (v , µ, ϑ) is denoted by and the Lagrangian for the single guiding center is given by Here, the unit vector parallel to the magnetic field is written as where the field strength is given by and the ith contravariant component B i of the magnetic field is expressed in Eq. (11). The Lagrangian L GC shown in Eq. (83) represents Littlejohn's guiding-center Lagrangian 4 written using the general spatial coordinates and the Eulerian picture. We now describe the particle's motion in the Lagrangian picture by representing the guiding center position coordinates, parallel velocity, magnetic moment, and gyrophase at time t as the functions x i L (x n 0 , v 0 , µ 0 , ϑ 0 , t 0 ; t), v L (x n 0 , v 0 , µ 0 , ϑ 0 , t 0 ; t), µ L (x n 0 , v 0 , µ 0 , ϑ 0 , t 0 ; t), and ϑ L (x n 0 , v 0 , µ 0 , ϑ 0 , t 0 ; t), respectively, where x n 0 , v 0 , µ 0 , and ϑ 0 denote their values at the initial time t 0 . Then, the distribution function F (x i , v , µ, ϑ, t) at time t is related to that at time t 0 by where In the Eulerian picture, the temporal change rates of the guiding center position, parallel velocity, magnetic moment, and gyrophase are denoted by the functions u i x (x m , v , µ, ϑ, t), u v (x m , v , µ, ϑ, t), u µ (x m , v , µ, ϑ, t), and u ϑ (x m , v , µ, ϑ, t), respectively, and they are related to those in the Lagrangian picture by whereḟ = ∂f (x n 0 , v 0 , µ 0 , ϑ 0 , t)/∂t represents the time derivative of an arbitrary function f (x n 0 , v 0 , µ 0 , ϑ 0 , t) with (x n 0 , v 0 , µ 0 , ϑ 0 ) kept fixed. It can be shown from Eqs. (86) and (87) that F satisfies (88) The virtual displacement of the particle's trajectory in the (x i , v , µ, ϑ) space is represented by the variations of the Lagrangian representations of the particle's motion as δx i The variations in the guiding center position, parallel velocity, magnetic moment, and gyrophase are represented in the Eulerian picture by δx i E (x n , v , µ, ϑ, t), δv E (x n , v , µ, ϑ, t), δµ E (x n , v , µ, ϑ, t), and δϑ E (x n , v , µ, ϑ, t), which are related to those in the Lagrangian picture by Using Eqs. (87) and (89), the variations in the functional forms of u i x , u v , u µ , and u ϑ due to the virtual displace-ment of the particle's trajectory are given by The variation in the distribution function due to the virtual displacement of the particle's trajectory is written by using Eqs. (86) and (89) as (91) Using Eqs. (88), (90), and (91), we find that the variation in the action integral I DK due to the virtual displacement of particle's trajectory is written as where (∂L GC /∂x i ) u , (∂L GC /∂v ) u , (∂L GC /∂µ) u , and (∂L GC /∂ϑ) u denote the derivatives of L GC in x i , v , µ, and ϑ, respectively, with (u i x , u ϑ ) kept fixed in L GC , and the time derivative along the particle's trajectory is represented by We now use the Eulerian variation principle which implies that the collisionless drift kinetic equation for the distribution function F can be derived from the condition that δI DK = 0 for arbitrary variations δx i E , δv E , δµ E , and δϑ E which vanish on the boundaries of the integral region. We first use δI DK /δx i E = 0 to obtain where p i represents the covariant vector component of the canonical momentum defined by (95) We should note that the distribution function F is included as a factor in δI DK /δx i E = 0 although it is omitted from Eq. (94) for simplicity in the same way as done in Sec. II.A. This omission of F is also done in the other equations obtained below from δI DK = 0 although it does not make a difference in deriving the resultant collisionless drift kinetic equation in Eq. (109). We can rewrite Eq. (94) as where the modified electric and magnetic fields are defined by and respectively. Next, δI DK /δv E = 0 is used to obtain from which we have Furthermore, δI DK /δµ E = 0 and δI DK /δϑ E = 0 yield and d dt respectively. Equations (96), (100), (101), and (102) are rewritten as and where Equations (103) and (104) are obtained by taking the vector and scalar products between the magnetic field and Eq. (96), respectively. Also, using Eq. (93), we can write Then, with the help of Eq. (108), it is clearly confirmed that Eqs. (103)-(106) represent the same guiding center motion equations as derived by Littlejohn from the guiding center Lagrangian. We can verify that the right-hand sides of Eqs. (103)-(106) are all independent of ϑ and that the magnetic moment µ is an invariant of motion as seen from Eq. (105). Substituting Eqs. (103)-(106) into Eqs. (88) and taking its average with respect to the gyrophase ϑ, the collisionless drift kinetic equation is derived as where F denotes the gyrophase-averaged distribution function,

B. Transformation of spatial coordinates
Here, in the same way as in Sec. II.B, the infinitesimal transformation of the spatial coordinates is given by Eq. (43) and the infinitesimal variation ξ i in the spatial coordinate x i is again regarded as an arbitrary function of only the spatial coordinates. However, it should be noted that the other variables (v , µ, ϑ) are independent of the choice of the spatial coordinates because they are defined from the relation of the velocity vector to the direction of the local magnetic field. This is in contrast to the case of Sec. II.B where the velocity components (v i ) i=1,2,3 are transformed as the contravariant vector components under the spatial coordinate transformation.
The spatial coordinate transformation given by Eq. (43) changes the Lagrangian representation of the guiding center position as represents the guiding center position at time t 0 in the transformed spatial coordinates. The distribution function is written in the transformed coordinates as (112) The initial distribution functions F ′ (x ′n 0 , v 0 , µ 0 , ϑ 0 , t 0 ) and F (x n 0 , v 0 , µ 0 , ϑ 0 , t 0 ) in the transformed and original coordinate systems are related to each other by (113) The variation δF in the functional form of the distribution function F due to the spatial coordinate transformation is defined by (114) Then, it is shown by using Eqs. (111), (112), (113), and (114) that δF can be represented by In the same way as seen in Sec. II.B, the variations in the functional forms of u i x , u v , u µ , u ϑ , φ, A i , and g ij due to the spatial coordinate transformation are denoted by δu i x , δu v , δu µ , δu ϑ , δφ, δA i , and δg ij , respectively. They can be represented by using the Lie derivative L ξ as The expressions of δφ, δA i , and δg ij in terms of the Lie derivative are shown in Eqs. (48) (52), and (54), respectively.

C. Derivation of the momentum balance equation
We can use Eqs. (80)-(83), (115), and (116) to derive the action integral I ′ DK = I DK +δI DK in the transformed spatial coordinates. Here, the variation δI DK in the action integral due to the spatial coordinate transformation is written as where the vector density J j DK and the tensor density T ij DK are defined by and respectively. In deriving Eq. (117), we also need to use Eqs. and are used. We see that the inertia term in the momentum balance equation, Eq. (120), contains only the parallel momentum component while the electric current eN V k in the Lorentz force term consists of the guiding-center current and the magnetization current 40 as shown in Eq. (122). The symmetric pressure tensor density P ij on the right-hand side of Eq. (120) is defined by where P ij CGL is given in the Chew-Goldberger-Low (CGL) form, 41 and π ij ∧ is the non-CGL part written as Here, the perpendicular component of the guiding center velocity is represented by The symmetric pressure tensor given by Eq. (123) with Eqs. (124) and (125) agrees with that given by Eq. (19) in Ref. 37 . The CGL pressure tensor shown in Eq. (124) contains the scalar (or isotropic) part, which represents background pressure, and the anisotropic part, the magnitude of which is considered to be smaller than the background pressure by the factor ∼ ρ/L in the neoclassical transport theory. Here, ρ and L represent the gyroradius and the equilibrium gradient scale length, respectively. On the other hand, the magnitude of the non-CGL pressure tensor defined in Eq. (125) is regarded as ∼ (ρ/L) 2 .
We next substitute Eq. (83) into Eq. (119). Then, the other condition, T ij DK = 0, derived from putting δI DK = 0 in Eq. (117) can be written as where P ij is the symmetric pressure tensor given by Eq. (123) and P ij c is define by Here, P ij c is an asymmetric tensor density representing the transport of the canonical momentum. The difference D ij between P ij and P ij c is written as

IV. DRIFT KINETIC SYSTEM WITH SELF-CONSISTENT FIELDS
In this section, not only the drift kinetic equations but also the equations for self-consistently generated electromagnetic fields are treated as constituents of the governing equations of the extended drift kinetic system. The Eulerian variational principle is used to present all the governing equations and to derive the momentum conservation law satisfied by them. The energy conservation law in the extended drift kinetic system is derived in Appendix C where the energy balance in the drift kinetic system considered in Sec. III is also obtained.

A. Quasineutrality and Ampère's law combined with drift kinetic equations
We here combine the quasineutrality condition and Ampère's law with the drift kinetic equations in order to simultaneously determine the electromagnetic fields and the distribution functions for all particle species. The action integral I DKF for deriving all the governing equations is written as where the Lagrangian density L DKF is given by Here, the subscript a represents the particle species. It is seen from Eq. (130) that L DKF contains the summation of the drift kinetic Lagrangian densities [see Eq. (81)] over all species and the magnetic energy density with the minus sign.
We now virtually let the trajectories of particles for all species, the electrostatic potential, and the vector potential vary infinitesimally. Then, the resulting variation δI DKF in the action integral I DKF is expressed as where (d/dt) a denotes the time derivative along the trajectory of the particle of species a [see Eq. (93)], δI DKF b represents the part which is written as the boundary integrals, and M k = g kl M l is the kth covariant component of the magnetization vector. The kth contravariant component of the magnetization vector is defined by The magnitude of the second term in the integrand on the right-hand side of Eq. (132) is smaller than that of the first term by the factor ∼ ρ/L. Except for this small correction, Eq. (132) agrees with the well-known expression of the magnetization vector. 40 For each particle species a, the same motion equations as shown in Eqs. (103)-(106) are derived from δI DKF /δx i aE = δI DKF /δv aE = δI DKF /δµ aE = δI DKF /δϑ aE = 0 and accordingly the same collisionless drift kinetic equation as Eq. (109) is obtained for the gyrophase-averaged distribution function F a ≡ F a dϑ/(2π).
The remaining governing equations of the system, namely, the quasineutrality condition and Ampère's law are derived from δI DKF /δφ = 0 and δI DKF /δA i = 0, respectively, as and where the ith contravariant component of the electric current vector density is defined by

B. The momentum conservation law
We now consider the transformation of the spatial coordinates given by Eq. (43) again. Under the spatial coordinate transformation, the variables (v , µ, ϑ) are kept fixed as noted in Sec. III.B. In the same way as in Eqs. (115) and (116), the spatial coordinate transformation causes the variations in the distribution function F a and the functional forms of (u i ax , u v a , u aµ , u aϑ ) which are written as and The variations in φ, A i , and g ij due to the spatial coordinate transformation are shown in Eqs. (48) (52), and (54), respectively. Using the expressions of these variations described above, we find that the variation δI DKF in the action integral I DKF caused by the spatial coordinate transformation is written in the form, (138) Here, J j DKF is given by where P j tot and Θ ij tot represent the total momentum vector density and the total symmetric pressure tensor density defined by and respectively. It should be noted that, in Eq. (140), the vector potential part of the canonical momentum does not contribute to the total momentum because of the quasineutrality condition, Eq. (133). The first term on the right-hand side of Eq. (141) is the particle part of the pressure tensor density defined by which consists of the CGL part, and the non-CGL part, The second term on the right-hand side of Eq. (141) is given by which represents the Maxwell stress tensor due to the magnetic field with the opposite sign. It is clear that Θ ij tot , Θ ij p , Θ ij f , P ij CGL , and π ij ∧ are all symmetric with respect to the interchange of the superscripts i and j.
The contravariant (i, j)-component T ij DKF of the tensor density appearing on the left-hand side of Eq. (138) is written as where the total asymmetric canonical pressure tensor density Π ij tot and the third-rank tensor density F ijk tot are defined by and respectively. We can immediately see that F ijk tot satisfies from which we have in the same way as in Eq. (73) Since the action integral I DKF is invariant under the spatial coordinate transformation, δI DKF written in Eq. (138) vanishes for any ξ j . Thus, the integrands at the interior and boundary points shown on the right-hand side of Eq. (138) should vanish separately so we obtain J j DKF = 0 and T ij DKF = 0. We find from Eqs. (139) and (146) that J j DKF = 0 represents the total momentum conservation law, and T ij DKF = 0 gives the relation of the total symmetric pressure tensor density Θ ij tot to the total asymmetric canonical pressure tensor density Π ij tot , Combining Eqs. (150) and (152) shows We clearly see that the relations between the two types of the pressure tensors shown in Eqs. (152) and (153) take the same forms as those given by Eqs. (77) and (78) in Sec. II.C, respectively. It is noted that the momentum conservation law similar to Eq. (151) was derived by Brizard and Tronci 36 for the guiding-center Vlasov-Maxwell system. In their model, the electromagnetic fields are determined by the full Maxwell equations including the Maxwell displacement current so that their system contains such rapid phenomena as the electromagnetic waves with the speed of light and the Maxwell stress due to the electric field. They also derived the symmetric pressure tensor including the same particle part as given by Eqs. (142)-(144) although they modified the magnetization term in the canonical momentum conservation law to transform the asymmetric pressure tensor to the symmetric one. Thus, their derivation is different from our direct derivation of the symmetric pressure tensor by taking the variation with respect to the metric tensor.
It is instructive here to consider momentum conservation law, Eq. (151), in the the equilibrium limit where the distribution functions are assumed to take the local Maxwellian form, F a = N a (m a /2πT a ) 3/2 exp[−( 1 2 m a v 2 + µB)/T a ] (note that, precisely speaking, this local Maxwellian distribution function is not the exact stationary solution but the zeroth-order one of the drift kinetic equation in the gyroradius ordering and that the deviation from the local Maxwellian appears in the firstorder solution). Then, it is found from Eqs. (142)-(144) that Θ ij p = P g ij where P ≡ a N a T a . We now use the conventional vector notation to rewrite Eq. (151) in the equilibrium state (∂/∂t = 0) as where Eqs.
where the current density is given by Eq. (135) as This formula is called the magnetization law 40 ; the current is represented by the sum of the flow of guiding centers and the curl of the magnetization M [Eq. (132)], which are given by the first and second terms on the right-hand side, respectively. As shown in Ref. 40 , it is found from using the local Maxwellian distribution functions that the sum of the perpendicular components of the first and second terms on the right-hand side of the above magnetization law gives the diamagnetic current, (c/B 2 )(B × ∇P ). Recall that the perpendicular component (u ax ) ⊥ of the guiding center velocity and the magnetization M are both produced from gyrations of particles around magnetic field lines. Even though finite gyroradius effects are not described by the guiding center distribution functions F a alone, such effects are partly included in (u ax ) ⊥ and M which help express the current properly and recover the familiar force balance equation in the MHD equilibrium as shown above.

V. EFFECTS OF COLLISIONS
We here investigate how collisions influence the momentum conservation laws and the momentum balance equation shown in Secs. II.C, III.C, and IV.B when the collision term is added into the right-hand sides of the Vlasov and drift kinetic equations. Effects of the collision term added into the right-hand side of Eq. (35) were already studied in Ref. 44 where it was shown how to evaluate the correction of the energy and momentum conservation laws due to the collision and other source terms. According to the prescription given in Ref. 44 , the modified conservation laws are obtained from the original ones with the time derivative of the distribution function being replaced as where K a is the term added into the the right-hand side of the kinetic equation to represent 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. When K a is added into the right-hand side of Eq. (35), Eq. (157) is applied to the momentum conservation law in Eq. (76), where the term ∂P j c /∂t contains a ∂F a /∂t as seen from Eq. (66). Then, we find that the resulting momentum balance equation is given by Eq. (76) with making the replacement, In the case where K a is given by the Coulomb collision operator (such as the Landau operator) which satisfies the conservation laws of the particles' number ( d 3 v K a = 0) and the momentum ( a d 3 v K a m a v j = 0), the velocity space integral vanishes in Eq. (158) and we have the momentum conservation law in the same form as that for the case of K a = 0. Also, it is noted in Ref. 44 that, even if K a contains some external source parts other than the collision term, the charge conservation law requires the condition a e a d 3 v K a = 0 which implies the correction term proportional to A j vanishes in Eq. (158). Next, let us consider the case where K is added into the right-hand side of the drift kinetic equation, Eq. (109), for a given particle species, in which the subscript representing the particle species is omitted. Here, K is regarded as gyrophase-averaged. Applying Eq. (157) to this case, we find the momentum balance equations is derived from Eq. (120) with the following replacement, (159) The parallel component of this derived momentum balance equation agrees with Eq. (18) in Ref. 37 where its perpendicular components are not derived. We see from Eq. (159) that the effect of K on the momentum balance equation for the single particle species is written as d 3 v Kmv b j . When K is given by the Coulomb collision operator, d 3 v Kmv b j represents the collisional transfer of the parallel momentum from the other particle species to the given species.
Since the momentum balance equation obtained by applying Eq. (159) to Eq. (120) is always valid for the distribution function which is the solution of the drift kinetic equation including K, it is also valid for each particle species even when the quasineutrality condition and Ampère's law are additionally imposed for the selfconsistent fields as in Sec. IV. Furthermore, we can use Eq. (157) in Eq. (151) to see how the total momentum conservation law for the drift kinetic system in the selfconsistent fields are modified by adding K a into the drift kinetic equation for the particle species a. The resultant momentum balance equation is given from Eq. (151) by putting This corresponds to the species summation of Eq.(159). When K a represents the Coulomb collision operator in the zero-gyroradius limit, it satisfies a d 3 v K a m a v = 0 and the momentum conservation law takes the same form as in Eq. (151). Note that the momentum conservation in Coulomb collisions is satisfied locally at the colliding particles' position which differs from the guidingcenter position. Therefore, if the finite gyroradius effect is taken into account, a d 3 v K a m a v does not generally vanish for the gyrophase-averaged collision operator K a at the fixed guiding-center position, which includes the classical transport processes. 44

VI. CONCLUSIONS
In this work, Eulerian variational formulations for kinetic plasma systems are presented. As examples, the Vlasov-Poisson-Ampère system and the drift kinetic systems are investigated. For the drift kinetic system, the additional case is also considered in which the quasineutrality condition and Ampère's law are included as supplementary governing equations to describe the selfconsistent fields.
For all cases treated here, general spatial coordinates are used to represent the action integrals and the governing equations which take the forms being invariant under an arbitrary (time-independent) transformation of spatial coordinates. Furthermore, the invariance of the action integral under the spatial coordinate transformation is made use of to derive the momentum conservation laws and/or the momentum balance in which the functional derivatives of the Lagrangians with respect to the metric tensor components yield the proper symmetric pressure tensors more directly than conventional techniques using translational and rotational symmetries or taking the moments of the kinetic equations.
It is also clarified how the momentum balances are influenced by adding the collision and/or external source terms into the kinetic equations. Since the invariance under the spatial coordinate transformations is valid independently whether the system has symmetric geometry or not, the present formulation can be applied to kinetic studies of plasmas confined in general magnetic configurations including nonaxisymmetric systems such as stellarators and heliotrons. 45 For example, the momentum balance equation derived here for the drift kinetic system are considered useful for verifications of accuracy of numerical simulations using Littlejohn's guiding center equations to study neoclassical transport processes in various magnetic geometries. The extension of the present study to the gyrokinetic system is now in progress and the results will be reported elsewhere.
is externally given as a time-independent one, B 0 (x) = ∇ × A 0 (x), and the electric field is written in terms of the electrostatic potential φ(x, t) as E(x, t) = −∇φ(x, t). The action integral I V P to describe the Vlasov-Poisson system is given by where the Lagrangian density L V P is written as Here, the single-particle Lagrangian L a for species a is defined by Eq. (9) where the covariant components A j (x n , t) of the vector potential are replaced with the time-independent ones A 0j (x n ). The field Lagrangian density L V P f is written in the general spatial coordinates In the same way as in Sec. II.A, we now consider the virtual displacement of the particle's trajectory, for which the variations in the particle's position and velocity are represented in the Eulerian picture as δx i aE and δv i aE , respectively [see Eq. (19)]. The electrostatic potential field φ is also virtually varied by δφ. However, since the vector potential A 0j is fixed, its virtual variation δA 0j does not appear. Then, the variation in the action integral I V P is given by where δI V P b represents the part which is determined from the values of δx i aE and δφ on the boundaries of the integral region. In deriving Eq. (A4), Eqs. (21) and (22) are used. Imposing the condition that δI V P = 0 for arbitrary variations δx i aE , δv i aE and δφ which vanish on the boundaries, the same equations as those in Eqs. (35) and (37) are obtained in the same manner as shown in Sec. II.A. Recalling that, in the present case, E i = −∂φ/∂x i because of ∂A 0i /∂t = 0, we confirm the fact that Eqs. (35) and (37) resulting from δI V P = 0 form the governing equations of the Vlasov-Poisson system.
To derive the momentum balance in the Vlasov-Poisson system, we next consider the infinitesimal spatial coordinate transformation as shown in Eq. (43) of Sec. II.B. In the same way as in Sec. II.B, the variations in v i , φ, F a , u i ax , and u i av due to the spatial coordinate transformation are denoted by δv i , δφ, δF a , δu i ax , and δu i av , respectively, which are defined by Eqs. (45), (48), (61), and (64). We should note that the spatial coordinate transformation also causes the variations in the metric tensor components [see Eq. (54)] as well as the variation δA 0i in the functional form of the contravariant component A 0i of the externally given vector potential where δA 0i is written in the same form as in Eq. (52), This is contrast to the case that δA 0i does not appear in considering the virtual variations to derive the governing equations of the Vlasov-Poisson system from δI V P = 0. Using Eqs. (48), (54), (61), (64), and (A5), it is found that the variation δI V P in the action integral I V P due to the spatial coordinate transformation is written as where J j V P and T ij V P are given by and respectively. The conditions, δI/δx i aE = δI/δv i aE = 0 and δI/δφ = 0, from which the Vlasov kinetic equation and Poisson's equation are derived, are also used in deriving Eq. (A6). In Eq. (A7), P j and P j c are the kinetic and canonical momentum densities which are defined by and (A10) respectively, and F j L represents the Lorentz force given by where B 0i ≡ g ij B j 0 , B i 0 ≡ (ǫ ijk / √ g)(∂A 0k /∂x j ), and ∂L a /∂A 0k = (e a /c)v k are used. Using the continuity equation derived from the Vlasov kinetic equation, we can confirm that the right-hand side of the first line in Eq. (A7) equals the last line. In addition, Eqs. (A7) and (A8) contain the symmetric pressure tensor Π ij and the canonical pressure tensor Π ij c which are defined by and respectively, where E i L ≡ g ij E Lj and E Li ≡ −∂φ/∂x i are used.
Because of the invariance of the action integral I V P under the general spatial coordinate transformation, δI V P vanishes for any ξ j , and accordingly, we have J j V P = 0 and T ij V P = 0 from Eq. (A6). The momentum balance in the Vlasov-Poisson system is obtained from J j V P = 0 as which agrees with that shown by Qin et al. 39 Another condition T ij V P = 0 gives the relation between the symmetric pressure tensor Π ij and the canonical pressure tensor Π ij c . The validity of T ij V P = 0 is also easily verified from Eqs. (A8), (A12), (A13), and ∂L a /∂A 0i = (e a /c)v i .

Appendix B: ENERGY CONSERVATION IN THE VLASOV-POISSON SYSTEM
In this Appendix, we consider the energy balance in the Vlasov-Poisson system. The energy conservation laws for the Vlasov-Poisson-Ampère system and the Boltzmann-Poisson-Ampère system are shown in Refs. 27 and 44 , respectively. In contrast to the case in Appendix A where the momentum balance in the Vlasov-Poisson system is derived, we do not need to use the general spatial coordinate system here. So we now use only the Cartesian coordinate system and represent three-dimensional vectors in terms of boldface letters. Either a Lagrangian or an Eulerian variational formulation can be used for the derivation of the energy balance although we here follow the Eulerian formulation to treat the variation of the action integral under translation in time. The infinitesimal time translation is represented by transforming the time coordinate as where ǫ is an infinitesimal constant. The time translation causes the variations δ t I V P in the action integral I V P , where I V P is defined in Eq. (A1) and δ t I V P is written as In this Appendix, we use δ t · · · to represent the variations associated with the time translation. The variations in u ax ≡ (u i ax ) i=1,2,3 , u av ≡ (u i av ) i=1,2,3 , φ, and F a due to the time translation are written as δ t u ax = −ǫ ∂u ax ∂t , δ t u av = −ǫ ∂u av ∂t , δ t φ = −ǫ ∂φ ∂t (B3) and respectively, where Eq. (6) is used. Then, substituting Eqs. (B3) and (B4) into Eq. (B2) and using δI V P /δx E = δI V P /δv E = 0 and δI V P /δφ = 0, we obtain where the canonical energy density E V P c and the canonical energy flux Q V P c are defined by and respectively. Here, the electrostatic electric field is represented by E L ≡ −∇φ.
Since the Lagrangian density L V P defined in Eq. (A2) for the Vlasov-Poisson system depends on time t only through the functions u ax , F a , and φ which are all determined by the variational principle (see Appendix A), the action integral I V P given in Eq. (A1) is invariant under the time translation. Therefore, δ t I V P vanishes for an arbitrarily chosen integral region [t 1 , t 2 ] × V and accordingly, the integrand in Eq. (B5) also vanishes. Thus, we obtain the local energy conservation law written as where the energy density E V P and the energy flux Q V P are defined by and respectively. Poisson's equation shown in Eq. (37) is also used for deriving Eq. (B8). The local energy conservation law shown in Eq. (B8) agrees with that obtained by Qin et al. 39