3.1. Transformation to a moving frame

Consider a velocity space coordinate transformation to a frame moving with velocity $\boldsymbol{u}$ :

(3.1) \begin{align} \boldsymbol{v} = \boldsymbol{v}' + \boldsymbol{u}(\boldsymbol{x},t). \end{align}

For this transformation, we have

(3.2) \begin{align} \frac {\partial \boldsymbol{v}'}{\partial t} &= -\frac {\partial \boldsymbol{u}}{\partial t},\\[-10pt]\nonumber\end{align}

(3.3) \begin{align} (\boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \breve {\boldsymbol{v}}')\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{v}'} &= -\left (\boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \breve {\boldsymbol{u}} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{v}'}. \end{align}

Using these in the Vlasov equation, we obtain

(3.4) \begin{align} \frac {\partial \bar {f}}{\partial t} -\frac {\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{v}'}\bar {f} + \boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot } [ (\boldsymbol{v}'+\boldsymbol{u})\bar {f} ] - [\left (\boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \breve {\boldsymbol{u}} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{v}'}] \boldsymbol{\cdot } [ (\boldsymbol{v}'{+}\boldsymbol{u})\bar {f} ] + \boldsymbol{\nabla} _{\boldsymbol{v}'}\boldsymbol{\cdot }(\boldsymbol{a}'\bar {f}) = 0, \end{align}

where $\bar {f} = \bar {f}(\boldsymbol{x}',\boldsymbol{v}',t')$ and

(3.5) \begin{align} \boldsymbol{a}' = \frac {q}{m}(\boldsymbol{E} + \boldsymbol{u}\times \boldsymbol{B}) + \frac {q}{m} \boldsymbol{v}'\times \boldsymbol{B}. \end{align}

Since $\boldsymbol{u}(\boldsymbol{x},t)$ does not depend on $\boldsymbol{v}'$ , we can rearrange various terms above to obtain

(3.6) \begin{align} \frac {\partial \bar {f}}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot }\left [ (\boldsymbol{v}'+\boldsymbol{u}){f} \right ]\nonumber\\& + \boldsymbol{\nabla} _{\boldsymbol{v}'}\boldsymbol{\cdot }\left [ \left ( -\frac {\partial \boldsymbol{u}}{\partial t} - \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} -\boldsymbol{v}'\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} + \frac {q}{m}\left (\boldsymbol{E}+\boldsymbol{u}\times \boldsymbol{B}\right ) + \frac {q}{m} \boldsymbol{v}'\times \boldsymbol{B} \right ) \bar {f} \right ] = 0. \end{align}

We can simplify this system by choosing $\boldsymbol{u}$ to be the fluid velocity

(3.7) \begin{align} \boldsymbol{u} = \frac {1}{n}\int \boldsymbol{v} f \thinspace {\rm d}^3\boldsymbol{v}, \end{align}

where $n$ is the number density. With this choice, we can use the equations for mass and momentum conservation,

(3.8) \begin{align} \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left [\rho \boldsymbol{u}\boldsymbol{u} + \boldsymbol{P} \right ] = \frac {q}{m} \rho \left [ \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} \right ] \quad \rightarrow \qquad\qquad\\[-10pt]\nonumber\end{align}

(3.9) \begin{align} \rho \left [ \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{u} \right ] + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{P} + \boldsymbol{u} \underbrace {\left [ \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left ( \rho \boldsymbol{u} \right ) \right ] }_{= 0} = \frac {q}{m} \rho \left [ \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} \right ], \end{align}

to obtain an equation for the velocity evolution

(3.10) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} + \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} = \frac {q}{m}\left (\boldsymbol{E}+\boldsymbol{u}\times \boldsymbol{B}\right ), \end{align}

where $\rho = m n$ is the mass density and $\boldsymbol{P}$ is the pressure tensor defined as

(3.11) \begin{align} \boldsymbol{P} = m \int (\boldsymbol{v}-\boldsymbol{u})\otimes (\boldsymbol{v}-\boldsymbol{u}) f \thinspace {\rm d}^3\boldsymbol{v} = m \int \boldsymbol{v}'\otimes \boldsymbol{v}' \bar {f} \thinspace {\rm d}^3\boldsymbol{v}'. \end{align}

We can then substitute (3.10) into (3.6) to simplify our transformed Vlasov equation to

(3.12) \begin{align} \frac {\partial \bar {f}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot }[ (\boldsymbol{v}'+\boldsymbol{u})\bar {f} ] + \boldsymbol{\nabla} _{\boldsymbol{v}'}\boldsymbol{\cdot }\left [ \left ( \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} -\boldsymbol{v}'\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} + \frac {q}{m} \boldsymbol{v}'\times \boldsymbol{B} \right ) \bar {f} \right ] = 0. \end{align}

We note that because the velocity coordinates of (3.12) have been transformed to move with the local fluid velocity, we must have the condition that

(3.13) \begin{align} \int \boldsymbol{v}' \bar {f} \thinspace {\rm d}^3\boldsymbol{v}' = 0. \end{align}

Multiplying (3.12) by $\boldsymbol{v}'$ and integrating over all velocity space shows that if (3.13) is satisfied at $t=0$ then it will be satisfied for all $t\gt 0$ , showing that this is an initial condition constraint on the distribution function.

In addition to solving the transformed Vlasov equation, to close the system of equations we also need to evolve the equation for the velocity (3.10), or momentum (3.8), with the pressure tensor determined by the the distribution function (3.11). The transformed Vlasov equation coupled to the momentum/velocity equation has the same physical content as the Vlasov equation in the laboratory frame.

For further reference, we derive the equation for the evolution of the pressure tensor by multiplying (3.12) with $m \boldsymbol{v}'\otimes \boldsymbol{v}'$ , integrating the velocity term by parts and then using the fact that $\boldsymbol{\nabla} _{\boldsymbol{v}'}\otimes \boldsymbol{v}'$ is the metric tensor in velocity space, to find that

(3.14) \begin{align} \frac {\partial \boldsymbol{P}}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }[ \boldsymbol{Q} + \boldsymbol{P}\otimes \breve {\boldsymbol{u}} ] + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }[ \boldsymbol{u}\otimes \boldsymbol{P}(\breve {\_},\_) + \boldsymbol{P}(\breve {\_},\_)\otimes \boldsymbol{u} ] - [ \boldsymbol{u}\otimes (\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P})\nonumber\\& + (\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P})\otimes \boldsymbol{u} )] =\frac {q}{m} [ \boldsymbol{B}\times \boldsymbol{P}(\breve {\_},\_) + \boldsymbol{B}\times \boldsymbol{P}(\_,\breve {\_}) ] . \end{align}

We can take the trace of this equation and note that $\textrm {Tr}(\boldsymbol{P}) = 3 p$ to derive an equation for the total scalar internal energy, $3p/2$ , as

(3.15) \begin{align} \frac {\partial }{\partial t} \frac {3}{2} p + \boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot }\left [ \boldsymbol{q} + \frac {3}{2} p \boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{P} \right ] = \boldsymbol{u}\boldsymbol{\cdot }[\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P}], \end{align}

where the heat-flux vector is defined as

(3.16) \begin{align} \boldsymbol{q} = \int \frac {1}{2} m v^\prime 2 \boldsymbol{v}' \bar {f} \thinspace {\rm d}^3\boldsymbol{v}' = \frac {1}{2} \textrm {Tr} ( \boldsymbol{Q} ), \end{align}

and the trace is over any two of the three slots of the fully symmetric $\boldsymbol{Q}$ . Here, $\boldsymbol{Q}$ is the third-order heat-flux tensor defined as

(3.17) \begin{align} \boldsymbol{Q} = m \int \boldsymbol{v}'\otimes \boldsymbol{v}'\otimes \boldsymbol{v}' \bar {f} \thinspace {\rm d}^3\boldsymbol{v}'. \end{align}

An alternate form of the pressure evolution equation can be derived by moving the $\boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot }[\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{P}]$ term to the right-hand side to write instead

(3.18) \begin{align} \frac {\partial }{\partial t} \frac {3}{2} p + \boldsymbol{\nabla} _{\boldsymbol{x}'}\boldsymbol{\cdot }\left [ \boldsymbol{q} + \frac {3}{2} p \boldsymbol{u} \right ] = -\boldsymbol{P} : \boldsymbol{\nabla} _{\boldsymbol{x}'}\otimes \boldsymbol{u}. \end{align}

In this form, we now recognise that the second moment of the Vlasov equation in the moving frame (3.12), gives the evolution of the internal energy of the plasma, and the total energy of the plasma can be reconstructed from the appropriate manipulation of the momentum equation (3.8) for the evolution of the total kinetic energy,

(3.19) \begin{align} \boldsymbol{u} \boldsymbol{\cdot }\left [ \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho \boldsymbol{u}\boldsymbol{u} + \boldsymbol{P} \right ) \right ] & = \boldsymbol{u} \boldsymbol{\cdot }\left [ \frac {q}{m} \rho \left ( \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} \right ) \right ] \notag \\[5pt] \frac {1}{2} \frac {\partial \rho |\boldsymbol{u}|^2}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left [\frac {1}{2} \rho |\boldsymbol{u}|^2 \boldsymbol{u} \right ] & = -\boldsymbol{u} \boldsymbol{\cdot }[\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{P}] + \frac {q}{m} \rho \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{E}, \end{align}

which when summed with (3.15) gives the identical total energy equation obtained from the untransformed Vlasov equation.

3.2. Transformation to CGL frame

We now consider a magnetised plasma. The pressure tensor $\boldsymbol{P}$ can be split into two parts, i.e.

(3.20) \begin{align} \boldsymbol{P} = \boldsymbol{P}^C + \boldsymbol{\varPi }^a, \end{align}

where $\boldsymbol{P}^C$ is the Chew-Goldberger-Low (CGL) pressure part of the tensor (Chew, Goldberger & Low Reference Chew, Goldberger and Low1956) given by

(3.21) \begin{align} \boldsymbol{P}^C = p_\parallel \boldsymbol{b}\otimes \boldsymbol{b} + p_\perp (\boldsymbol{g}-\boldsymbol{b}\otimes \boldsymbol{b}) \end{align}

and $\boldsymbol{\varPi }^a$ is the trace-free, agyrotropic part, with $\boldsymbol{g}$ the metric tensor in configuration space. We have defined $p_\parallel \equiv \boldsymbol{P} : \boldsymbol{b}\otimes \boldsymbol{b} = \boldsymbol{P}(\boldsymbol{b},\boldsymbol{b})$ , and because $\boldsymbol{\varPi }^a$ is trace free, we have $\textrm {Tr}(\boldsymbol{P}) = \textrm {Tr}(\boldsymbol{P}^C) = p_\parallel + 2 p_\perp \equiv 3p$ , where $p$ is the scalar pressure. Now, for any second-order tensor $\boldsymbol{A}$ , $\lambda$ and $\boldsymbol{r}$ are an eigenvalue/eigenvector pair if $\boldsymbol{A}(\_,\boldsymbol{r}) = \lambda \boldsymbol{r}$ . For the CGL part of the pressure tensor, it is straightforward to see that $\boldsymbol{P}^C(\_,\boldsymbol{b}) = p_\parallel \boldsymbol{b}$ and $\boldsymbol{P}^C(\_,\boldsymbol{\tau }) = p_\perp \boldsymbol{\tau }$ for all $\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\tau } = 0$ . Hence, in a frame with orthonormal unit vectors $(\boldsymbol{b}, \boldsymbol{\tau }_1, \boldsymbol{\tau }_2)$ , where $\boldsymbol{\tau }_{1,2}$ are orthogonal to $\boldsymbol{b}$ and $\boldsymbol{\tau }_{1}\times \boldsymbol{\tau }_{2} = \boldsymbol{b}$ , the CGL part of the pressure tensor will be diagonal.Footnote ³ However, in general, these vectors are not the eigenvectors of the full pressure tensor $\boldsymbol{P}$ . Nevertheless, we can always utilise this separation into gyrotropic and agyrotropic pressure tensors, and we call the frame in which the CGL part of the pressure tensor is diagonal the CGL frame.

We now transform the velocity coordinates to a local frame aligned with the $(\boldsymbol{b}, \boldsymbol{\tau }_1, \boldsymbol{\tau }_2)$ basis vectors. In this frame, we can write the transform

(3.22) \begin{align} \boldsymbol{v}' = \boldsymbol{v}'(\boldsymbol{v}'',\boldsymbol{x}'',t'') = v_\parallel \boldsymbol{b}(\boldsymbol{x}',t) + \boldsymbol{v}_\perp & = v_\parallel \boldsymbol{b}(\boldsymbol{x}',t) + v_\perp \cos \theta \boldsymbol{\tau }_1(\boldsymbol{x}',t)\nonumber\\[3pt]& \quad + v_\perp \sin \theta \boldsymbol{\tau }_2(\boldsymbol{x}',t), \end{align}

where $(w^1, w^2, w^3) \equiv (v_\parallel ,v_\perp ,\theta )$ form cylindrical velocity coordinates in the CGL frame. In general, the basis vectors depend both on position and time. For this transform, we can compute the tangent vectors in velocity space, $\tilde {\boldsymbol{e}}_{i}$ , as

(3.23) \begin{align} \tilde {\boldsymbol{e}}_{i} = \frac {\partial \boldsymbol{v}'}{\partial w^i}. \end{align}

These are $\tilde {\boldsymbol{e}}_{\parallel } = \boldsymbol{b}$ ,

(3.24) \begin{align} \tilde {\boldsymbol{e}}_{\perp } &= \cos \theta \boldsymbol{\tau }_1 + \sin \theta \boldsymbol{\tau }_2 ,\\[-10pt]\nonumber \end{align}

(3.25) \begin{align} \tilde {\boldsymbol{e}}_{\theta } &= -v_\perp \sin \theta \boldsymbol{\tau }_1 + v_\perp \cos \theta \boldsymbol{\tau }_2.\end{align}

The dual basis $\tilde {\boldsymbol{e}}^{i}$ are $\tilde {\boldsymbol{e}}^{\parallel } = \tilde {\boldsymbol{e}}_{\parallel }$ , $\tilde {\boldsymbol{e}}^{\perp } = \tilde {\boldsymbol{e}}_{\perp }$ and

(3.26) \begin{align} \tilde {\boldsymbol{e}}^{\theta } &= -\frac {1}{v_\perp } \sin \theta \boldsymbol{\tau }_1 + \frac {1}{v_\perp } \cos \theta \boldsymbol{\tau }_2. \end{align}

The Jacobian of the transform is simply $\mathcal{J}_v = v_\perp$ . In terms of $\tilde {\boldsymbol{e}}_{\perp }$ , we can write $\boldsymbol{v}_\perp = v_\perp \tilde {\boldsymbol{e}}_{\perp }$ .

We can again use (2.14) and (2.15) to derive the Vlasov equation in the CGL frame. First note that

(3.27) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{v}''}\boldsymbol{\cdot }\left ( \frac {\partial \boldsymbol{v}''}{\partial t} \right ) = \frac {\partial }{\partial t} (\boldsymbol{\nabla} _{\boldsymbol{v}''}\boldsymbol{\cdot }\boldsymbol{v}'') = 0 \end{align}

