2. Fluid equations of a multispecies plasma

We consider a magnetised plasma containing an arbitrary number of distinct species, indexed by the label $s$ . The fluid equations for a given species $s$ , of mass $m_s$ and electric charge $q_s$ , are obtained by taking velocity moments of the Boltzmann equation (Braginskii Reference Braginskii1965). The first three moment equations, evolving the density $n_s$ , the fluid momentum density $\boldsymbol{\mathcal{M}}_s = m_s n_s \boldsymbol{V}_s$ , with $\boldsymbol{V}_s$ the mean velocity of species $s$ , and the scalar pressure $p_s = n_s T_s$ , with $T_s$ the temperature of species $s$ , are given by

(2.1) \begin{align} &\frac {\partial n_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n_s \boldsymbol{V}_s) = \mathcal{S} _s,\\[-10pt]\nonumber \end{align}

(2.2) \begin{align} &\frac {\partial \boldsymbol{\mathcal{M}}_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{V}_s \boldsymbol{\mathcal{M}}_s) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s - q_s n_s\left (\boldsymbol{E} + \boldsymbol{V}_s \times \boldsymbol{B}\right ) = \boldsymbol{\mathcal{R}}_s ,\\[-10pt]\nonumber \end{align}

(2.3) \begin{align} &\frac {\partial }{\partial t}\left (\frac {3}{2} p_s \right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {3}{2}p_s \boldsymbol{V}_s + \boldsymbol{q}_s \right ) + \boldsymbol{P}_s : \boldsymbol{\nabla }\boldsymbol{V}_s= \mathcal{Q}_s + r_s \mathcal{K}_s , \end{align}

where $\boldsymbol{P}_s = p_s \mathbb{I} + \boldsymbol{\pi }_s$ is the pressure tensor of the species, $\boldsymbol{\pi }_s$ is the stress tensor, $\mathbb{I}$ is the identity matrix, $\boldsymbol{q}_s$ is the heat-flux density and $\mathcal{S}_s, \boldsymbol{\mathcal{R}}_s, \mathcal{Q}_s$ are source terms, arising from external drives and interactions between the species, e.g. via Coulomb collisions. In addition, we define $r_s \equiv \mathcal{S}_s/n_s$ as the rate at which particles of species $s$ enter or leave the system and $\mathcal{K}_s \equiv m_s n_s V_s^2/2$ is the fluid kinetic energy density of species $s$ . In (2.2) and (2.3), the notation $\boldsymbol{P} : \boldsymbol{\nabla }\boldsymbol{V} \equiv P^{ij} \boldsymbol{\nabla} _j V_{i}$ is introduced to represent the Frobenius product of the two tensors and the divergence of a tensor $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s \equiv \boldsymbol{\nabla} _i P^{ij}$ is defined as contracting over neighbouring indices. We express the application of a tensor on a vector with a similar abuse of notation and write $\boldsymbol{P} \boldsymbol{\cdot }\boldsymbol{V} \equiv P_{ij} V^j$ , with the tensorial property of the resulting quantity discernible from the context. For the sake of generality, we do not focus on the closure of the moment hierarchy in this work. We assume that an adequate closure is encoded in the heat-flux densities and stress tensors, as well as the contributions to the source terms arising from interparticle interactions. We note that the particle sources $\mathcal{S}_s$ cannot all be chosen independently, as charge neutrality imposes the constraint $\sum _s q_s \mathcal{S}_s = 0$ .

Since the turbulent phenomena we wish to model tend to evolve on scales much larger than the Debye length, $\lambda _D$ , and time scales much slower than the plasma frequency, $\omega _p$ , we assume

(2.4) \begin{align} \lambda _D \boldsymbol{\nabla }\sim \frac {1}{\omega _p}\frac {\partial }{\partial t} \equiv \epsilon _D \to 0 . \end{align}

In this limit, the plasma is instantaneously locally quasi-neutral and the electric field $\boldsymbol{E}$ is not determined by the Gauss law but rather by the quasi-neutrality condition $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ . Thus, Maxwell’s equations in the limit $\epsilon _D \to 0$ reduce to

(2.5) \begin{align} \boldsymbol{\nabla }\times \boldsymbol{B} &= \mu _0 \boldsymbol{J} = \mu _0 \sum _s q_s n_s \boldsymbol{V}_s ,\\[-10pt]\nonumber \end{align}

(2.6) \begin{align} \boldsymbol{\nabla }\times \boldsymbol{E} &= - \frac {\partial \boldsymbol{B}}{\partial t} ,\\[-10pt]\nonumber \end{align}

(2.7) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B} &= 0 ,\\[-10pt]\nonumber \end{align}

(2.8) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} &= 0 . \end{align}

Maxwell’s equations, coupled to the fluid equations, (2.1)–(2.3), yield a closed model which conserves energy, momentum, mass and charge.

The energy density $\mathcal{H}$ can be expressed as

(2.9) \begin{align} \mathcal{H} = \mathcal{H}_B + \sum _s \mathcal{H}_s = \frac {B^2}{2 \mu _0} + \sum _s \left (\mathcal{K}_s + \frac {3}{2}p_s\right )\! , \end{align}

where $\mathcal{H}_s = \mathcal{K}_s + \mathcal{U}_s$ is the energy of each species, composed of the kinetic energy, $\mathcal{K}_s$ , and the internal energy, $\mathcal{U}_s \equiv 3 p_s/2$ , while $\mathcal{H}_B \equiv B^2/2 \mu _0$ is the field energy, which contains only the magnetic energy density contribution in the quasi-neutral limit. To demonstrate the conservation of $\mathcal{H}$ , we take the scalar product of the momentum equation (2.2) with $\boldsymbol{V}_s$ and use the continuity equation (2.1). This leads to a transport equation for the fluid kinetic energy density $\mathcal{K}_s$ :

(2.10) \begin{align} \frac {\partial \mathcal{K}_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\mathcal{K}_s\boldsymbol{V}_s \right ) &= - \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s + q_s n_s \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{E} + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s - \mathcal{K}_s\frac {\mathcal{S}_s}{n_s} . \end{align}

Furthermore, we note that the electromagnetic field energy evolves according to Poynting’s theorem and, in a quasi-neutral magnetised plasma, where $\omega /\omega _p \to 0$ , the term due to the electric field energy, $\partial _t (\epsilon _0 E^2/2)$ , is negligible. Indeed, the displacement current scales as $\partial _t \boldsymbol{E}/c^2 \sim (\omega \varOmega _{ce}/\omega _{pe}^2) \mu _0 \boldsymbol{J}$ , where the electron time scales are the dominant ones, so in the neutral limit $\omega /\omega _{pe} \to 0$ , we have that $\partial _t \boldsymbol{E}/c^2 \to 0$ so long as $\varOmega _{ce}/\omega _{pe}$ does not diverge (which almost always holds in magnetised plasmas). In this limit, Poynting’s theorem is given by

(2.11) \begin{align} \frac {\partial \mathcal{H}_B}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} + \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{J} &= 0, \end{align}

with $\boldsymbol{S} \equiv \boldsymbol{E} \times \boldsymbol{B}/ \mu _0$ the Poynting vector, as can be observed by taking the scalar product of Faraday’s law with $\boldsymbol{B}$ and using the fact that $\boldsymbol{B} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{E}) = \mu _0 \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} + \boldsymbol{E} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{B})$ , leading to

(2.12) \begin{align} \frac {\partial }{\partial t}\left (\frac {B^2}{2 \mu _0}\right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} + \frac {1}{\mu _0}\boldsymbol{E} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{B}) = 0, \end{align}

and then imposing $\boldsymbol{\nabla }\times \boldsymbol{B} = \mu _0 \boldsymbol{J}$ . Thus, summing the transport equations for the internal energy, (2.3), the fluid kinetic energy, (2.10), and the magnetic field energy, (2.11), results in an equation for the total energy density:

(2.13) \begin{align} \frac {\partial \mathcal{H}}{\partial t} &+ \boldsymbol{\nabla }\boldsymbol{\cdot }\Big[\boldsymbol{S} + \sum _s (\mathcal{H}_s \boldsymbol{V}_s + \boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{V}_s + \boldsymbol{q}_s) \Big] = \sum _s \left (\mathcal{Q}_s + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s \right )\! , \end{align}

where a further summation over all the species $s$ has been performed and the symmetry property of the pressure tensor $\boldsymbol{P}_s = \boldsymbol{P}^T_s$ is used to write $\boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s + \boldsymbol{P}_s : \boldsymbol{\nabla }\boldsymbol{V}_s = \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{V}_s)$ in (2.13). The energy content of the system changes due to heat sources, $\mathcal{Q}_s$ , and forces, $\boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s$ . Defining the total power source $\mathcal{S}_{\mathcal{H}} = \sum _s(\mathcal{Q}_s + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s)$ and the energy flux $\boldsymbol{\varGamma }_{\mathcal{H}} \equiv \boldsymbol{S} + \sum _s (\mathcal{H}_s \boldsymbol{V}_s + \boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{V}_s + \boldsymbol{q}_s)$ , (2.13) can be written as a conservation law for $\mathcal{H}$ :

(2.14) \begin{align} \frac {\partial \mathcal{H}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varGamma }_{\mathcal{H}} = \mathcal{S}_{\mathcal{H}}. \end{align}

