Variational mean-fluctuation splitting and drift-fluid models

Within the Vlasov kinetic theory of charged particles, fluid models are usually obtained by well-established moment methods. In this process, it is convenient to split the dynamics of the mean flow velocity from the kinetics of velocity fluctuations. More specifically, one introduces the mean-flow velocity ${\boldsymbol{U}}({\bf{x}},t)$ defined by ${\boldsymbol{U}}({\bf{x}},t)\int F({\bf{x}},{\bf{v}},t){{\rm{d}}}^{3}v\,=\int {\bf{v}}F({\bf{x}},{\bf{v}},t){{\rm{d}}}^{3}v$ and rewrites the Vlasov kinetic equation in terms of the following system

$$\begin{eqnarray}&&{mn}\left(\displaystyle \frac{\partial {\boldsymbol{U}}}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{\boldsymbol{U}}\right)=-\mathrm{div}{\mathbb{P}}+{qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {\bf{B}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{l}\displaystyle \frac{{\rm{\partial }}\tilde{F}}{{\rm{\partial }}t}+\left({\bf{c}}+{\boldsymbol{U}}\right)\cdot \displaystyle \frac{{\rm{\partial }}\tilde{F}}{{\rm{\partial }}{\bf{x}}}+\left(\displaystyle \frac{q}{m}{\bf{E}}+\displaystyle \frac{q}{m}({\bf{c}}+{\boldsymbol{U}})\times {\bf{B}}\right.\\ \,\left.-({\bf{c}}+{\boldsymbol{U}})\cdot {\rm{\nabla }}{\boldsymbol{U}}-\displaystyle \frac{{\rm{\partial }}{\boldsymbol{U}}}{{\rm{\partial }}t}\right)\cdot \displaystyle \frac{{\rm{\partial }}\tilde{F}}{{\rm{\partial }}{\bf{c}}}=0,\end{array}\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&n=\int \widetilde{F}({\bf{x}},{\bf{c}},t){{\rm{d}}}^{3}c,\qquad {\mathbb{P}}=m\int {\bf{cc}}\widetilde{F}({\bf{x}},{\bf{c}},t){{\rm{d}}}^{3}c\end{eqnarray} \tag{ 3 }$$

are the zero-th and second kinetic moment, respectively, ${\bf{c}}={\bf{v}}-{\boldsymbol{U}}$ is the fluctuation velocity coordinate, and $\widetilde{F}({\bf{x}},{\bf{c}},t)=F({\bf{x}},{\bf{v}},t)$ is the probability density in the Eulerian frame moving with the mean-flow velocity U. The system (1)–(2) is the most fundamental example of a hybrid kinetic-fluid model, a class of models widely studied in the computational physics of magnetized plasmas [1]. A fluid model is obtained from (1)–(2) by making assumptions on the pressure tensor (density) ${\mathbb{P}}$. Most often, these assumptions make use of a closure, that is an ansatz on the form of the Vlasov distribution depending on the first three moments. The most common class of models is obtained by assuming a zero-mean Gaussian distribution in the relative velocity ${\bf{c}}$, with (usually diagonal) variance matrix given by ${\mathbb{P}}({\bf{x}},t)$ and zero-th velocity moment given by $n({\bf{x}},t)$.

The variational structure of the mean-fluctuation splitting (1)–(2) was exploited in [2] to derive a class of kinetic models for the description of magnetic reconnection in space plasmas. In turn, this variational structure is actually identical to that of hybrid MHD models in the pressure-coupling scheme [3, 4]. The present paper addresses the question of how the mean-fluctuation splitting applies to Littlejohn's kinetic equation for guiding-center motion [5] in order to shed a new light on the formulation of the so-called drift-ordered fluid models, or simply drift-fluids. In particular, while drift-fluid models may be obtained from different procedures, here we shall compare our results with certain classes of gyrofluid models [6–8] obtained within a similar variational approach. As we shall see, important features will emerge along the way, such as the role of the magnetization current and the parallel pressure as a source of magnetic reconnection. In addition, when the mean-fluctuation splitting is performed in the presence of the accompanying Maxwell's equations for the electromagnetic field, we shall see that the intricate structure of guiding-center theory leads to drift-ordered models much more involved than in the case of external electromagnetic fields. In order to simplify the treatment, an alternative guiding-center model is proposed along with its possible applications in hybrid kinetic-fluid models.

As anticipated above, the general variational structure of the mean-fluctuation splitting (1)–(2) is the same as that underlying certain hybrid MHD models for the description of energetic particle effects [3, 4]. This Euler–Poincaré variational structure [9–13] was also exploited in [2] to formulate fully kinetic reconnection models requiring lengthscales much smaller than the electron skin depth. These models were motivated by previously established moment truncations [14] and were recently discussed in [15]. Similar Euler–Poincaré variational structures have also been used recently to approach the kinetic-MHD model [16].

The Euler–Poincaré reduction [10] of Hamilton's variational principle for continuum dynamics is a powerful tool relating the celebrated Lagrangian and Eulerian pictures of fluid flows. In this context, while the Lagrangian paths possess arbitrary variations as in standard classical mechanics, the Eulerian quantities obtained upon reduction by the relabeling symmetry possess constrained variations as in Newcomb's early works [12]. For example, if ${\boldsymbol{\eta }}({{\boldsymbol{x}}}_{0},t)$ is a Lagrangian fluid path with arbitrary variation $\delta {\boldsymbol{\eta }}$, then the Eulerian velocity ${\boldsymbol{v}}({\boldsymbol{x}},t)$ such that $\dot{{\boldsymbol{\eta }}}({{\boldsymbol{x}}}_{0},t)\,={\boldsymbol{v}}({\boldsymbol{x}},t){| }_{{\boldsymbol{x}}={\boldsymbol{\eta }}({{\boldsymbol{x}}}_{0},t)}$ has variations $\delta {\boldsymbol{v}}=\delta (\dot{{\boldsymbol{\eta }}}\circ {{\boldsymbol{\eta }}}^{-1})$, where ◦ denotes composition of functions. A vector calculus exercise shows that $\delta {\boldsymbol{v}}={\partial }_{t}{\boldsymbol{\xi }}+{\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{\xi }}-{\boldsymbol{\xi }}\cdot {\rm{\nabla }}{\boldsymbol{v}}$, where ${\boldsymbol{\xi }}$ is the arbitrary infinitesimal displacement such that $\delta {\boldsymbol{\eta }}({{\boldsymbol{x}}}_{0},t)={\boldsymbol{\xi }}({\boldsymbol{x}},t){| }_{{\boldsymbol{x}}={\boldsymbol{\eta }}({{\boldsymbol{x}}}_{0},t)}$. In the specific case of the mean-fluctuation splitting (1)–(2), the geometric structure of its underlying Euler–Poincaré variational principle is more involved and it exploits the use of the so-called tangent lifts of the fluid path from physical space to phase space, so that Lagrangian fluid paths and particle paths on phase-space can be composed together to realize the change of frame [17] performed in the mean-fluctuation splitting. See [3, 4] for further discussions. Instead of reviewing the geometry underlying the Euler–Poincaré variational principle underlying (1)–(2), here we shall simply write the corresponding Euler–Poincaré equations as they arise from an arbitrary Lagrangian functional of the type $l({\boldsymbol{U}},{\bf{X}},F)$. Notice that we have conveniently dropped the tilde on the distribution function $\widetilde{F}({\bf{x}},{\bf{c}},t)$ from the previous section. Also, ${\bf{X}}({\bf{x}},{\bf{c}},t)=\left({\bf{w}}({\bf{x}},{\bf{c}},t),{\boldsymbol{a}}({\bf{x}},{\bf{c}},t\right)$ is a vector field on phase-space constructed in such a way that, by abusing notation, $\dot{{\bf{x}}}={\bf{w}}+{\boldsymbol{U}}$ (and also $\dot{{\bf{c}}}={\boldsymbol{a}}+{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}}$, by the same abuse of notation). Upon denoting ${\bf{z}}=({\bf{x}},{\bf{c}})$, the general form of the Euler–Poincaré equations is derived from the variational principle [4]

$$\begin{eqnarray*}&&\delta {\int }_{{t}_{1}}^{{t}_{2}}l({\boldsymbol{U}},{\bf{X}},F){\rm{d}}t=0\end{eqnarray*}$$

upon using the following variations

$$\begin{eqnarray*}\begin{array}{rcl} & & \delta {\boldsymbol{U}}={\partial }_{t}{\boldsymbol{\Xi }}+[{\boldsymbol{U}},{\boldsymbol{\Xi }}],\\ & & \delta {\bf{X}}={\partial }_{t}{\bf{Y}}+[{\bf{X}},{\bf{Y}}]+[{{\bf{X}}}_{{\boldsymbol{U}}},{\bf{Y}}]-[{{\bf{X}}}_{{\boldsymbol{\Xi }}},{\bf{X}}],\\ & & \delta F=-{{\rm{\nabla }}}_{{\bf{z}}}\cdot (F\,{\bf{Y}})-{{\rm{\nabla }}}_{{\bf{z}}}\cdot (F\,{{\bf{X}}}_{{\boldsymbol{\Xi }}}).\end{array}\end{eqnarray*}$$

Here, $[{\bf{P}},{\bf{R}}]=({\bf{P}}\cdot {\rm{\nabla }}){\bf{R}}-({\bf{R}}\cdot {\rm{\nabla }}){\bf{P}}$ is the vector field commutator, ${\boldsymbol{\Xi }}({\bf{x}},t)$ and ${\bf{Y}}({\bf{x}},{\bf{c}},t)$ are arbitrary displacement vector fields vanishing at the endpoints t₁ and t₂, while ${{\bf{X}}}_{{\boldsymbol{U}}}=({\boldsymbol{U}},{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}})$ is the phase-space vector field generated by the tangent lift of the fluid velocity vector field ${\boldsymbol{U}}$ onto phase-space (analogously for ${{\bf{X}}}_{{\boldsymbol{\Xi }}}$). Eventually, one is left with the following system in terms of functional derivatives:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{\delta l}{\delta {\boldsymbol{U}}}+{\unicode{x000A3} }_{{\boldsymbol{U}}}\,\displaystyle \frac{\delta l}{\delta {\boldsymbol{U}}}=-\displaystyle \int {\left({\unicode{x000A3} }_{{\bf{X}}}\displaystyle \frac{\delta l}{\delta {\bf{X}}}-F{{\rm{\nabla }}}_{{\bf{z}}}\displaystyle \frac{\delta l}{\delta F}\right)}_{{\bf{x}}}{{\rm{d}}}^{3}c\\ \quad +\,{\rm{\nabla }}\cdot \displaystyle \int {\bf{c}}{\left({\unicode{x000A3} }_{{\bf{X}}}\displaystyle \frac{\delta l}{\delta {\bf{X}}}-F{{\rm{\nabla }}}_{{\bf{z}}}\displaystyle \frac{\delta l}{\delta F}\right)}_{{\bf{c}}}{{\rm{d}}}^{3}c,\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{\delta l}{\delta {\bf{X}}}+{\unicode{x000A3} }_{{\bf{X}}+{{\bf{X}}}_{{\boldsymbol{U}}}}\,\displaystyle \frac{\delta l}{\delta {\bf{X}}}=F\,{{\rm{\nabla }}}_{{\bf{z}}}\displaystyle \frac{\delta l}{\delta F}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial F}{\partial t}+{\rm{\nabla }}\cdot \left[({\bf{X}}+{{\bf{X}}}_{{\boldsymbol{U}}})F\right]=0.\end{eqnarray} \tag{ 6 }$$

Here, ${\unicode{x000A3} }_{{\boldsymbol{W}}}{\boldsymbol{\sigma }}={\rm{\nabla }}\cdot ({\boldsymbol{W}}{\boldsymbol{\sigma }})+{\rm{\nabla }}{\boldsymbol{W}}\cdot {\boldsymbol{\sigma }}$ is the Lie derivative of a one-form density and the subscripts ${\bf{x}}$ and ${\bf{c}}$ denote the space and velocity components of vector-valued quantities on phase-space. As shown in [4], a remarkable property of the system (4)–(6) is that the equation (4) can be rewritten equivalently as a Lie transport equation as follows:

$$\begin{eqnarray}&&\left(\displaystyle \frac{\partial }{\partial t}+{\unicode{x000A3} }_{{\boldsymbol{U}}}\right)\,\left(\displaystyle \frac{\delta l}{\delta {\boldsymbol{U}}}-\int \displaystyle \frac{\delta l}{\delta {\bf{w}}}\,{{\rm{d}}}^{3}c+{{\rm{\nabla }}}_{{\bf{x}}}\cdot \int {\bf{c}}\,\displaystyle \frac{\delta l}{\delta {\boldsymbol{a}}}\,{{\rm{d}}}^{3}c\right)=0.\end{eqnarray} \tag{ 7 }$$

As we shall see, this relation has important consequences as it allows enforcing the zero-average condition on velocity fluctuations.