and

(3.28) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{v}''}\boldsymbol{\cdot }[ \boldsymbol{\nabla} _{\boldsymbol{x}'}\otimes \breve {\boldsymbol{v}}'' ] = \boldsymbol{\nabla} _{\boldsymbol{x}'} (\boldsymbol{\nabla} _{\boldsymbol{v}''}\boldsymbol{\cdot }\boldsymbol{v}'') = 0, \end{align}

because $\boldsymbol{\nabla} _{\boldsymbol{v}''}\boldsymbol{\cdot }\boldsymbol{v}'' = 3$ . With these identities, we can write the Vlasov equation in the CGL frame as

(3.29) \begin{align} \frac {\partial \bar {f}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}''}\boldsymbol{\cdot }[ (v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp +\boldsymbol{u}){f} ] + \frac {1}{v_\perp }\frac {\partial }{\partial w^i} \left ( v_\perp a^i \bar {f} \right ) = 0, \end{align}

where now $\bar {f} = \bar {f}(\boldsymbol{x}'',v_\parallel ,v_\perp ,\theta ,t)$ and $a^i = \tilde {\boldsymbol{e}}^{i}\boldsymbol{\cdot }\boldsymbol{a}''$ are the components of the acceleration in the CGL frame

(3.30) \begin{align} \boldsymbol{a}'' = \frac {\partial \boldsymbol{v}''}{\partial t} + (\boldsymbol{v}''+\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{v}'' + \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} -\boldsymbol{v}''\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} + \frac {q}{m} \boldsymbol{v}''\times \boldsymbol{B}, \end{align}

and $\boldsymbol{v}'' = v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp$ .

In the CGL frame, the phase-space incompressibility can be read off from (3.29) as

(3.31) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}''}\boldsymbol{\cdot }[ v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp +\boldsymbol{u}] + \frac {1}{v_\perp }\frac {\partial }{\partial w^i} ( v_\perp a^i) = 0. \end{align}

At this point, the physics content of (3.29) combined with the momentum equation (and Maxwell equations) is the same as the Vlasov--Maxwell system. No information has been lost in the process of the transformations to the local flow and CGL frames.

4.1. Fourier expansion in the CGL frame

We first expand the distribution function in velocity gyroangle using the Fourier series (dropping primes on $\boldsymbol{x}$ )

(4.1) \begin{align} \bar {f}(\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) & = \bar {f}_0(\boldsymbol{x},v_\parallel ,v_\perp ,t)\nonumber\\[3pt]& \quad + \sum _{n=1}^\infty \left [ \bar {f}_n^c(\boldsymbol{x},v_\parallel ,v_\perp ,t)\cos n\theta + \bar {f}_n^s(\boldsymbol{x},v_\parallel ,v_\perp ,t)\sin n\theta \right ]. \end{align}

Here $\bar {f}_0$ is the gyrotropic part and $\bar {f}_n^{c,s}$ are the agyrotropic parts of the distribution function. These are given by

(4.2) \begin{align} \bar {f}_0 = \frac {1}{2\pi } \int _0^{2\pi } \bar {f}(\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) \thinspace {\rm d}\theta \end{align}

and

(4.3) \begin{align} \bar {f}_n^c(\boldsymbol{x},v_\parallel ,v_\perp ,t) &= \frac {1}{\pi } \int _0^{2\pi } \bar {f}(\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) \cos n\theta \thinspace {\rm d}\theta ,\\[-10pt]\nonumber\end{align}

(4.4) \begin{align} \bar {f}_n^s(\boldsymbol{x},v_\parallel ,v_\perp ,t) &= \frac {1}{\pi } \int _0^{2\pi } \bar {f}(\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) \sin n\theta \thinspace {\rm d}\theta ,\end{align}

for $n\gt 0$ .

We can derive equations for each of the Fourier harmonics. However, to illustrate the approach, we focus on the gyrotropic distribution function $\bar {f}_0$ . We can derive the equation for the zeroth harmonic by integrating (3.29) over gyroangles to get

(4.5) \begin{align} \frac {\partial \bar {f}_0}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }[ (v_\parallel \boldsymbol{b} + \boldsymbol{u})\bar {f}_0 ] + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{M}_\perp \notag \\ & + \frac {1}{v_\perp } \frac {\partial }{\partial v_\parallel } \left ( v_\perp \frac {1}{2\pi } \int _0^{2\pi } a^\parallel \bar {f} {\rm d}{\theta } \right ) + \frac {1}{v_\perp } \frac {\partial }{\partial v_\perp } \left ( v_\perp \frac {1}{2\pi } \int _0^{2\pi } a^\perp \bar {f} {\rm d}{\theta } \right ) = 0. \end{align}

Here,

(4.6) \begin{align} \boldsymbol{M}_\perp (\boldsymbol{x},v_\parallel ,v_\perp ,t) & = \frac {1}{2\pi } \int _0^{2\pi } \boldsymbol{v}_\perp \bar {f} (\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) \thinspace {\rm d}{\theta }\nonumber\\& = \frac {v_\perp }{2} \left [ \bar {f}_1^c(\boldsymbol{x},v_\parallel ,v_\perp ,t) \boldsymbol{\tau }_1 + \bar {f}_1^s(\boldsymbol{x},v_\parallel ,v_\perp ,t) \boldsymbol{\tau }_2 \right ] \end{align}

is the contribution to the streaming term from agyrotropy of the distribution function.

We can compute $a^\parallel = \tilde {\boldsymbol{e}}^{\parallel }\boldsymbol{\cdot }\boldsymbol{a}'' = \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{a}''$ as

(4.7) \begin{align} a^\parallel = \boldsymbol{v}_\perp \boldsymbol{\cdot }\left [ \frac {\partial \boldsymbol{b}}{\partial t} + (v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{b} \right ] + \frac {1}{\rho } \boldsymbol{b}\boldsymbol{\cdot }\left [\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P}\right ] - \boldsymbol{b}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ], \end{align}

from which we find that

(4.8) \begin{align} \frac {1}{2\pi } \int _0^{2\pi } a^\parallel \bar {f} {\rm d}{\theta } & = \boldsymbol{b}\boldsymbol{\cdot }\left [ \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ] \bar {f}_0 + \boldsymbol{M}_\perp \boldsymbol{\cdot }\left [ \frac {\partial \boldsymbol{b}}{\partial t} + (v_\parallel \boldsymbol{b} + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{b} \right ] \notag \\[5pt]& \quad - \boldsymbol{b}\boldsymbol{\cdot }\left [\boldsymbol{M}_\perp \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}\right ] + \boldsymbol{F}_{\perp \perp } : \boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \boldsymbol{b}, \end{align}

where

(4.9) \begin{align} \boldsymbol{F}_{\perp \perp }(\boldsymbol{x},v_\parallel ,v_\perp ,t) \equiv \frac {1}{2\pi } \int _0^{2\pi } (\boldsymbol{v}_\perp \otimes \boldsymbol{v}_\perp ) \bar {f} (\boldsymbol{x},v_\parallel ,v_\perp ,\theta ,t) \thinspace {\rm d}{\theta }\end{align}

is a second-order symmetric tensor. Using the Fourier expansion of the distribution function, we get

(4.10) \begin{align} \boldsymbol{F}_{\perp \perp }(\boldsymbol{x},v_\parallel ,v_\perp ,t) & = \frac {v_\perp ^2\bar {f}_0}{2} \left ( \boldsymbol{g} - \boldsymbol{b}\otimes \boldsymbol{b} \right ) + \frac {v_\perp ^2\bar {f}_2^c}{4} \left ( \boldsymbol{\tau }_1\otimes \boldsymbol{\tau }_1 - \boldsymbol{\tau }_2\otimes \boldsymbol{\tau }_2 \right )\nonumber\\[5pt]& \quad + \frac {v_\perp ^2\bar {f}_2^s}{4} \left ( \boldsymbol{\tau }_1\otimes \boldsymbol{\tau }_2 + \boldsymbol{\tau }_2\otimes \boldsymbol{\tau }_1 \right ), \end{align}

where we used $\boldsymbol{g} = \boldsymbol{\tau }_1\otimes \boldsymbol{\tau }_1 + \boldsymbol{\tau }_2\otimes \boldsymbol{\tau }_2 + \boldsymbol{b}\otimes \boldsymbol{b}$ to rewrite the first term. Using this result, we can compute

(4.11) \begin{align} \boldsymbol{F}_{\perp \perp } : \boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \boldsymbol{b} & = \frac {v_\perp ^2\bar {f}_0}{2} \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} - \frac {v_\perp ^2\bar {f}_2^c}{4} \boldsymbol{b}\boldsymbol{\cdot }\left ( \boldsymbol{\tau }_1\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }_1 - \boldsymbol{\tau }_2\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }_2 \right )\nonumber\\[5pt]& \quad - \frac {v_\perp ^2\bar {f}_2^s}{4} \boldsymbol{b}\boldsymbol{\cdot }\left ( \boldsymbol{\tau }_1\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }_2 + \boldsymbol{\tau }_2\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }_1 \right ). \end{align}

The first term in this expression is the contribution to the parallel acceleration from the gyrotropic part of the distribution function, and the remaining two terms are the contributions from the agyrotropy of the distribution function. An alternate form of the first term can be obtained, since we can replace

(4.12) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} = \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\frac {\boldsymbol{B}}{B} = \frac {1}{B}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{B} + \boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}{\frac {1}{B}} = -\frac {1}{B} \boldsymbol{\nabla} _\parallel B, \end{align}

where $\boldsymbol{\nabla} _\parallel \equiv \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}$ .

Similarly, we can compute $a^\perp = \tilde {\boldsymbol{e}}^{\perp }\boldsymbol{\cdot }\boldsymbol{a}''$ as

(4.13) \begin{align} a^\perp &= \tilde {\boldsymbol{e}}^{\perp } \boldsymbol{\cdot }\left [ -v_\parallel \frac {\partial \boldsymbol{b}}{\partial t} -v_\parallel (v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{b} + \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - (v_\parallel \boldsymbol{b} + \boldsymbol{v}_\perp )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ].\nonumber\\[5pt] \end{align}

Since $v_\perp \tilde {\boldsymbol{e}}^{\perp } = \boldsymbol{v}_\perp$ , we find that

(4.14) \begin{align} v_\perp \frac {1}{2\pi } \int _0^{2\pi } a^\perp \bar {f} {\rm d}{\theta } & = \boldsymbol{M}_\perp \boldsymbol{\cdot }\left [ -v_\parallel \frac {\partial \boldsymbol{b}}{\partial t} -v_\parallel (v_\parallel \boldsymbol{b} + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{b} + \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ] \notag \\[5pt]& \quad - \boldsymbol{F}_{\perp \perp } : (v_\parallel \boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \boldsymbol{b} + \boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \boldsymbol{u}). \end{align}

Though the first term of the velocity gyroangle integrated perpendicular acceleration has only agyrotropic contributions arising from the coupling to the first Fourier harmonic, $\boldsymbol{M}_\perp$ , the $\boldsymbol{F}_{\perp \perp }$ term in $a^\perp$ does have a gyrotropic contribution as well -- see (4.10) and the term proportional to $\bar {f}_0$ . Putting all these expressions together, we finally obtain the equation for the gyrotropic distribution function:

(4.15) \begin{align} \frac {\partial \bar {f}_0}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b} + \boldsymbol{u})\bar {f}_0 \right ] + \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right )\bar {f}_0 + \frac {v_\perp ^2}{2}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \bar {f}_0 \right ] \notag \\[5pt]& + \frac {1}{v_\perp } \frac {\partial }{\partial v_\perp } \left [ \frac {v_\perp ^2}{2} \left ( \boldsymbol{b}\boldsymbol{\cdot }\{\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}\} - \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\{v_\parallel \boldsymbol{b}+ \boldsymbol{u}\} \right )\bar {f}_0 \right ] + \mathrm{AG} = 0. \end{align}

Here $\mathrm{AG}$ are the agyrotropic terms, which involve the first and second Fourier harmonics of the distribution function. The expressions for these terms are given in Appendix C. We note that to determine the coupling of $\bar {f}_0$ to the agyrotropic terms, further equations would be needed; for example, evolution equations for $\boldsymbol{M}_\perp$ and $\boldsymbol{F}_{\perp \perp }$ can be derived by multiplying the Vlasov equation in our transformed coordinates, (3.29), by $\boldsymbol{v}_\perp$ and $\boldsymbol{v}_\perp \otimes \boldsymbol{v}_\perp$ and integrating over the velocity gyroangle. In fact, this procedure of determining equations for, e.g. $\boldsymbol{M}_\perp$ and $\boldsymbol{F}_{\perp \perp }$ , provides a natural procedure for the determination of higher Fourier harmonics, at least through the second Fourier harmonic, because the coupling to the $f^{c,s}_{1,2}$ coefficients of the Fourier expansion appears only through these vector and tensor combinations of Fourier coefficients and plane perpendicular to the magnetic field, $\boldsymbol{\tau }_{1,2}$ .

In terms of the Fourier expansion, various velocity moments can be determined through integration over $v_\parallel$ and $v_\perp$ . For example, the number density is

(4.16) \begin{align} n &= 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \bar {f}_0, \end{align}

and the CGL pressure tensor components can be computed as

(4.17) \begin{align} p_\parallel &= 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } m v_\parallel ^2 \bar {f}_0,\\[-10pt]\nonumber \end{align}

(4.18) \begin{align} p_\perp &= 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \frac {1}{2} m v_\perp ^2 \bar {f}_0. \end{align}

Thus, the gyrotropic part of the distribution function contains the total density and total internal energy. The agyrotropic part of the pressure tensor is

(4.19) \begin{align} \boldsymbol{\varPi }^a &= 2\pi m \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \times \bigg[ v_\parallel (\boldsymbol{b}\otimes \boldsymbol{M}_\perp + \boldsymbol{M}_\perp \otimes \boldsymbol{b})\nonumber\\& \quad + \frac {v_\perp ^2\bar {f}_2^c}{4} \left ( \boldsymbol{\tau }_1\otimes \boldsymbol{\tau }_1 - \boldsymbol{\tau }_2\otimes \boldsymbol{\tau }_2 \right ) + \frac {v_\perp ^2\bar {f}_2^s}{4} \left ( \boldsymbol{\tau }_1\otimes \boldsymbol{\tau }_2 + \boldsymbol{\tau }_2\otimes \boldsymbol{\tau }_1 \right ) \bigg]. \end{align}