The source term $\mathcal{S}_{\mathcal{H}}$ vanishes if the system is not externally driven and particle interactions are elastic, as is the case for the Coulomb interaction (Braginskii Reference Braginskii1965). A transport equation can be similarly derived for the momentum density $\boldsymbol{\mathcal{M}} \equiv \sum _s \boldsymbol{\mathcal{M}}_s$ :

(2.15) \begin{align} \frac {\partial \boldsymbol{\mathcal{ M}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\sum _s (\boldsymbol{V}_s \boldsymbol{\mathcal{ M}}_s + \boldsymbol{P}_s) &= \boldsymbol{J} \times \boldsymbol{B} + \sum _s \boldsymbol{\mathcal{R}}_s , \end{align}

where, in the presence of elastic interactions and in the absence of external momentum input, the term $\sum _s \boldsymbol{\mathcal{R}}_s$ vanishes. Total mass is conserved by virtue of the continuity equation and, finally, charge is conserved by the requirement that $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ .

3. Drift approximation

We focus on low-frequency dynamics, therefore assuming that the physics of interest evolves on a scale of the order of the drift scale, $\partial _t \sim \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\sim \omega _{* s} \equiv \boldsymbol{k}_\perp \boldsymbol{\cdot }\boldsymbol{v}_{* s}$ , where $\boldsymbol{v}_{* s} = \boldsymbol{B} \times \boldsymbol{\nabla }p_s/(q_s n_s B^2)$ is the diamagnetic drift velocity of species $s$ and $\boldsymbol{k}_\perp$ is the wavenumber in the plane orthogonal to $\boldsymbol{B}$ (Zeiler et al. Reference Zeiler, Drake and Rogers1997). The drift scale is assumed much slower than the cyclotron frequency $\omega _{* s} \ll \varOmega _{cs}$ , thus $\partial _t/\varOmega _{cs} \sim \boldsymbol{ V}_s \boldsymbol{\cdot }\boldsymbol{\nabla }/\varOmega _{cs} \equiv \epsilon$ , with $\epsilon \ll 1$ a small parameter. This implies that the fluid velocities and length scales are ordered as $V_{\perp s} k_\perp \sim V_{\parallel s} k_\parallel \sim \epsilon \varOmega _{cs}$ , where $k_\perp = |\boldsymbol{k}_\perp |$ and $k_\parallel = \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{b}$ denote the perpendicular and parallel wavenumbers, and $V_{\perp s}$ and $V_{\parallel s}$ represent the corresponding perpendicular and parallel components of the fluid velocity. Our ordering choice is consistent with both long-wavelength/large-flow drift-reduced assumptions, $k_\perp \rho _{Ls} \sim \sqrt {\epsilon },\ V_{\perp s} \sim \sqrt {\epsilon } v_{Ts}$ , and short-wavelength/small-flow gyrokinetics, $k_\perp \rho _{Ls} \sim 1,\ V_{\perp s} \sim \epsilon v_{Ts}$ , where $v_{Ts} = \sqrt {T_s/m_s}$ is the thermal velocity of species $s$ and $\rho _{Ls} = v_{Ts}/\varOmega _{cs}$ its Larmor radius (Brizard & Hahm Reference Brizard and Hahm2007). Indeed, here we adopt a more general ordering than the one typically assumed in drift-reduced models, in which only large-scale dynamics are retained $k_\perp \rho _{Ls} \sim \sqrt {\epsilon }$ (Drake & Antonsen Reference Drake and Antonsen1984; Zeiler et al. Reference Zeiler, Drake and Rogers1997; Simakov & Catto Reference Simakov and Catto2003), as the large-scale assumption is not required to derive the results presented in this work.

In the $\epsilon = 0$ limit, key physical processes that give rise to electrostatic $E \times B$ turbulence are absent (Zeiler et al. Reference Zeiler, Drake and Rogers1997). In order to retain them in our model, we expand the fluid equations, (2.1)–(2.3), in powers of $\epsilon$ (Drake & Antonsen Reference Drake and Antonsen1984). Decomposing the fluid velocity $\boldsymbol{V}_s$ of species $s$ in terms of its parallel and perpendicular components with respect to the magnetic field $\boldsymbol{B}$ ,

(3.1) \begin{align} \boldsymbol{V}_s = \boldsymbol{v}_{\parallel s} + \boldsymbol{v}_{\perp s} , \end{align}

and taking the cross-product of the momentum equation (2.2) with $\boldsymbol{B}$ , we find that $\boldsymbol{v}_{\perp s}$ decomposes into the sum of drift contributions (Drake & Antonsen Reference Drake and Antonsen1984):

(3.2) \begin{align} \boldsymbol{v}_{\perp s} = \boldsymbol{v}_{E} + \boldsymbol{v}_{* s} + \boldsymbol{v}_{\boldsymbol{\pi }s} + \boldsymbol{v}_{\boldsymbol{\mathcal{R}} s} + \boldsymbol{v}_{p s}, \end{align}

that is, an $E\times B$ drift term, $\boldsymbol{v}_{E} = \boldsymbol{E} \times \boldsymbol{B}/B^2$ , a diamagnetic drift related to isotropic pressure, $\boldsymbol{v}_{* s}$ , a diamagnetic drift due to pressure anisotropies, $\boldsymbol{v}_{\boldsymbol{\pi }s} = \boldsymbol{B} \times \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\pi }_s/(q_s n_s B^2)$ , a drift due to momentum source/sinks, $\boldsymbol{v}_{\boldsymbol{\mathcal{R}} s} = \boldsymbol{\mathcal{R}}_s \times \boldsymbol{B}/(q_s n_s B^2)$ , and a drift related to the rate of change of the bulk flow, the polarisation velocity $\boldsymbol{v}_{p s}$ , here defined as

(3.3) \begin{align} \boldsymbol{v}_{p s} \equiv \frac {\boldsymbol{b}}{n_s \varOmega _{cs}} \times \left [\frac {\partial (n_s\boldsymbol{V}_s)}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(n_s \boldsymbol{V}_s \boldsymbol{V}_s) \right ]\! , \end{align}

where $\boldsymbol{b} = \boldsymbol{B}/B$ is the unit vector in the direction of $\boldsymbol{B}$ . Using the continuity equation (2.1), the polarisation velocity can also be expressed as

(3.4) \begin{align} \boldsymbol{v}_{p s} = \frac {\boldsymbol{b}}{ \varOmega _{cs}} \times \left (\frac {\text{d}_s \boldsymbol{V}_s}{\text{d}t} + r_s \boldsymbol{V}_s\right )\!, \end{align}

where ${\rm d}_s/{\rm d}t \equiv \partial _t + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\nabla}$ is the material derivative of species $s$ . The polarisation drift results from the particle’s finite inertia and is subleading in the imposed drift ordering, that is, $\boldsymbol{v}_{p s} \sim (d_t/\varOmega _{cs}) \boldsymbol{V}_s = O(\epsilon \boldsymbol{V}_s)$ , if the rate of particle injection is of the order of the drift frequency $r_s/\varOmega _{cs} \sim \epsilon$ , which we assume to hold. We can therefore write the total fluid velocity,

(3.5) \begin{align} \boldsymbol{V}_s = \boldsymbol{\bar v}_s + \boldsymbol{v}_{p s} = (\boldsymbol{v}_{\parallel s} + \boldsymbol{\bar v}_{\perp s}) + \boldsymbol{v}_{p s}, \end{align}

as the sum of a leading-order flow, $\boldsymbol{\bar v}_s$ , and the subleading correction, $\boldsymbol{v}_{p s}$ , with $\boldsymbol{v}_{\parallel s}, \boldsymbol{\bar v}_s$ and $ \boldsymbol{\bar v}_{\perp s}$ scaling as $O(\boldsymbol{V}_s)$ . Explicitly, the leading-order perpendicular flow is given by

(3.6) \begin{align} \boldsymbol{\bar v}_{\perp s} = \boldsymbol{v}_{E} + \boldsymbol{v}_{* s} + \boldsymbol{v}_{\boldsymbol{\pi }s} + \boldsymbol{v}_{\boldsymbol{\mathcal{R}} s} \equiv \boldsymbol{v}_{E} + \boldsymbol{w}_{\perp s}, \end{align}

where $\boldsymbol{w}_{\perp s}$ is the leading-order drift velocity of species $s$ in the co-moving $E \times B$ frame. Separating the fluid velocity $\boldsymbol{V}_s$ appearing in (3.4) into its leading and subleading components, one obtains

(3.7) \begin{align} \boldsymbol{v}_{p s} & = \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \bigg(\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s + \boldsymbol{v}_{p s} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + \frac {\partial \boldsymbol{v}_{p s}}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s}\nonumber\\& \quad + r_s \boldsymbol{v}_{p s} + \boldsymbol{v}_{p s} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s} \bigg). \end{align}

Equation (3.7) yields an implicit equation for $\boldsymbol{v}_{p s}$ , which is

(3.8) \begin{equation} \boldsymbol{v}_{p s} = F(\boldsymbol{\bar v}_s) + G(\boldsymbol{\bar v}_s, \boldsymbol{v}_{p s}) + K(\boldsymbol{v}_{p s}), \end{equation}

where

(3.9) \begin{align} F(\boldsymbol{\bar v}_s) &= \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s \right ) = O(\epsilon \boldsymbol{\bar v}_s),\\[-10pt]\nonumber \end{align}