In this section, we shall show how the mean-fluctuation splitting (1)–(2) arises from equations (4)–(6) by adopting the following Lagrangian, essentially given by the sum of the particle phase-space Lagrangian and the fluid Lagrangian for the mean-flow:

$$\begin{eqnarray}\begin{array}{rcl}l & = & \displaystyle \int F[\left(m{\bf{c}}+m{\boldsymbol{U}}+q{\bf{A}}+{\boldsymbol{\gamma }}\right)\cdot {\bf{w}}\\ & & -\displaystyle \frac{m}{2}| {\bf{c}}{| }^{2}-q{\rm{\Phi }}+\displaystyle \frac{m}{2}| {\boldsymbol{U}}{| }^{2}+q{\boldsymbol{U}}\cdot {\bf{A}}]{{\rm{d}}}^{3}x\,{{\rm{d}}}^{3}c.\end{array}\end{eqnarray} \tag{ 8 }$$

Here, $({\rm{\Phi }},{\bf{A}})$ are the potentials of the external (static) electromagnetic field, while ${\boldsymbol{\gamma }}({\bf{x}},t)$ is an extra dynamical variable playing the role of a Lagrange multiplier enforcing the necessary zero-average condition on velocity fluctuations, that is $\int F{\bf{w}}\,{{\rm{d}}}^{3}{\bf{c}}=0$. The mean-fluctuation splitting system (1)–(2) arises by replacing the functional derivatives

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{\delta l}{\delta {\bf{X}}} & = & F\left(m{\bf{c}}+m{\boldsymbol{U}}+q{\bf{A}}+{\boldsymbol{\gamma }},0\right),\\ \displaystyle \frac{\delta l}{\delta {\boldsymbol{U}}} & = & {mn}{\boldsymbol{U}}+{qn}{\bf{A}},\\ \displaystyle \frac{\delta l}{\delta F} & = & \left(m{\bf{c}}+m{\boldsymbol{U}}+q{\bf{A}}+{\boldsymbol{\gamma }}\right)\cdot {\bf{w}}-\displaystyle \frac{m}{2}| {\bf{c}}{| }^{2}\\ & & -q{\rm{\Phi }}+\displaystyle \frac{m}{2}| {\boldsymbol{U}}{| }^{2}+q{\boldsymbol{U}}\cdot {\bf{A}}\end{array}\end{eqnarray*}$$

in (4)–(6) and by making use of the constraint $\int F{\bf{w}}\,{{\rm{d}}}^{3}c=0$. As a first step, we notice that, the Lagrange multiplier ${\boldsymbol{\gamma }}$ can be set to zero after taking variations, since equation (7) implies

$$\begin{eqnarray*}&&\left(\displaystyle \frac{\partial }{\partial t}+{\unicode{x000A3} }_{{\boldsymbol{U}}}\right)(n{\boldsymbol{\gamma }})=0,\end{eqnarray*}$$

which indeed possesses the trivial solution ${\boldsymbol{\gamma }}({\bf{x}},t)=0$. In what follows we shall always select this particular solution after taking variations. With this in mind, after some vector calculus, equation (5) leads to

$$\begin{eqnarray*}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{lll}{\bf{w}}({\bf{x}},{\bf{c}},t) & = & {\bf{c}},\\ {\boldsymbol{a}}({\bf{x}},{\bf{c}},t) & = & \displaystyle \frac{q}{m}{\bf{E}}+\,\displaystyle \frac{q}{m}({\bf{c}}+{\boldsymbol{U}})\times {\bf{B}}\\ & & -{{\rm{\partial }}}_{t}{\boldsymbol{U}}-(2{\bf{c}}+{\boldsymbol{U}})\cdot {\rm{\nabla }}{\boldsymbol{U}}.\end{array}\end{array}\end{array}\end{array}\end{array}\end{eqnarray*}$$

Then, making use of this result in equation (4) yields (1) and (2).

The advantage of the treatment presented here lies in the fact that, at this point, the Lagrangian (8) may be modified at will, depending on the desired approximation. For example, the assumption of a negligible mean-flow inertia leads to discarding the first two ${\boldsymbol{U}}-$ terms in (8) and this is exactly the approach followed in [2]. In this paper, we shall present a similar application of this approach in the case when the fundamental kinetic equation (before the mean-fluctuation splitting is performed) is given by Littlejohn's drift-kinetic equation of guiding-center motion [5, 18]. This is precisely the topic of the next section and, as we shall see, the resulting system sheds a new light on different features of guiding-center dynamics.

Before we enter the discussion on the role of the the mean-fluctuation splitting in drift-kinetic theory, it is useful to introduce an alternative variational formulation of guiding-center theory, which differs from Littlejohn's well-known Lagrangian

$$\begin{eqnarray}&&L({\bf{x}},\dot{{\bf{x}}},{v}_{\parallel })=({{mv}}_{\parallel }+q{\bf{A}}({\bf{x}}))\cdot \dot{{\bf{x}}}-\displaystyle \frac{m}{2}{v}_{\parallel }^{2}-\mu B({\bf{x}})-q{\rm{\Phi }}({\bf{x}}).\end{eqnarray} \tag{ 9 }$$

Notice that here, we simply treat the magnetic moment μ as a physical constant thereby ignoring the gyrophase. While Littlejohn's Lagrangian (9) is entirely written on the guiding-center phase-space, an alternative Lagrangian is also available on physical space and this reads as follows:

$$\begin{eqnarray}&&{ \mathcal L }({\bf{x}},\dot{{\bf{x}}})=\displaystyle \frac{m}{2}{\left({\boldsymbol{b}}({\bf{x}})\cdot \dot{{\bf{x}}}\right)}^{2}+q\dot{{\bf{x}}}\cdot {\bf{A}}({\bf{x}})-\mu B({\bf{x}})-q{\rm{\Phi }}({\bf{x}}).\end{eqnarray} \tag{ 10 }$$

Then, upon denoting

$$\begin{eqnarray}&&{{\bf{E}}}^{* }=-{\rm{\nabla }}{\rm{\Phi }}-\displaystyle \frac{\mu }{q}{\rm{\nabla }}B,\qquad {{\bf{B}}}^{* }={\bf{B}}+\displaystyle \frac{m}{q}({\boldsymbol{b}}\cdot \dot{{\bf{x}}}){\rm{\nabla }}\times {\boldsymbol{b}},\end{eqnarray} \tag{ 11 }$$

the guiding-center equation emerges as a standard Euler–Lagrange equation on physical space, that is

$$\begin{eqnarray*}&&q{{\bf{E}}}^{* }+q\dot{{\bf{x}}}\times {{\bf{B}}}^{* }=m({\boldsymbol{b}}\cdot \ddot{{\bf{x}}}){\boldsymbol{b}}.\end{eqnarray*}$$

The general idea of carrying the guiding-center position along with the 3-dimensional particle velocity goes back a long time and we address the reader to [19]. The Lagrangian (10) may be helpful to identify conserved quantities associated to certain symmetries. However, in the context of the present paper, the Lagrangian (10) is particularly useful to understand the construction underlying the mean-fluctuation splitting. More specifically, since the mean-flow equation is a fluid equation in physical space, then a description of guiding-center motion in the configuration space becomes advantageous. Indeed, Lagrangian fluid paths are always defined on configuration space, rather than phase-space.

The Lagrangian (10) offers the possibility of lifting guiding-center motion on ${{\mathbb{R}}}^{4}$ to the standard six-dimensional phase-space with corresponding phase-space Lagrangian given by

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{{ \mathcal L }}^{{\rm{{\prime} }}}({\bf{x}},\dot{{\bf{x}}},{\bf{v}})=\left(m{\boldsymbol{b}}({\bf{x}}){\boldsymbol{b}}({\bf{x}})\cdot {\bf{v}}+q{\bf{A}}({\bf{x}}\right))\cdot \dot{{\bf{x}}}\\ \,-\,\displaystyle \frac{m}{2}{\left({\boldsymbol{b}}({\bf{x}})\cdot {\bf{v}}\right)}^{2}-\mu B({\bf{x}})-q{\rm{\Phi }}({\bf{x}}).\end{array}\end{array}\end{eqnarray} \tag{ 12 }$$

This Lagrangian is highly degenerate and obviously no information can be obtained on the perpendicular components of the velocity. Nevertheless, lifting guiding-center motion to ${{\mathbb{R}}}^{6}$ is particularly advantageous as a preliminary step for the formulation of hybrid kinetic-fluid descriptions. As we remarked above, this is due to the fact that fluid and particle paths can be composed with each other, thereby restoring the fundamental duality between Lagrangian and Eulerian pictures in continuum mechanics. As an example, this approach based on lifting guiding-center motion to ${{\mathbb{R}}}^{6}$ was recently exploited in [3] for the formulation of hybrid MHD models in the pressure-coupling scheme. This is precisely the strategy adopted later on in this paper.

In preparation for the sections below, here we shall also present the exact evolution equation for the guiding-center momentum density, which is computed as the first-order moment of Littlejohn's drift-kinetic equation for a static external electromagnetic field. Littlejohn's equation for the guiding-center distribution $f({\bf{x}},{v}_{\parallel },t)$ is considered here in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+{\rm{\nabla }}\cdot (f{\bf{u}})+\displaystyle \frac{\partial f}{\partial {v}_{\parallel }}(f{\alpha }_{\parallel })=0.\end{eqnarray} \tag{ 13 }$$

where $({\bf{u}},{\alpha }_{\parallel })$ is the phase-space vector field defined by

$$\begin{eqnarray}&&q{{\bf{E}}}^{* }+q{\bf{u}}\times {{\bf{B}}}^{* }=m{\alpha }_{\parallel }{\boldsymbol{b}},\end{eqnarray} \tag{ 14 }$$

so that, upon denoting ${B}_{\parallel }^{* }={\boldsymbol{b}}\cdot {{\bf{B}}}^{* }=B+({{mv}}_{\parallel }/q){\boldsymbol{b}}\,\cdot {\rm{\nabla }}\times {\boldsymbol{b}}$, we have

$$\begin{eqnarray}\begin{array}{rcl}{\bf{u}}({\bf{x}},{v}_{\parallel },t) & = & \displaystyle \frac{1}{{B}_{\parallel }^{* }}\left({v}_{\parallel }{{\bf{B}}}^{* }-{\boldsymbol{b}}\times {{\bf{E}}}^{* }\right),\\ {\alpha }_{\parallel }({\bf{x}},{v}_{\parallel },t) & = & \displaystyle \frac{q}{{{mB}}_{\parallel }^{* }}\,{{\bf{B}}}^{* }\cdot {{\bf{E}}}^{* }.\end{array}\end{eqnarray} \tag{ 15 }$$

As we shall see, the algebraic relation (14) may be extremely useful when computing moments of the drift-kinetic equation (13). Indeed, equation (14) is much more fundamental than its algebraic solution (15) for $({\bf{u}},{\alpha }_{\parallel })$.

At this point, we introduce the notation ${\bf{K}}=m\int {{fv}}_{\parallel }{\boldsymbol{b}}{\rm{d}}{v}_{\parallel }{\rm{d}}\mu$ for the momentum density and we are ready to compute its evolution equation. To this purpose, we write

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{m}\displaystyle \frac{\partial {\bf{K}}}{\partial t} & = & -{\rm{\nabla }}\cdot \displaystyle \int {{fv}}_{\parallel }{\bf{u}}{\boldsymbol{b}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu +\displaystyle \int {{fv}}_{\parallel }{\bf{u}}\cdot {\rm{\nabla }}{\boldsymbol{b}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu \\ & & +\displaystyle \int f{\alpha }_{\parallel }{\boldsymbol{b}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu \\ & = & -{\rm{\nabla }}\cdot \displaystyle \int {{fv}}_{\parallel }{\bf{u}}{\boldsymbol{b}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu +\displaystyle \int {{fv}}_{\parallel }{\rm{\nabla }}{\boldsymbol{b}}\cdot {\bf{u}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu \\ & & +\displaystyle \int f\left(q{\bf{u}}\times {\bf{B}}-\mu {\rm{\nabla }}B+q{\bf{E}}\right){\rm{d}}{v}_{\parallel }{\rm{d}}\mu .\end{array}\end{eqnarray} \tag{ 16 }$$

Eventually, by using the relation

$$\begin{eqnarray*}\begin{array}{c}\begin{array}{rcl} & & \displaystyle \int f({v}_{\parallel }{\rm{\nabla }}{\boldsymbol{b}}\cdot {\bf{u}}-\mu {\rm{\nabla }}B){\rm{d}}{v}_{\parallel }{\rm{d}}\mu ={\rm{\nabla }}{\bf{B}}\cdot {\bf{M}}\\ & & \,=\,({\rm{\nabla }}\times {\bf{M}})\times {\bf{B}}+{\rm{\nabla }}\cdot \left[({\bf{B}}\cdot {\bf{M}}){\bf{1}}-{\bf{BM}}\right],\end{array}\end{array}\end{eqnarray*}$$

we are left with the following evolution equation for the guiding-center momentum density:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\bf{K}}}{\partial t}=-{\rm{\nabla }}\cdot {\mathbb{T}}+{qn}{\bf{E}}+\left({\bf{J}}+{\rm{\nabla }}\times {\bf{M}}\right)\times {\bf{B}},\end{eqnarray} \tag{ 17 }$$

where we have defined the stress tensor, the current, and the magnetization respectively as

$$\begin{eqnarray}\begin{array}{rcl}{\mathbb{T}}({\bf{x}},t) & = & \displaystyle \int f\left[{{mv}}_{\parallel }^{2}{\boldsymbol{b}}{\boldsymbol{b}}+{{mv}}_{\parallel }{{\bf{u}}}_{\perp }{\boldsymbol{b}}\right.\\ & & \left.+{{mv}}_{\parallel }{\boldsymbol{b}}{{\bf{u}}}_{\perp }+\mu B({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}})\right]\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu ,\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\bf{M}}({\bf{x}},t)=-\int f\mu \left[{\boldsymbol{b}}-\displaystyle \frac{m}{\mu B}{v}_{\parallel }{{\bf{u}}}_{\perp }\right]\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu ,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\bf{J}}({\bf{x}},t)=q\int f{\bf{u}}\,{\rm{d}}{v}_{\parallel }{\rm{d}}\mu .\end{eqnarray} \tag{ 20 }$$

