Angular momentum and rotational energy of meanflows in toroidal magnetic fields

We derive the balance equation for the Favre averaged angular momentum in toroidal not necessarily axisymmetric magnetic field equilibria. These include tokamaks, stellarators, the reversed field pinch and field-reversed configurations. We find that the components of angular momentum are given by the covariant poloidal and toroidal components of \ExB and parallel flow velocities and we separately identify all relevant stress tensors, torques and source terms for each of these components. Our results feature the Favre stress generalisations of previously found Reynolds stresses like the diamagnetic or parallel \ExB stress. Further, we identify the magnetic shear as a source of poloidal \ExB angular momentum and discuss the mirror and the Lorentz force. Here, we find that the geodesic transfer term, the Stringer-Winsor spin-up term and the ion-orbit loss term are all part of the Lorentz force and are in fact one and the same term. Discussing the relation to angular velocity we build the inertia tensor with the help of the first fundamental form of a flux-surface. In turn, the inertia tensor is used to construct a flux-surface averaged rotational energy for \ExB surface flows of the plasma. The evolution of this rotational energy features a correction of previous results due to the inertia tensor. In particular, this correction suggests that density sources on the high-field side contribute much more to zonal flow energy generation than on the low field side.


Introduction
The double periodicity of a toroidal magnetic field configuration can be associated with two rotational degrees of freedom: toroidal and poloidal rotation. In a toroidally confined plasma both toroidal and poloidal rotation are observed and subject to intensive research.
Studies of toroidal rotation favour the toroidally symmetric tokamak case, where the symmetry leads to the exact conservation of the collective ‡ canonical angular momentum [1,2,3]. Of particular interest is the so-called intrinsic rotation, which refers to the ability of the plasma to spontaneously rotate without application of an external torque like neutral beam injection [4,5,6]. This is an important topic because toroidal rotation stabilizes the plasma against instabilities like the resistive wall mode.
The ideal toroidal symmetry of a tokamak is broken in stellarators and in reality also in tokamaks due to magnetic ripple effects from external field coils spacing. In fact, stellarator physics is different from tokamaks in some important aspects [7]. Neoclassical transport levels are much higher in a stellarator than in a tokamak even though stellarator optimization aims at reducing these levels down or below turbulent transport levels. More importantly however, the exact invariance of toroidal angular momentum is lost in a stellarator due to the lack of axial symmetry §. It is argued that in this case it is impossible for the plasma to rotate as fast as in (quasi-)axisymmetric devices [8] but that zonal flows may still develop.
Poloidal angular momentum, just as toroidal angular momentum, has two components in a general magnetic field, one stemming from the parallel velocity projected to the poloidal direction u ∥b ⋅ e ϑ , the other from the drifts perpendicular to the magnetic field u ⊥ ⋅ e ϑ (toroidal momentum analogously with e ϕ ). Here, u ∥ ≡ u ⋅b is the parallel flow velocity, u ⊥ ≡ b × (u ×b) is the perpendicular flow velocity,b is the magnetic unit vector and e ϑ and e ϕ are the covariant poloidal and toroidal base vectors. In reverse this means that both parallel velocity as well as the perpendicular drifts contribute to both toroidal as well as poloidal rotation. This is simply the geometrical observation that parallel and perpendicular directions versus poloidal and toroidal directions are different basis vectors for a flux-surface. This being said, the poloidal component of E × B velocity u E,ϑ ≡ u E ⋅ e ϑ gains significant interest because of its role in the formation of a transport barrier during the L-H transition [9,10,11]. The high confinement mode is accompanied by a narrow potential well just inside the ‡ after species and particle summation -individual particles exchange momentum through fluctuating electromagnetic fields § Axisymmery, axial symmetry and toroidal symmetry are used interchangeably throughout this manuscript.
separatrix of a diverted magnetic field geometry. The associated radial electric field drives a strongly sheared and flux-aligned E × B mean flow, which suppresses turbulence and thus reduces the radial flow of particles and heat out of the confined region. This E × B shear flow is believed to emerge out of turbulent fluctuations via the Reynolds stress, yet other mechanisms like the ion-orbit loss mechanism [12,13,14] or the Favre stress and background density gradient drive [15] are currently under discussion as well. Recent results suggest that the latter significantly alter the generation mechanism of E × B zonal flows for high density fluctuation amplitudes and steep density gradients [15,16]. The derivation of the separate evolution equations for the E × B and parallel velocity parts of both poloidal and toroidal angular momentum density is an open problem. Previous work is mostly restricted to toroidal symmetry [1,2,3], simplified magnetic field geometry [17,18,19,15] or delta-f modelling [17,18,19]. In this work we are interested in how the angular momentum anchors to the background magnetic field in the absence of a symmetry, what components appear in the complete stress tensor beside the ever present Reynolds stress and the impact of high fluctuation amplitudes and small gradient length scales.
It is instructive to introduce rotation also from a purely mechanical perspective. Consider a particle of mass m confined to a toroidal surface. Its Lagrangian reads L p = m(R 2φ2 + a 2θ2 ) 2 with the geometrical toroidal angle ϕ and poloidal angle ϑ. In an ideal torus the distance from the major axis R(ϑ) = R 0 + a cos ϑ, with R 0 the major radius, is independent of the geometric toroidal angle ϕ. The distance from the minor axis a = a 0 remains the minor radius a 0 . The Euler-Lagrange equations directly yield the conservation of toroidal angular momentumL ϕ = 0 with L ϕ = mR 2φ . This is a consequence of the independence of R and a of the toroidal angle ϕ. We then haveφ = L ϕ m(R 0 + a cos ϑ) 2 . Notice that the angular frequencyφ is higher on the torus inside ϑ = π than on the outside ϑ = 0, which we intuitively expect. In contrast, the equation for the poloidal angle is given by the nonlinear differential equation ϑ = −L 2 ϕ sin ϑ m 2 a 4 (R 0 a + cos ϑ) 3 . We observe that the poloidal angular momentum is not a conserved quantity in general. On a generally shaped toroidal flux-surface like that of a stellarator R as well as a depend on both ϕ and θ. Thus, neither toroidal nor poloidal angular momenta are conserved and ϑ and ϕ obey a coupled set of nonlinear differential equations.
In this contribution we calculate the toroidal and poloidal angular momentum balance separately for both the E × B and the parallel velocity part. The derivation rests upon a full-F gyro-kinetic formalism, where finite Larmor radius and polarization effects are taken in the long-wavelength limit. The main effect of the full-F formalism is the appearance of the Favre stress that replaces the conventional Reynolds stress and is a natural consequence of the density weighted flux-surface average -the Favre average [15]. The magnetic field geometry is arbitrary and we in particular do not invoke a toroidal symmetry. Thus, our results are applicable to fusion devices, such as tokamaks and stellarators. At the same time we make no assumptions on the form of the gyro-kinetic distribution function and our results thus apply and gyro-kinetic models as well. On the other side we allow a direct comparison to drift-reduced fluid equations due to the applied long-wavelength limit. We carefully recall the definition of angular momentum from the underlying particle Lagrangian in suitable coordinates and construct the inertia tensor with the help of the first fundamental form of general flux-surfaces. This enables us to then construct a rotational energy balance.
This manuscript is divided into the following parts. In Section 2 we review the magnetic field representation via flux-coordinates in order to setup suitable poloidal and toroidal angle coordinates. Our main derivation then rests on two pillars: (i) the definition of the gyro-kinetic action in Section 3, which encompasses our assumptions on the model and (ii) an ordering scheme presented in Section 4, which determines the long-wavelength limit. This enables us to then derive the poloidal and toroidal angular momentum balance. In Section 5 we apply the previously proposed Favre decomposition [15] in order to identify the signature of relative density fluctuations in both known and novel components of the stress tensor. In Section 6 we derive the relation between angular momentum and angular velocity and identify the inertia tensor. Furthermore, we find the time evolution for the rotational energy using the previously derived momentum balance. Finally, we discuss the significance of our results on various topics discussed in the literature in Section 7, including the electromagnetic field momentum, the ion orbit loss mechanism and the transition to simplified geometries. Appendix A provides a formulary intended as a quick reference list of the most often used relations and notations.