(3.10) \begin{align} G(\boldsymbol{\bar v}_s, \boldsymbol{v}_{p s}) &= \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\boldsymbol{v}_{p s} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + \frac {\partial \boldsymbol{v}_{p s}}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s} + r_s \boldsymbol{v}_{p s} \right ) = O(\epsilon \boldsymbol{v}_{p s} \boldsymbol{\bar v}_s),\\[-10pt]\nonumber \end{align}

(3.11) \begin{align} K(\boldsymbol{v}_{p s}) &= \frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{v}_{p s} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s}) = O(\epsilon \boldsymbol{v}_{p s} \boldsymbol{v}_{p s}). \end{align}

Traditionally, (3.8) is treated with ordinary perturbation theory: $\boldsymbol{v}_{p s} = \boldsymbol{v}_{p s}^{(1)} + \boldsymbol{v}_{p s}^{(2)} + \cdots$ , with $\boldsymbol{v}_{p s}^{(N)} \equiv O(\epsilon ^N \boldsymbol{\bar v}_s)$ . For example, solving (3.8) order by order, up to $O(\epsilon ^3 \boldsymbol{\bar v}_s)$ , yields

(3.12) \begin{align} \boldsymbol{v}_{p s}^{(1)} &= F(\boldsymbol{\bar v}_s) = \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s\boldsymbol{\bar v}_s\right )\!,\\[-10pt]\nonumber \end{align}

(3.13) \begin{align} \boldsymbol{v}_{p s}^{(2)} &= G(\boldsymbol{\bar v}_s, \boldsymbol{v}_{p s}^{(1)}) = \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\boldsymbol{v}_{p s}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + \frac {\partial \boldsymbol{v}_{p s}^{(1)}}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s}^{(1)} + r_s \boldsymbol{v}_{p s}^{(1)} \right )\! ,\\[-10pt]\nonumber \end{align}

(3.14) \begin{align} \boldsymbol{v}_{p s}^{(3)} &= G\big(\boldsymbol{\bar v}_s, \boldsymbol{v}_{p s}^{(2)}\big) + K\big(\boldsymbol{v}_{p s}^{(1)}\big) = G\big(\boldsymbol{\bar v}_s, \boldsymbol{v}_{p s}^{(2)}\big) + \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \big(\boldsymbol{v}_{p s}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_{p s}^{(1)}\big). \end{align}

The resulting drift-reduced model is presented in detail in Gath & Wiesenberger (Reference Gath and Wiesenberger2019). Truncating the expansion (3.12) to $O(\epsilon \boldsymbol{\bar v}_s)$ yields the models currently implemented in drift-reduced fluid codes (Dudson et al. Reference Dudson2015; Stegmeir et al. Reference Stegmeir2019; Giacomin et al. Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ; Düll et al. Reference Düll, Bufferand, Serre, Ciraolo, Quadri, Rivals and Tamain2024).

It is a known result that the models, based on the perturbative expansion in (3.12), are not conservative for both energy and momentum, with spurious source terms of $O(\epsilon )$ appearing in their time evolution (Reiser Reference Reiser2012; Halpern et al. Reference Halpern, Waltz and Bernard2023). The reason is that the expression for $\boldsymbol{v}_{p s}$ in (3.12) encodes the transport of the leading-order component of perpendicular momentum, which lacks advection by the polarisation drift (a subleading term; see (3.13)). On the other hand, the polarisation advection term is retained in other parts of the model, for instance in the continuity equation, which is expanded to $O(\epsilon )$ . The resulting inconsistency between the transport of momentum and density breaks the energy and momentum conservation properties of the model (Reiser Reference Reiser2012). We further note that, given the recursive nature of the perturbative expansion method, exact conservation is not retrieved by expanding the model to higher order in $\epsilon$ . Indeed, for the drift-reduced model expanded to $O(\epsilon ^N)$ with $N \geqslant 1$ , energy and momentum conservation is broken by terms of the same order, namely $O(\epsilon ^{N})$ .

The solution to the conservation problem of the drift approximation has been known since the original derivation of the model (Drake & Antonsen Reference Drake and Antonsen1984; Zeiler et al. Reference Zeiler, Drake and Rogers1997; Scott Reference Scott2003) and consists of retaining the polarisation advection term $\boldsymbol{b} \times (\boldsymbol{v}_{p s}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s)/\varOmega _{cs}$ in the perpendicular momentum transport equation (3.12). More precisely, the polarisation drift $\boldsymbol{v}_{p s}^{(1)}$ should be defined implicitly as

(3.15) \begin{align} \boldsymbol{v}_{p s}^{(1)} = \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + \boldsymbol{v}_{p s}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s \right ) \equiv \frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\frac {\text{d}_s^{(1)} \boldsymbol{\bar v}_s}{\text{d}t} + r_s \boldsymbol{\bar v}_s \right )\!, \end{align}

instead of the perturbative definition in (3.12). The expression for $\boldsymbol{v}_{p s}$ in (3.15) results in conservation laws for the leading-order contributions to energy and momentum, valid to all orders in $\epsilon$ . In the following section, we first invert (3.15) to express $\boldsymbol{v}_{p s}^{(1)}$ as a function of only the leading-order flow $\boldsymbol{\bar v}_s$ (§ 4). This was the main obstacle to overcome in obtaining an explicit and exactly conservative drift-reduced model. We then derive in § 5 a conservative drift-reduced model. Finally, mirroring previous work on the topic (Drake & Antonsen Reference Drake and Antonsen1984; Zeiler et al. Reference Zeiler, Drake and Rogers1997; Reiser Reference Reiser2012), we demonstrate the conservation properties of the obtained drift-reduced system in § 6. Since we consider the expansion of the model only to first order in $\epsilon$ , we suppress henceforth the superscript in $\boldsymbol{v}_{p s}^{(1)}$ to lighten the notation in the rest of the paper.

4. Solution of the implicit equation for $\boldsymbol{v}_{p s}$

To the best of our knowledge, the expression of $\boldsymbol{v}_{p s}$ , given by the solution of (3.15), has never been used in simulation codes, except for the case of an electrostatic cold-ion plasma in linear geometry (Reiser Reference Reiser2012). In that case, the symmetry properties of the magnetic field were exploited to obtain a vorticity equation and close the system equations. In general, however, the lack of an explicit form for the polarisation velocity as a function of the electric potential’s rate of change precludes the direct use of (3.15). In fact, in the absence of an explicit relation $\boldsymbol{v}_{p s} = \boldsymbol{v}_{p s}(\partial _t \phi )$ , imposing quasi-neutrality through $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ does not yield a scalar equation for $\phi$ . As a result, closing the drift-reduced fluid system requires explicit time evolution of the polarisation drift (on the cyclotron time scale), thereby defeating the purpose of the drift-reduced approximation (Halpern et al. Reference Halpern, Waltz and Bernard2023). In Scott (Reference Scott2003), the closure issue was avoided by assuming that the polarisation drift can be expressed in terms of the gradient of a scalar, but this relation does not hold in general.

In this section we seek to invert the expression for $\boldsymbol{v}_{p s}$ , given by (3.15), as a function of the leading-order velocity $\boldsymbol{\bar v}_s$ . Reordering terms in (3.15), we note that $\boldsymbol{v}_{p s}$ satisfies

(4.1) \begin{align} \boldsymbol{v}_{p s} - \frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{v}_{p s} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\bar v}_s = \frac {\boldsymbol{b}}{\varOmega _{cs}}\times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s \right )\! . \end{align}

We now show that the solution of (4.1) yields

(4.2) \begin{align}& \boldsymbol{v}_{p s} = \frac {1 + (\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s)/\varOmega _{cs} }{1 + \boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{\bar v}_s)/\varOmega _{cs} + \text{det}(\boldsymbol{\nabla }\boldsymbol{\bar v}_s)_\perp /\varOmega _{cs}^2} \boldsymbol{\cdot }\left [ \frac {\boldsymbol{b}}{\varOmega _{cs}}\times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s \right ) \right ]\!.\nonumber \\& \end{align}

For convenience, we define the inertial drift-velocity term, $\boldsymbol{U}_s$ , to represent the right-hand side of (4.1), that is,

(4.3) \begin{align} \boldsymbol{U}_s \equiv \frac {\boldsymbol{b}}{\varOmega _{cs}}\times \left (\frac {\partial \boldsymbol{\bar v}_s}{\partial t} + \boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s + r_s \boldsymbol{\bar v}_s \right )\!, \end{align}

which contains only leading-order quantities related to the velocity $\boldsymbol{\bar v}_s$ . We note that $\boldsymbol{U}_s \boldsymbol{\cdot }\boldsymbol{b} =0$ . Equation (4.1) can be written in component form as

(4.4) \begin{align} \left (\delta ^i_l - \frac {1}{\varOmega _{cs}}\epsilon ^{i j}_{\,\, \,\,k} b_j \boldsymbol{\nabla} _l \bar v_{s}^k \right )v_{ps}^l \equiv (Q^{-1}_s)^i_{\, \, l}v_{ps}^l = U^i_s, \end{align}