Here, we have denoted ${{\bf{u}}}_{\perp }=-{\boldsymbol{b}}\times {\boldsymbol{b}}\times {\bf{u}}$. As we can see, the CGL stress tensor in (18) is completed by two nongyrotropic terms, which have already appeared in recent work [20, 21]. Similarly, the definition (19) of the magnetization requires the addition of a nongyrotropic moment appeared in [20, 22].

We notice that the momentum density equation (17) obtained here is exact in the sense that no approximations were used in its derivation. This momentum density equation stands as a fundamental point of departure especially for the formulation of drift-fluid models, which should always arise from (17) upon adopting an adequate moment closure or truncation for the the stress tensor evolution. Indeed, unless a moment closure or truncation is performed, the whole drift-kinetic equation (13)–(14) must be carried over to close the system. Equivalently, in the next section we shall perform the mean-fluctuation splitting on Littlejohn's drift-kinetic equation by applying the variational approach presented earlier in section 1.2.

Starting from the mean-fluctuation system, second-order moment truncations may be obtained upon recalling the form (26)–(27) of the drift-kinetic equation (25) and by enforcing the gyrotropic condition

$$\begin{eqnarray}&&\int \,\tilde{f}\,{c}_{\parallel }{{\boldsymbol{w}}}_{\perp }\,{\rm{d}}{c}_{\parallel }\,{\rm{d}}\mu =0\end{eqnarray} \tag{ 33 }$$

within the moment equations, so that the pressure tensor (29) reduces to the usual CGL form. Also, in this case the magnetization (28) depends uniquely on the magnetic moment density thereby recovering the familiar-looking relation ${\boldsymbol{ \mathcal M }}=-\bar{\mu }{\boldsymbol{b}}$. The next step is the characterization of the magnetic moment density in terms of its equation of motion; it is easy to see that, upon neglecting the magnetic moment flux $\int \tilde{f}{c}_{\parallel }\mu \,{\rm{d}}{c}_{\parallel }\,{\rm{d}}\mu$, the magnetic moment density satisfies the continuity equation

$$\begin{eqnarray}&&{\partial }_{t}\bar{\mu }+{\rm{\nabla }}\cdot (\bar{\mu }{\boldsymbol{U}})=0.\end{eqnarray} \tag{ 34 }$$

Notice that this is the same equation satisfied by the particle number, that is the zero-th moment equation

$$\begin{eqnarray}&&{\partial }_{t}n+{\rm{\nabla }}\cdot (n{\boldsymbol{U}})=0.\end{eqnarray} \tag{ 35 }$$

At this point, one is led to dealing with the longitudinal pressure ${p}_{\parallel }=\,m\int \tilde{f}{c}_{\parallel }^{2}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu$. In the simplest case, the pressure is ignored by the cold-plasma closure, which resorts to the solution ansatz

$$\begin{eqnarray}&&\tilde{f}({\bf{x}},{c}_{\parallel },\mu ,t)=n({\bf{x}},t)\delta ({c}_{\parallel })\delta \left(\mu -\bar{\mu }({\bf{x}},t)/n({\bf{x}},t)\right).\end{eqnarray} \tag{ 36 }$$

Notice that this ansatz represents a fluid that is 'cold' only in the longitudinal direction, while a finite perpendicular temperature is incorporated in the magnetic moment density. For example, this approach was followed in [7] in the context of gyrofluids, although the equations presented therein are different from those obtained here. The cold drift-fluid ansatz (36) produces the following mean-flow equation

$$\begin{eqnarray}&&{mn}{\partial }_{t}{U}_{\parallel }{\boldsymbol{b}}+{mn}({\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}={qn}{\bf{E}}-\bar{\mu }{\rm{\nabla }}B+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }.\end{eqnarray} \tag{ 37 }$$

In [7], the above force $-\bar{\mu }{\rm{\nabla }}B$ associated to the magnetic drift is replaced by a gradient pressure term, thereby making the resulting drift-fluid model similar to those previously proposed in [6]. Also, we notice that the cold drift-fluid equation (37) possesses a variational structure similar to that in [8]. More specifically, equation (37) arises as the Euler–Poincaré equation for the variational principle

$$\begin{eqnarray}&&\delta {\int }_{{t}_{1}}^{{t}_{2}}\int \left(\displaystyle \frac{{mn}}{2}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}+{qn}{\boldsymbol{U}}\cdot {\bf{A}}-\bar{\mu }B-{qn}{\rm{\Phi }}\right){\rm{d}}t=0\end{eqnarray} \tag{ 38 }$$

with the variations $\delta {\boldsymbol{U}}={\partial }_{t}{\boldsymbol{\Xi }}+[{\boldsymbol{U}},{\boldsymbol{\Xi }}]$, $\delta n=-{\rm{\nabla }}\cdot (n{\boldsymbol{\Xi }})$ and $\delta \bar{\mu }=-{\rm{\nabla }}\cdot (\bar{\mu }{\boldsymbol{\Xi }})$. A similar variational structure will be used in later sections. Before proceeding further, it may be useful to remind the reader that the cold-plasma model stands as a moment closure as opposed to a moment truncation. Indeed, equation (37) follows directly from the solution ansatz (36) without any further assumptions. On the contrary, moment truncations are obtained by simply discarding higher-order moment contributions as shown below.

A more sophisticated way of dealing with the longitudinal pressure is by applying the moment method, that is by computing the evolution equation for ${p}_{\parallel }$ from (25) and then truncating to discard heat-flux contributions. To this purpose, we multiply (25) by ${c}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}$ and integrate over the parallel velocity before taking the trace. Upon denoting ${q}_{\parallel }=m\int \tilde{f}{c}_{\parallel }^{3}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu$, this proceeds as follows:

$$\begin{eqnarray*}&&\begin{array}{l}{\boldsymbol{b}}{\boldsymbol{b}}\displaystyle \frac{\partial }{\partial t}\displaystyle \int {c}_{\parallel }^{2}\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu =-{\boldsymbol{b}}{\boldsymbol{b}}{\rm{\nabla }}\cdot \displaystyle \int {c}_{\parallel }^{2}({\boldsymbol{w}}+{\boldsymbol{U}})\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu \\ \ +\,2{\boldsymbol{b}}\displaystyle \int {c}_{\parallel }\left[\displaystyle \frac{q}{m}({\boldsymbol{w}}+{\boldsymbol{U}})\times {{\bf{B}}}^{* }-{\boldsymbol{b}}{\partial }_{t}{U}_{\parallel }\right.\\ \ \left.-\,\displaystyle \frac{\mu }{m}{\rm{\nabla }}B-{c}_{\parallel }{\rm{\nabla }}{U}_{\parallel }-\displaystyle \frac{1}{2}{\rm{\nabla }}{U}_{\parallel }^{2}\right]\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu \\ =\,-{\boldsymbol{b}}{\boldsymbol{b}}{\rm{\nabla }}\cdot \left({p}_{\parallel }{\boldsymbol{U}}+{q}_{\parallel }{\boldsymbol{b}}\right)-{\boldsymbol{b}}{\boldsymbol{b}}{\rm{\nabla }}\cdot \displaystyle \int {c}_{\parallel }^{2}{{\boldsymbol{w}}}_{\perp }\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu \\ \ +\,2{\boldsymbol{b}}\displaystyle \int {c}_{\parallel }\left[{{\boldsymbol{w}}}_{\perp }\times \left(q{\bf{B}}+{\rm{\nabla }}\times ({c}_{\parallel }{\boldsymbol{b}}+{U}_{\parallel }{\boldsymbol{b}}\right)\right.\\ \ \left.-\,\displaystyle \frac{\mu }{m}{\rm{\nabla }}B\right]\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu \\ \ +\,\displaystyle \frac{2}{m}{\boldsymbol{b}}\left[{q}_{\parallel }{\boldsymbol{b}}\times {\rm{\nabla }}\times {\boldsymbol{b}}+{p}_{\parallel }\left({\boldsymbol{U}}\times {\rm{\nabla }}\times {\boldsymbol{b}}+{\boldsymbol{b}}\right.\right.\\ \ \left.\left.\times \,{\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}})-{\rm{\nabla }}{U}_{\parallel }\right)\right],\end{array}\end{eqnarray*}$$

where the first equality follows from writing (25) as in (26)–(27), while the second equality follows after writing ${\boldsymbol{w}}={{\boldsymbol{w}}}_{\perp }+{c}_{\parallel }{\boldsymbol{b}}$ and expanding. Then, taking the trace and invoking (33) leads to

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {p}_{\parallel }}{\partial t} & = & -{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{U}}+{q}_{\parallel }{\boldsymbol{b}}+{{\boldsymbol{q}}}_{\perp })\\ & & -2{\boldsymbol{b}}{\boldsymbol{b}}:({p}_{\parallel }{\rm{\nabla }}{\boldsymbol{U}}+{\rm{\nabla }}{{\boldsymbol{q}}}_{\perp })+2{q}_{\mu }{\rm{\nabla }}\cdot {\boldsymbol{b}},\end{array}\end{eqnarray} \tag{ 39 }$$

where ${{\boldsymbol{q}}}_{\perp }=m\int \tilde{f}{c}_{\parallel }^{2}{{\bf{w}}}_{\perp }{\rm{d}}{c}_{\parallel }{\rm{d}}\mu$ and ${q}_{\mu }=B\int \tilde{f}{c}_{\parallel }\mu \,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu$. Here, the special case ${{\boldsymbol{q}}}_{\perp }=0$ consistently recovers equation (16) in [23]. In turn, we notice that equation (39) differs quite substantially from the celebrated parallel CGL equation [24, 25]. Indeed, while (39) has been obtained here exclusively from guiding-center theory, the original CGL derivation [24, 25] relies on special solutions of the truncated equation for the full-orbit pressure tensor dynamics. Nevertheless, the parallel CGL equation coincides with (39) when heat flux contributions are neglected, so that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}_{\parallel }}{\partial t}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{p}_{\parallel }=-{p}_{\parallel }{\rm{\nabla }}\cdot {\boldsymbol{U}}-2{p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}:{\rm{\nabla }}{\boldsymbol{U}}.\end{eqnarray} \tag{ 40 }$$

In the context of MHD models, due to the frozen-in condition on the magnetic field, equation (40) also implies an adiabatic invariant and it can be obtained from a thermodynamic equation of state [26, 27]. As pointed out in [28], the frozen-in condition on the magnetic field is strictly necessary for the emergence of an adiabatic invariant associated to the parallel pressure. Here, we notice that (40) can be written in terms of the Lie derivative of the tensor density ${p}_{\parallel }{b}_{i}{b}_{j}\,{\rm{d}}{x}^{i}{\rm{d}}{x}^{j}\otimes {{\rm{d}}}^{3}x$ as follows:

$$\begin{eqnarray*}&&\mathrm{Tr}\left[\left(\displaystyle \frac{\partial }{\partial t}+{\unicode{x000A3} }_{{\boldsymbol{U}}}\right)({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\right]=0,\end{eqnarray*}$$

where $\mathrm{Tr}$ denotes the matrix trace and in this case the Lie derivative reads ${({\unicode{x000A3} }_{{\boldsymbol{U}}}S)}_{{ij}}={\boldsymbol{U}}\cdot {\rm{\nabla }}{S}_{{ij}}+({\rm{\nabla }}\cdot {\boldsymbol{U}}){S}_{{ij}}\,+{S}_{{kj}}{\partial }_{i}{U}^{k}+{S}_{{ki}}{\partial }_{j}{U}^{k}$. Equation (40) accompanies the drift-fluid equation, which extends the cold-fluid model (37) as follows:

$$\begin{eqnarray}&&\begin{array}{l}{mn}{\partial }_{t}{U}_{\parallel }{\boldsymbol{b}}+{mn}({\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}\\ \quad =\,-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})+{qn}{{\boldsymbol{E}}}^{* }+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }\\ \quad =\,-{\rm{\nabla }}\cdot [{p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}+\bar{\mu }B({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}})]\\ \qquad +\,{qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }-{\bf{B}}\times {\rm{\nabla }}\times {\boldsymbol{M}},\end{array}\end{eqnarray} \tag{ 41 }$$