Clearly, both the first and second Fourier harmonics contribute to $\boldsymbol{\varPi }^a$ . In the gyrotropic limit, it follows that we must have $\boldsymbol{\varPi }^a = 0$ .

The heat-flux vector defined in (3.16) is

(4.20) \begin{align} \boldsymbol{q} = (q_\parallel + q_\perp ) \boldsymbol{b} + 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \thinspace \frac {1}{2} m \big(v_\parallel ^2 + v_\perp ^2\big) \boldsymbol{M}_\perp , \end{align}

where

(4.21) \begin{align} q_\parallel &= 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \thinspace \frac {1}{2} m v_\parallel ^3 \bar {f}_0,\\[-10pt]\nonumber \end{align}

(4.22) \begin{align} q_\perp &= 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \thinspace \frac {1}{2} m v_\parallel v_\perp ^2 \bar {f}_0 \end{align}

are the parallel components of the parallel and perpendicular heat fluxes, respectively. These components are determined fully from the gyrotropic distribution function. On the other hand, the components of the heat-flux vector perpendicular to the magnetic field are determined from the first Fourier harmonic.

The characteristic velocities for the evolution of the gyrotropic distribution function can be read off from (4.15) with $\mathrm{AG}=0$ as

(4.23) \begin{align} \dot {\boldsymbol{x}} &= v_\parallel \boldsymbol{b} + \boldsymbol{u}, \end{align}

(4.24) \begin{align} \dot {v}_\parallel &= \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) + \frac {v_\perp ^2}{2}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b}, \end{align}

(4.25) \begin{align} \dot {v}_\perp &= \frac {v_\perp }{2} \left [ \boldsymbol{b}\boldsymbol{\cdot }\{\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}\} - \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\{v_\parallel \boldsymbol{b}+ \boldsymbol{u}\} \right ]. \end{align}

Furthermore, integrating the phase-space incompressibility condition (3.31) over all angles gives

(4.26) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\dot {\boldsymbol{x}} + \frac {\partial \dot {v}_\parallel }{\partial v_\parallel } + \frac {1}{v_\perp } \frac {\partial }{\partial v_\perp }(v_\perp \dot {v}_\perp ) = 0. \end{align}

Hence, in the gyrotropic limit, the evolution of the gyrotropic distribution function does not violate phase-space incompressibility. For reference, the portion of the parallel acceleration coming from the pressure tensor can be rewritten in a more familiar form as

(4.27) \begin{align} \frac {1}{\rho } \boldsymbol{b}\boldsymbol{\cdot }[\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P}] & = \frac {1}{\rho }\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} p_{\parallel } + \frac {(p_{\parallel }-p_{\perp })}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} + \frac {1}{\rho } \boldsymbol{b}\boldsymbol{\cdot }\left [ \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{\varPi }^a \right ], \notag \\[3pt] & = \frac {1}{\rho } \left [ \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left ( p_{\parallel } \boldsymbol{b} \right ) - p_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] + \frac {1}{\rho } \boldsymbol{b}\boldsymbol{\cdot }\left [ \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{\varPi }^a \right ]. \end{align}

4.2. Reduction to the parallel-kinetic-perpendicular-moment system

We now expand the gyrotropic distribution function in a Laguerre series in $v_\perp$ as

(4.28) \begin{align} \bar {f}_0(\boldsymbol{x},v_\parallel ,v_\perp ,t) = \sum _{n=0}^\infty F_{n}(\boldsymbol{x},v_\parallel ,t) G_M(v_\perp ,T_\perp ) L_n\left (\frac {m v_\perp ^2}{2 T_\perp }\right ), \end{align}

where $L_n(x)$ are Laguerre polynomials of order $n$ and

(4.29) \begin{align} G_M(v_\perp ,T_\perp ) = \frac {m}{2\pi T_\perp }e^{-m v_\perp ^2/2 T_\perp }. \end{align}

In these expressions, $T_\perp = T_\perp (\boldsymbol{x},t)$ is the perpendicular temperature, in terms of which the perpendicular pressure is $p_\perp = n T_\perp$ . The derivation of the equations for the Laguerre coefficients $F_k$ is complicated by the fact that $T_\perp$ depends on position as well as time. Typically, in this expansion, $T_\perp$ would need to be determined from another equation, such as an equation for the perpendicular pressure.

The $F_k(\boldsymbol{x},v_\parallel ,t)$ coefficient of this series is

(4.30) \begin{align} F_{k}(\boldsymbol{x},v_\parallel ,t) = 2\pi \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } L_k\left (\frac {m v_\perp ^2}{2 T_\perp }\right ) \bar {f}_0(\boldsymbol{x},v_\parallel ,v_\perp ,t). \end{align}

The number density and the parallel pressure (see (4.16) and (4.17)) are completely determined from the zeroth Laguerre coefficient of the gyrotropic distribution function:

(4.31) \begin{align} n &= \int _{-\infty }^\infty F_0 \thinspace {\rm d}{v_\parallel },\\[-10pt]\nonumber \end{align}

(4.32) \begin{align} p_\parallel &= \int _{-\infty }^\infty m v_\parallel ^2 F_0 \thinspace {\rm d}{v_\parallel }. \end{align}

The parallel components of the heat flux are

(4.33) \begin{align} q_\parallel &= \int _{-\infty }^\infty \frac {1}{2} m v_\parallel ^3 F_0 \thinspace {\rm d}{v_\parallel },\\[-10pt]\nonumber \end{align}

(4.34) \begin{align} q_\perp &= T_\perp \int _{-\infty }^\infty v_\parallel (F_0-F_1) \thinspace {\rm d}{v_\parallel }, \end{align}

which are determined from just the first two Laguerre coefficients of the gyrotropic distribution function. Finally, the constraint (3.13) on the $\boldsymbol{v}'$ moment of the distribution function shows that we must have

(4.35) \begin{align} \int _{-\infty }^\infty v_\parallel F_0 \thinspace {\rm d}{v_\parallel } &= 0,\\[-10pt]\nonumber \end{align}

(4.36) \begin{align} \int _0^\infty v_\perp \thinspace {\rm d}{v_\perp } \int _{-\infty }^\infty \thinspace {\rm d}{v_\parallel } \thinspace \boldsymbol{M}_\perp &= 0. \end{align}

4.2.1. Some preliminaries

To derive the equations for $F_{k}$ for $k\gt 0$ (recall that $L_0 = 1$ ), we need to account for the fact that $T_\perp$ depends both on space and time. To do this derivation then, we write

(4.37) \begin{align} L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \frac {\partial \bar {f}_0}{\partial t} = \frac {\partial }{\partial t} \left [ L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \bar {f}_0 \right ] - \bar {f}_0 \frac {\partial }{\partial t} \left [ L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \right ]. \end{align}

Now, using (B.5) we get

(4.38) \begin{align} \frac {\partial }{\partial t} L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) = -\frac {1}{T_\perp }\frac {\partial T_\perp }{\partial t} \left [ k L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) - k L_{k-1}\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \right ]. \end{align}

Also,

(4.39) \begin{align} L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})\bar {f}_0 \right ] & = \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u}) L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \bar {f}_0 \right ]\nonumber\\[5pt]& \quad - \bar {f}_0 \left \{ (v_\parallel \boldsymbol{b}+\boldsymbol{u}) \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \left [ L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \right ] \right \}, \end{align}

where, as above, we can write

(4.40) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}} \left [ L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \right ] = -\frac {1}{T_\perp }\boldsymbol{\nabla} _{\boldsymbol{x}} T_\perp \left [ k L_k\left (\frac {m v_\perp ^2}{2T_\perp }\right ) - k L_{k-1}\left (\frac {m v_\perp ^2}{2T_\perp }\right ) \right ]. \end{align}

4.2.2. The equation for $F_k$

We can multiply (4.15) by $L_k(m v_\perp ^2/2 T_\perp )$ and integrate over all $v_\perp$ to obtain

(4.41) \begin{align} \frac {\partial F_{k}}{\partial t} &+ \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})F_{k} \right ] + \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) F_{k}\right.\nonumber\\[5pt]& \quad +\left. [ (2k+1) F_k - k F_{k-1} - (k+1) F_{k+1} ]\frac {T_\perp }{m}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] \notag \\ &\quad = S_k^T(\boldsymbol{x},v_\parallel ) + S_k(\boldsymbol{x},v_\parallel ). \end{align}

The source $S^T_k$ appears because $T_\perp$ depends both on space and time. Using the expressions derived in the previous section, we can write this source term as

(4.42) \begin{align} S^T_k(\boldsymbol{x},v_\parallel ) = -\frac {1}{T_\perp } \left [ \frac {\partial T_\perp }{\partial t} + (v_\parallel \boldsymbol{b} + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} T_\perp \right ] (k F_{k} - k F_{k-1}). \end{align}

The source $S_k$ results from the integration by parts of the perpendicular velocity derivative. Calculating this term, we get

(4.43) \begin{align} S_k(\boldsymbol{x},v_\parallel ) = \left [ \boldsymbol{b}\boldsymbol{\cdot }(\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}) - \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }(v_\parallel \boldsymbol{b} + \boldsymbol{u}) \right ] (k F_{k} - k F_{k-1}). \end{align}

At this point we have all the equations we need to evolve the Laguerre coefficients of the gyrotropic distribution function. The number of Laguerre coefficients one will need will depend on the structure of the distribution function in the perpendicular coordinate. One yet undetermined quantity in these equations is $T_\perp$ . We can derive an explicit equation for this quantity; however, we can also pursue an alternate approach that will completely eliminate the need for an auxiliary equation for $T_\perp$ .

4.3. The first two Laguerre coefficients and perpendicular temperature

The evolution of the first two Laguerre coefficients of the gyrotropic distribution can be determined by first setting $k=0$ in (4.41). This substitution gives

(4.44) \begin{align} \frac {\partial F_{0}}{\partial t} &+ \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})F_{0} \right ] \notag \\ &+ \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) F_{0} + \mathcal{G}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] = 0, \end{align}

where we have defined

(4.45) \begin{align} \mathcal{G} \equiv \frac {T_\perp }{m}(F_0 - F_1). \end{align}

Furthermore, setting $k=1$ in (4.41) gives

(4.46) \begin{align} &\frac {\partial F_{1}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})F_{1} \right ] \notag \\[3pt] &\qquad + \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) F_{1} + \left ( 3 F_1 - F_{0} - 2 F_{2} \right )\frac {T_\perp }{m}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] \notag \\[6pt] & \quad = -\frac {1}{T_\perp } \left [ \frac {\partial T_\perp }{\partial t} + (v_\parallel \boldsymbol{b} + \boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} T_\perp \right ] (F_{1} - F_{0}) + S_1(\boldsymbol{x},v_\parallel ). \end{align}

Subtracting the second equation from the first, and multiplying throughout by $T_\perp /m$ , we get

(4.47) \begin{align} &\frac {\partial \mathcal{G}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})\mathcal{G} \right ] \notag \\[4pt] & \quad \quad + \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) \mathcal{G} + \left ( 4\mathcal{G} + \frac {2 T_\perp }{m}\{F_2-F_0\} \right )\frac {T_\perp }{m}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] \notag \\[5pt] &\quad = \left [ \boldsymbol{b}\boldsymbol{\cdot }(\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}) - \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }(v_\parallel \boldsymbol{b} + \boldsymbol{u}) \right ] \mathcal{G}. \end{align}

Instead of solving the equation for $F_1$ , we can solve this equation for $\mathcal{G}$ . The advantage of this approach is that the source $S^T_1(\boldsymbol{x},v_\parallel )$ due to the temporal and spatial variation of $T_\perp$ does not appear, avoiding the need to compute the time and spatial derivatives of $T_\perp$ . Furthermore, we claim that

(4.48) \begin{align} p_\perp = n T_\perp = T_\perp \int _{-\infty }^\infty F_0 \thinspace {\rm d}{v_\parallel } = \int _{-\infty }^\infty T_\perp (F_0-F_1) \thinspace {\rm d}{v_\parallel } = m \int _{-\infty }^\infty \mathcal{G} \thinspace {\rm d}{v_\parallel }, \end{align}

thus eliminating the need to evolve an explicit equation for the perpendicular temperature.

As a consistency check, we integrate (4.47) over all velocities to obtain

(4.49) \begin{align} \frac {\partial p_\perp }{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ q_\perp \boldsymbol{b} + p_\perp \boldsymbol{u} \right ] = p_\perp \boldsymbol{b}\boldsymbol{\cdot }[\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}] - q_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} - p_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u}. \end{align}

Likewise, multiplying (4.44) by $1/2 \thinspace m v_\parallel ^2$ and integrating over all velocities we get an evolution equation for the parallel pressure:

(4.50) \begin{align} \frac {1}{2} \frac {\partial p_\parallel }{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ q_\parallel \boldsymbol{b} + \frac {1}{2} p_\parallel \boldsymbol{u} \right ] = -p_\parallel \boldsymbol{b}\boldsymbol{\cdot }[\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}] + q_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b}. \end{align}

If we add (4.49) and (4.50), we obtain

(4.51) \begin{align} \frac {\partial }{\partial t} \left (\frac {1}{2} p_\parallel + p_\perp \right ) & + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ \left (q_\parallel + q_\perp \right ) \boldsymbol{b} + \left (\frac {1}{2} p_\parallel + p_\perp \right ) \boldsymbol{u} \right ]\nonumber\\[5pt]& = (p_\perp - p_\parallel ) \boldsymbol{b}\boldsymbol{\cdot }[\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}] - p_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u}. \end{align}

Recall that in the gyrotropic limit,

(4.52) \begin{align} \boldsymbol{P} = \boldsymbol{P}^C = p_\parallel \boldsymbol{b}\otimes \boldsymbol{b} + p_\perp (\boldsymbol{g}-\boldsymbol{b}\otimes \boldsymbol{b}), \end{align}

which implies that

(4.53) \begin{align} (p_\perp - p_\parallel ) \boldsymbol{b}\boldsymbol{\cdot }[\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}] - p_\perp \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{u} = -\boldsymbol{P} : \boldsymbol{\nabla} _{\boldsymbol{x}}\otimes \boldsymbol{u}. \end{align}