where $v_{ps}^l$ are the contravariant components of $\boldsymbol{v}_{p s}$ in an arbitrary coordinate system, $\boldsymbol{Q}^{-1}_s \equiv (Q^{-1}_s)^i_{\, \, l}$ is the linear operator we seek to invert, $\epsilon ^{ijk}$ is the Levi-Civita tensor in $\mathbb{R}^3$ and $\delta _{ij} = \delta ^{ij} = \delta ^i_j$ is the Kronecker delta, used to raise/lower indices in $\mathbb{R}^3$ . Einstein’s convention of summation over repeated indices is used throughout the paper. We simplify the calculation by considering a basis composed of $\boldsymbol{\hat w}$ , an arbitrary unit vector such that $\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{b} = 0$ , $ \boldsymbol{b} \times \boldsymbol{\hat w}$ and $\boldsymbol{b}$ . To lighten the notation we define $\boldsymbol{b} \times \boldsymbol{\hat w} \equiv \boldsymbol{\hat w}^*$ , and note that $\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{\hat w}^* = 0$ , $ \hat w^{* 2} = \hat w^2 = 1$ and $(\boldsymbol{\hat w}^*)^* = \boldsymbol{b} \times (\boldsymbol{b} \times \boldsymbol{\hat w}) = - \boldsymbol{\hat w}$ . Since the polarisation velocity lies in the plane orthogonal to $\boldsymbol{b}$ , we have

(4.5) \begin{align} \boldsymbol{v}_{p s} = \alpha _s \boldsymbol{\hat w} + \beta _s \boldsymbol{\hat w}^*, \end{align}

where $\alpha _s$ and $ \beta _s$ are two scalar functions of $\boldsymbol{U}_s$ to be determined. Substituting this expression for $\boldsymbol{v}_{p}$ in (4.1), we have

(4.6) \begin{align} \alpha _s \boldsymbol{\hat w} + \beta _s \boldsymbol{\hat w}^* - \alpha _s \frac {\boldsymbol{b} }{\varOmega _{cs}}\times (\boldsymbol{\hat w}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) - \beta _s \frac {\boldsymbol{b}}{\varOmega _{cs}}\times (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) = \boldsymbol{U}_s. \end{align}

Projecting (4.6) onto the basis elements $\boldsymbol{\hat w}$ and $ \boldsymbol{\hat w}^*$ , we obtain

(4.7) \begin{align} \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s &= \alpha _s \left (1- \boldsymbol{\hat w} \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v})\right ]\right ) - \beta _s \boldsymbol{\hat w} \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \left (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s\right )\right ]\! ,\\[-10pt]\nonumber \end{align}

(4.8) \begin{align} \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s &= \beta _s\left (1 -\boldsymbol{\hat w}^* \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s)\right ]\right ) - \alpha _s \boldsymbol{\hat w}^* \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s)\right ]\! . \end{align}

Equations (4.7) and (4.8) constitute a linear system for the unknowns $\alpha _s$ and $\beta _s$ . Solving for $\alpha _s$ and $\beta _s$ yields

(4.9) \begin{align} \alpha _s &= \frac {(1-\tilde \delta _s) \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s + \gamma _s \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s}{\varDelta _s},\\[-10pt]\nonumber \end{align}

(4.10) \begin{align} \beta _s &= \frac {\tilde \gamma _s \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s +(1 - \delta _s) \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s}{\varDelta _s} , \end{align}

where

(4.11) \begin{align} \delta _s &= \boldsymbol{\hat w} \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) \right ]\! ,\\[-7pt]\nonumber \end{align}

(4.12) \begin{align} \tilde \delta _s &= \boldsymbol{\hat w}^* \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) \right ]\! ,\\[-7pt]\nonumber \end{align}

(4.13) \begin{align} \gamma _s &= \boldsymbol{\hat w} \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) \right ]\! ,\\[-7pt]\nonumber \end{align}

(4.14) \begin{align} \tilde \gamma _s &= \boldsymbol{\hat w}^* \boldsymbol{\cdot }\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\bar v}_s) \right ]\! , \end{align}

and the determinant $\varDelta _s$ is given by

(4.15) \begin{align} \varDelta _s = 1 - (\delta _s + \tilde \delta _s) +(\delta _s \tilde \delta _s - \gamma _s \tilde \gamma _s) . \end{align}

Using the cyclical property of the triple product $\boldsymbol{a}_1 \boldsymbol{\cdot }(\boldsymbol{a}_2 \times \boldsymbol{a}_3) = \boldsymbol{a}_2 \boldsymbol{\cdot }(\boldsymbol{a}_3 \times \boldsymbol{a}_1)$ and $(\boldsymbol{\hat w}^*)^* = - \boldsymbol{\hat w}$ , the coefficients in (4.11)–(4.14) can be written in more transparent form as

(4.16) \begin{align} \delta _s &= - \boldsymbol{\hat w}^* \boldsymbol{\cdot }\left (\frac {\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s}{\varOmega _{cs}} \right ) \boldsymbol{\cdot }\boldsymbol{\hat w}^* ,\\[-10pt]\nonumber \end{align}

(4.17) \begin{align} \tilde \delta _s &= - \boldsymbol{\hat w} \boldsymbol{\cdot }\left (\frac {\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s}{\varOmega _{cs}} \right ) \boldsymbol{\cdot }\boldsymbol{\hat w} ,\\[-10pt]\nonumber \end{align}

(4.18) \begin{align} \gamma _s &= \boldsymbol{\hat w} \boldsymbol{\cdot }\left (\frac {\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s}{\varOmega _{cs}} \right ) \boldsymbol{\cdot }\boldsymbol{\hat w}^* ,\\[-10pt]\nonumber \end{align}

(4.19) \begin{align} \tilde \gamma _s &= \boldsymbol{\hat w}^* \boldsymbol{\cdot }\left (\frac {\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s}{\varOmega _{cs}} \right ) \boldsymbol{\cdot }\boldsymbol{\hat w} . \end{align}

Substituting the above expressions for the coefficients $\delta _s, \tilde \delta _s, \gamma _s$ and $ \tilde \gamma _s$ into (4.9), (4.10) and (4.15), we obtain

(4.20) \begin{align} \alpha _s &= \frac {1}{\varDelta _s}\left ( \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s + \left (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}\right ) \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s + \left (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}^*\right ) \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s \right )\! ,\\[-10pt]\nonumber \end{align}

(4.21) \begin{align} \beta _s &= \frac {1}{\varDelta _s}\left ( \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s + \left (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}^*\right ) \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s + \left (\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}\right ) \boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s\right )\! ,\\[-10pt]\nonumber \end{align}

(4.22) \begin{align} \varDelta _s &=1 + (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w} + \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}^*) \nonumber \\[4pt] &\quad + [(\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w})( \boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}^*) - (\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}^*)(\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{\hat w}) ] , \end{align}

with $\boldsymbol{T}_s \equiv (\boldsymbol{b}/\varOmega _{cs} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s )_\perp : W_\perp \to W_\perp$ the linear operator $\boldsymbol{b}/\varOmega _{cs} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s$ restricted to act on the subspace perpendicular to $\boldsymbol{b}(\boldsymbol{x},t)$ at position $\boldsymbol{x}$ and time $t$ , defined by $W_\perp \equiv \{\boldsymbol{y} \in \mathbb{R}^3 | \boldsymbol{b}(\boldsymbol{x},t) \boldsymbol{\cdot }\boldsymbol{y} = 0\}$ . Substituting (4.20) and (4.21) into (4.5), and recalling that $\boldsymbol{U}_s = \boldsymbol{\hat w}(\boldsymbol{\hat w} \boldsymbol{\cdot }\boldsymbol{U}_s) + \boldsymbol{w}^*(\boldsymbol{\hat w}^* \boldsymbol{\cdot }\boldsymbol{U}_s)$ , we obtain an expression for $\boldsymbol{v}_{p s}$ :

(4.23) \begin{align} \boldsymbol{v}_{p s} = \frac {1 }{\varDelta _s}(\boldsymbol{U}_s + [\boldsymbol{e}^i (\boldsymbol{e}_i \boldsymbol{\cdot }\boldsymbol{T}_s \boldsymbol{\cdot }\boldsymbol{e}^j) \boldsymbol{e}_j] \boldsymbol{\cdot }\boldsymbol{U}_s) , \end{align}

with $\boldsymbol{e}^i = \{\boldsymbol{\hat w}, \boldsymbol{\hat w}^* \}$ the orthonormal basis elements of $W_\perp$ . Since $\boldsymbol{\hat w}$ has arbitrary direction in the plane perpendicular to $\boldsymbol{b}$ , we can express the result in a basis-independent form as

(4.24) \begin{align} \boldsymbol{v}_{p s} = \frac {1 +\boldsymbol{T}_s }{\varDelta _s} \boldsymbol{\cdot }\boldsymbol{U}_s = \frac {1 +(\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s)_\perp / \varOmega _{cs}}{\varDelta _s} \boldsymbol{\cdot }\boldsymbol{U}_s . \end{align}

We note that $\boldsymbol{U}_s$ in (4.3) is an element of the orthogonal subspace $W_\perp$ and, therefore, $(\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s)_\perp \boldsymbol{\cdot }\boldsymbol{U}_s = (\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s) \boldsymbol{\cdot }\boldsymbol{U}_s$ . The determinant $\varDelta _s$ in basis-independent form is constructed from (4.22) and can be written as

(4.25) \begin{align} \varDelta _s = 1 + \text{tr}(\boldsymbol{T}_s) + \text{det}(\boldsymbol{T}_s) , \end{align}

or, explicitly, using the definition of $\boldsymbol{T}_s$ , as

(4.26) \begin{align} \varDelta _s = 1 + \text{tr}\left ( \left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s\right ]_\perp \right ) + \text{det}\left (\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s\right ]_\perp \right )\!. \end{align}

The trace term is given by