where ${\boldsymbol{M}}=-\bar{\mu }{\boldsymbol{b}}$. The moment model (41)–(40) emerges here as a slight extension of a model previously appeared in [8] and lacking the magnetization term. In turn, this magnetization term appears to be necessary for momentum conservation. The dynamics of the total drift-fluid momentum is actually the topic of the next section.

A third way to deal with pressure in drift-fluid models was first proposed by Pfirsch and Correa-Restrepo (PCR) [6] and it consists in replacing the kinetic energy term ${{mc}}_{\parallel }^{2}/2+\mu B$ in (21) by a fluid internal energy. The resulting fluid model can be obtained from the variational principle (38) upon subtracting the internal energy term ${mn}{ \mathcal U }(n)$. Here we shall follow Pfirsch's approach by adopting a barotropic fluid closure, while adiabatic extensions can be easily formulated to include entropy transport. Then, the following variant of PCR variational principle

$$\begin{eqnarray}&&\delta {\int }_{{t}_{1}}^{{t}_{2}}\int \left(\displaystyle \frac{{mn}}{2}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}+{qn}{\boldsymbol{U}}\cdot {\bf{A}}-{qn}{\rm{\Phi }}-{mn}U(n)\right){\rm{d}}t=0\end{eqnarray} \tag{ 42 }$$

produces the fluid equation

$$\begin{eqnarray}&&{mn}{\partial }_{t}{U}_{\parallel }{\boldsymbol{b}}+{mn}({\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}={qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }-{\rm{\nabla }}p,\end{eqnarray} \tag{ 43 }$$

where $p={mn}{ \mathcal U }^{\prime} (n)$. While Pfirsch and Correa-Restrepo obtained an analogous result by including the polarization drift in the variational principle (42), equation (43) coincides with PCR when those polarization effects are ignored (see equation (44) in [6]). Also, notice that the magnetic drift arising from the presence of the magnetic moment density in (37) is actually absent in (43). Instead, the magnetic drift has been replaced by the diamagnetic drift associated to the isotropic pressure; see section 7.2 in [19] for a discussion on how the coexistence of magnetic and diamagnetic drift is actually excluded by physical arguments.

In this section, we want to study the implications of the drift-fluid momentum equation in different cases. We begin the discussion by considering equation (41) in the equivalent Lie-derivative form

$$\begin{eqnarray*}&&\begin{array}{l}\left(\displaystyle \frac{\partial }{\partial t}+{\unicode{x000A3} }_{{\boldsymbol{U}}}\right)({qn}{\bf{A}}+{{mnU}}_{\parallel }{\boldsymbol{b}})=-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\\ \quad -\,\bar{\mu }{\rm{\nabla }}B+\displaystyle \frac{{mn}}{2}{\rm{\nabla }}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}+{qn}{\rm{\nabla }}({\boldsymbol{U}}\cdot {\bf{A}}),\end{array}\end{eqnarray*}$$

thereby leading to the following explicit equation for circulation dynamics:

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{c}\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\oint }_{{\rm{\Gamma }}(t)}\left[\displaystyle \frac{n}{\bar{\mu }}(q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\right]\cdot {\rm{d}}{\bf{x}}={\oint }_{{\rm{\Gamma }}(t)}\displaystyle \frac{1}{\bar{\mu }}\left[\displaystyle \frac{{mn}}{2}{\rm{\nabla }}{U}_{\parallel }^{2}\right.\\ \,\left.+\,{qn}{\rm{\nabla }}({\boldsymbol{U}}\cdot {\bf{A}})-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\right]\cdot {\rm{d}}{\bf{x}},\end{array}\end{array}\end{array}\end{eqnarray*}$$

for an arbitrary loop ${\rm{\Gamma }}(t)$ moving with velocity ${\boldsymbol{U}}$. Since ${\rm{\Gamma }}(t)$ is arbitrary, the above circulation equation leads to the vorticity evolution by applying Stokes theorem. Thus, the terms generating circulation are also terms generating vorticity.

Let us first consider the case of a cold drift-fluid, that is ${p}_{\parallel }\,=\,0$. Unlike the full-orbit version of the cold-fluid equation, we notice that circulation is generated by the gradients of the ratio $n/\bar{\mu }$, which then represent a source of magnetic reconnection that is instead absent in standard full-orbit fluid models. If we insist on full consistency with the full-orbit case, then reconnection sources can be eliminated by restricting to consider a constant ratio $n/\bar{\mu }$, which indeed is allowed by the equation ${\partial }_{t}(n/\bar{\mu })+{\boldsymbol{U}}\cdot {\rm{\nabla }}(n/\bar{\mu })=0$. Notice that this is in agreement with the barotropic (isentropic) limit of equation (47) in [26], where the authors propose defining a characteristic perpendicular length such that ${d}_{\perp }=\bar{\mu }/n={const}.$

Now let us consider the case in which $\bar{\mu }/n={d}_{\perp }$ and the longitudinal pressure is retained, so that the circulation law becomes

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\oint }_{{\rm{\Gamma }}(t)}(q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\cdot {\rm{d}}{\bf{x}}=-{d}_{\perp }{\oint }_{{\rm{\Gamma }}(t)}\left[\displaystyle \frac{1}{\bar{\mu }}{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\right]\cdot {\rm{d}}{\bf{x}}.\end{eqnarray*}$$

As we can see, circulation is now generated by the longitudinal pressure. For the sake of completeness, we should remark that in this case (41) produces the following equation of motion for the hydrodynamic helicity:

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\displaystyle \int (q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\cdot (q{\bf{B}}+m{\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}})){{\rm{d}}}^{3}x\\ =\,-2{d}_{\perp }\displaystyle \int (q{\bf{B}}+m{\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}}))\cdot \left[\displaystyle \frac{1}{\bar{\mu }}{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\right]{{\rm{d}}}^{3}x.\end{array}\end{eqnarray*}$$

Thus, the longitudinal pressure stands as another source of vorticity. We notice that this is very different from the case of isentropic full-orbit fluid models, whose pressure forces do not generate vorticity. Now, since the circulation generated by the longitudinal pressure also represents a source of magnetic reconnection, one has to be careful about the physical effects incorporated in the model. For example, the Hall-MHD model does not comprise magnetic reconnection in the non-resistive case.

A possible way to formulate drift-fluid models whose pressure forces mimic the full-orbit case is found in the PCR approach [6]. Indeed, in this case, equation (43) preserves the circulation integral ${\oint }_{{\rm{\Gamma }}}(q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\cdot {\rm{d}}{\bf{x}}$, thereby leading to helicity conservation

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\displaystyle \int (q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\cdot (q{\bf{B}}+m{\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}})){{\rm{d}}}^{3}x\\ \quad =\,{qm}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\displaystyle \int (B+{B}_{\parallel }^{* }){U}_{\parallel }\,{{\rm{d}}}^{3}x=0\end{array}\end{eqnarray*}$$

and reproducing the analogous result from full-orbit models. Thus, in this sense, the PCR drift-fluid method leads to drift-fluid models that successfully reproduce guiding-center drifts while retaining the invariants from the full-orbit theory.

As pointed out in section 3.1, the existence of a longitudinal adiabatic invariant ${\partial }_{t}({p}_{\parallel }{B}^{2}/{n}^{3})+{\boldsymbol{U}}\cdot {\rm{\nabla }}({p}_{\parallel }{B}^{2}/{n}^{3})=0$ arises in CGL fluid theory [24] since in that particular case the magnetic field is frozen-in. In the general case, no longitudinal adiabatic invariant can be shown to emerge from standard guiding-center theory. On the other hand, one may be tempted to enforce the existence of such an invariant by invoking a CGL equation of state of the type

$$\begin{eqnarray}&&{p}_{\parallel }\,=\,\displaystyle \frac{{m}^{3}{n}^{3}}{{B}^{2}}h(s),\end{eqnarray} \tag{ 44 }$$

where h is some fixed positive function and $s({\bf{x}},t)$ is a transported scalar satisfying ${\partial }_{t}s+{\boldsymbol{U}}\cdot {\rm{\nabla }}s=0$. For example, this scalar can be the specific entropy [27], although here we shall chose the specific magnetic moment $s=\bar{\mu }/n$. Then, since ${p}_{\parallel }\,=\,m\int \tilde{f}{c}_{\parallel }^{2}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu$, the guiding-center kinetic energy associated to velocity fluctuations can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int \tilde{f}\left(\displaystyle \frac{m}{2}{c}_{\parallel }^{2}+\mu B\right){\rm{d}}{c}_{\parallel }\,{\rm{d}}\mu ={mn}\,{ \mathcal U }(n,s)\\ \quad =\,{mn}\left(\displaystyle \frac{{m}^{2}{n}^{2}}{2{B}^{2}}h(s)+\displaystyle \frac{s}{m}B\right),\end{array}\end{eqnarray} \tag{ 45 }$$

thereby leading to the following anisotropic modification of the PCR variational principle (42):

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{c}\delta {\int }_{{t}_{1}}^{{t}_{2}}\int \left(\displaystyle \frac{{mn}}{2}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}+{qn}{\boldsymbol{U}}\cdot {\bf{A}}\right.\\ \left.-{qn}{\rm{\Phi }}-{mn}{ \mathcal U }(n,s)\right){\rm{d}}t=0.\end{array}\end{array}\end{eqnarray} \tag{ 46 }$$

Here, the variations are again given by $\delta {\boldsymbol{U}}={\partial }_{t}{\boldsymbol{\Xi }}+[{\boldsymbol{U}},{\boldsymbol{\Xi }}]$, $\delta n=-{\rm{\nabla }}\cdot (n{\boldsymbol{\Xi }})$, and $\delta \bar{\mu }=-{\rm{\nabla }}\cdot (\bar{\mu }{\boldsymbol{\Xi }})$. Thus, we compute

$$\begin{eqnarray*}&&\begin{array}{c}\begin{array}{l}\delta {p}_{\parallel }\,=\,\delta \left(\displaystyle \frac{{m}^{3}{n}^{3}}{{B}^{2}}h(s)\right)=-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{\Xi }})\\ \,-\,2{p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}:{\rm{\nabla }}{\boldsymbol{\Xi }}+2\displaystyle \frac{{p}_{\parallel }}{B}{\boldsymbol{b}}\cdot {\rm{\nabla }}\times ({\boldsymbol{\Xi }}\times {\bf{B}}),\end{array}\end{array}\end{eqnarray*}$$

where the last term appears by writing $0=\delta {\bf{B}}={\rm{\nabla }}\,\times ({\boldsymbol{\Xi }}\,\times {\bf{B}})-{\rm{\nabla }}\times ({\boldsymbol{\Xi }}\times {\bf{B}})$, since we are dealing with an external electromagnetic field. The mean flow equation reads as

$$\begin{eqnarray}&&\begin{array}{l}{mn}{\partial }_{t}{U}_{\parallel }{\boldsymbol{b}}+{mn}({\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}={qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }\\ -\,{\rm{\nabla }}\cdot [{p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}+\bar{\mu }B({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}})]+{\bf{B}}\times {\rm{\nabla }}\times ((\bar{\mu }+{p}_{\parallel }/B){\boldsymbol{b}}).\end{array}\end{eqnarray} \tag{ 47 }$$

Thus, the last term shows that enforcing the longitudinal adiabatic invariant in guiding-center motion produces an extra magnetization current, since the magnetization density is now redefined to include parallel pressure contributions, i.e. ${\boldsymbol{M}}=-({p}_{\parallel }/B+\bar{\mu }){\boldsymbol{b}}$. In addition, we observe that if $\bar{\mu }={d}_{\perp }n$ as in section 3.2, then the drift-fluid model (35)–(47) preserves Kelvin's circulation ${\oint }_{{\rm{\Gamma }}(t)}(q{\bf{A}}+{{mU}}_{\parallel }{\boldsymbol{b}})\cdot {\rm{d}}{\bf{x}}$ as well as the hydrodynamic helicity $\int (B+{B}_{\parallel }^{* }){U}_{\parallel }\,{{\rm{d}}}^{3}x$. Thus, while invoking extra magnetization effects, the CGL equation of state recovers consistency with full-orbit models in terms of circulation conservation and so it does not comprise possible sources of magnetic reconnection. The same arguments hold for possible extensions of the CGL equation of state that have previously been proposed in [29]. Therein, the parallel pressure ${p}_{\parallel }$ is expressed as a more general function ${{ \mathcal F }}_{\parallel }$ of the ratio ${n}^{3}/{B}^{2}$, so that ${p}_{\parallel }\,=\,{{ \mathcal F }}_{\parallel }({n}^{3}/{B}^{2})$. Analogously, one defines two functions ${{ \mathcal F }}_{\perp }^{(1)}$ and ${{ \mathcal F }}_{\perp }^{(2)}$ so that the magnetic moment density is expressed as $\bar{\mu }=n{{ \mathcal F }}_{\perp }^{(1)}({n}^{3}/{B}^{2})\,+{{nB}}^{-1}{{ \mathcal F }}_{\perp }^{(2)}({n}^{3}/{B}^{2})$. Then, one can still define a (barotropic) fluid internal energy ${ \mathcal U }(n)$ by an appropriate modification of (45), thereby preserving both circulation and hydrodynamic helicity.