Thus, utilising the fact that $3 p/2 = p_\parallel /2 + p_\perp$ , we obtain the gyrotropic limit of (3.18) from the second moment of the $F_0$ equation and the zeroth moment of the $\mathcal{G}$ equation.

The statement that we evolve the perpendicular pressure $p_\perp$ by evolving $\mathcal{G}$ is non-trivial. In effect, we have absorbed the $T_\perp$ normalisation of the Laguerre expansion that appears in our fundamental ansatz (4.28) into the evolution of $\mathcal{G}$ , even though this normalisation must typically be determined independently and not depend on Laguerre coefficients. Thus, our spectral expansion in Laguerres can fully leverage the optimised spatially and temporarily dependent $T_\perp$ normalisation without the need for auxiliary equations. Note this formulation of the equation for $\mathcal{G}$ indicates that we must have the constraint

(4.54) \begin{align} \int _{-\infty }^\infty F_1 \thinspace {\rm d}{v_\parallel } = 0. \end{align}

From (4.46), can show that this is an initial-value constraint: that we must ensure (4.46) is true at $t=0$ for it to be satisfied for all $t\gt 0$ .

4.4. The lowest-order PKPM system of equations

We define the lowest-order PKPM system as the zeroth Fourier harmonic and the first two Laguerre coefficients, with the truncation that all higher Fourier harmonics and all higher Laguerre coefficients are zero. With these approximations, we are left with the following two coupled four-dimensional kinetic equations:

(4.55) \begin{align} \frac {\partial F_{0}}{\partial t} &+ \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})F_{0} \right ] \notag \\ &+ \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) F_{0} + \mathcal{G}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] = 0, \end{align}

(4.56) \begin{align} \frac {\partial \mathcal{G}}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [ (v_\parallel \boldsymbol{b}+\boldsymbol{u})\mathcal{G} \right ] \notag \\[5pt] & + \frac {\partial }{\partial v_\parallel } \left [ \boldsymbol{b}\boldsymbol{\cdot }\left ( \frac {1}{\rho }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{P} - v_\parallel \boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u} \right ) \mathcal{G} + \left [ 4\mathcal{G} - 2 \frac {T_\perp }{m} F_0 \right ]\frac {T_\perp }{m}\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{b} \right ] \notag \\[5pt] & = \left [ \boldsymbol{b}\boldsymbol{\cdot }(\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{u}) - \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }(v_\parallel \boldsymbol{b} + \boldsymbol{u}) \right ] \mathcal{G}. \end{align}

The system is closed via the solution of an equation for the flow velocity, $\boldsymbol{u}$ , which comes from conservation of momentum,

(4.57) \begin{align} & \frac {\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\left [\rho \boldsymbol{u} \otimes \boldsymbol{u} + \boldsymbol{P} \right ] = \rho \frac {q}{m} \left [ \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} \right ],\\[-10pt]\nonumber \end{align}

(4.58) \begin{align} & \boldsymbol{P} = p_\parallel \boldsymbol{b}\otimes \boldsymbol{b} + p_\perp (\boldsymbol{g}-\boldsymbol{b}\otimes \boldsymbol{b}),\\[-10pt]\nonumber \end{align}

(4.59) \begin{align} & p_\parallel = m \int v_\parallel ^2 F_0 \thinspace {\rm d}v_\parallel ,\\[-10pt]\nonumber \end{align}

(4.60) \begin{align} & p_\perp = m \int \mathcal{G} \thinspace {\rm d}v_\parallel , \end{align}

plus Maxwell’s equations for the evolution of the electromagnetic fields,

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

(4.62) \begin{align} \epsilon _0\mu _0\frac {\partial \boldsymbol{E}}{\partial t} - \boldsymbol{\nabla} _{\boldsymbol{x}}\times \boldsymbol{B} &= -\mu _0\boldsymbol{J},\\[-10pt]\nonumber \end{align}

(4.63) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{E}&= \frac {\varrho _c}{\epsilon _0},\\[-10pt]\nonumber \end{align}

(4.64) \begin{align} \boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{B}&= 0, \end{align}

with the coupling to the electromagnetic fields arising directly from the momentum equation through the plasma currents,

(4.65) \begin{align} \boldsymbol{J} = \sum _s \frac {q_s}{m_s} \rho _s \boldsymbol{u}_s. \end{align}

5.1. Connection to Kulsrud’s KMHD and Ramos’ FLR kinetic theory

One of the most foundational models in the theory of magnetised plasmas is KMHD, often referred to as Kulsrud’s KMHD (Kulsrud Reference Krommes1964, Reference Kuldinow, Yamashita, Mansour and Hara1983). Kinetic magnetohydrodynamics is a hybrid fluid-kinetic model that transforms the kinetic equation in an identical fashion as the PKPM derivation outlined here, first transforming velocity coordinates to a bulk fluid frame, specifically the $\boldsymbol{E} \times \boldsymbol{B}$ frame, and then transforming the velocity coordinates to a field-aligned coordinate system $(v_\parallel , v_\perp , \theta )$ . However, following these transformations, an asymptotic expansion is performed with small parameter $\epsilon = m/e$ or simply $\epsilon = 1/e$ , where $m$ is the species mass and $e$ is the elementary charge. Because the bulk fluid flow being utilised is only the $\boldsymbol{E} \times \boldsymbol{B}$ velocity, this asymptotic expansion yields an immediate dynamical consequence: the parallel electric field is also $\mathcal{O}(\epsilon )$ smaller than the perpendicular electric field, and to lowest order the distribution function is gyrotropic.

In this limit, the gyrotropic distribution function evolves as

(5.1) \begin{align} \frac {\partial \bar {f}_0}{\partial t} & + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left [ \left (v_\parallel \boldsymbol{b} + \boldsymbol{U}_E \right ) \bar {f}_0 \right ] \notag \\[4pt] & + \frac {\partial }{\partial v_\parallel } \left [ \left (-\boldsymbol{b} \boldsymbol{\cdot }\frac {\partial \boldsymbol{U}_e}{\partial t} - \boldsymbol{b} \boldsymbol{\cdot }\left \{\boldsymbol{U}_e \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{U}_e \right \} + \frac {q_s}{m_s} E_\parallel\right.\right.\nonumber\\[4pt] & \qquad\qquad\left.\left. - \boldsymbol{b} \boldsymbol{\cdot }\left \{v_\parallel \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{U}_e \right \} + \frac {v_\perp ^2}{2} \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{b} \right ) \bar {f}_0 \right ] \notag \\[4pt] & + \frac {\partial }{\partial v_\perp ^2/2} \left [ \frac {v_\perp ^2}{2} \left ( \boldsymbol{b} \boldsymbol{\cdot }\left \{ \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}\boldsymbol{U}_e \right \} - \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left \{ v_\parallel \boldsymbol{b} + \boldsymbol{U}_e \right \} \right ) \bar {f}_0 \right ] = 0, \end{align}

where

(5.2) \begin{align} \boldsymbol{U}_E = \frac {\boldsymbol{E} \times \boldsymbol{B}}{|\boldsymbol{B}|^2} \end{align}

is the $\boldsymbol{E} \times \boldsymbol{B}$ velocity, and we have transformed the perpendicular velocity derivative from $1/v_\perp \partial _{v_\perp }$ to a derivative on the variable $v_\perp ^2/2$ for convenience since everywhere the perpendicular velocity appears in this equation it appears as $v_\perp ^2/2$ . Note that typically, a further coordinate transformation is performed from $v_\perp ^2/2$ to the magnetic moment $\mu = v_\perp ^2/(2B)$ , where $B = |\boldsymbol{B}|$ is the magnitude of the magnetic field, to further simplify the kinetic equation by eliminating the derivative in $v_\perp ^2/2$ since the magnetic moment is conserved by the dynamics of this equation. Irrespective of whether one utilises $v_\perp$ , $v_\perp ^2/2$ or $\mu$ as a velocity coordinate, (5.1) is often referred to as the drift-kinetic equation (Frieman et al. Reference Frei, Jorge and Ricci1966; Hinton & Wong Reference Hénon1985).

At first glance, it would seem that the only difference between this approach and the PKPM approach, at least when keeping only the zeroth Fourier harmonic, i.e. the gyrotropic component of the distribution function, is which particular hydrodynamic velocity is utilised for the initial velocity coordinate transformation. The characteristics of the particles in all the dimensions are identical save for a variable substitution of $\boldsymbol{U}_E \rightarrow \boldsymbol{u}$ and a subsequent manipulation of the parallel forces to utilise conservation of momentum,

(5.3) \begin{align} -\boldsymbol{b} \boldsymbol{\cdot }\frac {\partial \boldsymbol{u}}{\partial t} - \boldsymbol{b} \boldsymbol{\cdot }\left [\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{u} \right ] + \frac {q}{m} E_\parallel = \boldsymbol{b} \boldsymbol{\cdot }\left [ \frac {1}{\rho } \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{P} \right ]. \end{align}

Indeed, we expect on scales larger than the gyroradius and in plasmas where the bulk velocity is dominated by the $\boldsymbol{E} \times \boldsymbol{B}$ velocity that the PKPM model contains all the physics of KMHD.

The differences between the PKPM approach and the classical approach of Kulsrud are made more salient when considering how the system is closed. To close the KMHD system, the kinetic equation in (5.1), for each species, is coupled to a set of bulk fluid equations for the evolution of the total mass density and total momentum,

(5.4) \begin{align} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left [ \rho \boldsymbol{U} \right ] = 0, \qquad\qquad\\[-10pt]\nonumber\end{align}

(5.5) \begin{align} \frac {\partial \rho \boldsymbol{U}}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\Big[ \rho \boldsymbol{U}\otimes \boldsymbol{U} + \sum _s \boldsymbol{\mathcal{P}}_s \Big] = \boldsymbol{J} \times \boldsymbol{B}, \end{align}

where $\rho$ is the total mass density, often approximated as simply the ion mass density, and $\boldsymbol{U}$ is the centre-of-mass velocity, often approximated as simply the ion bulk velocity. Here, $\boldsymbol{\mathcal{P}}_s$ is the pressure tensor of species $s$ ,

(5.6) \begin{align} \boldsymbol{\mathcal{P}}_s = p_{\parallel _s} \boldsymbol{b}\otimes \boldsymbol{b} + p_{\perp _s} (\boldsymbol{g}-\boldsymbol{b}\otimes \boldsymbol{b}),\\[-10pt]\nonumber \end{align}

(5.7) \begin{align} p_{\parallel _s} = \int m_s (v_\parallel - u_{\parallel _s})^2 f_s {\rm d}v_\parallel d \left (\frac {v_\perp ^2}{2} \right ),\\[-13pt]\nonumber \end{align}

(5.8) \begin{align} p_{\perp _s} = \int m_s \frac {v_\perp ^2}{2} f_s {\rm d}v_\parallel d \left (\frac {v_\perp ^2}{2} \right ),\qquad \end{align}

where we have reintroduced the species index to the mass, parallel velocity and distribution function for clarity. We can see that $\boldsymbol{\mathcal{P}}_s$ is the same gyrotropic pressure tensor we defined earlier for the lowest-order PKPM system, but with a different means of computing the parallel pressure and assuming the only bulk perpendicular velocity is $\boldsymbol{E} \times \boldsymbol{B}$ ; in other words, there is an implicit ordering that the remaining guiding centre drifts, such as $\boldsymbol{\nabla} _{\boldsymbol{x}} B$ and curvature drifts, are all smaller than $\boldsymbol{E} \times \boldsymbol{B}$ . The forces on the bulk motion after summing over all species only require $\boldsymbol{J}$ , the current density of the plasma, usually given simply by

(5.9) \begin{align} \boldsymbol{J} = \frac {1}{\mu _0} \boldsymbol{\nabla} _{\boldsymbol{x}} \times \boldsymbol{B}, \end{align}

if one also assumes non-relativistic flows, though generalisations of the KMHD approach to relativity are possible via the relativistic generalisations of guiding centre theory (Vandervoort Reference Vandervoort1960; Ripperda et al. Reference Ramos2018; Bacchini et al. 2020; Trent et al. Reference Trent, Christian, Chan, Psaltis and Özel2023, Reference Trent, Roley, Golden, Psaltis and Özel2024). Finally, the time evolution of the magnetic field is given by Faraday’s equation

(5.10) \begin{align} \frac {\partial \boldsymbol{B}}{\partial t} = - \boldsymbol{\nabla} _{\boldsymbol{x}} \times \boldsymbol{E}, \end{align}

and in the commonly employed non-relativistic limit an Ohm’s law for the electric field gives the electric field

(5.11) \begin{align} \boldsymbol{E} = \boldsymbol{U} \times \boldsymbol{B}. \end{align}

So, up to the coupling to the pressure tensor of each evolved species, the fluid equations, Faraday’s law and the Ohm’s law for the electric field are exactly the same equations as those of ideal MHD. It is only the closure -- how the pressure is determined for the evolution of the momentum -- that is affected by the kinetic response of the plasma.

Immediately we note two subtleties to the evolution of the KMHD system: the total momentum equation permits the development of perpendicular motions that are not $\boldsymbol{E} \times \boldsymbol{B}$ , but these bulk motions do not explicitly feedback on the kinetic response of the plasma, and the parallel electric field is as yet unspecified by the Ohm’s law. Firstly, the bulk motions drive perpendicular currents such as

(5.12) \begin{align} \boldsymbol{J}_{\perp _{\boldsymbol{\mathcal{P}}}} & = -\frac {1}{|\boldsymbol{B}|^2} \left [ \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left ( \sum _s \boldsymbol{\mathcal{P}}_s \right ) \times \boldsymbol{B} \right ] \notag \\[8pt] & = -\underbrace {\frac {\boldsymbol{\nabla} _{\boldsymbol{x}} \left (\sum _s p_{\perp _s} \right ) \times \boldsymbol{B}}{|\boldsymbol{B}|^2}}_{\textrm {diamagnetic}} - \underbrace {\sum _s \left (p_{\perp _s} - p_{\parallel _s} \right ) \frac {\left (\boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{b} \right ) \times \boldsymbol{B}}{|\boldsymbol{B}|^2}}_{\textrm {curvature}}, \end{align}

which we have labelled as the diamagnetic and curvature-driven perpendicular currents, respectively. However, as we noted earlier, the form of the kinetic equation in KMHD assumes that there are no further contributions of perpendicular bulk flows to the computation of the perpendicular pressure. Thus, for example, any diamagnetic currents that should modify the equilibrium pressure profile are not correctly captured by KMHD in its current form, at least dynamically, because the response of the kinetic equation does not include how the particles react to the presence of other perpendicular flows. These perpendicular currents do implicitly modify the local magnetic field, which itself may modify the plasma pressure through, e.g. adiabatic heating and cooling, but the interplay between the development of perpendicular flows and heating of the plasma is not contained in the KMHD system as written. Importantly, these sorts of time-dependent interplays between different components of the system, while subdominant in the asymptotics, may be a critical component of well-behaved numerical discretisations of time-dependent partial differential equations, as including the interplay allows the system to relax to the desired equilibrium instead of trying to satisfy potentially stiff constraint equations.