(4.27) \begin{align} \text{tr}\left ( \left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s\right ]_\perp \right ) = \frac {1}{\varOmega _{cs}}\epsilon _i^{\,\, jk}b_k \boldsymbol{\nabla} _j \bar v_s^i = \frac {1}{\varOmega _{cs}}\boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{\bar v}_s). \end{align}

Meanwhile, for any linear application $A : \mathbb{R}^3 \to \mathbb{R}^3$ , the determinant of its restriction to the $W_\perp$ subspace can be expressed as $\text{det}(A_\perp ) = \epsilon _{\perp }^{kl} \epsilon _{\perp \,\, ij} A^i_{\,\, k} A^{j}_{\,\,l}/2$ , where $\epsilon _\perp ^{ij} \equiv \epsilon ^{ijk}b_k$ denotes the induced Levi-Civita symbol on $W_\perp$ . The determinant term in (4.26) becomes

(4.28) \begin{align} \text{det}\left (\left [\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s\right ]_\perp \right ) &= \frac {1}{2} \epsilon _{\perp }^{kl}\epsilon _{\perp ij}\left (\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s \right )^i_{\,\, k} \left (\frac {\boldsymbol{b}}{\varOmega _{cs}} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s \right )^j_{\,\, l} \nonumber \\[5pt] &= \frac {1}{2 \varOmega _{cs}^2} \epsilon _{\perp }^{kl} \epsilon _{\perp \, i j} (\boldsymbol{\nabla }\boldsymbol{\bar v}_s)^i_{\,\, k} (\boldsymbol{\nabla }\boldsymbol{\bar v}_s)^j_{\,\, l} = \frac {1}{\varOmega _{cs}^2}\text{det} ( \boldsymbol{\nabla }\boldsymbol{\bar v}_s )_\perp , \end{align}

where we used the fact that $(\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s)^i_{\,\, k} = \epsilon _\perp ^{ij} \boldsymbol{\nabla} _j \bar v_k$ and $\epsilon _{\perp }^{ki} \epsilon _{\perp \, kj} = (\varPi _\perp )^i_{\,\, j}$ . Finally, the determinant $\varDelta _s$ in (4.25) becomes

(4.29) \begin{equation} \varDelta _s = 1 + \frac {\boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{\bar v}_s)}{\varOmega _{cs}} + \frac {\text{det} ( \boldsymbol{\nabla }\boldsymbol{\bar v}_s )_\perp }{\varOmega _{cs}^2} , \end{equation}

which contains a correction due to the parallel component of the vorticity $\boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{\bar v}_s)/\varOmega _{cs} \sim O(\epsilon )$ and a nonlinear correction related to flow shear $\text{det} ( \boldsymbol{\nabla }\boldsymbol{\bar v}_s )_\perp /\varOmega _{cs}^2 \sim O(\epsilon ^2)$ . For a linear geometry in the cold-ion approximation, and for the case where the leading-order ion perpendicular flow is purely the $E \times B$ drift $\boldsymbol{\bar v}_{\perp } = \boldsymbol{v}_{E}$ , (4.29) reduces to equation (A6) in Reiser (Reference Reiser2012). The final solution for the polarisation velocity is therefore given by (3.15), or, in more succinct form, by

(4.30) \begin{align} \boldsymbol{v}_{p s} = \boldsymbol{Q}_s(\boldsymbol{\bar v}_s) \boldsymbol{\cdot }\boldsymbol{U}_s, \end{align}

where

(4.31) \begin{align} \boldsymbol{Q}_s(\boldsymbol{\bar v}_s) = \frac {1 + (\boldsymbol{b} \times \boldsymbol{\nabla }\boldsymbol{\bar v}_s)/\varOmega _{cs} }{1 + \boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{\bar v}_s)/\varOmega _{cs} + \text{det}(\boldsymbol{\nabla }\boldsymbol{\bar v}_s)_\perp /\varOmega _{cs}^2}. \end{align}

Within the drift-ordering $\boldsymbol{\bar v}_s \boldsymbol{\cdot }\boldsymbol{\nabla }/\varOmega _{cs} \sim \epsilon$ , the application $\boldsymbol{Q}_s^{-1}$ is invertible, as the determinant $\varDelta _s = 1 + O(\epsilon )$ does not vanish. In Appendix A, we verify that the expression of $\boldsymbol{v}_{p s}$ given in (4.30) satisfies (4.1).

We note that (4.30) represents a non-perturbative result, containing terms of all orders in $\epsilon$ , as can be seen by Taylor-expanding the denominator in (4.31). All of the subleading terms contained in $\boldsymbol{Q}_s$ are required for the existence of exact conservation laws.

5. Conservative drift-reduced model

Based on (3.15) we construct now the drift-reduced model, valid to $O(\epsilon )$ , which admits energy and momentum conservation laws satisfied exactly, i.e. to all orders in $\epsilon$ . The model evolves the density, parallel velocity and pressure of each species and is closed by a vorticity equation, an associated Poisson equation for $\phi$ and Ampère’s equation for the vector potential $\boldsymbol{A}$ .

Expressing the total fluid velocity as $\boldsymbol{V}_s = \boldsymbol{\bar v}_s + \boldsymbol{v}_{p s}$ in (2.1) and (2.3), we obtain

(5.1) \begin{align} &\frac {\partial n_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }[n_s (\boldsymbol{\bar v}_s + \boldsymbol{v}_{p s})] = \mathcal{S}_s ,\\[-10pt]\nonumber \end{align}

(5.2) \begin{align} &\frac {\partial }{\partial t}\left ( \frac {3}{2} p_s \right ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [\frac {3}{2}p_s (\boldsymbol{\bar v}_s + \boldsymbol{v}_{p s}) + \boldsymbol{q}_s \right ] + \boldsymbol{P}_s : \boldsymbol{\nabla }(\boldsymbol{\bar v}_s + \boldsymbol{v}_{p s}) = \mathcal{Q}_s +r_s \mathcal{\overline { K}}_s , \end{align}

where $\boldsymbol{v}_{p s}$ is given either implicitly by (3.15) or explicitly via (4.2), and we define $\mathcal{\overline { K}}_s \equiv m_s n_s \bar v_s^2/2$ to be the leading-order component of the fluid kinetic energy. The equation for the velocity component parallel to the magnetic field $v_{\parallel s} = \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{b}$ of species $s$ is obtained by projecting the $O(\epsilon )$ expanded momentum (2.2),

(5.3) \begin{equation} \frac {\partial (m_s n_s \boldsymbol{\bar v}_s)}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(m_s n_s \boldsymbol{V}_s \boldsymbol{\bar v}_s) + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s - q_s n_s(\boldsymbol{E} + \boldsymbol{V}_s \times \boldsymbol{B}) = \boldsymbol{\mathcal{R}}_s, \end{equation}

onto $\boldsymbol{b}$ , obtaining

(5.4) \begin{equation} \frac {\partial (m_s n_s v_{\parallel s})}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(m_s n_s v_{\parallel s} \boldsymbol{V}_s) - m_s n_s \boldsymbol{\bar v}_s \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t } + (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\parallel = \mathcal{R}_{\parallel s} + q_s n_s E_{\parallel }, \end{equation}

with $(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\parallel \equiv \boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)$ and, similarly, for other quantities.

To derive the vorticity equation, we impose the quasi-neutrality constraint $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ . We begin with the force-balance equation for the leading-order component of the total momentum density, $\boldsymbol{\mathcal{\overline M}} = \sum _s m_s n_s \boldsymbol{\bar v}_s$ . As shown in § 6, this is given by

(5.5) \begin{align} \frac {\partial \boldsymbol{\mathcal{\overline M}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\sum _s ( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s + \boldsymbol{P}_s) &= \boldsymbol{J} \times \boldsymbol{B} + \sum _s \boldsymbol{\mathcal{R}}_s , \end{align}

where the total current $\boldsymbol{J} = \boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp + \boldsymbol{J}_{p}$ is the sum of the leading-order components, $\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp = \sum _s q_s n_s \boldsymbol{v}_{\parallel s} + \sum _s q_s n_s \boldsymbol{\bar v}_{\perp s}$ , and the subleading polarisation term, $\boldsymbol{J}_{p} = \sum _s q_s n_s \boldsymbol{v}_{p s}$ . We note the similarity between (5.5) and the non-drift-reduced form in (2.15). Taking the cross-product of (5.5) with $\boldsymbol{B}$ leads to

(5.6) \begin{equation} \boldsymbol{J} = \boldsymbol{J}_\parallel + \frac {\boldsymbol{b}}{B}\times \frac {\partial \boldsymbol{\mathcal{\overline M}}}{\partial t} + \frac {\boldsymbol{b}}{B}\times \left (\sum _s[\boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s + \boldsymbol{P}_s) - \boldsymbol{\mathcal{R}}_s ] \right )\!, \end{equation}

which can be rewritten as

(5.7) \begin{equation} \boldsymbol{J} = \boldsymbol{J}_\parallel + \frac {\partial }{\partial t}\left (\frac {\boldsymbol{b} \times \boldsymbol{\mathcal{\overline M}}}{B}\right ) + \frac {\boldsymbol{b}}{B}\times \left (\sum _s[\boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s + \boldsymbol{P}_s) - \boldsymbol{\mathcal{R}}_s ] \right ) - \frac {\partial }{\partial t}\left (\frac {\boldsymbol{b}}{B}\right )\times \boldsymbol{\mathcal{\overline M}}. \end{equation}

Imposing the quasi-neutrality constraint $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ , we obtain the vorticity equation for $\varpi \equiv -\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{b} \times \boldsymbol{\mathcal{\overline M}}/B)$ :