In the first model which we consider, the fluctuation terms in the Lagrangian (21) are replaced by the full-orbit fluctuation terms appearing in (8). This step produces the following modified Lagrangian

$$\begin{eqnarray}\begin{array}{rcl}l & = & \displaystyle \int F\left[\left(m{\bf{c}}+m{\boldsymbol{b}}{\boldsymbol{b}}\cdot {\boldsymbol{U}}+q{\bf{A}}+{\boldsymbol{\gamma }}\right)\cdot {\bf{w}}\right.\\ & & \left.-\displaystyle \frac{m}{2}{c}^{2}-q{\rm{\Phi }}+\displaystyle \frac{m}{2}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}+q{\boldsymbol{U}}\cdot {\bf{A}}\right]{{\rm{d}}}^{3}x\,{{\rm{d}}}^{3}c,\end{array}\end{eqnarray} \tag{ 48 }$$

where the distribution function is $F=F({\bf{x}},{\bf{c}},t)$. The resulting equations of motion are again obtained by specializing (4)–(6). In this case, we obtain ${\bf{w}}={\bf{c}}$ and

$$\begin{eqnarray*}&&\begin{array}{l}m({\boldsymbol{a}}+{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}})-({\bf{c}}+{\boldsymbol{U}})\times \left[q{\bf{B}}+m{\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}})\right]\\ \quad =\,-m{\boldsymbol{b}}{\partial }_{t}{U}_{\parallel }+q{\bf{E}}-m{\rm{\nabla }}{\boldsymbol{U}}\cdot {\bf{c}}-\displaystyle \frac{m}{2}{\rm{\nabla }}{U}_{\parallel }^{2}.\end{array}\end{eqnarray*}$$

Then, upon denoting

$$\begin{eqnarray*}&&{{\bf{E}}}^{\star }={\bf{E}}-(m/q){\rm{\nabla }}{\boldsymbol{U}}\cdot {\bf{c}},\qquad {{\bf{B}}}^{\star }={\bf{B}}+(m/q){\rm{\nabla }}\times ({U}_{\parallel }{\boldsymbol{b}}),\end{eqnarray*}$$

and after noticing that ${{\rm{\nabla }}}_{{\bf{z}}}\cdot ({\bf{X}}+{{\bf{X}}}_{{\boldsymbol{U}}})=0$ (see the notation in section 1.2), the kinetic equation reads

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{\partial F}{\partial t}+({\bf{c}}+{\boldsymbol{U}})\cdot \displaystyle \frac{\partial F}{\partial {\bf{x}}}+\left[(q/m)({\bf{c}}+{\boldsymbol{U}})\times {{\bf{B}}}^{\star }\right.\\ \quad \left.+\,(q/m){{\bf{E}}}^{\star }-{\boldsymbol{b}}{\partial }_{t}{U}_{\parallel }-{U}_{\parallel }{\rm{\nabla }}{U}_{\parallel }\right]\cdot \displaystyle \frac{\partial F}{\partial {\bf{c}}}=0.\end{array}\end{eqnarray*}$$

In addition, since ${\left({\unicode{x000A3} }_{{\bf{X}}}(\delta l/\delta {\bf{X}})-F{{\rm{\nabla }}}_{{\bf{z}}}(\delta l/\delta F)\right)}_{{\bf{c}}}=0$ and

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \int {\left({\unicode{x000A3} }_{{\bf{X}}}\displaystyle \frac{\delta l}{\delta {\bf{X}}}-F{{\rm{\nabla }}}_{{\bf{z}}}\displaystyle \frac{\delta l}{\delta F}\right)}_{{\bf{x}}}{{\rm{d}}}^{3}c\\ \quad =\,-{\rm{\nabla }}\cdot {\mathbb{P}}-\displaystyle \frac{{mn}}{2}{\rm{\nabla }}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2}-{qn}{\rm{\nabla }}({\boldsymbol{U}}\cdot {\bf{A}})-{qn}{\bf{E}},\end{array}\end{eqnarray*}$$

the fluid equation (4) reads

$$\begin{eqnarray}&&{mn}{\partial }_{t}{U}_{\parallel }{\boldsymbol{b}}+{mn}({\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}=-{\rm{\nabla }}\cdot {\mathbb{P}}+{qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* },\end{eqnarray} \tag{ 49 }$$

with ${{\boldsymbol{B}}}^{* }={\bf{B}}+(m/q){U}_{\parallel }{\rm{\nabla }}\times {\boldsymbol{b}}$ and we recall ${\mathbb{P}}\,=\int F{\bf{cc}}\,{{\rm{d}}}^{3}c$. We notice that, upon introducing the hat map notation ${\widehat{A}}_{{ij}}=-{\epsilon }_{{ijk}}{A}_{k}$, the pressure tensor evolution reads

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\mathbb{P}}}{\partial t}+({\boldsymbol{U}}\cdot {\rm{\nabla }}){\mathbb{P}}+({\rm{\nabla }}\cdot {\boldsymbol{U}}){\mathbb{P}}+{\mathbb{P}}\cdot {\rm{\nabla }}{\boldsymbol{U}}\\ \quad +\,{\left({\mathbb{P}}\cdot {\rm{\nabla }}{\boldsymbol{U}}\right)}^{T}+\left[(q/m){\widehat{B}}^{* }-\widehat{\omega },{\mathbb{P}}\right]=-{\rm{\nabla }}\cdot {\bf{Q}},\end{array}\end{eqnarray} \tag{ 50 }$$

where $[\cdot ,\cdot ]$ denotes the matrix commutator, ${\boldsymbol{\omega }}={\rm{\nabla }}\times {\boldsymbol{U}}$ is the mean flow vorticity, and ${\bf{Q}}=\int F{\bf{ccc}}\,{{\rm{d}}}^{3}c$ is the heat flux tensor. It is easy to see that this equation possesses gyrotropic solutions recovering the CGL evolution equations [24, 25] and thus this approach recovers the gyrofluid model in [8]. In addition, at this point one can adopt an isotropic equation of state as in the standard approach to fluid dynamics; this leads naturally to equation (43) as in the PCR model [6] discussed previously. Alternatively, retaining non-gyrotropic components of the full pressure tensor can be advantageous, for example, in capturing magnetic reconnection effects while still exploiting the simplifications offered by the use of the guiding-center approximation to describe the mean-fluid flow fluid.

In this section, we shall consider the alternative approach to that described above. That is, here we shall treat the mean fluid flow as a standard full-orbit charged fluid, while we shall retain the guiding-center description to describe fluctuation dynamics. This leads immediately to the following Lagrangian:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{rcl}l & = & \displaystyle \int F\left[\left(m{\boldsymbol{b}}{\boldsymbol{b}}\cdot {\bf{c}}+m{\boldsymbol{U}}+q{\bf{A}}+{\boldsymbol{\gamma }}\right)\cdot {\bf{w}}\right.\\ & & -\displaystyle \frac{m}{2}{\left({\boldsymbol{b}}\cdot {\bf{c}}\right)}^{2}-\mu B-q{\rm{\Phi }}\\ & & \left.+\displaystyle \frac{m}{2}{U}^{2}+q{\boldsymbol{U}}\cdot {\bf{A}}\right]{{\rm{d}}}^{3}x\,{{\rm{d}}}^{3}c\,{\rm{d}}\mu .\end{array}\end{array}\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 51 }$$

Again, the equations of motion are obtained by specializing (4)–(6). In this case, we obtain ${\boldsymbol{b}}\cdot {\bf{w}}={\boldsymbol{b}}\cdot {\bf{c}}$ as well as, after a vector calculus exercise,

$$\begin{eqnarray*}&&\begin{array}{l}m\left[{\boldsymbol{b}}\cdot ({\boldsymbol{a}}+{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}})+({\bf{w}}+{\boldsymbol{U}})\cdot {\rm{\nabla }}{\boldsymbol{b}}\cdot {\bf{c}}\right]{\boldsymbol{b}}\\ \qquad -\,q({\bf{w}}+{\boldsymbol{U}})\times {{\boldsymbol{ \mathcal B }}}^{\star }\\ \quad =\,-m{\partial }_{t}{\boldsymbol{U}}+q{{\bf{E}}}^{* }-m({\boldsymbol{b}}\cdot {\bf{c}}){\rm{\nabla }}({\boldsymbol{b}}\cdot {\boldsymbol{U}})-\displaystyle \frac{m}{2}{\rm{\nabla }}{U}^{2},\end{array}\end{eqnarray*}$$

where we have used

$$\begin{eqnarray*}&&{{\boldsymbol{ \mathcal B }}}^{\star }={\bf{B}}+(m/q)({\boldsymbol{b}}\cdot {\bf{c}}){\rm{\nabla }}\times {\boldsymbol{b}}+(m/q){\rm{\nabla }}\times {\boldsymbol{U}}.\end{eqnarray*}$$

Then, upon using (24), we have the kinetic equation

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rm{\partial }}\tilde{f}}{{\rm{\partial }}t}+{\rm{\nabla }}\cdot \left[\left({c}_{\parallel }+{U}_{\parallel })\displaystyle \frac{\tilde{f}{{\boldsymbol{ \mathcal B }}}^{\star }}{{{ \mathcal B }}_{\parallel }^{\star }}-\displaystyle \frac{m}{q}\displaystyle \frac{\tilde{f}{\boldsymbol{b}}}{{{ \mathcal B }}_{\parallel }^{\star }}\right.\right.\\ \,\left.\times \,\left(\displaystyle \frac{q{{\bf{E}}}^{\ast }}{m}-{c}_{\parallel }{\rm{\nabla }}{U}_{\parallel }-\frac{{\rm{\nabla }}{U}^{2}}{2}-\displaystyle \frac{{\rm{\partial }}{\boldsymbol{U}}}{{\rm{\partial }}t}\right)\right]\\ \,+\,\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}{c}_{\parallel }}\left[\displaystyle \frac{\tilde{f}{{\boldsymbol{ \mathcal B }}}^{\star }}{{{ \mathcal B }}_{\parallel }^{\star }}\cdot \left(\displaystyle \frac{q{{\bf{E}}}^{\ast }}{m}-{c}_{\parallel }{\rm{\nabla }}{U}_{\parallel }-\frac{{\rm{\nabla }}{U}^{2}}{2}-\displaystyle \frac{{\rm{\partial }}{\boldsymbol{U}}}{{\rm{\partial }}t}\right)\right]=0.\end{array}\end{eqnarray} \tag{ 52 }$$

This can be equivalently written in the form (26), where the vector field $({\boldsymbol{w}},{a}_{\parallel })$ is now given by the solution of the algebraic equation

$$\begin{eqnarray}&&\begin{array}{l}{{ma}}_{\parallel }{\boldsymbol{b}}-q({\boldsymbol{w}}+{\boldsymbol{U}})\times {{\boldsymbol{ \mathcal B }}}^{\star }\\ \quad =\,-m{\partial }_{t}{\boldsymbol{U}}+q{{\bf{E}}}^{* }-{{mc}}_{\parallel }{\rm{\nabla }}{U}_{\parallel }-\displaystyle \frac{m}{2}{\rm{\nabla }}{U}^{2}.\end{array}\end{eqnarray} \tag{ 53 }$$

Upon taking the second moment and by neglecting the heat flux contributions, it is easy to see that this formulation of the kinetic equation (52) recovers the CGL equation (40) for the parallel pressure ${p}_{\parallel }$.