Furthermore, to obtain an equation for the parallel electric field, we must go to higher order in our expansion to find the parallel force balance equation from the evolution of the electron momentum

(5.13) \begin{align} E_\parallel = -\frac {m_e}{e n_e} \boldsymbol{b} \boldsymbol{\cdot }\left [ \frac {\partial n_e \boldsymbol{u}_e}{\partial t} + \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (n_e \boldsymbol{u}_e\otimes \boldsymbol{u}_e \right ) + \frac {1}{m_e} \boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{\mathcal{P}}_e \right ] \approx -\frac {1}{e n_e} \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{\mathcal{P}}_e, \end{align}

assuming that the electron mass is small and, thus, the electron inertia is ignorable. Unfortunately, (5.13) is only a steady-state equation, and thus, its application to solving a time-dependent partial differential equation poses a non-trivial numerical difficulty. It need not be the case that (5.13) holds instantaneously, and enforcing that it does may be impossible without some sort of numerical operator that diffuses spurious fluctuations in the parallel electric field.

Similar to the feedback of perpendicular currents on the particle dynamics, while the equilibrium is contained at the desired order in the derivation of the equation system, the relaxation to this equilibrium is not. In the case of solving this constraint equation for the parallel electric field, in the physical system there will be the development of parallel currents due to a parallel electric field, which then must relax to a state of $J_\parallel \approx 0$ and a parallel electric field supported by the parallel pressure force. But then in a time-dependent partial differential equation, we may develop spurious fluctuations in the parallel electric field caused by non-zero parallel currents since this equation system contains a kinetic equation for every evolved species. We are thus not guaranteed to satisfy the identity $J_\parallel = 0$ and (5.13) instantaneously in time as each species independently develops parallel flows.

These two subtleties, the lack of feedback of the other guiding centre drifts on the kinetic response of the plasma and the parallel electric field equation only arising due to a steady-state parallel force balance equation, were principal motivations for the work of Ramos (Reference Quataert, Peterson, Pogge and Polidan2008) and Ramos (Reference Ramos2016) and the derivation of what Ramos (Reference Quataert, Peterson, Pogge and Polidan2008) referred to as FLR kinetic theory. By transforming to the total bulk flow frame, irrespective of the make up of that bulk flow, whether large parallel flows exist or all the guiding centre drifts are comparable in magnitude, all of the benefits of the KMHD formalism can be realised while eliminating the ambiguities of how to couple the fluid-kinetic system of equations. Not only the perpendicular electric field, but the parallel electric field too, are both eliminated from the kinetic equation via conservation of momentum. Furthermore, with a Fourier expansion in the velocity gyroangle, we can clearly identify both the evolution of the gyrotropic component of the distribution function and how the gyrotropic component couples to each subsequent Fourier harmonic as well as the impact of FLR effects on the plasma’s evolution.

In this regard, the PKPM formalism is the realisation of the formalism outlined in Ramos (Reference Quataert, Peterson, Pogge and Polidan2008) to transform the velocity coordinates to the total flow frame but, to our knowledge, never numerically implemented anywhere. By transforming the velocity coordinates of each distribution function to their particular bulk flow frame and determining that bulk flow from conservation of momentum of that particular plasma species, the difficulties in handling the macroscopic parallel electric fields, parallel currents and the evolution of the equilibrium magnetic field due to, e.g. diamagnetic effects, are all self-consistently captured by the interplay and coupling of Maxwell’s equations with each species’ momentum equations. Where the PKPM approach diverges from the approach outlined in Ramos (Reference Quataert, Peterson, Pogge and Polidan2008) is the avoidance of any asymptotic ordering of our transformed Vlasov equation, and we instead reduce the number of velocity degrees of freedom compared with the full Vlasov--Maxwell system through the spectral expansions in velocity gyroangle and perpendicular velocity. In lieu of expanding the transformed Vlasov equation in powers of $\epsilon = \rho _s/L$ , where $\rho _s$ is the plasma species’ gyroradius and $L$ is the gradient scale length, to include FLR effects to some order in $\epsilon$ , the PKPM model can simply add Fourier harmonics to some desired order, and we consider it likely that even the first Fourier harmonic is thus a super-set of the physics contained in the approach outlined in Ramos (Reference Quataert, Peterson, Pogge and Polidan2008). This argument is not to say that the PKPM model cannot benefit from asymptotic reductions; for example, if desired and appropriate, further numerical savings can be straightforwardly attained by employing reductions of Maxwell’s equations, such as the Darwin approximation to eliminate the speed of light as the fastest time scale in the problem (Schmitz & Grauer Reference Schekochihin and Cowley2006; Pezzi et al. Reference Parker and Dellar2019), because the field-particle couplings are entirely contained within the momentum equation.

The fact that the PKPM model as constructed can handle the interplay between all of the guiding centre drifts and develop macroscopic parallel electric fields from the self-consistent evolution of the kinetic equation and corresponding momentum equation can be contrasted with other KMHD-like models in the literature. For example, KMHD-like models such as the kglobal model (Arnold et al. Reference Arnold, Drake, Swisdak and Dahlin2019; Drake et al. Reference Dong, Wang, Hakim, Bhattacharjee, Slavin, DiBraccio and Germaschewski2019) avoid the numerical difficulties associated with the macroscopic parallel electric field by leveraging physics intuition about how the system attempts to restore parallel pressure balance and zero net parallel current, introducing a ‘cold’ electron fluid to instantaneously maintain quasi-neutrality and provide the expected return current that balances the local parallel flow of the ‘hot’, kinetic electrons. A total parallel pressure balance equation can then be defined to give the instantaneous parallel electric field similar to an Ohm’s law prescription for the perpendicular electric field. In reality of course, because one is solving a kinetic equation this return current should be self-consistently contained as a sub-population of electrons in the electron distribution function. So, while kglobal does successfully apply the KMHD formalism, this formalism is leveraged only under a prescribed model for how the system maintains net zero parallel current.

The approach in models like kglobal has been highly successful (Arnold et al. Reference Arnold, Drake, Swisdak, Guo, Dahlin, Chen, Fleishman, Glesener, Kontar, Phan and Shen2021, Reference Arnold, Drake, Swisdak, Guo, Dahlin and Zhang2022), a clear demonstration for why it has arguably been frustrating historically that KMHD resisted such easy discretisation. We thus argue the PKPM formalism presents a step forward in applying the same intuition that guided Kulsrud (Reference Krommes1964) and later Ramos (Reference Quataert, Peterson, Pogge and Polidan2008) to a diverse array of plasma systems. The dynamical interplay of the guiding centre drifts with the kinetic response of the plasma, the evolution of the macroscopic electromagnetic fields including macroscopic parallel electric fields, and the self-consistent treatment of the distinct particle populations that may arise to, e.g. drive the parallel currents to zero, are all handled by the PKPM formalism at reduced computational cost compared with solving the Vlasov equation in full generality.

5.2. Comparison to other (hybrid) spectral method approaches

Given the widespread popularity of spectral methods for velocity space discretisations of kinetic equations, from fully kinetic (Holloway Reference Hoffmann, Balestri and Ricci1996; Delzanno Reference Conley, Juno, TenBarge, Barbhuiya, Cassak, Howes and Lichko2015; Parker & Dellar Reference Ohia, Egedal, Lukin, Daughton and Le2015; Vencels et al. Reference Vencels, Delzanno, Manzini, Markidis, Peng and Roytershteyn2016; Roytershteyn & Delzanno Reference Rosin, Schekochihin, Rincon and Cowley2018; Koshkarov et al. Reference Kempski, Quataert, Squire and Kunz2021; Pagliantini et al. Reference Northrop2023; Issan et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2024; Schween & Reville Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2024; Schween et al. Reference Schekochihin, Cowley, Kulsrud, Rosin and Heinemann2025) to gyrokinetics (Mandell et al. Reference Mallet, Eriksson, Swisdak and Juno2018; Frei et al. Reference Francisquez, Bernard, Mandell, Hammett and Hakim2020; Hoffmann et al. Reference Hoffmann, Frei and Ricci2023b ; Frei et al. Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2024) to drift kinetics (Parker et al. Reference Ohia, Egedal, Lukin, Daughton and Le2016) and other reduced kinetic models (Zocco & Schekochihin Reference Zocco and Schekochihin2011; Loureiro et al. Reference Liu, Hesse, Guo, Daughton, Li, Cassak and Shay2013; Zocco et al. Reference Zocco, Loureiro, Dickinson, Numata and Roach2015; Loureiro et al. Reference Liu, Cassak, Li, Hesse, Lin and Genestreti2016), it is worth discussing how the PKPM approach compares to these other techniques for discretising velocity space. The key differences are three-fold: the direct optimisation of the spectral basis with our coordinate transformations that avoids the difficulty in handling time and spatially dependent shift and normalisation factors, the hybrid nature of the PKPM approach, which does not perform a spectral expansion in all the velocity degrees of freedom, and the lack of transformation of configuration space coordinates that distinguishes the gyroangle Fourier harmonic expansion from the gyroaveraging procedure in spectral gyrokinetic codes. All of these differences arise due to specific goals of the PKPM model; for example, as discussed in the beginning, we anticipate a spectral basis in $v_\parallel$ can perform quite poorly on the phase-space structures that typically arise in magnetised plasma dynamics, such as field-aligned beams or trapped particle distributions.

Likewise, we focused in § 2 on why we pursued a path of transforming the velocity coordinates of the Vlasov equation to optimise the spectral basis. We reiterate that not only do we avoid any assumptions on the temporal or spatial variation of the flow velocity with our coordinate transformation to move with the local flow velocity, but that the linear combination of Laguerre coefficients to produce the second kinetic equation for $\mathcal{G}$ , (4.47), in § 4 eliminates the need for an auxiliary equation for the time and spatially dependent Laguerre normalisation. Thus, we avoid the restriction in other spectral approaches that commonly assume a fixed and/or uniform shift and normalisation in the spectral expansion (Vencels et al. Reference Vencels, Delzanno, Manzini, Markidis, Peng and Roytershteyn2016; Koshkarov et al. Reference Kempski, Quataert, Squire and Kunz2021; Frei et al. Reference Francisquez, Juno, Hakim, Hammett and Ernst2023b , Reference Frei, Hoffmann, Ricci, Brunner and Tecchioll2024).Footnote ⁴ Furthermore, we show in § 6 test cases that illustrate the utility of not performing the spectral expansion in $v_\parallel$ due to the non-trivial structure that can arise in a variety of magnetised systems, such as magnetic reconnection, similar to other groups use of a hybrid approach mixing spherical harmonics with finite element methods to optimise their simulations of cosmic ray transport (Schween & Reville Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2024; Schween et al. Reference Schekochihin, Cowley, Kulsrud, Rosin and Heinemann2025).

We thus use this section to draw a particular contrast between the spectral gyrokinetic approach and the PKPM model for the handling of the Laguerre expansion. There is an important consequence from the transformation from particle position to gyrocentre position in spectral velocity representations of the gyrokinetic equation: the Laguerre representation of the Bessel function utilised to perform the gyroaveraging of various quantities, such as in the electrostatic potential, requires a sum over Laguerres (Zocco et al. Reference Zocco, Loureiro, Dickinson, Numata and Roach2015):

(5.14) \begin{align} J_0\left (k_\perp \rho _i v_\perp \right ) = e^{-\frac {1}{4}k_\perp ^2 \rho _i^2} \sum _{n=0}^{\infty } \frac {\left (k_\perp ^2 \rho _i^2 / 4\right )^n}{n!} L_n\left (v_\perp ^2\right ). \end{align}

Here, $J_0$ is the zeroth-order Bessel function of the first kind, $k_\perp$ is the wavenumber perpendicular to the magnetic field and $\rho _i$ is the ion gyroradius. An accurate treatment of gyroaveraging for large $k_\perp \rho _i$ thus requires an increasing number of Laguerre coefficients; see Appendix B of Mandell et al. (Reference Mallet, Eriksson, Swisdak and Juno2018). While accuracy improvements for low Laguerres have been utilised in the gyrofluid approach (Dorland & Hammett Reference Decker1993), and a number of nonlinear calculations have shown reasonable results in the fusion context, with works such as Hoffmann et al. (Reference Forslund and Shonk2023a ) and Mandell et al. (Reference Malkov and Drury2024) utilising as few as 3–4 Laguerres, in other cases such as the unstable entropy mode in Hoffmann et al. (Reference Francisquez, Juno, Hakim, Hammett and Ernst2023b ), more Laguerre resolution is needed to properly account for the FLR effects impact on the transport.

No such coupling of Laguerres occurs in the PKPM model because we have kept the configuration space coordinates untransformed; we obtain the distribution function at the particle position, not the gyrocentre position. In the PKPM formalism, each successive Fourier harmonic in velocity gyroangle will have its own Laguerre expansion. We can thus identify the physics of each Fourier harmonic Laguerre coefficient by Laguerre coefficient. For example, the first Fourier harmonic gives the evolution of $\boldsymbol{M}_\perp$ and the zeroth Laguerre coefficient of $\boldsymbol{M}_\perp$ can be utilised to obtain the component of the agyrotropic pressure tensor proportional to $\boldsymbol{M}_\perp$ in (4.19). Likewise, the first Laguerre coefficient of $\boldsymbol{M}_\perp$ can be utilised to obtain the heat flux of the perpendicular temperature perpendicular to the magnetic field in (4.20), the $v_\perp ^2 \boldsymbol{M}_\perp$ term. We defer a systematic comparison of the accuracy of FLR effects in the PKPM approach compared with a spectral gyrokinetic code to a future publication, but we note for now that because the PKPM model does not transform configuration space coordinates to gyrocentre coordinates or perform any gyroaveraging, the Laguerre couplings remain local at every $k_\perp$ , and thus, there is no obvious need for large Laguerre resolution at large $k_\perp$ in the PKPM approach. What specific impact this representation of the distribution function in perpendicular velocity and gyroangle, instead of gyrocentre, perpendicular velocity and gyrophase, has on the physics of, e.g. the entropy cascade (Schekochihin et al. Reference Rowan, Sironi and Narayan2009; Tatsuno et al. Reference St-Onge, Kunz, Squire and Schekochihin2009) is as yet undetermined, but at an initial glance the couplings that can require high Laguerre resolution are not present in the PKPM approach.