(5.8) \begin{align} \frac {\partial \varpi }{\partial t} & = \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\Bigg[\sum _s \frac {\boldsymbol{b}}{B}\times \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s) \Bigg]\nonumber\\[5pt]& \quad + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \left ( \frac {(\boldsymbol{\varPi }_\perp - \boldsymbol{b} \boldsymbol{b})\boldsymbol{\cdot }\boldsymbol{\nabla }\times \boldsymbol{E}}{B^2}\right ) \times \boldsymbol{\mathcal{\overline M}} \right ]\!, \end{align}

where we use Faraday’s law and, since $\boldsymbol{b} \boldsymbol{\cdot }\partial _t \boldsymbol{b} = 0$ , $\partial _t(\boldsymbol{b}/B) = 2 \boldsymbol{b} \boldsymbol{b} \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{E}) /B^2 - \boldsymbol{\nabla }\times \boldsymbol{E}/B^2 = (\boldsymbol{b} \boldsymbol{b} - \boldsymbol{\varPi }_\perp ) \boldsymbol{\cdot }(\boldsymbol{\nabla }\times \boldsymbol{E}) /B^2$ , with $\boldsymbol{\varPi }_\perp = \mathbb{I} - \boldsymbol{b} \boldsymbol{b}$ the perpendicular projection operator. In (5.8), we also express the pressure and momentum drive terms as $ \sum _s \boldsymbol{b} \times (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s - \boldsymbol{\mathcal{R}}_s)/B = \sum _s q_s n_s \boldsymbol{\bar v}_{\perp s} = \boldsymbol{\bar J}_\perp$ . Using $\boldsymbol{V}_s = \boldsymbol{\bar v}_s + \boldsymbol{v}_{p s}$ in (5.8), we obtain

(5.9) \begin{align} \frac {\partial \varpi }{\partial t} &= \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\Bigg[\sum _s \frac {\boldsymbol{b}}{B}\times \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\bar v}_s \boldsymbol{\mathcal{\overline M}}_s) \Bigg]\nonumber\\[6pt]& \quad + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \left ( \frac {(\boldsymbol{\varPi }_\perp - \boldsymbol{b} \boldsymbol{b})\boldsymbol{\cdot }\boldsymbol{\nabla }\times \boldsymbol{E}}{B^2}\right ) \times \boldsymbol{\mathcal{\overline M}} \right ] \nonumber \\[6pt] & \quad + \boldsymbol{\nabla }\boldsymbol{\cdot }\Bigg[\sum _s \frac {\boldsymbol{b}}{B}\times \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{Q}_s(\boldsymbol{\bar v}_s) \boldsymbol{\cdot }\boldsymbol{U}_s \boldsymbol{\mathcal{\overline M}}_s) \Bigg], \end{align}

where $\boldsymbol{v}_{p s}$ is expressed in terms of leading-order quantities via (4.30), with $\boldsymbol{U}_s$ given in (4.3). The last term in (5.9) is the conservative correction to the usual vorticity equation found in the literature, accounting for polarisation advection of the leading-order momentum $\boldsymbol{\mathcal{\overline M}}$ . Having evolved the scalar vorticity in time, the Poisson equation, defined via $\varpi = -\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{b} \times \boldsymbol{\mathcal{\overline M}}/B)$ , takes its usual form (see e.g. Giacomin et al. Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ):

(5.10) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\rho }{B^2}\boldsymbol{\nabla} _\perp \phi \right ) = \varpi - \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\sum _s \frac {(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\perp - \boldsymbol{\mathcal{R}}_{\perp s} }{B \varOmega _{cs}}\right ) - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\rho }{B^2}\left .\frac {\partial \boldsymbol{A}}{\partial t}\right |_\perp \right )\!, \end{equation}

where the terms due to $\boldsymbol{\mathcal{R}}_{\perp s}$ and $\partial _t \boldsymbol{A}|_\perp$ are usually neglected because of additional ordering assumptions on the size of perpendicular friction terms and strength of electromagnetic effects (Zeiler et al. Reference Zeiler, Drake and Rogers1997). In (5.10), $\rho \equiv \sum _s m_s n_s$ is the total mass density of the plasma and the perpendicular component of a given vector $\boldsymbol{W}$ is defined using the projection operator, $\boldsymbol{W}_\perp \equiv \boldsymbol{\varPi }_\perp \boldsymbol{\cdot }\boldsymbol{W}$ .

To summarise, the drift-reduced fluid model we propose, obeying exact conservation laws for the leading-order energy $\mathcal{\overline { H}}$ , leading-order momentum $\boldsymbol{\mathcal{\overline M}}$ , mass and charge (as shown in § 6), is given by the moment equations (5.1), (5.2) and (5.4), the vorticity equation (5.9) for $\varpi$ , the Poisson equation (5.10) for $\phi$ , the polarisation drift expression (4.2) for $\boldsymbol{v}_{p s}$ and the Ampère equation for $\boldsymbol{A}$ , which in the Coulomb gauge, $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{A} = 0$ , is given by

(5.11) \begin{align} \boldsymbol{\nabla} ^2 \boldsymbol{A} = - \mu _0 (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp + \boldsymbol{J}_{p}). \end{align}

Computing the polarisation velocity requires evaluating $\partial _t \boldsymbol{\bar v}_s$ in (4.2), which is itself a function of $\boldsymbol{v}_{p s}$ . Solution of this coupled system can be implemented either implicitly via an iterative procedure or explicitly by using $\boldsymbol{v}_{p s}$ from the previous time step to first evolve the dynamical fields, then calculating $\partial _t \boldsymbol{\bar v}_s$ to update $\boldsymbol{v}_{p s}$ . Devising a numerical implementation which preserves the conservation properties of the system is, however, left for future work. The commonly used non-conservative drift-reduced vorticity equation is obtained from (5.9) by keeping only the leading-order terms. This amounts to neglecting the last term in (5.9), yielding

(5.12) \begin{align} \frac {\partial \varpi }{\partial t} & \simeq \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [\sum _s \frac {\boldsymbol{b}}{B}\times \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{\bar v}_s \boldsymbol{\mathcal{\overline M}}_s) \right ]\nonumber\\[9pt]&\quad + \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \left ( \frac {(\boldsymbol{\varPi }_\perp - \boldsymbol{b} \boldsymbol{b})\boldsymbol{\cdot }\boldsymbol{\nabla }\times \boldsymbol{E}}{B^2}\right ) \times \boldsymbol{\mathcal{\overline M}} \right ]\!, \end{align}

which is the form found in the literature and implemented in drift-reduced turbulence codes (Gath & Wiesenberger Reference Gath and Wiesenberger2019; Stegmeir et al. Reference Stegmeir2019; Giacomin et al. Reference Giacomin, Ricci, Coroado, Fourestey, Galassi, Lanti, Mancini, Richart, Stenger and Varini2022b ; Düll et al. Reference Düll, Bufferand, Serre, Ciraolo, Quadri, Rivals and Tamain2024). For the case considered in Reiser (Reference Reiser2012), the use of a conservative formulation had a small influence on the overall turbulent dynamics. In fact, the corrections needed to ensure energy–momentum consistency are subdominant in the drift-ordering expansion. However, in regimes where flow shear becomes large, and the tensor $\boldsymbol{Q}_s$ in (4.31) deviates sufficiently from unity, a self-consistent treatment of the vorticity equation, as given in (5.9), could become necessary to accurately describe the nonlinear dynamics.

6. Conservation properties

To prove that the system of equations (5.1), (5.2), (5.4), (5.9) and (5.10) conserves both energy and momentum, we follow the procedure in Halpern et al. (Reference Halpern, Waltz and Bernard2023). To obtain the transport equation for the parallel kinetic energy density $\mathcal{K}_{\parallel s} \equiv m_s n_s v_{\parallel s}^2/2$ of species $s$ , we multiply (5.4) by $v_{\parallel s}$ . After straightforward algebra, this leads to

(6.1) \begin{align} \frac {\partial \mathcal{K}_{\parallel s}}{\partial t} & + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\mathcal{K}_{\parallel s}\boldsymbol{V}_s\right ) - m_s n_s v_{\parallel s} \boldsymbol{\bar v}_s \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} + \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s) = \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s\nonumber\\& \quad + q_s n_s \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{E} - r_s \mathcal{K}_{\parallel s}. \end{align}

We note that $\boldsymbol{b} \boldsymbol{\cdot }\text{d}_s \boldsymbol{b}/\text{d}t = \text{d}_s (b^2/2)/\text{d}t= 0$ and, as a consequence, $\boldsymbol{\bar v} \boldsymbol{\cdot }d_s \boldsymbol{b} /dt = \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }d_s\boldsymbol{b}/dt$ . Adding the internal energy $\mathcal{U}_s = 3 p_s/2$ contribution in (5.2)–(6.1), we have

(6.2) \begin{align} \frac {\partial (\mathcal{K}_{\parallel s} + \mathcal{U}_s) }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\!\left [(\mathcal{K}_{\parallel s} + \mathcal{U}_s)\boldsymbol{V}_s + \boldsymbol{q}_s\!\right ] &- m_s n_s v_{\parallel s} \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} + \boldsymbol{P}_s: \boldsymbol{\nabla }\boldsymbol{V}_s + \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\!(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s) \nonumber \\ &= \mathcal{Q}_s + \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s + q_s n_s \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{E} + r_s (\mathcal{\overline { K}}_{s} -\mathcal{K}_{\parallel s}). \end{align}