Preliminary: the magnetic field in flux-coordinates
A toroidal magnetic field equilibrium can be represented by so-called flux-coordinates {ρ, ϑ, ϕ} (Reference [20] calls them magnetic coordinates) where the magnetic field lines appear straight Here, ψ p (ρ) is the poloidal flux and ψ t (ρ) is the toroidal flux and we have dψ p = ι(ρ)dψ t where we introduced the rotational transform ι. Further, ρ is any radially increasing flux label, ϑ is the poloidal flux angle and ϕ is the toroidal flux angle coordinate. Note that ϑ increases in the counter-clockwise direction in the poloidal plane while ϕ increases clockwise if viewed from above to get a right-handed coordinate system. We emphasize that in general ϕ and ϑ are different from the geometric angles. In this manuscript we always refer to flux angles when speaking of the toroidal and poloidal angles or directions and will highlight when these angles coincide with the geometric angles.
There are many different toroidal flux coordinate systems, notably Hamada and Boozer coordinates [21,20].
In Fig. 1 we show an example of a numerically integrated [22] flux-coordinate system for an axisymmetric tokamak magnetic field. Here, we show the lines of constant ψ p in colour and the lines of constant poloidal flux angle ϑ in white. The toroidal flux angle ϕ coincides with the geometric angle. The magnetic field B 2 = dA 1 can be written as a total differential of the magnetic potential which notably identifies A ϕ = ψ p (ρ) and A ϑ = ψ t (ρ). At the same time dB 2 = d ○ dA 1 = 0 immediately as d ○ d = 0 for the exterior derivative d. This is the coordinate-free expression of vanishing divergence. We formulate Eqs. (1) and (2) in terms of differential forms, which we here introduce because the gyro-kinetic theory heavily relies on them (for an excellent introduction to differential geometry for physicists see Frankel's text [23]). An interesting (if somewhat aloof) property of using differential forms is that they (and therefore the magnetic field) can be defined without the existence of a metric tensor. Recall for example that the 1-form dϑ symbolizes the planes that are constructed by keeping ϑ constant and varying ρ and ϕ, which is a purely topological operation. In contrast, the gradient basis vector ∇ϑ is the vector that is perpendicular to the planes of constant ϑ, which requires a metric to define.
We are of course aware of the practicality that the physicist's notation of Eq. (1) provides We are here able to identify the poloidal B p ∶= ∇ψ p ×∇ϕ and toroidal B t ∶= ∇ψ t × ∇ϑ parts of the magnetic field vector B. With the choice of signs in Eq. (3) and assuming ∇ψ p points radially outwards, we get a lefthanded field-line winding when going in the positive ϕ direction since B p ∼ −e ϑ . Furthermore, notice the useful properties where e ϕ and e ϑ are the covariant basis vectors, that is the vectors that generate the directional derivatives along ϕ and ϑ, or in other words, e ϕ is the tangent vector to the line that we get when keeping ρ and ϑ constant and varying ϕ (e ϑ analogous). We emphasize that we mean these two vectors when we speak of toroidal e ϕ and poloidal e ϑ directions in contrast to the ∇ϕ and ∇ϑ directions. For example, in Fig. 1 e ϕ points perpendicularly out of the plane while e ϑ is tangent to the contours of ψ p (!) and in particular does not point in the same direction as ∇ϑ, which has component out of the flux-surface as well.
When we deal with a symmetric field independent of the geometric toroidal angle, we will choose ϕ as the geometric toroidal angle and keep ϑ as a fluxcoordinate with ∇ϑ ⋅ ∇ϕ = ∇ρ ⋅ ∇ϕ = 0 as we do in Fig. 1. This type of coordinates is known as symmetry flux or PEST coordinates [24]. Notice that we do not use the geometric poloidal angle since we want to keep the form Eq. (1). A useful property of this type of coordinate is that qR 2 √ g ≡ I(ρ) is a flux function, which allows us to write Last, note that all flux coordinates are problematic when an X-point with ∇ψ p = 0 is present in or close to the domain of interest. The poloidal flux ψ p is continuous and well-defined across the separatrix. However, the toroidal flux ψ t as well as the poloidal flux angle ϑ are only well-defined up to but not including or across the separatrix and furthermore ι −1 diverges on the separatrix. This is expected since the poloidal component of B vanishes at the X-point. In practice, the divergence manifests for example in Fig. 1 where the coordinate lines for ϑ are distorted when getting close to the separatrix on the low field side of the tokamak. We remark here that in fact any coordinate system with a flux label as the first coordinate is problematic when an X-point is present. In Reference [25] we find that a vanishing Laplacian of ψ p at the X-point is sufficient for the existence of a coordinate system but the volume element necessarily diverges at the X-point.
Last, we introduce the flux surface average (see for example [21]) as an average over a small volume -a differential shell centered around the flux-surface. We define where we define v(ψ p ) ∶= ∫ ψp 0 dV as the volume flux label. In flux coordinates we have dA = √ g ∇ρ dϑdϕ.
The average fulfills the identity Also note that for any divergence free vector field ∇ ⋅ j = 0 and a flux function f (ψ p ) we have which is proven straightforwardly. In summary, using flux coordinates for the following derivation defines suitable angle coordinates as well as poloidal and toroidal directions. We expect the resulting expressions to be valid for any flux coordinate system within the closed field-line region up to the separatrix. We remark that the numerical issues of flux coordinates close to the separatrix do not affect the theoretical results presented here.

Model definition
In this section we define our gyro-kinetic model and discuss the approximations that go into it. Our goal is to set up a model suitable for edge and scrape-off layer conditions. Literature on the derivation of gyro-kinetic models based on Lie-transform perturbation theory include the rather technical review [26] and references therein. A friendlier tutorial can be found in [27] or the more recent [28]. Here, we start directly with the gyrocentre Poincaré 1-form expressed in the transformed phase-space coordinates Z ∶= {X, w ∥ , µ, θ}, with gyrocentre coordinate X, parallel canonical moment w ∥ , magnetic moment µ, gyro-angle θ γ ∶= qA + mw ∥b ⋅ dX + m q µdθ (10) with species mass m and charge q and we omit the species label. We have the magnetic background potential A⋅dX ≡ A 1 from Eq. (2) and the background magnetic field unit vectorb ∶= B B.
In fluxcoordinates Eq. (10) explicitly reads We remark that this 1-form is already enlightening because it immediately identifies as the toroidal angular momentum and as the poloidal angular momentum. Recall here that angular momentum is defined as the canonically conjugate momentum to the angle coordinate. In anticipation of the following discussion we here remark that qψ p and qψ t will lead to the toroidal and poloidal components of the E × B velocity contribution. The parallel velocity contribution is given by the two components of the magnetic field unit vector b ϑ and b ϕ as expected. Unfortunately however, the definitions for toroidal and poloidal angular momentum in Eqs. (12) and (13) are not coordinate invariant and therefore care must be taken when comparing results from different coordinate systems. This is evident since the value of b ϕ and b ϑ depend on the choice of coordinates. Physically, we attribute this to different reference points/axes for the rotation that different angle coordinates entail. The symplectic 2-form, defined by the Poincaré 1-form, w ∶= dγ, defines the geometry of phase-space much the same way the metric tensor g defines the geometry of ordinary space. The difference is that ω defines areas instead of distances and is skewsymmetric instead of symmetric (see [23]). In 6dimensional phase-space coordinates we have Note the covariant vector components b i (withb T ∶= (b 1 , b 2 , b 3 )) and the appearance of the determinant of the metric tensor g in the definition of the cross- (17) Notice that the volume form is proportional to B * ∥ not just B * ∥ as often noted since it needs to remain positive. More importantly, it is apparent that the coordinate system possesses a (coordinate) singularity at mw ∥ = −qB (∇ ×b) ∥ , where B * ∥ = 0 and thus det(ω) = 0. This destroys the symplecticity of the 2-form ω, the volume form Eq (17) vanishes and the inverse of ω diverges (and thus the equations of motion). It is questionable how we can deal with this singularity especially when we later integrate over the phase-space volume to form the field equations. Furthermore, when deriving gyro-fluid models terms ∝ (B * ∥ ) −1 prevent identifying velocity space moments that involve B * ∥ in the volume element. This problem is often ignored in the literature or circumvented by requiring (∇ ×b) ∥ = 0 and we will follow this approach in this work. For a low-β stellarator ∇ × B = 0, however for general tokamak magnetic fields the requirement is only approximately fulfilled. As we will show in Section 7.1 the problem is also resolved by simplifying the magnetic field to purely toroidal or poloidal. Interestingly, the requirement (∇ ×b) ∥ = 0 relates to the integrability condition for vector fields perpendicular to the magnetic field. The Frobenius theorem [23] states that planes perpendicular tob(x) everywhere exist in the sense that there exist functions λ(x) and f (x) such that λ(x)b(x) = ∇f if and only if (∇ ×b) ⋅b = 0. In other words we surmise that the existence of drift-planes is a prerequisite for gyrokinetic and -fluid models.

Our Hamiltonian reads
with the effective gyro-centre potentials where we define the field Hamiltonian H f ∶= qΨ − qw ∥ A 1,∥ + q 2 A 2 1,∥ 2m to contain all terms dependent on the electromagnetic field perturbations φ and A 1,∥ . The potential φ is in fact a first order term where the zeroth order φ 0 has been neglected. The first order perturbation A 1,∥ is not to be confused with the zeroth order magnetic field potential A 1 . Finally, see Table A1 in Appendix A for definitions of ∇ ⊥ and ∆ ⊥ . Here, we follow [1,2] and use the Hamiltonian formulation with mw ∥ ∶= mv ∥ +qA 1,∥ such that the electromagnetic field variations appear in the Hamiltonian only and do not disturb the symplectic geometry (10). We note that we (i) neglect all terms k 3 ⊥ ρ 3 0 and higher in the Hamiltonian (this especially neglects the second order guiding centre contributions, which according to [29] leads to guiding centre drifts in the polarization equation). In particular, both the polarization contribution (the last term in Eq. (20)) as well as the finite Larmor radius effects are taken in the long-wavelength limit [30].
(ii) neglect compressional Alfvén waves entering through A 1,⊥ [31] (iii) neglect all terms non-linear in the magnetic potential A 1,∥ (except in the parallel kinetic energy). This approximation implies the absence of A 1,∥ terms in the polarization and of φ terms in the parallel Ampère law [31] equation and vice versa φ terms in the parallel Ampère law and therefore decouples the two equations, which is numerically desirable Our model is comparable to Reference [32] with the difference that we additionally take the longwavelength limit in the gyro-average operator. We note that with our approximations the Hamiltonian formulation with w ∥ is entirely equivalent to the symplectic formulation using v ∥ in the sense that the resulting equations are the same. The Hamiltonian Desirable might be an understatement. We are not aware of any successful attempts to numerically solve the completely coupled set of equations in a turbulence simulation.
formulation is more convenient here since γ is timeindependent. We also remark that the gyro-average and polarization corrections in our gyro-kinetic model Eq. (18) resemble the second order guiding centre transformation terms in guiding-centre models [33]. However, since we logically start with and approximate a gyro-kinetic model we will keep referring to our model as gyro-kinetic.
We introduce the gyro-kinetic particle distribution function F (Z, t) ≡ F (X, w ∥ , µ, t) (independent of gyro-angle θ, which is averaged out). The Vlasov equation states Here, S is a general kinetic source term on gyro-centre phase-space S(X, w ∥ , µ, t). With the kinetic source function S we formally model effects like for example non-elastic collisions, collisions with neutrals, heating of the plasma, or plasma fuelling. We call S a source understanding that it can act as a sink as well. Next, with the 1-form γ in Eq. (10) and the Hamiltonian H in Eq. (18) we can define a particle Lagrangian Together with the volume form in Eq. (17) and the phase space distribution function F we can then define the system Lagrangian L p ∶= ∑ s ∫ vol(Z)F (Z, t)L p (Z,Ż, t), where we sum over species. Finally, we close the system with a field Lagrangian and propose the action integral where dV ∶= √ gd 3 X is the spatial volume form. The action in Eq. (23) plus the Vlasov equation (21) are the central relations in every gyro-kinetic model. They completely define the system that we investigate. In particular this means that S contains all approximations to our model and that the following calculations are exact.
We remark that (i) the neglect of the electric energy E 2 in the field part of Eq. (23) leads to quasineutrality (that is a vanishing right hand side in Eq. (42)) (ii) the magnetic field energy in (23) neglects the A 1,∥ (∇ × B) ∥ contribution from the background magnetic field. This leads to the omission of the background equilibrium current in the Ampère equation (43). This approximation is in line with (∇ ×b) ∥ = 0.

The Vlasov-Maxwell equations
In the Lagrangian picture [34] the equations of motion can be retrieved from the action Eq. (23) by expressing Z = Z(Z 0 , t 0 ; t), using F (Z, t) = F 0 (Z 0 , t 0 ) by the Vlasov equation (21), taking the integration to the initial positions and time ¶ and then varying δS δZ = 0. This indeed recovers the Euler-Lagrange equations d dt The application of the Euler Lagrange equations (24) yields the Hamilton equations of motion (using where we define Z i H ≡Ż i as the components of the Hamiltonian vector field on phase space. Here, i Z is the inner product with the vector field Z H and d is the total differential. The particle trajectories are given by the streamlines of Z H (with J ∶= ω −1 ) The time-derivative of any phase-space function along the trajectory is then given by Here and in the following we useŻ synonymously with Z H . In particular, the derivative of the Hamiltonian gives where we use Eq. (26) and the antisymmetry of J. Explicit expressions for the inverse of the symplectic 2-form Eq. (14) and the gradient of the Hamiltonian (18) are The µ component of ∂H contains corrections due to ¶ Technically, here we also need to know that the volume form is conserved in time B(z)d 6 Z = B(z 0 )d 6 z 0 , something that we will need to show explicitly.
the fluctuating electric field B # ∶= B + m∆ ⊥ φ 2q 2 B. An explicit expression for the components of Z H (oṙ Z) can now be formed given Eqs. (29) and (30) The phase space volume vol = ω ∧ ω ∧ ω is conserved along the particle trajectories where the Lie derivative on differential forms is given by Cartan's formula [23] The conservation of volume thus translates to a vanishing divergence of the Hamiltonian vector field in phase space Notice that volume conservation does not mean that B is conserved along particle trajectories, we rather havė B =Ẋ ⋅ ∇B. The conservation of the particle distribution function F (X, v ∥ , µ, t) is expressed by the gyro-kinetic Vlasov equation dF dt = S, which together with phase space volume conservation (37) reads in conservative form The Vlasov-equation Eq. (38) together with the equations of motion Eq. (31)-(33) forms the first half of the Vlasov-Maxwell system.
In order to derive the Maxwell equations we first define the velocity space moment operator [35] ζ ∶= dw ∥ dµdθm 2 BF ζ (39) where ζ(X, w ∥ , µ, t) is any function defined on phasespace and the integration encompasses the entire velocity space. Notice that we name the first few fluid moments N ∶= 1 , N W ∥ ∶= w ∥ and P ⊥ ∶= µB and give a comprehensive list in Appendix A.2.
We also define the moment operator for the source function S analogous to the velocity space moment operator for the gyro-kinetic distribution function F in Eq. (39) Analogous to the moments of F we name the source moments S N ∶= 1 S , S P⊥ = µB S , etc. Using Eq. (38) together with the fact that ∂ ∂t and ∇ commute with the velocity integral and F vanishes for w ∥ = ±∞ we find the important identity [35] The variation of the action (23) with respect to φ(x) yields the quasi-neutrality equation and with respect to A 1,∥ the parallel Ampère law where we used that γ does not depend on either φ or Notice the appearance of F B inside the divergence/Laplace operators. Carrying out the variations with the help of Eq. (44) in the polarization and Ampère equations (42) and (43) and identifying the velocity space moments (39) yields + ⋅ M gy and the gyro-kinetic polarization and magnetization densities Note that the parallel component of the polarization current j pol ⋅b = ∂P gy ⋅b ∂t = 0 vanishes in Eq. (46). Also, the parallel part of the magnetization density M gy ∥ ∶= − µ b ≡ −P ⊥b B does not contribute to the parallel magnetization current.
In total, we now explicitly derived the equations of the Vlasov-Maxwell system in Eq. (38), (45) and (46) together with the equations of motion in (31)-(33).

Interlude: relation between gyro-fluid and fluid moments
Gyro-fluid quantities like 1 = N (X, t) or v ∥ = U ∥ (X, t) are given in gyro-centre coordinates X and are thus not directly comparable to the physical fluid quantities, which we denote with lower case letters is the distribution function in particle phase-space (and we here overburden the use of v as the velocity on top of the volume flux-label).
We need to use the gyrokinetic phase-space coordinate transformations to transform between particle and gyro-kinetic phasespace moments. Helpfully, Reference [26] relates the coordinate transformation to the variational derivative of the action. With our action (23) we obtain where ξ is ζ transformed to particle coordinates and ξ v ∶= ∫ d 3 vξf is the particle phase-space moment operator. Thus, ξ v is the physical fluid moment corresponding to ζ . In Eq. (49) we immediately see + The attentive reader will notice that Eqs. (45) and (46) are only semi-elliptic since the projection tensor h is only positive semi-definite. Concerns about existence and uniqueness of solutions are dealt with under "degenerate partial differential equations" in the mathematical literature. In particular, the field of stochastic differential equations contains a solution to the Dirichlet problem, see for example Reference [36]. that the actual fluid moment ξ v equals the gyro-fluid moment ζ up to an order O(ρ 2 0 k 2 ⊥ ) correction. For example the density transforms as The right hand side terms appear exactly in the polarization equation (45), which we obtained from the variational principle. This shows that Eq. (45) is the gyrofluid version of quasineutrality ∑ s qn = 0. It is possible to invert the relation between gyrofluid and fluid quantities. We follow [37,38] and explicitly express the first two gyro-fluid quantities N and U ∥ in terms of the true fluid quantities n and u ∥ in the long-wavelength limit up to order (ρ 0 k ⊥ ) 2 .
Note that we neglect the potential part in Eq. (52) since we miss the corresponding term in the Hamiltonian. The moments of S transform back to particle phase space analogous to Eq. (49). This is because the coordinate transformation works for any phasespace function, not just the distribution function F . For example, we have where S n is the true fluid particle source term. We are now able to formulate the only constraint we have for the source term namely that it should conserve the total electric charge via where we define the polarization source Note that Eq. (54) is completely analogous to Eq. (45).

Poloidal and toroidal E × B momentum
With the model developed in Sections 2 and 3 we are now ready to start the derivation of the balance equations for the angular momentum density. We begin by computing the time derivative of qA ϕ = qψ p , which is the first part of the toroidal angular momentum (12) where we separated the field Hamiltonian H f . Now, to simplify the right hand side of Eq. (56) we need to relate the variational derivative to ordinary derivatives. Consider a generic Hamiltonian In order to proceed we need to assume that η ⋅ ∇ commutes with ∇ ⊥ and ∆ ⊥ . This is not true in general which is why we now revert to a long-wavelength or drift-kinetic ordering [19,39], where we order (i) the frequency of turbulent fluctuations compared to the ion gyro-frequency as small This in particular orders the E × B velocity compared to the ion thermal velocity as u E c s,i ∼ δ where c s,i = T i m i (iii) all derivatives on the magnetic field (vectors) as that is parallel derivatives on the magnetic field variation scale. This orders k ∥ k ⊥ ∼ δ 2 Note that Reference [28] orders ρ i L B ∼ δ 4 . However, this would completely neglect all curvature terms in our scheme. In our ordering the Hamiltonian H f (18) appears to be second order and we now neglect all terms of order δ 4 on the right hand side of Eq. (56).
With the above orderings we directly have that This equation is a useful identity and in fact holds for any vector field that commutes with ∇ ⊥ and ∆ ⊥ . It links the ordinary derivative on H to the variational derivative and correction terms that appear as exact divergences.
Summing over all species, integrating over velocity space and inserting our Hamiltonian from Eq. (18) we get Now, we focus on the term qψ p = qẊ ⋅ ∇ψ p on the left hand side of Eq. (56). First we insert q into the velocity space moment equation (41). We find ∂ t q + ∇ ⋅ qẊ = q S . Under species summation we see that we can identify the polarization equation (45) ∑ s qN = ∇ ⋅ P gy and analogously the quasineutrality for the sources Eq. (54) ∑ s qS N = ∇ ⋅ S P . The next step is to apply the flux-surface average Eq. (7) to obtain ∂ v ∂ t ⟨P gy ⋅ ∇v⟩ + ∑ s ⟨ qẊ ⋅ ∇v ⟩ − ⟨S P ⋅ ∇v⟩ = 0. After volume integration ∫ v 0 dv (the inner integration boundary vanishes) and multiplying with dψ p dv we obtain which recovers the radial part of the polarization current j pol ≡ ∂P gy ∂t. We stress that Eq. (59) is an important identity [1]. It links the derivative of the poloidal flux or in fact the first part of the toroidal angular momentum of particles to the polarization current and sources. Now, we further investigate the terms appearing from Eq. (59) by explicitly inserting our polarization density (47) The second term can be simplified using the dynamical One key ingredient for Eq.
⟩ in the longwavelength limit. Now, we add the terms Eq. (58) and (60) and use our ordering to eliminate the magnetic field derivatives to get where we use the abbreviation ∇ v = ∇v ⋅ ∇ and ∇ η = η ⋅ ∇ and imply species summation to present this intermediate result. We also used that η commutes with ∇ ⊥ in the long-wavelength limit and that the fluxsurface average of ∇ ⋅ (ηh) vanishes exactly. Furthermore, we replace the gyro-centre quantities by their particle analogons, which is possible in the long-wavelength limit (see Eqs. (51)). Finally, the term A 1,∥ v ∥ K κ ⋅ ∇ψ p in Eq. (56) vanishes under species summation and the parallel Ampère law in the longwavelength limit. With the help of Eq. (4) we then finally arrive at Here, we define the E × B drift u E , the grad-B drift u ∇B , the diamagnetic drift u D , the curvature drift u κ , the first order magnetic fluctuations B 1,⊥ and the electromagnetic magnetization density M em ⊥ (different from M gy by a factor 2, fluid instead of gyro-fluid quantities and the long-wavelength limit) and b 1,⊥ ∶= B 1,⊥ B. Equation (61) describes the evolution of the toroidal E × B angular momentum density and is the first result of this paper. The second term on the left side is the average over the convective acceleration term with radial velocity u v E +u v D , the sum of E ×B and diamagnetic velocity. In Section 5 we will show that this term can be split into an advective part and components of the turbulent stress tensor. The remaining terms on the left hand side are two stress terms stemming from magnetic fluctuations. On the right hand side the Lorentz force originating from the "free" current j f ∶= ∑ s qn(u κ + u ∇B ) appears. The source term signifies that particles are created with instant E × B velocity, which adds to the angular momentum.
Another point we note is that the poloidal analogue of Eq. (61) follows immediately. Recall Eq. (5) together with ι∇ψ t = ∇ψ p in flux coordinates. This yields u E,ϑ = ∇φ⋅∇ψ t B 2 = ι −1 u E,ϕ and thus from Eq. (61) directly follows the equation for the poloidal E × B angular momentum density Equation (65) exhibits a similar structure as Eq. (61) with the additional appearance of the magnetic shear σ ∶= ∂ v ln ι [21]. Depending on its sign the shear term can both dampen and generate poloidal E × B angular momentum. Before we continue with the identification of the various stress terms in Section 5 and a more detailed interpretation of our results, we first derive the equations for the remaining angular momentum components, namely the poloidal and toroidal angular momentum components stemming from u ∥ . As it turns out we will get the full parallel momentum balance as a by-product. Finally, recall that both u E,ϕ = ∇φ ⋅ ∇ψ p B 2 and u E,ϑ = ∇φ ⋅ ∇ψ t B 2 are related to the radial electric field, a fact that will lead to the identification of the electromagnetic field angular momentum density in Section 7.2.

Parallel (angular) momentum
We now turn to the parallel terms in the toroidal canonical momentum γ ϕ = qψ p + mw ∥ b ϕ as well as the poloidal canonical momentum γ ϑ = qψ t + mw ∥ b ϑ . Repeating the ordering scheme from the previous section one could assume that ∂u ∥ ∂t ∼ δ 2 and argue that therefore only O(δ 2 ) terms should be kept in our ordering. However, we note that the ions accelerate very slowly. This is because nu ∥,i ∼ ∇ ∥ p i ∼ δ 3 . Note that this requires ∇ ⊥ u ∥,i to be small as well. In contrast, the electron velocity is mainly determined by with η ∥ being the parallel resistivity. In order to reflect these considerations we order This ordering mandates that the terms ∂ t m i u ∥,i ∼ ∂ t m i u E,ϕ ∼ δ 3 are similar in size. The parallel ion velocity itself is larger than the E × B velocity but we order its time derivative smaller by the same factor. In total, we again do not assume toroidal symmetry but we do assume a long-wavelength limit and keep terms up to O(δ 3 ). We start with (for η ∈ {ϕ, ϑ}) With the vector triple product rule applied to (b × ∇H) × B * we see mẇ ∥b = q(Ẋ × B) + mw ∥Ẋ × (∇ × b) − µB∇ ln B − ∇H f . Next, we noteẊ ⋅ ∇b η = −(Ẋ × (∇ ×b)) η +Ẋ i ∂ η b i (Notice that we do not use the covariant derivative here since b η ≡b ⋅ê η is a scalar quantity and (a ⋅ ∇b) η ≠ a ⋅ ∇b η ; the first is the component of a covariant derivative while the second is the directional derivative of b η ). Finally, we have The second identity follows immediately from the equations of motion (32). In the long-wavelength limit (and under species summation to make qv ∥ A 1,∥ vanish) we can write mw ∥Ẋ we interpret as a generalized curvature contribution.
With similar arguments as in the previous section we can recover the variational derivatives in ∑ s b η ∇ ∥ H f using Eq. (57). However, the remaining terms are all O(δ 5 ) and can be safely neglected in our ordering.
Taking the velocity space moment we arrive at We note that m ∂ Note here that the curvature terms vanish under the divergence in the long-wavelength limit. As a final step we again apply the flux-surface average and note that with implied species summation the term ∂ qA 1,∥ ∂t vanishes using the polarization equation. Then we have while the average parallel momentum reads The two components of Eq. (72) complement the previously derived E × B velocity components in Eqs. (65) and (61).
In total, Eqs (61), (65), (72) and (73) form the basis of our discussion for the remainder of the manuscript.

Favre averaged momentum equations
In order to distinguish the effects of turbulent fluctuations and flux-surface averaged quantities a Reynolds decomposition is traditionally used to rewrite the averaged evolution equations. For any function Unfortunately, as we point out in Reference [15] the Reynolds decomposition technique does not lead to well-behaved terms when the absolute density n appears in the nonlinear terms in the sense that (i) absolute density fluctuationsñ appear instead of relative density fluctuationsñ ⟨n⟩, (ii) the radial advective part is not correctly recovered and (iii) effects from the density gradient ∂ v ln ⟨n⟩ are not evident. We will thus follow [15] and introduce the so-called Favre decomposition. Consider a term of the form ⟨nh⟩. If we multiply and divide by ⟨n⟩, we can write ⟨nh⟩ ≡ ⟨n⟩ h . Here, we introduce the so-called Favre average h ∶= ⟨nh⟩ ⟨n⟩ which can be understood as a density weighted Reynolds average. We note that this definition is species dependent through the dependence on the species density n. The Favre average then allows the definition of the Favre decomposition The Favre average reduces to the Reynolds average for small fluctuation amplitudes or if the density is a fluxfunction h = ⟨h⟩ + ⟨ñh⟩ ⟨n⟩ ≈ ⟨h⟩. Reference [15] also reported u ϑ ≈ ⟨u ϑ ⟩ within a few percent since u ϑ ≈ 0 in gyro-fluid simulations.

Favre averaged covariant E × B velocity
We first apply the Favre average technique to the continuity equation ∂ t n + ∇ ⋅ (nu) = S n to get where we define the average radial velocity and we use If we now replace all terms of the form ⟨nh⟩ with ⟨n⟩ h in Eq. (61) and use gh = g h + ĝĥ (Eq. (A.10)) and u v D = 0 (long-wavelength limit), we get and similarly in Eq. (65) we get Note that the density ⟨n⟩ is species dependent and therefore we cannot divide Eqs. (79) and (80) by ⟨n⟩.
What is usually possible is to neglect the electron mass, which reduces the sum Eqs. (79) and (80) to a sum over all ion species. Equations (79) and (80) describe the evolution of the Favre averaged covariant components of the E ×B velocity in general, not necessarily axisymmetric magnetic field geometry in the long-wavelength limit. On the left hand side we find a radial advection term proportional to U v [15]. The first term on the right hand side is the total perpendicular stress T v ⊥,η , which consists of the perpendicular Favre stress F v ⊥,η and the Maxwell stress M v η . We note here that we define the Favre stress as a kinematic stress ( "stress divided by mass density") with units m 2 /s 2 as opposed to the Maxwell stress which has units of stress N/m 2 .
The kinematic Favre stress F v ⊥,η contains the E × B Favre stress F v E,η . As Reference [15] points out the E × B Favre stress F v E,η can be written as 40] and the often neglected [5] triple term appears on the right-hand side of Eq. (85). An advantage of the Favre decomposition is that the density fluctuations are automatically contained as relative fluctuation levels as is evident in the triple term in Eq. (85). An analogous identity to Eq. (85) holds for the diamagnetic Favre stress F v which encompasses the diamagnetic Reynolds stress R v D,η ∶= ⟨ u E,η u v D ⟩ [19]. Note that the diamagnetic Favre stress is asymmetric in contrast to E × B Favre stress. In this form the diamagnetic stress consist of the radial component of the diamagnetic velocity together with the η component of the E × B velocity. This is a consequence of using the pressure equation to evaluate the time-derivative of the diamagnetic velocity [19], which we have done using Eq. (60). Otherwise the transpose of the diamagnetic stress consisting of the radial E × B component and the η component of the diamagnetic velocity appears [41]. We elaborate further on different interpretations of the E×B angular momentum density in Section 7.2. In addition to F v E,η and F v D,η we find the stress term F v F,η that appears for fluctuating magnetic field b v 1,⊥ . This term is in fact a remainder of the actual magnetic flutter Favre stress term m u E,η u ∥ b v 1,⊥ that would appear, had we not neglected the A 1,∥ nonlinearities in the Hamiltonian (18) (through the variaton in Eq. (57)). We expect F v F,η to vanish for small relative density fluctuations and to only play a role for O(ñ ⟨n⟩) ∼ 1 fluctuation amplitudes, due to the similar dependence as the second term in the Favre stresses [15].
The Maxwell stress M v η consists of the symmetric vacuum field contribution M v B,η and the asymmetric magnetization stress term M v M,η . The role of the vacuum Maxwell stress M v B,η on the generation of sheared E × B flows was highlighted previously in for example [42,43]. The novel asymmetric magnetization stress M v M,η appears in its present form analogously to the diamagnetic stress as a consequence of using the pressure equation (60). In Section 7.2 we will encounter its transpose in the full electromagnetic field stress tensor. It notably contains a contribution from the parallel heat flux q ∥ + p ⊥ u ∥ and physically originates in the magnetization term in parallel Ampère's law Eq. (46).
As was highlighted in [15] the density gradient ∂ v ln ⟨n⟩ contributes to the evolution of E × B shear flow. Consider We emphasize that both E × B and diamagnetic Favre stresses appear in the term m ⟨n⟩ F v ⊥,η ∂ v ln ⟨n⟩ and that this term is non-zero even if ∂ v F v ⊥η vanishes. This is particularly interesting for the steep density gradient that develops during the transition to H-mode.
Contrary to the toroidal angular momentum density, the poloidal E × B angular momentum density in Eq. (80) is influenced by a gradient in the rotational transform profile or magnetic shear σ = ∂ v ln ι. The magnetic shear is known to influence the E × B shear flow evolution [44,45]. In particular the shear dampens drift-wave turbulence and leads to narrow zonal flows [45]. Furthermore, it dampens the Kelvin-Helmholtz instability, which would otherwise be driven by the E × B velocity shear [44]. In Eq. (80), we explicitly identify two magnetic shear contributions. The first shear term m ⟨n⟩ U v u E,ϑ σ corresponds to roughly exponential growth or damping of poloidal flows, assuming that the average radial velocity is constant (which is a good estimate since it is approximately the E × B radial particle transport). The second shear term appears analogous to the density gradient ∂ v ln n term and contributes even if ⟨n⟩ F v ⊥,ϑ and M v ϑ are "radially" homogeneous (no volume derivative).
On the right hand side of Eq. (79) and (80) we further find the components of the Lorentz force originating from the radial curvature drift current j f defined in Eq. (A.21). The grad-B induced current part of this term is the Stringer-Winsor spin-up term [46,47,48]. In order to see this recall that (j f × B) ϕ = j f ⋅ ∇ψ p ∼ p ⊥ K(ψ p ) that is the radial component of the free current (see A1 for the definition of the curvature operator K). The same term was found in δF drift-fluid models [17,43,18] and was there called the geodesic transfer term. In any case the term is known to excite geodesic acoustic modes and to both dampen or drive zonal flows depending on the parameter regime [48,17]. We further discuss this term in relation to the ion orbit loss mechanism in Section 7.3.
Finally, on the right side of Eq. (79) and (80) we find source related terms contained in S u E,η defined in Eq. (84). The term ⟨ S n u E,ϑ ⟩ in Eq. (80) describes the poloidal spin-up mechanism for poloidally asymmetric particle sources described in [47]. In Eq. (79) we find an equivalent term also for the toroidal E × B velocity. A toroidally asymmetric particle source can generate or dampen toroidal E × B velocity. The second source term is proportional to the difference between Reynolds and Favre averaged E × B velocity ⟨u E,η ⟩ − u E,η . For small density fluctuations we thus expect this term to vanish and only contribute for large fluctuation amplitudes.

Favre averaged parallel velocity
For the parallel angular momentum components (72) we have where U v is given in Eq. (78) and we identify With the Favre average we re-write Eq. (73) into where U v is given in Eq. (78) and we identify Analogous to Eq. (79) in Eqs (93) and (88) we find a radial advection term of momentum by U v followed by various stress terms contained in T v ∥ respectively T v ∥,η . Again, we define the Favre stress as a kinematic stress and analogous relations to Eq. (85) hold for the parallel Favre stress components. The parallel E × B Favre stress F v E,∥ respectively F v E,∥,η now depends on fluctuations in the parallel velocity instead of E × B velocity. The Reynolds stress analogue of F v E,∥ is well-known in the literature on intrinsic toroidal rotation (see e.g. [5]), however we point out here that F v E,∥,η is the actual component that drives angular momentum u ∥ b η instead of just parallel momentum u ∥ . The parallel magnetic flutter Favre stress term F v F,∥ respectively F v F,∥,η is a transfer term appearing for magnetic fluctuations b 1,⊥ . The kinetic stress term K v ∥ respectively K v ∥,η is related to the kinetic dynamo mechanism as for example discussed for the reversed field pinch in Reference [49,50]. On the right hand side we find the mirror force term − ⟨p ⊥ ∇ ∥ ln B⟩ respectively − ⟨b η p ⊥ ∇ ∥ ln B⟩.
In the equation for the parallel angular momentum Eq. (88) we find an additional geometrical correction to the mirror force. Finally, the momentum source term S u ∥ respectively S u ∥ bη represents angular momentum generation by external sources. Note that with the definition of a velocity source S u ∥ via S nu ∥ ∶= nS u ∥ + u ∥ S n we can write and analogous for S u ∥ bη . Eq. (98) now consists of the Favre averaged velocity source plus a contribution from a poloidally asymmetric source term analogous to Eq. (84). We comment here on the appearance of the Lorentz force in the equation for the parallel angular momentum Eq. (88).
The Lorentz force acts perpendicularly to the magnetic field line and should not contribute to the parallel momentum at all. Indeed, we can further simplify the right hand side of Eq. (88) to where we defineκ η ∶= e η ⋅ κ − b i ∂ η b i with the curvature κ ∶=b ⋅ ∇b. Now, only the component of the mirror force −p ⊥ b η ∇ ∥ ln B and a geometric correction term appear. To see the mirror force recall the sign of the magnetic moment vector µ = −µb and the guiding centre parallel magnetization density [26,29] M gy ∥ = µ = − µ b = −P ⊥b B. The force acting on magnetic dipoles is [51] which is what appears in Eq. (73), while F d,η = −P ⊥ ∂ η ln B appears in (72). The perpendicular part gives rise to the ∇B drift. Last, notice that ∇ ∥ ln B = −∇ ⋅b such that Pressure fluctuations are required to affect the the angular momentum generation via the mirror force.

Favre averaged total angular momentum density
The velocity equations (79)/(80) and (88) can be easily cast back into conservative form using the continuity equation (77). Summing up the results, we finally find the evolution of the total average poloidal and toroidal angular momentum density where δ is the Kronecker delta. The magnetic shear term only contributes to the poloidal angular momentum. The convective term proportional to U v vanishes under volume integration up to a surface contribution as does the total stress term T v ⊥,η + T v ∥η . In Eq. (102) we further find that the momentum transfer to the background magnetic field is mediated by the mirror force and the generalized curvature force term on the right hand side. Clearly, the Lorentz force term cancels in the total angular momentum density evolution. Finally, we recover the external source terms on the right hand side.
We see that the total angular momentum in Eq. (102) is given by the covariant components of the E × B and parallel velocities. Comparing this to the total advection velocity u ∶= Ẋ = u E + u κ + u ∇B + u ∥b + u ∥ b 1,⊥ that appears in the continuity equation ∂ t n+∇⋅(nu) = S n we see that the curvature, grad-B and magnetic flutter velocities do not appear in the angular momentum (102) even though we at least expected the magnetic flutter term as an order O(δ) term. At this point recall Eq. (59), which identifies the radial polarization current with the macroscopic expression for the angular momentum density (except u ∥ b η ). The polarization density P gy is directly connected to the definition of the Hamiltonian (18) through the variational principle. Since we neglected the second order guiding center corrections we accordingly miss the guiding center polarization density [2,29] and thus the corresponding curvature terms in our angular momentum density. On the other hand we also neglected the nonlinear terms in A 1,∥ in the Hamiltonian, which accounts for the missing magnetic flutter velocity mv ∥ b 1,⊥ in the polarization [26] and thus angular momentum density (102).

Angular momentum and angular velocity
In Section 5 we have derived equations for the covariant components of the E × B and parallel velocity, which add up to the total angular momentum density in Eq. (102). We now focus on the angular momentum as a vector quantity. We define We are now interested only in the part of the flow that stays within a given flux-surface, because this flow can be constructed from the covariant ϕ and ϑ components of u L that we have available. To see this, we formulate the projection tensor onto the flux surfaces with the contravariant radial unit vectorρ ∶= ∇v ∇v . With this we can split the flow velocity according to u L = u L ψp + uρ Lρ where we define the surface or rotational velocity where we follow [24] and introduce the surface operator ∇ S ∶= h S ⋅ ∇. We thus have L i = u E,i + u ∥ b i for i ∈ {ϕ, ϑ, ρ}. As expected we do not need the radial component of u L to construct the surface flow in Eq. (105).
It is now important to see that ∇ S ϕ and ∇ S ϑ form the contravariant basis of the flux surface as a stand-alone manifold and analogous e ϕ and e ϑ are its covariant basis vectors In fact, explicitly writing h S into components we realize that all components h S,ρk vanish for k ∈ {ϕ, ϑ, ρ}. We thus define I as the two-dimensional tensor consisting of the non-zero components of h S , that is The interested reader will recognize I as the the first fundamental form of flux surfaces parameterized with ϑ and ϕ. The first fundamental form I can be interpreted as the two-dimensional metric tensor of the flux-surface thought as a standalone structure and is thus an intrinsic structure of the magnetic flux surfaces (and in particular has a well-defined expression in every coordinate system). Unfortunately, the flux-surface average is not an intrinsic surface operation since it requires the knowledge of the volume form √ g to compute. Also, note that the components of I and its inverse are given by I ϕϕ = e ϕ ⋅ e ϕ , I ϑϕ = e ϑ ⋅ e ϕ , I ϑϑ = e ϑ ⋅ e ϑ and I ϕϕ = ∇ S ϕ ⋅ ∇ S ϕ, I ϑϕ = ∇ S ϑ ⋅ ∇ S ϕ, I ϑϑ = ∇ S ϑ ⋅ ∇ S ϑ respectively. Now, the fundamental form I has another interpretation, namely as the inertia tensor of rotations in ϑ and ϕ. To see this recall that the contravariant components of the surface velocity L, L ϕ and L ϑ are actually the angular velocities with units s −1 , because the particle trajectory is given byφ = L ϕ andθ = L ϑ . In contrast, the covariant components L i = I ij L j for i, j ∈ {ϑ, ϕ} form the angular momentum as it results in Eq. (102) that is mL i has units kgm 2 s −1 . This leaves mI as the (kinematic) inertia tensor that connects the angular velocity and angular momentum of a fluid element rotating on a flux-surface.

Mean and fluctuating angular momentum
Consider now the mean surface velocity field generated by Favre averaged covariant ϕ and ϑ velocity components The time evolution of m ⟨n⟩ L m is directly given by Eq. (102). First, we emphasize that the corresponding angular velocity components of L m , L i = I ij L i are not flux functions since the inertia tensor does not commute with the flux-surface average and thus L m ≠ L ϑ e ϑ + L ϕ e ϕ or in other words, if angular momentum is a flux-function then angular velocity cannot be at the same time. In fact, we perform the splitting L i = L i +L i expecting that the relative fluctuationsL i are small and that u i is well-described by its Favre average L i . A priori, these arguments of course also hold the other way, if angular velocity were a flux-function then angular momentum cannot be at the same time and we should split the angular velocities. At this point recall the discussion in the introduction. When angular momentum is conserved, a particle moves faster closer to the axis (for example on the high field side in Fig. 1). We take this as an indication that angular velocities are not well-described by flux-surface averages, while angular momenta are. Furthermore, in the equations in Section 5 (for example Eq. (79)) we see that the average angular momentum is fed by turbulent fluctuations through the stress tensor, which we interpret as an indication that fluctuationsL i and notL i become small.

Total energy evolution
Before we construct a zonal or mean flow rotational energy we first focus on the total energy evolution of our system. We follow Reference [37] and derive the pressure equations (the thermal energy) for p ⊥ and p ∥ directly from the moment evolution equation (41). We point out that we need to keep terms one order higher in the energy conservation law than in the momentum conservation law, that is O(δ 4 ) in our ordering. This is due to the fundamental property of the gyro-kinetic system [26] that a higher order Hamiltonian needs to be kept in the system to obtain polarization effects and an exact energy invariant. If we thus neglect all terms of order O(δ 5 ), use parallel Ampère's law (46) and apply the species summation we get We formally summarize all total divergences into the term j v E,p . An interesting side-remark here is to view the energy conservation Eq. (108) to lowest order, which leaves Bernoulli's identity ⟨p ⊥ + p ∥ 2 + mu 2 ∥ 2⟩ = const along fluid trajectories. On the right side of Eq. (108) appears the energy exchange term ⟨j f ⋅ E ⊥ + j ∥ E ∥ ⟩ as well as the pressure source terms (heating).
On the other side using the definition of Ψ in Eq. (20) and the polarization equation (45) we find which recovers the E × B kinetic energy density in the last term on the right hand side. Interestingly, a completely analogous relation holds for the term qΨ S (by replacing µB with µB S and N with S N in Eq. (109)) since we require the sources to preserve quasineutrality in Eq. (54). Applying Eq. (41) to qΨ = q∂ t Ψ + qẊ ⋅ ∇Ψ and using (109) and (58) for Ψ under species summation and neglecting again terms of order O(δ 5 ) the result is given by where we identify the total mass density ρ M ∶= ∑ s mn since the E × B velocity is species independent and again summarize all divergences into the formal j v E,ψ term. The density source S n either generates or destroys kinetic E × B energy depending on its sign. The sum of Eqs. (108) and (110) recovers the conservation of the flux-surface averaged total energy of our model since the energy exchange term ⟨j f ⋅ E ⊥ + j ∥ E ∥ ⟩ cancels.

Mean rotational energy evolution
The direct approach to a rotational energy density is the kinetic energy of the surface flow velocity L This energy is equivalent to subtracting the radial E × B energy ⟨ρ M u E,v u v E 2⟩ from the total kinetic energy density ∑ s ⟨mn(u 2 E + u 2 ∥ ) 2⟩. It is now important to realize that contrary to the parallel kinetic energy the E × B rotational energy density can be related to the (species summed) angular momentum evolution. This is because the E × B drift velocity is equal for all species. We can write where we define and here introduce the total mass density in the Favre averages for any (possibly species dependent) function h. If h is species independent Eq. (115) simplifies to h M = ⟨ρ M h⟩ ⟨ρ M ⟩. With u E,ϑ = ι −1 u E,ϕ we can simplify further Here, we introduce the inertia factor I 0 ∶= (ι −1 , 1)I −1 (ι −1 , 1) T . For a purely toroidal magnetic field we have I 0 = R −2 as expected. For symmetry flux coordinates we have g ϑϑ = R 2 (∇ψ p ) 2 I 2 ι 2 , g ϕϑ = 0 and g ϕϕ = R 2 and thus I 0 = R −2 (1+I 2 ∇ψ p 2 ) = B 2 ∇ψ p 2 . The inertia factor vanishes for a slab magnetic field. In this case our zonal flow energy agrees with [15] and in the case of small density fluctuations also with its δF analogue [43,18]. Since I 0 is time-independent we can use the evolution equations for the density Eq. (77) and angular momentum (79) to get where we neglected the term ⟨nu ⋅ ∇I 0 ⟩ in the continuity equation as small in our ordering and we have In Eq. (117) we find the term u v M as the convective velocity for the zonal flow energy. On the right side the derivative of the total perpendicular stress T v ⊥,ϕ given by Eq.
(87) appears. Thus, the Favre and Maxwell stress given by fluctuating velocities and the fluctuating magnetic field in Eqs. (82) and (83) respectively together with a gradient in the density ∂ v ln ⟨n⟩ can appear as sources for zonal flow energy. The E × B Favre stress F v E,ϕ was already identified as a source for zonal flow energy in a slab geometry in [15]. The vacuum field Maxwell stress M v B,ϕ and the E × B Reynolds stress R v E,ϕ (contained in our F v ⊥,ϕ according to Eq. (85)) appear in similar form in δF models [43,18]. Compared to these previous findings we find the additional appearance of the diamagnetic Favre stress F v D,ϕ contained in F v ⊥,η and the magnetization stress M v M,ϕ contained in M v ϕ . In addition, we find the inertia correction factor I 0 M that vanishes only in the simple slab geometry. On the right hand side of Eq. (117) we further find the Lorentz force term, which includes the geodesic transfer term. This term represents an energy transfer to the internal energy density Eq. (108) since we know the Lorentz force to transfer angular momentum to the parallel angular momentum density. Disregarding the inertia correction factor I 0 this term was also identified earlier to transfer energy to the zonal flow [48,43,18].
The second term on the right hand side is a novel term that appears for fluctuating radial velocitŷ u v and the inertia factorÎ 0 . In order to estimate the importance of the inertia factor we plot I 0 for an exemplary tokamak equilibrium in Fig. 2. We immediately see that the inertia factor is not a flux function and is much smaller on the low-field side than on the high-field side. Furthermore, it diverges at the X-point and the O-point. Note that at the X-and O-point also the toroidal component of the E ×B velocity u E,ϕ vanishes since the magnetic field is purely toroidal (and thus the zonal flow energy remains finite. Further, the divergence at the X-point is an integrable singularity as shown in Fig. 3, where we plot the flux-surface average ⟨I 0 ⟩. Here, we mainly see that there appear gradients close to the separatrix and in the core of the domain. Finally, on the right hand side of Eq. (116) we find the source term S zonal . This term contains a contribution from the density source S n proportional to the inertia factor and the square toroidal E ×  B velocity. The sign of this contribution depends only on the sign of S n . Comparing to Fig. 2 we see that the inertia factor is almost 2 orders of magnitude higher on the high field side than on the low field side. A particle source on the tokamak high field side is a far more effective source for zonal flow energy than on the low field side. A second contributor is the angular momentum source defined in Eq. (84), which we already discussed to be pronounced for poloidally asymmetric particle sources.

Simplified magnetic field geometries
It is common in the existing literature to reduce the full three-dimensional magnetic field geometry to simplify expressions. The general magnetic field in Eq. (3) with both toroidal and poloidal components reduces to a purely toroidal magnetic field for ψ p = 0 and the purely poloidal field for ψ t = 0. All our results so far hold for the general magnetic field without axisymmetry. We thus first discuss the poloidal and toroidal fields without assuming axisymmetry. A glance at the gyrokinetic 1-form Eq. (11) convinces us that in each of these cases both the poloidal and toroidal angular momentum have a single component. In a poloidal field the poloidal angular momentum contains only the parallel velocity u ∥ b ϑ while the toroidal angular momentum consists only of the E × B flow u E,ϕ and vice versa for the purely toroidal magnetic field geometry.
For the poloidal field the resulting evolution equations are actually already available and we do not need to compute anything further. The relevant equations are Eq. (79) and the ϑ component of (72). For the purely toroidal magnetic field the parallel momentum balance is given by the ϕ component of (72), however the E × B momentum is problematic since ι is zero and thus Eq. (80) does not hold. Furthermore, since ψ p vanishes the flux-surface average needs to be redefined with the help of ψ t .
In the following we will discuss the axisymmetric case for the general, the purely toroidal and the purely poloidal magnetic fields, which allows further simplifications.

7.1.1.
General axisymmetric magnetic field An axisymmetric magnetic field can be written as in Eq. (6) and is a general feature of the tokamak configuration. It is well known that in this case the toroidal angular momentum density is a conserved quantity [1,2]. In our derivation axisymmetry leads to the full toroidal angular momentum conservation (up to external sources) in the ϕ component of Eq. (102). The ϕ derivatives in the first two terms on the right hand side vanish and the magnetic shear does not contribute. Note that compared to Ref. [1] we cast the time derivative of half the diamagnetic drift on the right hand side of our equation. We comment more on this feature in Sec. 7.2.

Purely toroidal, axisymmetric magnetic field
In the axisymmetric case we discuss here we can write (with cylindrical coordinates R, Z and toroidal angle ϕ).
The gyro-kinetic 1-form Eq. (11) becomes γ = (qψ t dZ + mw ∥ b ϕ dϕ + mµdθ q and now has symmetry in both the R and Z-direction, which makes both qψ t (R) and mw ∥ R conserved quantities separately. This is in fact an important point to emphasize. The purely toroidal magnetic field has two symmetries and thus two exactly conserved quantities instead of just one in the general axisymmetric geometry. In the derivation of the poloidal E × B momentum in Section 4, all we have to do is replace ∇ψ p with ∇ψ t (R) = B(R)ê R , which defines η ∶=ê ϕ × ∇ψ t B(R) ≡ê Z . Equation (80) thus reads (with zero magnetic shear) In the limit B 0 (R) = B 0 and without A 1,∥ and finite Larmor radius effects this equation agrees with [15]. The parallel angular momentum balance Eq. (88) now reduces to Due to the symmetry in R and Z neither the Lorentz force, nor the mirror force appears in Eqs. (122) and (123).

Purely poloidal, axisymmetric magnetic field
The poloidal field approximation with ψ t = 0 is particularly interesting for the field-reversed configuration [52]. We will here investigate the axisymmetric case since, as discussed before, the non-axisymmetric case is already covered. The Poincaré 1-form Eq. (11) reads γ = qψ p dϕ + mw ∥ √ g ϑϑ −1 dϑ + m q µdθ. This results in B * ∥ = B * ⋅ê ϑ = B ⋅ê ϑ ≡ B p withê ϑ ∶= e ϑ e ϑ . The approximation clearly breaks at the X-point where B p = 0, however this point might be redundant since flux coordinates themselves do not exist on the last closed flux-surface where ι −1 diverges as we discussed in Section 2.
It is interesting to note that toroidal symmetry now leads to the exact conservation of γ ϕ = qA ϕ = qψ p sinceê ϑ has no component in dϕ in a symmetric situation. The toroidal angular momentum conservation in the poloidal field approximation thus contains only the toroidal component of the E × B motion. In this case we can write (note that Eq. (4) still holds) η ∶=ê ϑ × ∇ψ p B p = e ϕ which is possible with Eq. (4), B = B pêϑ and e ϑ ⋅ e ϕ = g ϑϕ = 0. The vector η thus points in the actual toroidal direction and does not have a poloidal component. We further have u v E = − dv dψp ∂φ ∂ϕ B p In comparison, we have that u E,ϑ = 0, that is in the poloidal field approximation u E =ê ϑ × ∇φ B p has no poloidal component. The non-zero part of the momentum fluxes is thus where we used that the ϕ component of j f ×B vanishes with ((ê ϑ × ∇ ln B p ) ×ê ϑ ) ϕ = ((ê ϑ × κ) ×ê ϑ ) ϕ = 0 due to the symmetry. This means that in the poloidal field approximation there is no transfer term between E × B motion and parallel momentum just as in the purely toroidal magnetic field in Eq. (122). In contrast the equation for the parallel momentum in toroidally symmetric cases becomes (with b ϕ = 0 In contrast to the purely toroidal magnetic field here we find the mirror force and the geometric correction in the poloidal direction on the right hand side. This means that the background magnetic field acts as a source/sink of parallel momentum.

The momentum of electromagnetic fields in matter
We now note that we have the possibility to rewrite Eq. (61) using identity Eq. (60) to cast the diamagnetic drift under the time derivative (using with Equation (126) is the evolution equation for the electromagnetic momentum flux g ∶= D × B. The electric part in the displacement field D ∶= 0 E + P em vanishes because we neglected the corresponding field part of the action (23) and have quasineutrality. The momentum tensor has the form T with the magnetizing field H ∶= B 1,⊥ µ 0 − M em ⊥ . The momentum flux g and tensor T correspond to the ones given in Reference [53]. With the identification of the Lorentz force density f L = j f × B on the right hand side we can write Eq. (126) as Notice the minus in the Lorentz force, which is a signature that g ϕ is indeed the momentum flux for the electromagnetic field rather than for the plasma itself. Furthermore, the form of the Lorentz force motivates the identification of j f as the free current as opposed to the bound polarization current. The ϑ component of the momentum flux follows by multiplying Eq. (131) with ι −1 Here, notably a contribution from the magnetic shear appears on the right hand side as a coupling term to the external magnetic field. In Eq. (128) we define the electromagnetic polarization charge P em analogous to the magnetization M em ⊥ (129) (which we repeat here for convenience) and different from the gyro-centre polarization charge P gy by half the diamagnetic drift. We remark that neither of these quantities is uniquely defined. The form P em and M em ⊥ highlights the physical origin of polarization and magnetization in gyro-kinetic models. Here, we can view the plasma as a collection of charged discs that can be magnetized and polarized. The disc polarization π ∶= mb ×Ẋ B stems from the drift velocities and reflects that due to the drifts the gyro-orbits are no longer closed [26,29]. Macroscopically, in our model we have P em = mnb × (u E + u D ) B. On the other side, the magnetization M em ⊥ contains the moving electric dipole contribution. An electric dipole π that moves with velocity v ∥b along the magnetic field lines induces a magnetic moment µ = π × v ∥b . However, we only find the diamagnetic part to the moving dipole contribution. We are missing the contribution mnu ∥b × ∇φ B 2 since we neglected the corresponding nonlinear coupling terms in the Hamiltonian (18).
Finally, we note that Eq. (126) can also be viewed as a relation for the radial force density ⟨(P em × B) ϕ ⟩ = ∑ s ⟨ m qB 2 (qn∇ ⊥ φ + ∇ ⊥ p) ⋅ ∇ψ p ⟩ where the force density −qnE ⊥ + ∇ ⊥ p appears inside the bracket on the right hand side. If the right hand side of Eq. (126) is zero, the radial pressure gradient and the radial electric field strength balance each other. Alternatively, we can rewrite Eq. (126) as We emphasize that the evolution equation for the E × B flow Eq. (61), the evolution for the electromagnetic field momentum Eq. (126) and the interpretation as a radial force density or the sum of E × B and diamagnetic drifts in Eq. (133) are completely equivalent views of the same result. In particular, physical arguments made with one of the three equations immediately translate into the other two. The form closest to the drift-fluid vorticity equation [41,54] is Eq. (133). To compare one needs to take the flux-surface average over the generalized vorticity equation and then integrate over the volume. This immediately allows the interpretation of Eq. (126) as the volume integrated equation for a divergence free current or a closed current loop.

Relation to the ion orbit loss mechanism
The ion orbit loss mechanism [13,14,12,55,3,6] refers to the idea that ion orbits close to the Xpoint end on the divertor target or the vessel wall and are thus lost to the confined plasma region. It is thought that the poloidal magnetic field close to the X-point is small such that the grad-B curvature drift velocity dominates over the parallel velocity making ions drift across the separatrix. This then generates a net flux of positive charge out of the confined region. In particle phase space the ions that are on a loss orbit are situated on a "loss-cone" encompassing ions with small parallel velocity and large perpendicular velocity / magnetic moment. It is reported that the perpendicular kinetic energy of the loss cone reaches down to thermal energies [13].
The ion orbit loss is often invoked in models explaining the L-H transition [12,13,14], where it is thought that the outward current leaves a small region inside the separatrix negatively charged, which generates a strong radial electric field. This field in turn drives a strong poloidal shear E × B flow that then forms the transport barrier typical for the high confinement mode. On the other side the same idea is used to explain intrinsic toroidal rotation [55,3,6], the observation that the plasma rotates toroidally without controlled external sources like NBI. The main ingredient here is to assume that the rate by which ions enter loss orbits depends on the direction of their parallel velocity. This then generates an asymmetry between losses of so-called co-and counter-current ions. Since ions carry toroidal momentum, the preferential loss in one direction accelerates the plasma in the other.
Since our derivation of poloidal and toroidal angular momentum balance makes no assumption on the form of the distribution function F (in particular it does not assume that F is Maxwellian) and the particle orbits are retained via Eqs. (31) the ion orbit loss mechanism must consequentially be contained in our results. Here, we want to identify the relevant terms for both poloidal and toroidal rotation.
The net surface integrated current ∫ ψp j ⋅ dA = ⟨j ⋅ ∇v⟩ * flowing through a flux-surface ψ p , in particular the separatrix, by magnetic drifts is given * Recall the definition of the flux-surface average Eq. (7) to see that this is indeed the area integral by s ⟨qn(u ∇B + u κ ) ⋅ ∇v⟩ = ⟨j f ⋅ ∇v⟩ where we inserted the definition of curvature and grad-B drifts Eqs. (A. 19) and (A.18) and the velocity space moments to emphasize the origin of j f as particle drifts. At this point recall again that (4) and (5). The term described in Eq. (134) is nothing but the Lorentz force term that appears in our momentum equations in Section 5 and which we already identified as the Stringer-Winsor spin-up or geodesic transfer term. The ion orbit loss contribution must be contained in the first term on the right side of Eq. (134) since it was argued that ions with large µ and small v ∥ fall on loss orbits. A signature of ion orbit loss would be if the ion term in Eq. (134) is larger than the electron contribution at or close to the separatrix. At this point we notice that for favourable curvature drift direction the curvature vectors counteralign with ∇ψ p (K(ψ p ) < 0, decelerate) on the top and align (K(ψ p ) > 0, accelerate) on the bottom of the tokamak. In order for the flux-surface average in Eq. (134) to yield a non-vanishing result we therefore need an up-down asymmetry of the pressure in the fluxsurface. Furthermore, we notice that for our example tokamak equilibrium in Fig. 1 we have ⟨K κ ⋅ ∇ψ p ⟩ = 0. Indeed, more generally we find ∇ ⋅ K ∇B = −∇ ⋅ K κ = −K κ ⋅ ∇ ln B ∼ O(δ 6 ), which results in ∂ v ⟨K κ ⟩ = ∂ v (⟨K κ ⋅ ∇ψ p ⟩ dv dψ p ) ∼ O(δ 6 ) and thus ⟨K κ ⋅ ∇ψ p ⟩ ≈ 0. This means that only the fluctuations in p ⊥ , p ∥ and nu 2 ∥ contribute and we can write♯ 7.3.1. Poloidal E × B flow Even though, as argued in Eq. (102) in Section 5, the Lorentz force does not generate net poloidal momentum, it does generate E × B momentum, respectively a radial electric field u E,ϑ ∼ ∇φ ⋅ ∇ψ t . We thus conclude that the ion-orbit loss mechanism may indeed contribute to the radial electric field through the Lorentz force.
♯ If we assume p ∥ = p⊥ = p, we can further simplify ⟨j f ⋅ ∇ψp⟩ = ∑ s ⟨b ×∇p B ⋅ ∇ψp⟩ + ⟨m nu 2 ∥ Kκ ⋅ ∇ψp⟩ where we use that K ∇B + Kκ = K and ∇ ⋅ K = 0 (see A1). Then we find the radial component of the diamagnetic drift ⟨u v D ⟩ dψp dv in the first term on the right hand side.
On the other hand, we emphasize that the Lorentz force is not the only candidate that contributes to the poloidal E × B flow generation. Any other term in Eq. (80) could be equally important. Besides the E × B Favre stress we identified for example the diamagnetic stress F v D,η or the density gradient and magnetic shear related terms as additional candidates that may be equally relevant for the L-H transition.

Intrinsic toroidal rotation
The ion loss mechanism is through the Lorentz force indeed contained in the toroidal angular momentum conservation for E × B (79) and parallel (88) angular momentum. However, as we discussed in Eq. (102) in Section 5 the Lorentz force does not actually generate net angular momentum, neither poloidal nor toroidal. A loss of ions through the separatrix does thus not generate toroidal angular momentum. As is shown in Eq. (102) for an axisymmetric equilibrium the only sources for toroidal angular momentum are the actual source terms S nu ∥ and S n on the right hand side. We need additional symmetry breaking mechanisms like for example the asymmetric turbulent diffusion proposed by Reference [3]. Here, a preferential loss of coor counter-current ions through the separatrix generates a net momentum gain for the remaining plasma inside the separatrix. Its signature must be found in the surface contributions of the volume integrated Eq. (102).

Comparison to parallel acceleration
The argument was made [56,57,58,59] that in experimental measurements the parallel velocity u ∥ is measured and not the parallel momentum density nu ∥ . It was concluded that therefore u ∥ respectively ⟨u ∥ ⟩ should be the quantity that theoretical work should focus on when discussing intrinsic rotation. In our view, neither premise nor conclusion of this hypothesis holds. First, the velocity can be measured at the same position and time as the density with for example velocity space tomography [60] (and it should be noted that it is the velocity with respect to the line of sight rather than the parallel velocity that is actually measured in charge exchange diagnostics). Second, u ∥ is not the angular momentum; u ∥ b ϕ ≈ u ∥ R is and only part of it at that. Also, recall that even though it is not technically wrong to compute ⟨u ∥ ⟩ the fluxsurface average (7) is a volume average and should be taken over density like quantities (like nu ∥ ). Finally, what comes out of a gyro-kinetic moment expansion (as performed in [56,57,58,59]) is the gyro-fluid parallel velocity U ∥ , not the actually measured fluid velocity u ∥ . As we discuss in Sec. 3.3 care must be taken when comparing gyro-fluid quantities like U ∥ to the actually physically measured fluid quantity u ∥ due to the involved coordinate transformation of Eq. (49), which for U ∥ is given in Eq. (52). The time evolution equation for u ∥ reads in our ordering (keeping terms up to O(δ 3 )) The terms that appear beside the time derivative are in order the parallel advection term, the E × B advection term, the parallel pressure gradient term, the mirror force term and the last two terms form the parallel electric field. In Eq. (136) we see the local parallel acceleration of a single (ion) species. However, working with accelerations instead of force densities as in Eq. (93) does not reveal that after species summation and flux-surface averaging all net internal forces vanish and only external forces remain. As collectively generated, internal forces neither the pressure gradient nor the electric field can be the source of an intrinsic rotation profile. We point out here that the only external force that is able to make a contribution, the mirror force term ⟨p ⊥ ∇ ∥ ln B⟩, was neglected in [56,57,58,59].

Conclusions
Our main results are the Favre averaged covariant poloidal and toroidal velocity evolution equations (79), (80), and (88) applicable in arbitrary magnetic field geometry including tokamaks, the reversed field pinch, the field-reversed configuration and stellarators. The E ×B equations (79), (80) and the parallel components in Eq. (88) sum to the total angular momentum in Eq. (102). In our full-F gyro-kinetic formulation the perpendicular Favre stress F v ⊥,η appears in the E × B part of the angular momentum as a natural extension of the Reynolds stress through the density weighted fluxsurface average -the Favre average [15]. The perpendicular Favre stress consists of the previously found E × B contribution F v E,η , but also of the novel diamagnetic F v D,η and magnetic flutter F v F,η contributions defined in Eq. (82). Besides the Favre stress, the vacuum Maxwell stress M v B,η and magnetization stresses M v M,η defined in Eq. (83) appear. Furthermore, the Lorentz force originating from the curvature and grad-B drift induced currents represents a source for E × B angular momentum density. Finally, poloidally asymmetric density sources Eq. (84) contribute to angular momentum generation.
Analogous to the E × B part, the parallel component of the angular momentum density Eq. (88) is generated by the parallel Favre stress F v ∥,η in Eq. (91) as well as the kinetic stress K v ∥,η Eq. (90) stemming from magnetic fluctuations. The Lorentz force appears with an opposite sign as in the E × B equation thus vanishing in the summed total angular momentum density in Eq. (102), both toroidally as well as poloidally. In addition, in Eq. (88) the mirror force appears as a source of parallel momentum density.
We construct the inertia tensor from the first fundamental form in Eq. (106).
The relevant discussion is based on the mean flow generated by the covariant, Favre averaged velocity components that we investigate in the first part of the paper. From there we construct the rotational energy in Eq. (111). The E × B part of this energy can be split into a mean "zonal" and fluctuating part and we present the evolution of the mean in Eq. (117) using the previously derived evolution equations for angular momentum. The main finding compared to a simplified geometry is the appearance of a correction factor due to the inertia tensor, which in particular modifies the effect of the density source on the right hand side. A density source on the high field side is a more effective source of zonal flow energy than on the low field side.
We show that we recover previous results obtained in simplified geometries. Interestingly, the purely toroidal magnetic field leads to the exact conservation of both the poloidal E × B velocity as well as the parallel angular momentum density. This is because an additional symmetry is introduced by this geometry. We also point out that the ion orbit loss mechanism as outlined in the literature is identical to the "geodesic transfer term" and the "Stringer-Winsor spin-up mechanism" and is contained in our results in the Lorentz force term on the right hand side of the poloidal E × B angular momentum equation (80). Finally, we clarify several misconceptions in connection with "parallel acceleration" relating previous findings to our results.
The main drawback of our derivation is the longwavelength limit, which effectively reduces our model to a drift-kinetic model and misses higher order finite Larmor radius and polarization effects that could play a role for the L-H transition. We mainly perform this limit in order to avoid tedious geometrical correction factors stemming from for example perpendicular derivatives on the magnetic field unit vector in Section 4 and to avoid the introduction of a fluid closure of the infinite expansions in the polarisation and gyro-averages [30]. On the other side this enables to easily recover the fluid moments in our equations and to compare to existing drift-kinetic models via Eq. (126).
Our results can be used to verify simulation results. The application of these results within full-F gyro-fluid models is subject of ongoing research. However, as previously stated, the available equations in this work are by no means restricted to gyro-fluid models since the derivation contains no assumption on the form of the distribution function. Thus the presented results are relevant also for other frameworks beyond gyro-fluid models like for example gyro-kinetic or drift-fluid models.
The experimental validation of our results may be challenging due to the number of terms that appear in the evolution equations (79), (80), and (88) that in particular require the measurement of plasma potential, parallel velocity, density, pressure and possibly magnetic field fluctuations at the same time and positions. Further, a problematic operation is the flux-surface average. The argument that a time average over the measurement interval equates the flux-surface average only holds if the measured quantity is a flux-function in the first place. On the other hand we provide the theoretical foundation for a discussion of the dominant physical mechanisms that generate poloidal and toroidal angular momentum density and rotational energy in any toroidal magnetic field configuration. Table A1. Definitions of geometric operators with b i the contra-variant components ofb and g ij the contra-variant elements of the metric tensor. We assume (∇ ×b) ∥ = 0.

Name Symbol Definition
Projection Tensor ∥   Table A2. List of the first few gyro-fluid moments: gyro-fluid density N , parallel canonical velocity W ∥ the perpendicular and parallel pressure (P⊥ and P ∥ )/ temperature T⊥ and T ∥ as well the parallel heat flux Q ∥ .
The relation between gyro-fluid quantities N (X, t), U ∥ (X, t), ... given in gyro-centre coordinates X and the physical fluid quantities, which we denote with lower case letters n(x, t) ∶= ∫ d 3 vf (x, v, t), u ∥ (x, t) ∶= ∫ d 3 vv ∥ f (x, v, t) ..., where f (x, v, t) is the distribution function in particle phase-space (and we here overburden the use of v as the velocity instead of the volume flux-label) is given by Eq. (49) This relation can be inverted in the long-wavelength limit up to order δ 2 as for example in Eqs. (51) and (52) Analogous relations hold for the moments of the gyrokinetic source function S N and S N U ∥ .