6.1. Parallel electrostatic shock

Shock waves, specifically collisionless shock waves, are omnipresent in our universe. In astrophysical systems where the collisional mean free path is large, the dynamics of the nonlinear wave steepening and rapid conversion of bulk kinetic energy into other forms of energy such as particle acceleration and heating in these shock waves are not mediated by collisions, but a myriad of collisionless processes such as wave--particle interactions and kinetic instabilities. To demonstrate the PKPM model, we consider the case of a parallel shock, where the incoming supersonic flow is aligned with the local magnetic field. In this shock geometry, on short time scales, there is the potential for electrostatic shocks driven by nonlinear steepening of ion acoustic modes, but only up to a critical Mach number (Forslund & Shonk Reference Fermi1970; Sorasio et al. Reference Shu2006). If the incoming flows are sufficiently large, there will be no slow down of the supersonic flow, and the fast plasma can propagate freely, smoothly transitioning around obstacles or interpenetrating the ambient medium through which the plasma is propagating.

To demonstrate that the PKPM model reproduces this transition from shocked flows to interpenetrating flows and, thus, contains a complete description of the collisionless physics of electrostatic shocks, we perform two simulations of colliding plasma flows, similar to previous studies of electrostatic shocks. We initialise an electron--proton plasma in a box of length $L_x = 3584 \lambda _D$ , where $\lambda _D = v_{th_e}/\omega _{pe}$ is the electron Debye length, $v_{th_e} = \sqrt {T_e/m_e}$ is the electron thermal velocity and $\omega _{pe} = \sqrt {e^2 n_0/\epsilon _0 m_e}$ is the electron plasma frequency. Here, $T_e$ is the electron temperature, $m_e$ is the electron mass, $e$ is the elementary charge, $n_0$ is the reference electron density and $\epsilon _0$ is the permittivity of free space. For the purposes of these simulations, all of these quantities are normalised to values of 1.0 so that all length scales are normalised with respect to the electron Debye length, all velocities are normalised with respect to the electron thermal velocity, and all time scales are normalised with respect to the inverse electron plasma frequency. This electron--proton plasma is initialised with Maxwellian distributions and a supersonic flow for both species in the negative $x$ direction towards a reflecting wall at $x=0$ , which leads to a shock wave or interpenetrating plasma propagating in the positive $x$ direction. A continuous supply of plasma is provided with a copy boundary condition at $x=L_x$ Footnote ⁵ in exact analogy to the ‘reflecting wall’ set-up commonly employed in particle-in-cell simulations of collisionless shocks and identical to the initial conditions of previous continuum simulations of collisionless shocks.

Specific shock parameters are as follows: we utilise the real proton--electron mass ratio $m_p/m_e = 1836$ , a proton--electron temperature ratio $T_p/T_e = 0.25$ , a reference magnetic field strength $\boldsymbol{B} = B_0 \hat {\boldsymbol{x}} = \hat {\boldsymbol{x}}$ so the magnetic field points in the $x$ direction the entirety of the simulation on the time scale of the electrostatic dynamics, and an electron--electron collisionality $\nu _{ee} = 10^{-6} \omega _{pe}$ , with the ion--ion collisionality commensurately decreased by the square root of the mass ratio and temperature ratio to the 3/2 power. We consider two Mach numbers, one below and one above the critical Mach number, where $M_s = u_{x_0}/c_s$ is the Mach number, $u_{x_0}$ is the magnitude of the upstream flow velocity and $c_s = \sqrt {T_e/m_p}$ is the sound speed. The two simulations have upstream Mach numbers of $M_s = 3.0$ and $M_s = 5.0$ . We choose to utilise a colder proton population to reduce the ion acoustic damping and, thus, generate stronger shock waves in the cases when the plasma produces a shock. The critical Mach number is $M_s \sim 3.0$ , but we note that this critical Mach number was determined for plasmas with only one-velocity dimension (Forslund & Shonk Reference Fermi1970), and the PKPM model is constructed to be three dimensional in velocity space. As such, we expect with this definition of the sound speed, without any additional $\mathcal{O}(1)$ factors for the adiabatic index of the plasma, that the $M_s = 3.0$ simulation will shock, while the $M_s = 5.0$ simulation will not shock.Footnote ⁶ This subtlety of the conditions under which a one-velocity dimensional plasma shocks compared with a three-velocity dimensional plasma has already been proven to be a use case for the PKPM model, as an early version of the PKPM model was utilised to understand measurements of shock formation, and lack of shock formation compared with fluid model predictions, in the plasma-Jet driven magneto-inertial fusion experiment (Cagas et al. Reference Braginskii2023).

The final input file specifications are our velocity space extents, grid resolutions and polynomial orders for the DG finite element method we utilise. The velocity grid extents are $[-8 v_{th_e}, 8 v_{th_e}]$ for the electrons and $[-64 v_{th_p}, 64 v_{th_p}]$ for the protons for the $M_s = 3.0$ simulation and $[-128 v_{th_p}, 128 v_{th_p}]$ for the protons for the $M_s = 5.0$ simulation.Footnote ⁷ We utilise $N_x = 1792$ grid points in configuration space, $\Delta x = 2 \lambda _D$ , with linear polynomials in configuration space, and $\Delta v_\parallel = 0.0625 v_{th_s}$ for both species with quadratic polynomials in velocity space, so that the electrons have $N_{v_\parallel } = 256$ and the protons have $N_{v_\parallel } = 2048$ and $N_{v_\parallel } = 4096$ velocity grid points for the $M_s = 3.0$ and $M_s = 5.0$ simulations, respectively. These simulations utilised 32 56-core compute nodes on the Frontera cluster at the Texas Advanced Computing Center; the $M_s = 3.0$ ran for half an hour for a total of 16 node hours, while the $M_s = 5.0$ used an hour of computing time for a total of 32 node hours.

Figure 1. Electron (left column) and proton (right column) mass density (top row), momentum density (middle top row), total energy density (middle bottom row) and parallel pressure (bottom row) at $t=1500\omega _{pe}^{-1}$ for upstream Mach numbers $M_s = 3.0$ (black) and $M_s = 5.0$ (red). All quantities are normalised to their upstream values for ease of comparison between the $M_s = 3.0$ and $M_s = 5.0$ cases since their upstream flows and energies are different. The characteristics of a shock wave are clearly identifiable in the $M_s = 3.0$ simulation: a sharp pile-up of the density, a rapid stagnation of the flow, significant electron heating over the same length scale and a decrease in the ion energy from the rapid conversion of ion energy into both electron heating and electromagnetic energy. On the other hand, the $M_s = 5.0$ case exhibits no such sharp transitions, with a smooth gradient up to a total mass density $\rho \sim 2$ and total momentum $\rho u_x \sim 0.0$ for both the electrons and protons, corresponding to two interpenetrating beams of plasma.

Currently we utilise large velocity space extents for the protons in shock simulations due to the significant transient generated by the large $\partial u_x/\partial x$ at early times, since with this reflecting wall set-up $\partial u_x/\partial x$ is infinite at $t=0, x = 0$ . The result of this large $\partial u_x/\partial x$ is the generation of a small amount of reflected particles at very high velocity in the hydrodynamic frame that propagate upstream and have no bearing on the evolution of the plasma. Once the transient has relaxed, the vast majority of the plasma is contained in a velocity space volume a quarter to an eighth in size, and we consider it worthwhile future work to explore an optimised initial set-up of shock simulations in this model to reduce the need to resolve this transient high velocity reflection. Nevertheless, we wish to emphasise that the fact that this numerical implementation can successfully evolve this transient stably is evidence of the robustness of this implementation, and thus, suggests that the PKPM model is not by any means a numerically brittle model despite its complexity.

We show in figure 1 the evolution of the fluid moments at $t=1500 \omega _{pe}^{-1}$ for the two upstream Mach numbers. All fluid moments are normalised to their upstream values to facilitate cross-comparison between the $M_s = 3.0$ and $M_s = 5.0$ simulation since the upstream flows and energies are different for these two simulations. We note the key features of the $M_s = 3.0$ indicative of a shock: a sudden density pile-up, a sharp stagnation of the flow in conjunction with this density pile-up, rapid electron heating with the electron parallel pressure increasing by nearly a factor of 15 compared with the density increasing by a factor of 3, and a commensurate decrease in the total ion energy and ion pressure corresponding to the conversion of ion energy into electron heating and electromagnetic energy. In contrast, the $M_s = 5.0$ simulation shows a smoother transition to the downstream region that is more characteristic of two interpenetrating beams of plasma: the mass density and ion energy density downstream are $\rho _s \sim 2, \mathcal{E}_p \sim 2$ , and so the two colliding plasma beams are simply adding their density and energy.

We can see further evidence of the shocked versus unshocked flows examining the distribution functions for these two upstream Mach numbers, specifically the $F_0$ coefficient corresponding to integrating the distribution function over both $v_\perp$ and $\theta$ . We plot in figure 2 the electron and proton $F_0$ distribution functions in both the local fluid flow frame provided by the numerical solution of the PKPM model, and in the lab frame as would traditionally be obtained from the solution of the Vlasov--Maxwell system of equations; see, for example, Juno et al. (Reference Imbert-Gerard, Paul and Wright2018) for simulations of electrostatic shocks with Gkeyll’s Vlasov–Maxwell solver. Note that these distribution function plots are normalised to their respective maximum values on the grid, e.g. $F_{0_s} = F_{0_s}/\max (F_{0_s})$ . For the $M_s = 3.0$ simulation, we can more clearly identify the characteristics of a shock, specifically the trapping of both electrons and protons in the downstream region around $x \sim 25 \lambda _D$ . On the other hand, the $M_s = 5.0$ simulation shows no significant broadening of the electron distribution function. Furthermore, the proton distribution function in the $M_s = 5.0$ simulation is simply two interpenetrating beams of protons at close to the upstream flow velocity of $u_{x_0} = 10 v_{th_p}$ . We draw particular attention to the local fluid flow frame distribution functions for the protons as evidence of the successful non-trivial numerical implementation of this model. We observe that as the reflected population propagates upstream, the fluid flow frame distribution function naturally adjusts due to the pressure forces to produce an approximately even distribution function with only high odd moments, such as heat fluxes. The first moment of the $F_0$ distribution function thus remains 0 as we expect.

Figure 2. The $F_0$ distribution function in the local fluid flow frame for the electrons (left column) and protons (left middle column), and the $F_0$ distribution function in the lab frame for the electrons (right middle column) and protons (right column) for the $M_s = 3.0$ (top row) and $M_s = 5.0$ (bottom row) simulations. In the lab frame, the incoming proton beam is centred at the upstream velocity, $u_{x_0} = 6.0 v_{\mathrm{th}_i}$ and $u_{x_0} = 10.0$ , as we expect, and the characteristics of the shock with the trapped electron and ion populations are identifiable in the $M_s = 3.0$ simulation, while the $M_s = 5.0$ simulation shows only two distinct ion beams propagating through each other. We also draw attention to the form of the proton distribution function in the local fluid flow frame and emphasise that these are the distribution functions that are directly solved for by the numerical method. As expected, the distribution function adjusts in the local fluid flow frame to preserve the identity that the first moment is zero. Note that these distribution function plots are normalised to their respective maximum values on the grid, e.g. $F_{0_s} = F_{0_s}/\max (F_{0_s})$ .

This initial electrostatic evolution only couples the $F_0$ and $\mathcal{G}$ kinetic equations through the collision operator; and in both of these simulations, but especially the $M_s = 3.0$ case where a true shock develops, a significant temperature anisotropy develops, with $T_\parallel \gt T_\perp$ . Therefore, as the shock propagates, electromagnetic instabilities can be excited that will generate large transverse fluctuations and change the shock from being a purely parallel shock to a quasi-parallel, or even locally quasi-perpendicular, shock. The fastest growing modes in this case, such as the electron firehose and Alfvén ion cyclotron instabilities, require FLR effects to properly capture their growth and saturation, and thus, only retaining the zeroth Fourier harmonic is inadequate for accurately simulating the transition from an electrostatic to an electromagnetic parallel shock (Gary et al. Reference Frieman and Chen2001). Nevertheless, recent work using extended fluid models to model the proton parallel firehose instability suggests that accounting for the agyrotropic components of the pressure tensor provides an accurate model for these temperature anisotropy-driven instabilities via the FLR effects approximated by evolving the full pressure tensor (Walters et al. Reference Walters, Klein, Lichko, Juno and TenBarge2024). Thus, adding one, or at most two, Fourier harmonics will be sufficient to capture the necessary FLR effects to simulate both where in wavenumber space the fastest growing modes exist and their overall saturation due to the generated magnetic field structure, making the PKPM approach a potentially powerful tool for modelling of collisionless shocks in these particular parameter regimes.

We note that at higher Mach numbers where significant particle acceleration occurs due to processes such as shock-drift acceleration (Paschmann et al. Reference Pagliantini, Manzini, Koshkarov, Delzanno and Roytershteyn1982; Sckopke et al. Reference Schmitz and Grauer1983; Anagnostopoulos & Kaliabetsos Reference Anagnostopoulos and Kaliabetsos1994; Anagnostopoulos et al. Reference Anagnostopoulos, Rigas, Sarris and Krimigis1998; Ball & Melrose 2001; Anagnostopoulos, Tenentes & Vassiliadis Reference Anagnostopoulos, Tenentes and Vassiliadis2009) and diffusive shock acceleration (Fermi Reference Egedal, Le and Daughton1949, Reference Ellison1954; Blandford & Ostriker Reference Berlok, Pakmor and Pfrommer1978; Ellison Reference Drake, Arnold, Swisdak and Dahlin1983; Blandford & Eichler Reference Ball and Melrose1987; Decker Reference Cockburn and Shu1988; Malkov & Drury Reference Liu, Hesse, Genestreti, Nakamura, Burch, Cassak, Bessho, Eastwood, Phan, Swisdak, Toledo-Redondo, Hoshino, Norgren, Ji and Nakamura2001), the PKPM model is likely to be inefficient for representing the kinetic response of the plasma. The particle distribution functions that result from these shock acceleration processes are typically highly agyrotropic, and would thus necessitate a large number of Fourier harmonics and Laguerre coefficients to resolve. But, we emphasise that in the case considered here, at relatively low Mach number, the relative inexpensiveness of these simulations means that, once the necessary FLR effects are included for correctly modelling the saturation of the excited temperature anisotropy instabilities, the multiscale nature of these subcritical, collisionless shocks that transition from electrostatic to electromagnetic shocks will be possible with the PKPM model. Furthermore, because the PKPM model is ultimately defined per species, with the transformation of the velocity coordinates to move with the local flow velocity employing that particular species flow velocity, we can imagine unique hybrid modelling approaches where a fully kinetic proton treatment could be coupled to a PKPM electron treatment, thus reducing the computational expense of modelling electrons at modest Mach numbers where the electrons stay magnetised through the shock.