To find the transport equation for the leading-order perpendicular kinetic energy $\mathcal{\overline { K}}_{\perp s} \equiv m_s n_s \bar v_{\perp s}^2/2 = \mathcal{\overline { K}}_{s} -\mathcal{K}_{\parallel s}$ and the magnetic field energy $\mathcal{H}_B$ , we multiply the quasi-neutrality equation $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = 0$ by $\phi$ and express this as

(6.3) \begin{align} 0 = \boldsymbol{\nabla }\boldsymbol{\cdot }( \phi \boldsymbol{J} ) - (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp + \boldsymbol{J}_{p}) \boldsymbol{\cdot }\boldsymbol{\nabla }\phi , \end{align}

where we recall $\boldsymbol{J} \equiv \sum _s q_s n_s \boldsymbol{V}_s = \boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp + \boldsymbol{J}_{p}$ . The total polarisation current can be expressed as $\boldsymbol{J}_{p} \equiv \sum _s \boldsymbol{J}_{p s} = \sum _s q_s n_s \boldsymbol{v}_{p s}$ , with $\boldsymbol{v}_{p s}$ given by (3.15). The polarisation current density of species $s$ , $\boldsymbol{J}_{p s}$ , is therefore

(6.4) \begin{align} \boldsymbol{J}_{p s} = \frac {m_s n_s}{B} \boldsymbol{b}\times \left (\frac {\text{d}_s \boldsymbol{\bar v}_s}{\text{d} t} + r_s \boldsymbol{\bar v}_s \right )\!. \end{align}

As a first step, we note that, since the polarisation current is orthogonal to $\boldsymbol{b}$ , we can express

(6.5) \begin{align} \boldsymbol{J}_{p} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi = \boldsymbol{J}_{p} \boldsymbol{\cdot }\boldsymbol{\nabla} _\perp \phi = - \boldsymbol{J}_{p} \boldsymbol{\cdot }[\boldsymbol{b} \times (\boldsymbol{b} \times \boldsymbol{\nabla }\phi )], \end{align}

and, given that $ \boldsymbol{\nabla }\phi = -\boldsymbol{E} - \partial _t \boldsymbol{A}$ , we can rewrite (6.5) as

(6.6) \begin{align} \boldsymbol{J}_{p} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi = - \boldsymbol{J}_{p} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{v}_{E}) - \boldsymbol{J}_{p} \boldsymbol{\cdot }\partial _t \boldsymbol{A}. \end{align}

Moreover, since $\boldsymbol{v}_{E}$ is species-independent, we obtain

(6.7) \begin{equation} \boldsymbol{J}_{p} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi = - \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{v}_{E}) - \boldsymbol{J}_{p} \boldsymbol{\cdot }\partial _t \boldsymbol{A}. \end{equation}

Expressing the $E \times B$ drift in (6.7) as $\boldsymbol{v}_{E} = \boldsymbol{\bar v}_{\perp s} - \boldsymbol{w}_{\perp s}$ , with $\boldsymbol{w}_{\perp s}$ the drift velocity of species $s$ in the frame co-moving with the $E \times B$ velocity (cf. (3.6)), we have

(6.8) \begin{align} \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{v}_{E}) = \sum _s \left [m_s n_s \left (\frac {\text{d}_s \boldsymbol{\bar v}_s}{\text{d}t} + r_s \boldsymbol{\bar v}_s \right ) \boldsymbol{\cdot }\boldsymbol{\bar v}_{\perp s} \right ] - \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}). \end{align}

Decomposing $\boldsymbol{\bar v}_s = v_{\parallel s} \boldsymbol{b} + \boldsymbol{\bar v}_{\perp s}$ into its parallel and perpendicular components, the first term on the right-hand side of (6.8) becomes

(6.9) \begin{align} \sum _s \left [m_s n_s \left (\frac {\text{d}_s \boldsymbol{\bar v}_s}{\text{d}t} + r_s \boldsymbol{\bar v}_s \right ) \boldsymbol{\cdot }\boldsymbol{\bar v}_{\perp s} \right ] &= \sum _s m_s n_s \left (\boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{\bar v}_{\perp s}}{\text{d} t} + r_s \bar v_{\perp s}^2\right ) \nonumber \\[5pt] &+ \sum _s m_s n_s v_{\parallel s} \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} . \end{align}

The first term on the right-hand side of (6.9) leads to a transport law for the leading-order perpendicular kinetic energy of species $s$ , $ \mathcal{\overline { K}}_{\perp s}$ . Indeed, we have

(6.10) \begin{align} m_s n_s \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\left (\frac {\text{d}_s \boldsymbol{\bar v}_{\perp s}}{\text{d}t} + r_s \boldsymbol{\bar v}_{\perp s} \right ) &= \frac {\partial \mathcal{\overline { K}}_{\perp s}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\mathcal{\overline { K}}_{\perp s} \boldsymbol{V}_s) + r_s \mathcal{\overline { K}}_{\perp s}, \end{align}

where the continuity equation ${\rm d}_s n_s/ {\rm d}t = r_s n_s - n_s \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{V}_s$ is used. Using (6.7)–(6.10), we find

(6.11) \begin{align} \boldsymbol{J}_{p} \boldsymbol{\cdot }\boldsymbol{\nabla }\phi = &-\sum _s \left [\frac {\partial \mathcal{\overline { K}}_{\perp s}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\mathcal{\overline { K}}_{\perp s} \boldsymbol{V}_s) + r_s \mathcal{\overline { K}}_{\perp s} \right ] \nonumber \\[5pt] &- \sum _s m_s n_s v_{\parallel s} \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} + \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) - \boldsymbol{J}_{p} \boldsymbol{\cdot }\partial _t \boldsymbol{A}. \end{align}

Given (6.11), we can rewrite (6.3) as

(6.12) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }( \phi \boldsymbol{J} )& - (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) \boldsymbol{\cdot }\boldsymbol{\nabla }\phi - \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) + \boldsymbol{J}_{p} \boldsymbol{\cdot }\partial _t \boldsymbol{A} \nonumber \\ & +\sum _s \left [\frac {\partial \mathcal{\overline { K}}_{\perp s}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\mathcal{\overline { K}}_{\perp s} \boldsymbol{V}_s) + r_s \mathcal{\overline { K}}_{\perp s} \right ] + \sum _s m_s n_s v_{\parallel s} \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} = 0.\nonumber\\ \end{align}

The term $\boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{J})$ can be expressed in terms of the Poynting flux as

(6.13) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{J}) = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} + \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {\partial _t \boldsymbol{A} \times \boldsymbol{B}}{\mu _0} \right )\!, \end{align}

where $\boldsymbol{\nabla }\times \boldsymbol{B} = \mu _0 \boldsymbol{J}$ is used. Writing the term $\boldsymbol{\nabla }\boldsymbol{\cdot }(\partial _t \boldsymbol{A} \times \boldsymbol{B}/\mu _0)$ as

(6.14) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {\partial _t \boldsymbol{A} \times \boldsymbol{B}}{\mu _0} \right ) = \frac {\boldsymbol{B}}{\mu _0} \boldsymbol{\cdot }\partial _t \boldsymbol{\nabla }\times \boldsymbol{A} - \frac {\boldsymbol{\nabla }\times \boldsymbol{B}}{\mu _0} \boldsymbol{\cdot }\partial _t \boldsymbol{A}= \frac {\partial \mathcal{H}_B}{\partial t} - \boldsymbol{J} \boldsymbol{\cdot }\partial _t \boldsymbol{A}, \end{align}

where we recall that $\mathcal{H}_B = B^2/(2 \mu _0)$ is the magnetic energy density, we obtain

(6.15) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }(\phi \boldsymbol{J}) = \frac {\partial \mathcal{H}_B}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} - \boldsymbol{J} \boldsymbol{\cdot }\partial _t \boldsymbol{A} . \end{align}

Substituting (6.15) into (6.12), we find the transport equation for the sum of the field energy, $\mathcal{H}_B$ , and the total leading-order perpendicular kinetic energy, $\mathcal{\overline { K}}_\perp = \sum _s \mathcal{\overline { K}}_{\perp s}$ , given by

(6.16) \begin{align} & \frac {\partial (\mathcal{H}_B + \mathcal{\overline { K}}_\perp )}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\bigg(\boldsymbol{S} + \sum _s \mathcal{\overline { K}}_{\perp s} \boldsymbol{V}_s\bigg) + (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) \boldsymbol{\cdot }\boldsymbol{E} \nonumber \\ & + \sum _s r_s \mathcal{\overline { K}}_{\perp s} + \sum _s m_s n_s v_{\parallel s} \boldsymbol{\bar v}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} - \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = 0 . \end{align}

We now add the contribution from the parallel and internal energy densities in (6.2), having summed it over the species indices $s$ , to (6.16). We thus obtain that the total leading-order energy density, $\mathcal{\overline { H}} \equiv \mathcal{H}_B + \sum _s(\mathcal{\overline { K}}_s + \mathcal{U}_s)$ , evolves according to