In order to write the accompanying fluid equation, we first verify that ${\left({\unicode{x000A3} }_{{\bf{X}}}(\delta l/\delta {\bf{X}},)\right)}_{{\bf{c}}}=F{{\rm{\nabla }}}_{{\bf{c}}}(\delta l/\delta F)$ and then we compute

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \int {\left({\unicode{x000A3} }_{{\bf{X}}}\displaystyle \frac{\delta l}{\delta {\bf{X}}}-F{{\rm{\nabla }}}_{{\bf{z}}}\displaystyle \frac{\delta l}{\delta F}\right)}_{{\bf{x}}}{{\rm{d}}}^{3}c\,{\rm{d}}\mu ={\rm{\nabla }}\cdot { \mathcal P }\\ \quad +\,{\bf{B}}\times {\rm{\nabla }}\times {\boldsymbol{ \mathcal M }}-\displaystyle \frac{{mn}}{2}{\rm{\nabla }}{U}^{2}-{qn}{\rm{\nabla }}({\boldsymbol{U}}\cdot {\bf{A}})-{qn}{\bf{E}},\end{array}\end{eqnarray*}$$

with the same notation as in (28) and (29). Thus, equation (4) yields the mean-flow equation in the form

$$\begin{eqnarray}&&\begin{array}{l}{mn}\left({\partial }_{t}{\boldsymbol{U}}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{\boldsymbol{U}}\right)\\ \,\,\,\,\,\,=-{\rm{\nabla }}\cdot { \mathcal P }+{qn}{\bf{E}}+\left({qn}{\boldsymbol{U}}+{\rm{\nabla }}\times {\boldsymbol{ \mathcal M }}\right)\times {\bf{B}}.\end{array}\end{eqnarray} \tag{ 54 }$$

We recall that this equation is accompanied by the transport equations (34) and (35) for the magnetic moment density and the particle density, respectively. Thus, we notice that, even after enforcing a gyrotropic truncation so that $\int {c}_{\parallel }{{\boldsymbol{w}}}_{\perp }\tilde{f}\,{\rm{d}}{c}_{\parallel }{\rm{d}}\mu =0$, the guiding-center motion of fluctuation dynamics manifests in the mean-flow equation through the magnetic drift and the parallel pressure contribution. Indeed, in this case equation (54) leads to the gyrofluid model given by the equation

$$\begin{eqnarray}&&\begin{array}{l}{mn}\left({\partial }_{t}{\boldsymbol{U}}+{\boldsymbol{U}}\cdot {\rm{\nabla }}{\boldsymbol{U}}\right)\\ \quad =\,{qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {\bf{B}}-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})-\bar{\mu }{\rm{\nabla }}B,\end{array}\end{eqnarray} \tag{ 55 }$$

which is itself accompanied by the parallel CGL equation (40), as well as (34) and (35). This suggests that, when fluctuations are goverend by guiding-center motion, anisotropic fluid models based on the CGL pressure tensor ${p}_{\parallel }{\boldsymbol{b}}+{p}_{\perp }({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}})$ should actually be accompanied by a magnetization term appearing in the perpendicular pressure and arising from the relation

$$\begin{eqnarray*}\begin{array}{rcl}-\bar{\mu }{\rm{\nabla }}B & = & {\rm{\nabla }}{\bf{B}}\cdot {\boldsymbol{M}}=-{\bf{B}}\times {\rm{\nabla }}\times {\boldsymbol{M}}\\ & & +{\rm{\nabla }}\cdot \left[({\bf{B}}\cdot {\boldsymbol{M}}){\bf{1}}-{\bf{B}}{\boldsymbol{M}}\right]\\ & = & {\bf{B}}\times {\rm{\nabla }}\times (\bar{\mu }{\boldsymbol{b}})-{\rm{\nabla }}\cdot \left(\bar{\mu }B({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}}\right).\end{array}\end{eqnarray*}$$

However, the magnetization term cancels in the presence of a frozen-in condition for the magnetic field [26, 27] and this is the reason why magnetization effects do not appear explicitly in the original CGL model.

These last two sections showed how the mean-fluctuation splitting leads to different modeling options depending on the description adopted for the mean flow and the fluctuation dynamics. In the next sections, we shall go back to standard drift-fluid models in which both the mean flow and fluctuations are described by guiding-center theory and we shall illustrate the level of difficulty of these models when they are coupled to the Maxwell equations for the self-consistent evolution of the electromagnetic field.

In the presence of dynamic electromagnetic fields, the mean-fluctuation splitting Lagrangian (21) is added to the Maxwell Lagrangian

$$\begin{eqnarray}&&{l}_{\mathrm{Max}}=\displaystyle \frac{1}{2}\int \,\left({\varepsilon }_{0}{\left|{\rm{\nabla }}{\rm{\Phi }}+\displaystyle \frac{\partial {\bf{A}}}{\partial t}\right|}^{2}-\displaystyle \frac{1}{{\mu }_{0}}| {\rm{\nabla }}\times {\bf{A}}{| }^{2}\right){{\rm{d}}}^{3}x,\end{eqnarray} \tag{ 58 }$$

so that the new Lagrangian is

$$\begin{eqnarray}&&{\ell }=l+{l}_{\mathrm{Max}},\end{eqnarray} \tag{ 59 }$$

where l is given by (21). We take the occasion to emphasize that, consistently with Littlejohn's theory, the Lagrangian (59) assumes an E × B speed much smaller than the particle thermal speed [18]. Alternatively, should these speeds be comparable, we would simply subtract the E × B energy term ${m}_{h}{B}^{-2}{\left|{\bf{E}}\times {\boldsymbol{b}}\right|}^{2}/2$ from the guiding center kinetic energy ${m}_{h}{v}_{\parallel }^{2}/2+\mu B$, thereby producing extra magnetization terms along with polarization effects [32] which are instead ignored here. Then, as a result of the Euler–Poincaré variational principle $\delta {\int }_{{t}_{1}}^{{t}_{2}}{\ell }\,{\rm{d}}t=0$, the equations (4)–(6) do not change in form, although the time-dependent magnetic field now appears in (5) so that (23) becomes

$$\begin{eqnarray*}&&\begin{array}{l}m\left[{\boldsymbol{b}}\cdot ({\boldsymbol{a}}+{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}})+({\bf{u}}+{\boldsymbol{U}})\cdot {\rm{\nabla }}{\boldsymbol{b}}\cdot {\bf{c}}+\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}\cdot {\bf{c}}\right]{\boldsymbol{b}}\\ \,-\,q({\bf{w}}+{\boldsymbol{U}})\times {{\boldsymbol{ \mathcal B }}}^{* }=-m{\boldsymbol{b}}{\partial }_{t}({\boldsymbol{b}}\cdot {\boldsymbol{U}})\\ \,+\,q{{\bf{E}}}^{* }-m({\boldsymbol{b}}\cdot {\bf{c}}){\rm{\nabla }}({\boldsymbol{b}}\cdot {\boldsymbol{U}})-m\displaystyle \frac{1}{2}{\rm{\nabla }}{\left({\boldsymbol{b}}\cdot {\boldsymbol{U}}\right)}^{2},\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{E}}}^{* } & = & -\displaystyle \frac{\partial {\bf{A}}}{\partial t}-{\rm{\nabla }}{\rm{\Phi }}\\ & & -(\mu /q){\rm{\nabla }}B-(m/q)({\boldsymbol{b}}\cdot {\bf{c}}+{\boldsymbol{b}}\cdot {\boldsymbol{U}})\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}.\end{array}\end{eqnarray} \tag{ 60 }$$

Now, in the case of a time-dependent magnetic field, the relation (24) is modified as [3]

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial \tilde{f}}{\partial t}+{\rm{\nabla }}\cdot \displaystyle \int ({\bf{w}}+{\boldsymbol{U}})F\,{{\rm{d}}}^{2}{c}_{\perp }+\displaystyle \frac{\partial }{\partial {c}_{\parallel }}\displaystyle \int \left[{\boldsymbol{b}}\cdot ({\boldsymbol{a}}+{\bf{c}}\cdot {\rm{\nabla }}{\boldsymbol{U}})\right.\\ \quad \left.+\,({\bf{w}}+{\boldsymbol{U}})\cdot {\rm{\nabla }}{\boldsymbol{b}}\cdot {\bf{c}}+{\bf{c}}\cdot {\partial }_{t}{\boldsymbol{b}}\right]F\,{{\rm{d}}}^{2}{c}_{\perp }=0.\end{array}\end{eqnarray} \tag{ 61 }$$

and since (22) remains unchanged, one obtains formally the same drift-kinetic equation (25), although with the definition (60) in the form ${{\bf{E}}}^{* }={\bf{E}}-(\mu /q){\rm{\nabla }}B-(m/q)({c}_{\parallel }+{U}_{\parallel }){\partial }_{t}{\boldsymbol{b}}$. Notice that the drift-kinetic equation (61) can be again written as in (26)–(27).

As a next step, we turn our attention to the fluid equation (4). Since the functional derivatives of the Lagrangian ℓ in (59) with respect to $({\boldsymbol{U}},{\bf{X}},F)$ are the same as the functional derivatives of the Lagrangian l in (21) with respect to $({\boldsymbol{U}},{\bf{X}},F)$, equation (4) produces the following variant of the drift-fluid equation (30):

$$\begin{eqnarray}&&\begin{array}{l}{mn}({\partial }_{t}{U}_{\parallel }+{\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}=-{\rm{\nabla }}\cdot { \mathcal P }+{qn}{\bf{E}}\\ \quad -\,{{mnU}}_{\parallel }{\partial }_{t}{\boldsymbol{b}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }-{\bf{B}}\times {\rm{\nabla }}\times {\boldsymbol{ \mathcal M }},\end{array}\end{eqnarray} \tag{ 62 }$$

where ${\boldsymbol{ \mathcal M }}({\bf{x}},t)$ and ${ \mathcal P }({\bf{x}},t)$ are given by (28) and (29), respectively, and we have used (60) in the relation (27) defining the vector field ${\boldsymbol{w}}({\bf{x}},{c}_{\parallel },t)$. Notice that the electric field and its gradients appear on the right-hand side of (62) in three different instances: while the first appearance is explicit, the second is through ${\partial }_{t}{\boldsymbol{b}}$ (by Faraday's law, ${\partial }_{t}{\bf{B}}=-{\rm{\nabla }}\times {\bf{E}}$), and the third is through the vector field ${\boldsymbol{w}}$ itself, which is contained in the magnetization ${\boldsymbol{ \mathcal M }}$. Again, since the perpendicular component of the velocity ${{\boldsymbol{U}}}_{\perp }$ is absent in ${\boldsymbol{ \mathcal M }}({\bf{x}},t)$, ${ \mathcal P }({\bf{x}},t)$, and ${\partial }_{t}{\boldsymbol{b}}$, this allows finding an explicit algebraic expression for the velocity and equation (31) is modified as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{U}} & = & {U}_{\parallel }\displaystyle \frac{{{\boldsymbol{B}}}^{* }}{{B}_{\parallel }^{* }}+\displaystyle \frac{{\boldsymbol{b}}}{{{mnB}}_{\parallel }^{* }}\\ & & \times \left({\bf{B}}\times {\rm{\nabla }}\times {\boldsymbol{ \mathcal M }}+{\rm{\nabla }}\cdot {\mathbb{P}}-{qn}{\bf{E}}+{{mnU}}_{\parallel }{\partial }_{t}{\boldsymbol{b}}\right).\end{array}\end{eqnarray} \tag{ 63 }$$

At this point, we are ready to write down Maxwell's equations by taking variations of ℓ in (59) with respect to Φ and ${\bf{A}}$. While the first yield the usual Gauss law ${\rm{\nabla }}\cdot {\bf{E}}={qn}$, the second leads to the following Ampère's law

$$\begin{eqnarray}&&{\mu }_{0}^{-1}{\rm{\nabla }}\times {\bf{B}}-{\varepsilon }_{0}\displaystyle \frac{\partial {\bf{E}}}{\partial t}={qn}{\boldsymbol{U}}+{\rm{\nabla }}\times \left({\boldsymbol{ \mathcal M }}+\displaystyle \frac{{{mnU}}_{\parallel }}{B}\,{{\boldsymbol{U}}}_{\perp }\right).\end{eqnarray} \tag{ 64 }$$

We notice that the term in parenthesis coincides with the full magnetization (19) appearing in the momentum equation (17). Upon recalling the vector field ${\boldsymbol{w}}$ in (27), with the definition (60) in the form ${{\bf{E}}}^{* }={\bf{E}}-(\mu /q){\rm{\nabla }}B\,-(m/q)({c}_{\parallel }+{U}_{\parallel }){\partial }_{t}{\boldsymbol{b}}$, this fact is verified by substituting ${v}_{\parallel }={c}_{\parallel }+{U}_{\parallel }$ and ${\bf{u}}={\boldsymbol{w}}+{\boldsymbol{U}}$ in (19). Here, it is useful to explicitly adapt equation (28) to the present case, that is

$$\begin{eqnarray*}\begin{array}{rcl}{\boldsymbol{ \mathcal M }}({\bf{x}},t) & = & -\bar{\mu }{\boldsymbol{b}}+\displaystyle \frac{m}{B}\displaystyle \int {c}_{\parallel }[({c}_{\parallel }+{U}_{\parallel })\displaystyle \frac{f{{\bf{B}}}_{\perp }^{\ast }}{{B}_{\parallel }^{\ast }}-\displaystyle \frac{f{\boldsymbol{b}}}{{B}_{\parallel }^{\ast }}\\ & & \times \left({\bf{E}}-\mu {\rm{\nabla }}B-({c}_{\parallel }+{U}_{\parallel }){{\rm{\partial }}}_{t}{\boldsymbol{b}}\right)]\,{\rm{d}}{c}_{\parallel }\,{\rm{d}}\mu ,\end{array}\end{eqnarray*}$$

where ${{\bf{B}}}^{* }={\bf{B}}+(m/q)({c}_{\parallel }+{{\boldsymbol{U}}}_{\parallel }){\rm{\nabla }}\times {\boldsymbol{b}}$. This shows that gradients of the electric field up to third order appear in the right-hand side of Ampère's law (64) through ${\boldsymbol{ \mathcal M }}$ and ${{\boldsymbol{U}}}_{\perp }$, thereby making this equation extremely challenging in practice.

Finally, we conclude that applying the mean-fluctuation splitting method to the guiding-center-Maxwell system leads to the equations of motion (25), (62)–(63), and (64). These are accompanied by the definitions (28) and (29), where the vector field ${\boldsymbol{w}}$ satisfies (27) upon using ${{\bf{E}}}^{* }={\bf{E}}\,-(\mu /q){\rm{\nabla }}B\,-(m/q)({c}_{\parallel }+{U}_{\parallel }){\partial }_{t}{\boldsymbol{b}}$. However, as remarked above, this system is particularly challenging because of the third-order gradients of the electric field occurring in Ampère's law (64). This problem is partially addressed in the fluid approximation, for which only second order gradients appear in the current balance. This is presented in the next section.

Upon enforcing the gyrotropic approximation (33), the mean-flow equation (62) reduces to

$$\begin{eqnarray}&&\begin{array}{l}{mn}({\partial }_{t}{U}_{\parallel }+{\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}=-{\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})\\ \quad +\,{qn}{\bf{E}}-\bar{\mu }{\rm{\nabla }}B+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }-{{mnU}}_{\parallel }{\partial }_{t}{\boldsymbol{b}}.\end{array}\end{eqnarray} \tag{ 65 }$$

Already at this stage, upon replacing Faraday's law in the last term, we see that, unlike equation (49), only first-order gradients of the electric field appear on the right-hand side. Then, for consistency, one applies (33) to Ampère's law (64) thereby obtaining

$$\begin{eqnarray}&&{\mu }_{0}^{-1}{\rm{\nabla }}\times {\bf{B}}-{\varepsilon }_{0}\displaystyle \frac{\partial {\bf{E}}}{\partial t}={qn}{\boldsymbol{U}}-{\rm{\nabla }}\times \left(\bar{\mu }{\boldsymbol{b}}-\displaystyle \frac{{{mnU}}_{\parallel }}{B}\,{{\boldsymbol{U}}}_{\perp }\right),\end{eqnarray} \tag{ 66 }$$

where the fluid velocity is found from (65) as follows:

$$\begin{eqnarray*}\begin{array}{rcl}{\boldsymbol{U}} & = & {U}_{\parallel }\displaystyle \frac{{{\boldsymbol{B}}}^{* }}{{B}_{\parallel }^{* }}+\displaystyle \frac{{\boldsymbol{b}}}{{{mnB}}_{\parallel }^{* }}\\ & & \times \left({\rm{\nabla }}\cdot ({p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}})+\bar{\mu }{\rm{\nabla }}B-{qn}{\bf{E}}+\displaystyle \frac{{{mnU}}_{\parallel }}{B}{\left({\rm{\nabla }}\times {\bf{E}}\right)}_{\perp }\right).\end{array}\end{eqnarray*}$$