6.2. Moderate guide-field magnetic reconnection

Equally ubiquitous to collisionless shocks as a mechanism by which plasma’s rearrange their energy budget is the phenomenon of magnetic reconnection, whereby magnetic fields change their topology to a lower energy state and transfer this excess energy to the plasma. In collisionless plasma systems where the resistivity is very small, the onset of magnetic reconnection is considered to be a fundamentally kinetic process. The plasma can only break field lines at microscopic length scales where the plasma demagnetises and the plasma particles are no longer constrained to follow the magnetic field.

Importantly though, magnetic reconnection is commonly observed componentwise; many plasma systems are undergoing what is referred to as ‘guide-field’ reconnection. In the geometry of these systems, we can divide the dynamical fields into a set of planar reconnecting components and a component perpendicular to the plane known as the guide field that the plasma particles will still try to ‘stick’ to even as the in-plane reconnecting components vanish at the point of magnetic reconnection where the magnetic field is changing its topology. We can thus ask: even with retaining only the zeroth Fourier harmonic and thus constraining the plasma to be gyrotropic, how well does the PKPM model capture guide-field reconnection?

A number of studies have shown gyrokinetic models of reconnection compare favourably to fully kinetic simulations in the limit that the guide field is strong (TenBarge et al. Reference Stone, Tomida, White and Felker2014; Muñoz et al. Reference Melville, Schekochihin and Kunz2015), $B_g \gg B_0$ , where $B_g$ is the magnitude of the guide field and $B_0$ is the magnitude of the reconnecting field. Alternatively, one can think of this limit as the $\delta B/B \ll 1$ limit -- a natural limit for standard derivations of gyrokinetics that assume a strong background field and only weak perturbations. But, we have continually emphasised that the PKPM approach is not an asymptotic one, and thus, we need not restrict ourselves to a strong guide-field limit. So long as there is a guide field of at least modest strength, do the plasma particles stay magnetised through the layer? The answer to this question is yes: a number of studies (Swisdak et al. Reference Squire, Schekochihin, Quataert and Kunz2005; Le et al. Reference Kunz, Schekochihin and Stone2013; Egedal, Le & Daughton Reference Dougherty2013) have shown that even at $B_g \approx 0.1{-}0.2 B_0$ , or $\delta B/B \sim 5{-}10$ , the electrons are still magnetised through the reconnection layer, provided the electrons are sufficiently light and a realistic proton--electron mass ratio is utilised.

We consider two moderate guide-field cases: $B_g = B_0$ and $B_g = 0.5 B_0$ , or $\delta B/B_g = 1$ and $\delta B/B_g = 2$ . We initialise a force-free current sheet in 2X1V, $(x,y,v_\parallel )$ , geometry of the form

(6.1) \begin{align} J_{x_e} & = -\frac {B_0}{w_0} \frac {\tanh ({y}/{w_0}) \textrm {sech}^2 ({y}/{w_0})}{\sqrt { ({B_g}/{B_0})^2 + \textrm {sech}^2 ({y}/{w_0})}},\\[-8pt]\nonumber \end{align}

(6.2) \begin{align} J_{z_e} & = \frac {B_0}{w_0} \textrm {sech}^2 \left (\frac {y}{w_0} \right ),\\[-8pt]\nonumber \end{align}

(6.3) \begin{align} B_x & = -B_0 \tanh \left (\frac {y}{w_0} \right ),\\[-8pt]\nonumber \end{align}

(6.4) \begin{align} B_z & = B_0 \sqrt {\left (\frac {B_g}{B_0} \right )^2 + \textrm {sech}^2 \left (\frac {y}{w_0} \right )}, \end{align}

where current is given entirely to the electrons to satisfy Ampere’s law. A $\sim 1\%$ Geospace Environmental Modelling (GEM)-like perturbation is utilised (Birn et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2001), along with perturbations to the first 20 wave modes with $\sim 1\%$ noise in $B_x, B_y,$ and $J_z$ to break the symmetry of the GEM-like perturbation and accelerate the development of the reconnection, similar to TenBarge et al. (Reference Stone, Tomida, White and Felker2014). These simulations utilise the standard GEM reconnection challenge proton--electron temperature ratio $T_p/T_e = 5$ , a realistic proton--electron mass ratio, $m_p/m_e = 1836$ , and a similar box size as the original GEM reconnection challenge: $L_x = 8 \pi d_p, L_y = 4 \pi d_p$ , where $d_p = c/\omega _{pp}$ is the proton inertial length and $\omega _{pp} = \sqrt {e^2 n_0/\epsilon _0 m_p}$ is the proton plasma frequency.

Other parameters are as follows: $\beta _e = 2\mu n_0 T_e/B_0^2 = 1/6$ defined in terms of the upstream in-plane magnetic field, velocity space extents $[-8 v_{th_s}, 8 v_{th_s}]$ for both the electrons and protons, $v_{th_e}/c = 1/16$ , $\nu _{ee} = 10^{-2} \varOmega _{cp}$ , where $\varOmega _{cp} = eB_0/m_p$ is the proton cyclotron frequency defined in terms of the upstream in-plane magnetic field strength. $\nu _{pp}$ is commensurately smaller by the square root of the mass ratio.Footnote ⁸ We show the results of three different simulations, $B_g = B_0$ and $B_g = 0.5 B_0$ both with $N_x \times N_y \times N_{v_\parallel } = 896 \times 448 \times 32$ corresponding to $\Delta x \sim 1.2 d_e, \Delta v_\parallel = 0.5 v_{th_s}$ , and one higher resolution simulation with $B_g = 0.5 B_0$ and $N_x \times N_y \times N_{v_\parallel } = 1792 \times 896 \times 32$ , corresponding to $\Delta x \sim 0.6 d_e$ . We utilise periodic boundary conditions in $x$ , a reflecting wall for the plasma and conducting walls for the electromagnetic fields in $y$ , and zero-flux boundary conditions in $v_\parallel$ . Like the parallel electrostatic shock simulations, these simulations also utilise linear polynomials in configuration space and quadratic polynomials in velocity space. We employ a small amount of hyperdiffusion in the momentum equation on scales comparable to the electron inertial length, $\nu _{Hyp} = 10^{-2} d_e^4$ for the lower resolution simulations and $\nu _{Hyp} = 10^{-3} d_e^4$ for the higher resolution simulation. The total cost of these simulations is relatively modest; the $896\,\times\,448$ resolution simulations required 24 hours on 128 56-core nodes on the Frontera cluster at the Texas Advanced Computing Center, ${\sim}3000$ node hours in total, and the higher resolution simulation cost the expected ${\sim}25\,000$ Frontera node hours from the doubling of the resolution and halving of the size of the time step.

Figure 3. Evolution of the out-of-plane current density, $J_z$ , with contours of the in-plane magnetic field superimposed by computing $A_z$ , the out-of-plane vector potential, from the in-plane $B_x$ and $B_y$ for the $B_g = B_0$ (left) and $B_g = 0.5 B_0$ (right) lower resolution simulations. We observe morphologies of the current layer consistent with Le et al. (Reference Kunz, Schekochihin and Stone2013), which found at lower electron $\beta _e$ a transition from a regime at a lower guide field in which an extended current layer forms from the magnetised electrons developing strong anisotropy and driving a perpendicular current across field lines, to a regime in which the magnetic tension in the guide field causes the current and density to peak near the diagonally opposed separator field lines and negate the impact of the electron anisotropy on the magnetic field’s overall tension (see figure 5). This contrast is especially clear at $t=20 \varOmega _{ci}^{-1}$ and $30 \varOmega _{ci}^{-1}$ as the reconnection rate reaches its peak values (see figure 7) and we can see a more concentrated current layer in the $B_g = B_0$ simulation compared with the extended current layer in the $B_g = 0.5 B_0$ simulation.

We show in figure 3 the evolution of the out-of-plane current density $J_z$ with contours of the in-plane magnetic field superimposed by computing $A_z$ , the out-of-plane vector potential, from the in-plane $B_x$ and $B_y$ for the lower resolution $896 \times 448$ simulations for the two guide fields. Irrespective of the guide field, we observe the typical reconnection morphology from the thinning of the current layer: the magnetic field lines pinch towards a centrally located X point where the field is reconnecting, and the lower energy state of the field’s rearranged topology leads to an acceleration of the plasma on either side of the X point. In fact, the $B_g = 0.5 B_0$ simulation also seems to be developing a secondary instability in the current sheet, and we plot in figure 4 a zoom in of the current sheet at both the lower resolution, $896 \times 448$ , larger hyperdiffusion simulation and the higher resolution, $1792 \times 896$ , smaller hyperdiffusion simulation.

Figure 4. Zoom in of the $B_g = 0.5 B_0$ simulation with lower resolution and larger hyperdiffusion (left), and higher resolution and smaller hyperdiffusion (right). While the mode is identifiable in the lower resolution simulation, the secondary instability is especially prominent at increased resolution.

Figure 5. Comparison of the electron parallel temperature normalised to the initial electron temperature (top), electron perpendicular temperature normalised to the initial electron temperature (middle top), electron temperature anisotropy (middle bottom) and electron firehose criteria (bottom) at $t = 20 \varOmega _{ci}^{-1}$ for the $B_g = B_0$ simulation (left) and $B_g = 0.5 B_0$ simulation (right). In both cases, a significant electron anisotropy from an excess of parallel pressure develops in the current layer, but a depletion of electron perpendicular pressure in the $B_g = 0.5 B_0$ simulation further increases the electron anisotropy in the layer. Combined with the lower guide field and, thus, a weaker magnetic field at the X point, the electron firehose criteria is much closer to marginal stability $\varLambda _{\textit{firehose}} \sim 0$ for the $B_g = 0.5 B_0$ simulation. Thus, the electrons more significantly modify the tension in the magnetic field at the reconnecting X point compared with the higher guide-field simulation, driving a perpendicular current that spreads the current layer into the exhaust.

We draw particular attention to details of the reconnection morphology in figures 3 and 4, which are consistent with previous fully kinetic simulations of guide-field reconnection at realistic mass ratio: at a lower guide field the current layer is extended into the exhaust (Le et al. Reference Kunz, Schekochihin and Stone2013). The explanation for this extended exhaust in the $B_g = 0.5 B_0$ simulation compared with the more peaked current density in the $B_g = B_0$ simulation is identical to the physics of the fully kinetic simulations in Le et al. (Reference Kunz, Schekochihin and Stone2013), as we show in figure 5. While both simulations self-consistently develop a reasonably large temperature anisotropy in the layer, the temperature anisotropy in the lower guide field, $B_g = 0.5 B_0$ , simulation is large enough that, combined with the lower magnitude guide field, the magnetic tension at the X point is reduced, driving a large perpendicular current that broadens the layer into the exhaust, $\boldsymbol{J}_\perp \sim (p_\perp - p_\parallel ) \boldsymbol{\nabla} _{\parallel } \boldsymbol{b} \times \boldsymbol{B}/|\boldsymbol{B}|^2$ , where $\boldsymbol{\nabla} _\parallel = \boldsymbol{b} \boldsymbol{\cdot }\boldsymbol{\nabla} _{\boldsymbol{x}}$ . We observe that the electron parallel firehose criterion (Li & Habbal Reference Le, Egedal, Daughton, Fox and Katz2000)

(6.5) \begin{align} \varLambda _{firehose} = 1 + \frac {\beta _{\parallel _e}}{2} \left (\frac {T_{\perp _e}}{T_{\parallel _e}} - 1 \right ), \end{align}

where $\beta _{\parallel _e} = 2 \mu _0 n_e T_{\parallel _e}/|\boldsymbol{B}|^2 = 2 \mu _0 p_{\parallel _e}/|\boldsymbol{B}|^2$ is the electron parallel plasma $\beta$ , is much closer to marginal stability, $\varLambda _{firehose} \sim 0$ in the $B_g = 0.5 B_0$ simulation. On the other hand, the modest values observed in the $B_g = B_0$ simulation correspond to only minor modifications of the magnetic tension due to the electron anisotropy; the larger guide field is able to maintain the overall tension in the field in the exhaust.

Figure 6. Evolution of the different components of the energy normalised to the total energy at $t=0$ including the total kinetic energy, $\rho _s |\boldsymbol{u}_s|^2/2$ , for the electrons and protons, the total internal energy, $3p_s/2 = p_{\parallel _s}/2 + p_{\perp _s}$ , for the electrons and protons, the electric field energy, $\epsilon _0 |\boldsymbol{E}|^2/2$ , and the magnetic field energy, $|\boldsymbol{B}|^2/2\mu _0$ , comparing both different guide fields (left) and different resolutions for the $B_g = 0.5 B_0$ simulation (right). We observe a conversion of magnetic energy into initially proton kinetic energy at the onset of magnetic reconnection, followed by heating of the plasma as both the electron and proton internal energies increase. Consistent with Shay et al. (Reference Sckopke, Paschmann, Bame, Gosling and Russell2014), we observe that the overall electron internal energy increase is relatively insensitive to the guide-field strength, and consistent with Rowan et al. (Reference Roberg-Clark, Drake, Reynolds and Swisdak2019), we observe that the relative heating of the protons versus the electrons is reduced at a larger guide field, as less magnetic energy is converted to plasma energisation in a stronger guide field for the moderate plasma $\beta$ case considered here.

The exact transition to this extended current layer when the firehose criterion becomes sufficiently small, but not necessarily unstable, is consistent with other studies that have examined the impact of proton pressure anisotropy on the overall magnetic field tension and the propagation of Alfvén waves in anisotropic plasmas (Bott et al. Reference Blandford and Eichler2021, Reference Blandford and Ostriker2025). In these studies, an effective Alfvén speed must be defined, which decreases with increasing pressure anisotropy, corresponding to a reduced capability of the magnetic field to regulate the motion of a magnetised plasma and enhanced perpendicular transport, usually modelled as a Braginskii-like viscous stress (Squire et al. Reference Squire, Kunz, Arzamasskiy, Johnston, Quataert and Schekochihin2017b ). In fact, at a lower guide field, we would likely observe an analogous phenomena to the nonlinear interruption of Alfvén waves (Squire, Quataert & Schekochihin Reference Spitzer2016, Reference Sovinec, Glasser, Gianakon, Barnes, Nebel, Kruger, Schnack, Plimpton, Tarditi and Chu2017a ), but from the electrons destabilising themselves due to the enhanced anisotropy in the extended layer. Such simulations could potentially extend the results of previous studies that have found the firehose criterion to be a constraint on the outflows in fully kinetic simulations of low guide-field reconnection (Haggerty et al. Reference Greene2018) into regimes of even moderate guide field where, as we find, the electrons are still magnetised through the current layer. We emphasise though that, as with the results of the parallel electrostatic shock, at least the first Fourier harmonic would be required to approximate the FLR effects that govern the saturation of firehose modes.