(6.17) \begin{align} \frac {\partial \mathcal{\overline { H}} }{\partial t} &+ \boldsymbol{\nabla }\boldsymbol{\cdot }\bigg(\boldsymbol{S} + \sum _s [\mathcal{\overline { H}}_s\boldsymbol{V}_s + \boldsymbol{q}_s]\bigg) + \sum _s \boldsymbol{P}_s : \boldsymbol{\nabla }\boldsymbol{V}_s + \sum _s \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s) \nonumber \\ &+ (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) \boldsymbol{\cdot }\boldsymbol{E} - \sum _s q_s n_s \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{E} - \sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) \nonumber \\ &\qquad \qquad = \sum _s (\mathcal{Q}_s + \boldsymbol{v}_{\parallel s} \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s ) , \end{align}

where we introduce the species-specific energy density, $\mathcal{\overline { H}}_s \equiv \mathcal{\overline { K}}_s + \mathcal{U}_s$ . Finally, we rewrite the term $\sum _s \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s})$ by decomposing $\boldsymbol{v}_{p s} = \boldsymbol{v}_{\perp s} - \boldsymbol{\bar v}_{\perp s} = \boldsymbol{v}_{\perp s} - \boldsymbol{v}_{E} - \boldsymbol{w}_{\perp s}$ :

(6.18) \begin{align} &\boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = q_s n_s (\boldsymbol{v}_{\perp s} - \boldsymbol{v}_{E} - \boldsymbol{w}_{\perp s}) \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = q_s n_s (\boldsymbol{v}_{\perp s} - \boldsymbol{v}_{E}) \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}).\nonumber \\& \end{align}

Given that $\boldsymbol{v}_{E} = \boldsymbol{E} \times \boldsymbol{B}/B^2$ , one deduces that the term $ qn \boldsymbol{v}_{E} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = -q_s n_s \boldsymbol{w}_{\perp s} \boldsymbol{\cdot }\boldsymbol{E}$ . Furthermore, we have

(6.19) \begin{align} q_s n_s (\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = \boldsymbol{B} \times (\boldsymbol{v}_{* s} +\boldsymbol{v}_{\boldsymbol{\pi }s} + \boldsymbol{v}_{\boldsymbol{\mathcal{R}} s}) = - (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\perp + \boldsymbol{\mathcal{R}}_{\perp s}. \end{align}

Substituting the expression (6.19) into (6.18) we obtain

(6.20) \begin{align} \boldsymbol{J}_{p s} \boldsymbol{\cdot }(\boldsymbol{B} \times \boldsymbol{w}_{\perp s}) = - \boldsymbol{v}_{\perp s} \boldsymbol{\cdot }(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s) + \boldsymbol{v}_{\perp s} \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s + q_s n_s \boldsymbol{w}_{\perp s} \boldsymbol{\cdot }\boldsymbol{E}, \end{align}

which, substituted into (6.17), yields

(6.21) \begin{align} \frac {\partial \mathcal{\overline { H}} }{\partial t} &+ \boldsymbol{\nabla }\boldsymbol{\cdot }\bigg(\boldsymbol{S} + \sum _s [\mathcal{\overline { H}}_s \boldsymbol{V}_s + \boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{V}_s + \boldsymbol{q}_s]\bigg) \nonumber \\[4pt] &+ (\boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp ) \boldsymbol{\cdot }\boldsymbol{E} - \sum _s q_s n_s (\boldsymbol{v}_{\parallel s} + \boldsymbol{w}_{\perp s} ) \boldsymbol{\cdot }\boldsymbol{E} = \sum _s(\mathcal{Q}_s + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s ) . \end{align}

Recalling that

(6.22) \begin{equation} \boldsymbol{J}_\parallel + \boldsymbol{\bar J}_\perp = \sum _s q_s n_s (\boldsymbol{v}_{\parallel s} + \boldsymbol{w}_{\perp s}), \end{equation}

we find that the ohmic heating terms in (6.21) cancel, and the leading-order energy transport equation is given by

(6.23) \begin{align} \frac {\partial \mathcal{\overline { H}} }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\bigg(\boldsymbol{S} + \sum _s [\mathcal{\overline { H}}_s\boldsymbol{V}_s + \boldsymbol{P}_s \boldsymbol{\cdot }\boldsymbol{V}_s + \boldsymbol{q}_s]\bigg) &= \sum _s (\mathcal{Q}_s + \boldsymbol{V}_s \boldsymbol{\cdot }\boldsymbol{\mathcal{R}}_s ) . \end{align}

Comparing this result with the time evolution of the total energy density $\mathcal{H}$ , given by (2.13), we find that they have identical form. The model given by (5.1), (5.2) and (5.4), with $\boldsymbol{v}_{p s}$ given by (3.15), and coupled to the quasi-neutral Maxwell equations (2.5)–(2.8), therefore conserves the leading-order component $\mathcal{\overline { H}}$ of the total energy $\mathcal{H} = \mathcal{\overline { H}}[1 + O(\epsilon )]$ exactly, that is, to all orders in $\epsilon$ . This result can be extended to plasmas with finite Debye length, though for simplicity we have assumed $\partial _t / \omega _p = 0$ .

The leading-order component of the total fluid momentum density $\boldsymbol{\mathcal{\overline M}} \equiv \sum _s m_s n_s \boldsymbol{\bar v}_s = \boldsymbol{\mathcal{M}}_\parallel + \boldsymbol{\mathcal{\overline M}}_\perp$ is also a conserved quantity. The equation for the parallel momentum of species $s$ , $ \boldsymbol{\mathcal{M}}_{\parallel s} \equiv m_s n_s \boldsymbol{v}_{\parallel s}$ , is obtained by multiplying (5.4) by $\boldsymbol{b}$ :

(6.24) \begin{equation} \frac {\partial \boldsymbol{\mathcal{M}}_{\parallel s}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{V}_s \boldsymbol{\mathcal{M}}_{\parallel s}) - \mathcal{M}_{\parallel s} \frac {\text{d}_s \boldsymbol{b}}{\text{d}t} - \boldsymbol{\mathcal{\overline M}}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} \boldsymbol{b} + (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\parallel \boldsymbol{b} - q_s n_s \boldsymbol{E}_\parallel = \boldsymbol{\mathcal{R}}_{\parallel s}, \end{equation}

where $\boldsymbol{\mathcal{\overline M}}_{\perp s} \equiv m_s n_s \boldsymbol{\bar v}_{\perp s}$ is defined. The evolution of the perpendicular momentum $\boldsymbol{\mathcal{\overline M}}_{\perp s}$ is obtained from the force balance equation:

(6.25) \begin{equation} m_s n_s \left .\frac {\text{d}_s \boldsymbol{\bar v}_s}{\text{d}t}\right |_\perp + r_s m_s n_s \boldsymbol{\bar v}_{\perp s} - q_s n_s \boldsymbol{B} \times \boldsymbol{\bar v}_{\perp s} = q_s n_s \boldsymbol{V}_s \times \boldsymbol{B}, \end{equation}

which is a rewriting of the polarisation velocity definition (3.15). Using the continuity equation (5.1), (6.25) can be recast as

(6.26) \begin{equation} \left .\frac {\text{d}_s \boldsymbol{\mathcal{\overline M}}_s}{\text{d}t}\right |_\perp + \boldsymbol{\mathcal{\overline M}}_{\perp s} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{V}_s - q_s n_s \boldsymbol{B} \times \boldsymbol{\bar v}_{\perp s} = q_s n_s \boldsymbol{V}_s \times \boldsymbol{B}. \end{equation}

Using the fact that

(6.27) \begin{equation} \left .\frac {\text{d}_s \boldsymbol{\mathcal{\overline M}}_s}{\text{d}t}\right |_\perp = \frac {\text{d}_s \boldsymbol{\mathcal{\overline M}}_{\perp s}}{\text{d}t} + \boldsymbol{\mathcal{\overline M}}_{\perp s} \boldsymbol{\cdot }\frac {\text{d}_s \boldsymbol{b}}{\text{d}t} \boldsymbol{b} + \mathcal{M}_{\parallel s} \frac {\text{d}_s \boldsymbol{b}}{\text{d}t}, \end{equation}

we have, upon summing together (6.24) and (6.26), that

(6.28) \begin{equation} \frac {\partial \boldsymbol{\mathcal{\overline M}}_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s) + (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\parallel \boldsymbol{b} - q_s n_s \boldsymbol{B} \times \boldsymbol{\bar v}_{\perp s} = \boldsymbol{\mathcal{R}}_{\parallel s} + q_s n_s (\boldsymbol{E}_\parallel + \boldsymbol{V}_s \times \boldsymbol{B}). \end{equation}

Finally, given that $q_s n_s \boldsymbol{B} \times \boldsymbol{\bar v}_{\perp s} = q_s n_s \boldsymbol{E}_\perp - (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{P}_s)_\perp + \boldsymbol{\mathcal{R}}_{\perp s}$ , we have

(6.29) \begin{align} \frac {\partial \boldsymbol{\mathcal{\overline M}}_s}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s + \boldsymbol{P}_s) &= q_s n_s (\boldsymbol{E} + \boldsymbol{V}_s \times \boldsymbol{B}) + \boldsymbol{\mathcal{R}}_s , \end{align}

and, upon summing over all species $s$ , we obtain the analogue of the non-drift-reduced result, (2.15), that is,

(6.30) \begin{align} \frac {\partial \boldsymbol{\mathcal{\overline M}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\sum _s ( \boldsymbol{V}_s \boldsymbol{\mathcal{\overline M}}_s + \boldsymbol{P}_s) &= \boldsymbol{J} \times \boldsymbol{B} + \sum _s \boldsymbol{\mathcal{R}}_s , \end{align}

showing that the leading-order contribution to the total fluid momentum is conserved.