Then, we notice that Ampère's law (66) now involves up to second-order gradients of the magnetic field, thereby providing a simplification over the full mean-fluctuation system. It is easy to see that analogous second-order gradient terms also appear in the current balance (57) of the full drift-kinetic Maxwell system.

In the present setting, an auxiliary equation must be retained for the parallel pressure and this can be provided by the CGL equation (40), which is indeed unaffected by the coupling to Maxwell's equations. Unless one resorts to the cold plasma closure (36), a further possibility to deal with the pressure terms in (65) is provided by the PCR approach, which replaces the parallel pressure and magnetic moment terms by the gradient of an isotropic pressure associated to an internal energy. Alternatively, one can resort to the CGL equation of state (44) from section 3.3, so that the internal energy in (45) now depends also on the magnetic field, i.e. ${ \mathcal U }={ \mathcal U }(n,s,B)$. In this case, one considers the variational principle determined by the sum of the Maxwell Lagrangian ${l}_{\mathrm{Max}}$ and the Lagrangian in (46). Then equation (67) is modified as follows:

$$\begin{eqnarray}&&\begin{array}{l}{mn}({\partial }_{t}{U}_{\parallel }+{\boldsymbol{U}}\cdot {\rm{\nabla }}{U}_{\parallel }){\boldsymbol{b}}={qn}{\bf{E}}+{qn}{\boldsymbol{U}}\times {{\boldsymbol{B}}}^{* }-{{mnU}}_{\parallel }{\partial }_{t}{\boldsymbol{b}}\\ -\,{\rm{\nabla }}\cdot [{p}_{\parallel }{\boldsymbol{b}}{\boldsymbol{b}}+\bar{\mu }B({\bf{1}}-{\boldsymbol{b}}{\boldsymbol{b}})]+{\bf{B}}\times {\rm{\nabla }}\times ((\bar{\mu }+{p}_{\parallel }/B){\boldsymbol{b}}),\end{array}\end{eqnarray} \tag{ 67 }$$

while Ampère's law becomes

$$\begin{eqnarray}&&\begin{array}{l}{\mu }_{0}^{-1}{\rm{\nabla }}\times {\bf{B}}-{\varepsilon }_{0}\displaystyle \frac{\partial {\bf{E}}}{\partial t}\\ \quad =\,{qn}{\boldsymbol{U}}-{\rm{\nabla }}\times \left[\displaystyle \frac{1}{B}\left({p}_{\parallel }{\boldsymbol{b}}+\bar{\mu }{\bf{B}}-{{mnU}}_{\parallel }{{\boldsymbol{U}}}_{\perp }\right)\right].\end{array}\end{eqnarray} \tag{ 68 }$$

Once again, we observe how the enforcement of the longitudinal adiabatic invariant through the CGL equation of state (44) leads to redefining the magnetization current via the replacement $\bar{\mu }\to \bar{\mu }+{p}_{\parallel }/B$ within the magnetization density vector. Also, we emphasize that the coupling to Maxwell's equations leaves preservation of Kelvin's circulation entirely unaffected thereby leading again to the conservation of the hydrodynamic helicity.

In this section, we present a revisitation of guiding-center theory that averages out the fast gyromotion around a high-intensity background magnetic field. We emphasize that this is only one among several possible approaches and, in the present case, the main advantage is probably the level of simplicity of the resulting theory.

As a first step, we decompose the total magnetic field into a background and a time-dependent part, that is

$$\begin{eqnarray*}&&{\bf{B}}({\bf{x}},t)={{\bf{B}}}_{0}({\bf{x}})+\widetilde{{\bf{B}}}({\bf{x}},t).\end{eqnarray*}$$

Notice that we ${\rm{\nabla }}\cdot {\bf{B}}({\bf{x}},t)=0$ enforces ${\rm{\nabla }}\cdot \widetilde{{\bf{B}}}({\bf{x}},t)\,=-{\rm{\nabla }}\,\cdot {{\bf{B}}}_{0}({\bf{x}})$. Here, the modified guiding-center theory is obtained by resorting to a treatment recently appeared in [33], based on kinematic relations and averaged Lagrangians. Although closely related to Littlejohn's approach in [5], the present method does not invoke any previous knowledge on the gyroradius or the magnetic moment. As customary in guiding-center theory, we start the discussion with the phase-space Lagrangian in the form

$$\begin{eqnarray}&&L({\bf{r}},\dot{{\bf{r}}})=(\epsilon {\bf{v}}+{\bf{A}}({\bf{r}},t))\cdot \dot{{\bf{r}}}-\displaystyle \frac{\epsilon }{2}{v}^{2},\end{eqnarray} \tag{ 69 }$$

where ${{\bf{B}}}_{0}({\bf{x}})+\widetilde{{\bf{B}}}({\bf{x}},t)={\rm{\nabla }}\times {\bf{A}}({\bf{x}},t)$ and is the smallness parameter. Notice that here we have omitted the electric potential, which can always be restored without essential changes in the final result. Following the treatment in [33], here we decompose the particle position as follows

$$\begin{eqnarray}&&{\bf{r}}(t)={\bf{X}}(t)+\epsilon \rho (t){\boldsymbol{a}}({\bf{X}}(t),t),\end{eqnarray} \tag{ 70 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{a}}({\bf{x}},t)=R({\bf{x}},t){{\bf{e}}}_{1}({\bf{x}}).\end{eqnarray} \tag{ 71 }$$

While ${{\bf{e}}}_{1}({\bf{x}})$ is an arbitrary field such that ${{\boldsymbol{b}}}_{0}\times {{\bf{e}}}_{1}=0$, the rotation matrix R is constructed as a rotation around ${{\boldsymbol{b}}}_{0}$ by an angle θ and it is written explicitly as

$$\begin{eqnarray*}&&\begin{array}{l}R(\theta (t),{{\boldsymbol{b}}}_{0}({\bf{x}}))={\bf{1}}+\sin \theta (t){\widehat{{\boldsymbol{b}}}}_{0}({\bf{x}})\\ \quad +\,(1-\cos \theta (t)){\widehat{{\boldsymbol{b}}}}_{0}({\bf{x}}){\widehat{{\boldsymbol{b}}}}_{0}({\bf{x}}),\end{array}\end{eqnarray*}$$

where we write $\theta ={\epsilon }^{-1}{\rm{\Theta }}$ and we use the hat map introduced in (50). Further details on these definitions and the following calculations are found in [33]. Notice that here we are restricting to consider a gyroradius that is always perpendicular to the background magnetic field, while in the conventional approach the gyroradius is perpendicular to the total magnetic field. As customary in guiding-center motion, we expand the magnetic potential as ${\bf{A}}({\bf{r}},t)={\bf{A}}({\bf{X}},t)+\epsilon \rho {\boldsymbol{a}}({\bf{X}},t)\cdot {\rm{\nabla }}{\bf{A}}({\bf{X}},t)$ to write

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\bf{A}}({\bf{r}},t)\cdot \dot{{\bf{r}}} & = & {\bf{A}}({\bf{X}},t)\cdot \dot{{\bf{X}}}\\ & & +\displaystyle \frac{\epsilon }{2}{\rho }^{2}\dot{{\rm{\Theta }}}{{\boldsymbol{b}}}_{0}({\bf{X}})\cdot {\bf{B}}({\bf{X}},t)+{ \mathcal G }+O({\epsilon }^{2}).\end{array}\end{array}\end{eqnarray} \tag{ 72 }$$

Here, ${ \mathcal G }$ denotes the sum of terms that are either total time derivatives or linear in ${\boldsymbol{a}}({\bf{X}},t)$. The total time derivatives in the Lagrangian are well known to be irrelevant to the resulting dynamics. Likewise, the terms linear in ${\boldsymbol{a}}({\bf{X}},t)$ vanish under the final averaging process over θ and are therefore of no interest for the present purpose. As a further step, we use an orthonormal basis $({{\boldsymbol{b}}}_{0},{\boldsymbol{a}},{\boldsymbol{h}})$ to write ${\bf{v}}={v}_{\parallel }{{\boldsymbol{b}}}_{0}+n{\boldsymbol{a}}+w{\boldsymbol{h}}$, so that

$$\begin{eqnarray}\begin{array}{rcl}\epsilon {\bf{v}}\cdot \dot{{\bf{r}}} & = & \epsilon {\bf{v}}\cdot (\dot{{\bf{X}}}+\rho \dot{{\rm{\Theta }}}{{\boldsymbol{b}}}_{0}({\bf{X}})\times {\boldsymbol{a}}({\bf{X}},t))\\ & & +O(\epsilon )=\epsilon {\bf{v}}\cdot \dot{{\bf{X}}}-\epsilon \rho \dot{{\rm{\Theta }}}w+O(\epsilon ).\end{array}\end{eqnarray} \tag{ 73 }$$

Finally, replacing (72) and (73) in (69) and averaging over θ leads to the following averaged Lagrangian

$$\begin{eqnarray*}\begin{array}{rcl}\langle L\rangle & = & ({v}_{\parallel }{{\boldsymbol{b}}}_{0}+{\epsilon }^{-1}{\bf{A}})\cdot \dot{{\bf{X}}}+\left(\displaystyle \frac{1}{2}{\epsilon }^{-1}{\rho }_{L}{B}_{\parallel }-w\right){\rho }_{L}\dot{\theta }\\ & & -\displaystyle \frac{1}{2}({v}_{\parallel }^{2}+{w}^{2}+{n}^{2}),\end{array}\end{eqnarray*}$$