We next examine the macroscopic evolution of these simulations, plotting the overall energy evolution and reconnection rate in figures 6 and 7, respectively, for different guide fields and resolutions. The energy evolution of these moderate guide-field simulations is consistent with previous fully kinetic simulations: the overall heating of the electrons has only a weak dependence on guide-field strength (Shay et al. Reference Sckopke, Paschmann, Bame, Gosling and Russell2014), and as the guide-field strength increases, the amount of energy that the protons receive relative to the electrons decreases (Rowan, Sironi & Narayan Reference Roberg-Clark, Drake, Reynolds and Swisdak2019). The reconnection rate, computed as the time rate of change of the out-of-plane vector potential, ${\rm d}A_z/{\rm d}t$ , at the location of the maximum parallel electric field, $E_\parallel = \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{b}$ , peaks at the expected normalised value of $\sim 0.1$ (Shay et al. Reference Schween, Schulze and Reville1999; Liu et al. Reference Li and Habbal2017; Cassak, Liu & Shay Reference Cagas, Juno, Hakim, LaJoie, Chu, Langendorf and Srinivasan2017; Liu et al. Reference Le, Egedal, Ohia, Daughton, Karimabadi and Lukin2022), where the normalisation is defined in the standard way to the initial upstream, in-plane magnetic field strength multiplied by the initial upstream, in-plane Alfvén velocity. Additionally, outside of a slightly faster reconnection onset in the higher resolution $B_g = 0.5 B_0$ simulation, we find no sensitivity to our results with increasing resolution and lowering of the hyperdiffusion, suggesting that the kinetic response of the plasma is not modified by the hyperdiffusion model and the macroscopic dynamics of the reconnection are insensitive to this hyperdiffusion.

Figure 7. Reconnection rate as a function of time computed from the time rate of change of the out-of-plane vector potential, ${\rm d} A_z/{\rm d}t$ , at the location of maximum parallel electric field, $E_\parallel = \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{b}$ . Regardless of resolution or guide field, we observe a steady peak value of ${\sim}0.1$ in the standard normalised units dividing ${\rm d} A_z/{\rm d}t$ by the initial, upstream in-plane magnetic field strength multiplied by the initial, upstream in-plane Alfvén speed (Shay et al. Reference Schween, Schulze and Reville1999; Liu et al. Reference Li and Habbal2017; Cassak et al. Reference Cagas, Juno, Hakim, LaJoie, Chu, Langendorf and Srinivasan2017; Liu et al. Reference Le, Egedal, Ohia, Daughton, Karimabadi and Lukin2022).

Figure 8. Comparison of the reconnecting electric field, $E_z$ , and the individual components of the generalised Ohm’s law (6.6) for the $B_g = B_0$ (left) and $B_g = 0.5 B_0$ (right) simulations at $t = 20 \varOmega _{cp}^{-1}$ at a cut in $x$ through the current sheet approximately where the current density peaks in figure 3. As expected, the Hall term supports the electric field away from the current sheet, but the Hall term goes to zero neat the X point where both $u_{x_e}$ and $u_{y_e}$ go to zero from the stagnation of the flow. The reconnecting electric field’s dynamics is then governed by derivatives of the off-diagonal pressure tensor, which in this case is the gyrotropic electron pressure tensor. The combination of electron pressure anisotropy and the changing magnetic field geometry drive perpendicular currents that break the frozen-in condition as the electrons are no longer advecting with the magnetic field, thus allowing for the conversion of magnetic energy to plasma energy.

Figure 9. A zoom in of the $B_g = 0.5 B_0$ simulations with lower resolution and larger hyperdiffusion (left), and higher resolution and lower hyperdiffusion (right). We note that the overall layer width is only marginally affected by resolution and hyperdiffusion, $\Delta \sim 0.2 d_p \sim 8.5 d_e$ , where the proton and electron inertial lengths are defined with respect to the initial uniform density. The physics of the reconnecting electric field is completely unchanged qualitatively; the competition between $\partial _x P_{xz}$ and $\partial _y P_{yz}$ ultimately governs the reconnecting electric field evolution.

To further determine how accurately the PKPM model captures the kinetic response of the plasma at the X point, we can determine the physics of the out-of-plane electric field that governs the reconnection dynamics by rearranging the electron momentum equation:

(6.6) \begin{align} E_z = \frac {m_e}{q_e} \left (\underbrace {\frac {\partial u_{z_e}}{\partial t} + u_{x_e} \frac {\partial u_{z_e}}{\partial x} + u_{y_e} \frac {\partial u_{z_e}}{\partial y}}_{Inertia} + \frac {1}{\rho _e} \left [\frac {\partial P_{xz_e}}{\partial x} + \frac {\partial P_{yz_e}}{\partial y} \right ] \right ) - \underbrace {\left (u_{x_e} B_y - u_{y_e} B_x \right )}_{Hall}. \end{align}

Here, we have labelled components of this generalised Ohm’s law that correspond to the electron inertia and Hall terms (Vasyliunas Reference Vasyliunas1975). Note that because the reconnection is fairly steady at $t=20\varOmega _{cp}^{-1}$ (see figure 7) we will ignore $\partial u_{z_e}/\partial t$ in our analysis, as we expect this term to be small if the reconnection rate is relatively constant. We confirm in figures 8 and 9 the expected behaviour for this generalised Ohm’s law following the results of previous guide-field reconnection studies (Swisdak et al. Reference Squire, Schekochihin, Quataert and Kunz2005; Le et al. Reference Kunz, Schekochihin and Stone2013; Egedal et al. Reference Dougherty2013).

The electrons stay magnetised through the layer and, thus, provide the necessary support for the reconnection electric field through the electron’s off-diagonal pressure tensor components, $P_{xz} = (p_\parallel - p_\perp ) b_x b_z$ and $P_{yz} = (p_\parallel - p_\perp ) b_y b_z$ . In this case, the off-diagonal pressure tensor components that have historically been found to be the important components do not come from complex particle orbits but simply from the changing magnetic geometry and anisotropy of the electrons that develops as they stay magnetised through the layer.Footnote ⁹ Instead of the fundamental agyrotropy that develops in the zero or low guide-field case (Hesse et al. Reference Hammett and Perkins2018), it is the gyrotropic electron pressure tensor that can break the frozen-in condition (Egedal Reference Dorland and Hammett2002). In other words, the electrons no longer simply advect with the magnetic field, but directly change the magnetic field due the perpendicular currents driven by the gyrotropic electron pressure tensor. We note that this change in the field driven by perpendicular currents is not necessarily true field line ‘breaking’. figure 9 shows that while the overall current sheet width is relatively insensitive to resolution and hyperdiffusion, there is a region in which the sum of these various components of Ohm’s law and the reconnecting electric field diverge that does become smaller with increasing resolution. Thus, there is a region inside the current layer where the final field line topology changes, $b_x = b_y = 0$ , which is not captured in the lowest-order PKPM model with no Fourier harmonics.

Importantly, while a number of fluid models have been developed to exploit this fact that electron anisotropy can break the frozen-in condition and determine the current sheet morphology (Le et al. Reference Kunz, Jones and Zhuravleva2009; Ohia et al. Reference Nicastro, Mathur, Elvis, Drake, Fang, Fruscione, Krongold, Marshall, Williams and Zezas2012; Cassak et al. Reference Brizard and Hahm2015), we reiterate that this PKPM model employed here is fundamentally kinetic. The results of, e.g. how the electron anisotropy develops and affects the current sheet morphology are handled self-consistently without the need for any artificial limiters on the electron firehose condition (Ohia et al. Reference Northrop2015), and recent studies such as Walters et al. (Reference Walters, Klein, Lichko, Juno and TenBarge2024) suggest that with only one or two Fourier harmonics, a self-consistent saturation of firehose modes can be included with this approach. As it is, we can clearly identify the impact of the kinetic physics by examining higher velocity moments than typically included in fluid models, such as the total parallel heat flux plotted in figure 10 for electrons and protons, normalised to the initial thermal streaming values $\rho _s v_{th_s}^3$ , for the higher resolution, $1792 \times 896$ , $B_g = 0.5 B_0$ simulation at $t=20 \varOmega _{cp}^{-1}$ .

These heat fluxes are plotted componentwise so that we can separate the flow of $T_\parallel$ due to $q_\parallel$ and $T_\perp$ due to $q_\perp$ in the $x$ versus $y$ direction. These heat fluxes are not small compared with thermally streaming particles, especially $q_{\parallel }$ for both the electrons and ions. In other words, these collisionless heat fluxes are likely strongly influencing the dynamical evolution of $T_\parallel$ for each species. More importantly, $q_{\parallel _e}$ and $q_{\perp _e}$ are the opposite signs, towards the X point and away from the X point, respectively, in the extended current layer of this $B_g = 0.5 B_0$ simulation. These oppositely directed heat fluxes, in addition to the conservation of magnetic moment $\mu _e \propto T_{\perp _e}/B$ in the layer where the electrons are magnetised, further explain why there is a cooling in $T_{\perp _e}$ observed in figure 5 and that this cooling is more pronounced in the lower guide-field case where the magnetic field strength decreases more at the X point. Thus, the collisionless heat fluxes are enhancing the very temperature anisotropy in the layer that is modifying the tension in the magnetic field and driving the extended current layer.

Figure 10. Comparison of the different components of the parallel heat flux normalised to the initial values of a thermally streaming plasma, $\rho _e v_{th_e}^3$ , for both the electrons (left) and protons (right). We multiply both heat fluxes by the individual components of the magnetic field unit vector to separate the flow of $T_\parallel$ due to $q_\parallel$ and $T_\perp$ due to $q_\perp$ in the $x$ versus $y$ direction. We draw particular attention to how large the values of $q_{\parallel _e}$ are, suggesting that the exact evolution of $T_{\parallel _e}$ is strongly influenced by collisionless heat fluxes, and the fact that $q_{\parallel _e}$ and $q_{\perp _e}$ are the opposite sign in the current layer, which further explains the heightened temperature anisotropy in the $B_g = 0.5 B_0$ simulation. Furthermore, as large as the electron heat fluxes are, the ion heat fluxes compared with thermal streaming are not small either, with significant energy fluxes in the exhaust as the ions mix and heat in the outflows.

Figure 11. Zoomed-in out-of-plane current density for the $B_g = 0.5 B_0$ simulation with superimposed vectors denoting the direction of the in-plane flow direction and corresponding plots of the electron and proton distribution functions at the peak of the current density and inside one of the excited fluctuations from the secondary instability in the extended current sheet. Note that the distribution functions are plotted in their respective flow frames, so the $F_0$ distribution functions have a zero first moment, and any skewness or multi-modality to the $F_0$ distribution functions are manifestations of heat fluxes or counter-streaming field-aligned flows. While we need the sign of the magnetic field unit vector to determine which direction the vector heat flux points, we can clearly identify from the overall structure of the electron’s $F_0$ and $\mathcal{G}$ distribution functions that $q_{\parallel _e}$ and $q_{\perp _e}$ point in opposite directions, consistent with figure 10. In both cases, for the electrons, there is a depletion of $F_0$ relative to $\mathcal{G}$ for negative $v_\parallel '$ particles and an excess of $F_0$ relative to $\mathcal{G}$ for positive $v_\parallel '$ particles. Likewise, in the location of the excited secondary instability we have counter-streaming super-thermal beams along the field line, in addition to vortical structures in the in-plane flow velocity, which suggests that the electrons are unstable to a shearing instability. We identify this instability as the electron Kelvin–Helmholtz mode from its wavelength $k d_e \sim 1$ and its similarity to previous reconnection studies that excited secondary Kelvin–Helmholtz instabilities in the layer (Fermo et al. Reference Fabian2012).

In addition, at this same time of $t=20 \varOmega _{cp}^{-1}$ , we plot in figure 11 the zoom in of the out-of-plane current density, now with superimposed vectors denoting the direction of the in-plane flow, to examine the secondary instability that we observed to be developing in the extended current sheet in figure 4. We also show in figure 11 one-dimensional plots of the velocity distribution function within both the peak current density and one of the fluctuations of the secondary instability in the extended current layer. We identify this secondary instability now from its wavelength, $k d_e \sim 1$ , the vortical structure present in the in-plane velocity, and the counter-streaming superthermal electron beams as the electron Kelvin–Helmholtz instability. This instability has been previously measured in fully kinetic simulations at a larger guide field, $B_g = 2 B_0$ , and lower proton--electron mass ratio, $m_p/m_e = 25$ (Fermo, Drake & Swisdak Reference Fabian2012). The exact conditions under which this Kelvin–Helmholtz instability can be excited are likely a consequence of both mass ratio and guide field, as the excitement of this secondary instability seems most prevalent in these extended magnetised current layers and did not show up in the $B_g = B_0$ simulation with our chosen realistic proton--electron mass ratio.

The distribution function structure observed in this $B_g = 0.5 B_0$ simulation includes a number of other interesting features. While we need the sign of the magnetic field unit vector to obtain the sign of the vector heat flux, the depletion of negative velocity particles in $F_0$ compared with $\mathcal{G}$ , and likewise the excess of positive velocity particles in $F_0$ compared with $\mathcal{G}$ , easily explains the fact that $q_{\parallel _e}$ and $q_{\perp _e}$ have the opposite sign in figure 10. Furthermore, the broadened electron distribution function at the X point is consistent with trapped electron distribution functions in past studies of electron energisation in guide-field reconnection (Egedal et al. Reference Dougherty2013). Finally, the proton distribution functions also show clearly how the protons develop non-trivial heat fluxes, even if these heat fluxes are appreciably smaller than the electron heat fluxes.

Thus, even with a moderate guide field, because the electrons stay magnetised through the layer, the PKPM model provides a powerful cost-effective tool for kinetic simulations of magnetic reconnection. Self-consistent temperature anisotropies and parallel heat fluxes can develop that modify the current layer morphology in agreement with previous fully kinetic simulations. Secondary instabilities that may be excited via field-aligned beams are likewise included in this approach. Given the sensitivity of the kinetic evolution of the electrons to the proton--electron mass ratio, the PKPM model may provide a breakthrough in understanding the dynamics of guide-field reconnection, even in regimes beyond which magnetised models such as gyrokinetics can be applied.