where we have introduced the Larmor radius ${\rho }_{L}=\epsilon \rho$ and we have defined ${B}_{\parallel }={{\boldsymbol{b}}}_{0}\cdot {\bf{B}}$. Thus, taking variations yields n = 0 along with

$$\begin{eqnarray}&&{\rho }_{L}=\displaystyle \frac{\epsilon w}{{B}_{\parallel }},\qquad \quad \dot{\theta }=-{\epsilon }^{-1}{B}_{\parallel },\qquad \quad {\rho }_{L}^{2}{B}_{\parallel }={const}.\end{eqnarray} \tag{ 74 }$$

We notice that here the gyration frequency associated to the gyromotion around b₀ is actually time-dependent as it comprises contributions from the time-dependent magnetic field, since ${B}_{\parallel }={B}_{0}+({{\boldsymbol{b}}}_{0}\cdot \tilde{{\boldsymbol{b}}})\widetilde{B}$.

At this point, the above Lagrangian can be cast in a more familiar form upon using the first relation in (74) and introducing the parallel magnetic moment ${\mu }_{\parallel }={{mw}}^{2}/(2{B}_{\parallel })$. Upon restoring the electric potential, one obtains

$$\begin{eqnarray}&&L=({{mv}}_{\parallel }{{\boldsymbol{b}}}_{0}+q{\bf{A}})\cdot \dot{{\bf{X}}}-\epsilon {\mu }_{\parallel }\dot{\theta }-\displaystyle \frac{m}{2}{v}_{\parallel }^{2}-{\mu }_{\parallel }{{\boldsymbol{b}}}_{0}\cdot {\bf{B}}-q{\rm{\Phi }}.\end{eqnarray} \tag{ 75 }$$

This is a modification of Littlejohn's Lagrangian in the sense that the time-dependent magnetic potential ${\bf{A}}$ appears only in the magnetic moment term and in the minimal coupling term, while all other terms involve only the fixed background field ${{\bf{B}}}_{0}$. The resulting equations of motion read

$$\begin{eqnarray*}&&q{{\bf{E}}}^{\dagger }+q\dot{{\bf{X}}}\times {{\bf{B}}}^{\dagger }=m{\dot{v}}_{\parallel }{{\boldsymbol{b}}}_{0},\end{eqnarray*}$$

where we have defined

$$\begin{eqnarray*}&&{{\bf{B}}}^{\dagger }={\bf{B}}+(m/q){v}_{\parallel }{\rm{\nabla }}\times {{\boldsymbol{b}}}_{0},\qquad {{\bf{E}}}^{\dagger }={\bf{E}}-({\mu }_{\parallel }/q){\rm{\nabla }}{B}_{\parallel }.\end{eqnarray*}$$

The explicit equations of motion can then be derived as usual by using standard vector algebra. Then, upon defining ${B}_{\parallel }^{\dagger }={B}_{\parallel }+(m/q){v}_{\parallel }{{\boldsymbol{b}}}_{0}\cdot {\rm{\nabla }}\times {{\boldsymbol{b}}}_{0}$, the corresponding kinetic equation reads

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial F}{\partial t}+\displaystyle \frac{\partial }{\partial {\bf{x}}}\cdot \left[\displaystyle \frac{F}{{B}_{\parallel }^{\dagger }}\left({v}_{\parallel }{{\bf{B}}}^{\dagger }-{{\boldsymbol{b}}}_{0}\times {{\bf{E}}}^{\dagger }\right)\right]+\,\displaystyle \frac{q}{m}\displaystyle \frac{\partial }{\partial {v}_{\parallel }}\left(\displaystyle \frac{F}{{B}_{\parallel }^{\dagger }}\,{{\bf{B}}}^{\dagger }\cdot {{\bf{E}}}^{\dagger }\right)=0.\end{array}\end{eqnarray} \tag{ 76 }$$

At this point, the mean-fluctuation splitting may be applied thereby returning an alternative class of gyrofluid models to those presented previously. In this case, the gyrotropic condition (33) (here, gyrotropic means with respect to b₀), it is easy to see that, upon neglecting heat-flux contributions, the parallel pressure obeys equation (40) with ${\boldsymbol{b}}$ replaced by b₀. Here, one can resort to the CGL equation of state from section 3.3 (upon replacing B → B₀) or, alternatively, to a non-gyrotropic pressure by following the approach in section 4.1.

While the pressure tensor arising from the present theory is non-gyrotropic, one advantage of this model is the fact that magnetization currents do not lead to the difficulties mentioned in section 5.1. Indeed, upon denoting ${\bar{\mu }}_{\parallel }\,=\int F{\mu }_{\parallel }\,{\rm{d}}v{\rm{d}}{\mu }_{\parallel }$, Ampère's law corresponding to the present guiding-center model reads as follows:

$$\begin{eqnarray}&&{\mu }_{0}^{-1}{\rm{\nabla }}\times {\bf{B}}-{\varepsilon }_{0}\displaystyle \frac{\partial {\bf{E}}}{\partial t}={qn}{\boldsymbol{U}}-{\rm{\nabla }}\times ({\bar{\mu }}_{\parallel }{{\boldsymbol{b}}}_{0}),\end{eqnarray} \tag{ 77 }$$

which represents a substantial simplification to equation (64). A similar Ampère law is obtained by enforcing the CGL equation of state ${p}_{\parallel }\,=\,h(s){m}^{3}{n}^{3}/{B}_{0}^{2}$, except in this case one replaces ${\bar{\mu }}_{\parallel }\to {p}_{\parallel }/B+{\bar{\mu }}_{\parallel }$ in (77). Then, the present formulation offers the opportunity for drift-ordered models carrying the self consistent evolution of the electromagnetic field while avoiding the difficulties encountered in the standard guiding-center theory. In the next section, we shall present an application of this guiding-center model to a hybrid kinetic-fluid scheme for energetic particle effects.

This section presents an application of the modified guiding-center theory beyond drift-fluid models. More specifically, we shall follow the variational approach in [34] to formulate a hybrid kinetic-drift MHD model in the current-coupling scheme (CCS) aimed at capturing energetic particle effects in fusion devices. Further details are presented in [1, 4, 34–38], to which we refer the reader for more discussions on hybrid MHD models for energetic particles. Hybrid MHD models describe coupled systems comprising a MHD fluid bulk that interacts with an ensemble of energetic alpha particles obeying kinetic theory (drift-kinetic in this case). These energetic particles are naturally generated in fusion reactions and their effects have been studied over several decades [39]. Motivated by the variational principles presented in [4, 34], here we shall construct a hybrid model starting from the following variational principle

$$\begin{eqnarray}&&\begin{array}{l}\delta {\displaystyle \int }_{{t}_{1}}^{{t}_{2}}\left[\displaystyle \int F\left(\left({{mv}}_{\parallel }{{\boldsymbol{b}}}_{0}+q{\bf{A}}\right)\cdot {\boldsymbol{u}}\right.\right.\\ \quad \left.-\,\displaystyle \frac{m}{2}{v}_{\parallel }^{2}-{\mu }_{\parallel }{B}_{\parallel }-q{\boldsymbol{U}}\cdot {\bf{A}}\right){{\rm{d}}}^{3}x\,{\rm{d}}{v}_{\parallel }\,{\rm{d}}\mu \\ \quad \left.+\,\displaystyle \int \left(\displaystyle \frac{\rho }{2}{U}^{2}-\rho { \mathcal U }(\rho )-\displaystyle \frac{1}{2{\mu }_{0}}{B}^{2}\right){{\rm{d}}}^{3}x\right]\,{\rm{d}}t=0.\end{array}\end{eqnarray} \tag{ 78 }$$

Here, the Lagrangian is given by the sum of the guiding-center Lagrangian (first two lines) and the MHD Lagrangian (last line). The minimal coupling term at the end of the second line emerges as a gauge choice; see [34] for further discussions. The variational principle (78) involves the variations

$$\begin{eqnarray*}&&\delta {\boldsymbol{U}}={\partial }_{t}{\boldsymbol{\Xi }}+[{\boldsymbol{U}},{\boldsymbol{\Xi }}],\qquad \quad \delta \rho =-{\rm{\nabla }}\cdot (\rho {\boldsymbol{\Xi }})\end{eqnarray*}$$

for ρ and ${\boldsymbol{U}}$, while we also have

$$\begin{eqnarray*}\begin{array}{rcl}\delta {\bf{A}} & = & {\boldsymbol{\Xi }}\times {\bf{B}}-{\rm{\nabla }}({\boldsymbol{\Xi }}\cdot {\bf{A}}),\qquad \delta {\bf{X}}={\partial }_{t}{\bf{Y}}+[{\bf{X}},{\bf{Y}}],\\ \delta F & = & -{{\rm{\nabla }}}_{{\bf{z}}}(F{\bf{Y}}).\end{array}\end{eqnarray*}$$

Here, ${\bf{X}}({\bf{x}},{v}_{\parallel })=({\boldsymbol{u}}({\bf{x}},{v}_{\parallel }),a({\bf{x}},{v}_{\parallel }))$ is the guiding-center vector field such that ${\partial }_{t}F+{\rm{\nabla }}\cdot ({\boldsymbol{u}}F)+{\partial }_{{v}_{\parallel }}({aF})=0$ and ${\bf{Y}}$ is again a phase-space displacement vector field analogous to ${\boldsymbol{\Xi }}$. The resulting hybrid model is formed of equation (76) with ${{\bf{E}}}^{\dagger }=-{\boldsymbol{U}}\times {\bf{B}}-({\mu }_{\parallel }/q){\rm{\nabla }}{B}_{\parallel }$, and the two advection equations

$$\begin{eqnarray}&&{\partial }_{t}{\bf{B}}+{\rm{\nabla }}\times ({\bf{B}}\times {\boldsymbol{U}})=0,\qquad \qquad {\partial }_{t}\rho +{\rm{\nabla }}\cdot (\rho {\boldsymbol{U}})=0.\end{eqnarray} \tag{ 79 }$$

Also, we obtain the following MHD momentum equation

$$\begin{eqnarray}&&\begin{array}{l}\rho {\partial }_{t}{\boldsymbol{U}}+\rho \left({\boldsymbol{U}}\cdot {\rm{\nabla }}\right){\boldsymbol{U}}\\ =\,\left({\mu }_{0}^{-1}{\rm{\nabla }}\times {\bf{B}}+{qn}{\boldsymbol{U}}-{{\bf{J}}}_{\mathrm{gc}}-{\rm{\nabla }}\times {{\bf{M}}}_{\mathrm{gc}}\right)\times {\bf{B}}-{\rm{\nabla }}p,\end{array}\end{eqnarray} \tag{ 80 }$$

where ${{\bf{M}}}_{\mathrm{gc}}=-\int F{\mu }_{\parallel }{{\boldsymbol{b}}}_{0}\,{\rm{d}}{v}_{\parallel }{\rm{d}}{\mu }_{\parallel }$ and ${{\bf{J}}}_{\mathrm{gc}}=q\int F({v}_{\parallel }{{\bf{B}}}^{\dagger }-{{\boldsymbol{b}}}_{0}\,\times {{\bf{E}}}^{\dagger })/{B}_{\parallel }^{\dagger }\,{\rm{d}}{v}_{\parallel }{\rm{d}}{\mu }_{\parallel }$.

While this equation has a similar structure to the corresponding momentum equation in the Hamiltonian CCS from [34] (see equations (35) and (37) therein), we notice that substantial simplifications occur in the present case. In first instance, no time derivative of the magnetic field orientation occurs in the effective electric field ${{\bf{E}}}^{\dagger }$. Also, and perhaps more importantly, the magnetization term in (80) does not involve moments of the kinetic distribution other than the magnetization density. This last aspect establishes an analogy between the present (simplified) Hamiltonian CCS model and the non-Hamiltonian variant implemented in the MEGA code [36]. Indeed, one of the main differences between the Hamiltonian CCS in [34] and its non-Hamiltonian counterpart lies indeed in the fact that the former is substantially more involved due to the presence of extra magnetization terms involving higher velocity moments of the guiding-center distribution. This difficulty is removed in the present variant, which also conserves the standard cross helicity $\int {\bf{B}}\cdot {\boldsymbol{U}}\,{{\rm{d}}}^{3}x$. As noticed in [34], the momentum equation underlying the MEGA code breaks energy conservation and this is mainly due to the fact that the magnetization density defined in the MEGA code involves the time-dependent magnetic field direction ${\boldsymbol{b}}({\bf{x}},t)$. Indeed, the presence of the entire time-dependent magnetic field in the magnetization density requires the addition of the extra terms in (19), which instead are ignored in the MEGA code [36]. On the other hand, in the case of the model (76)–(80) energy conservation is still ensured by the variational structure (78), while the definition of the magnetization density is substantially simplified.