Hamiltonian kinetic-Hall magnetohydrodynamics with fluid and kinetic ions in the current and pressure coupling schemes

Abstract

We present two generalized hybrid kinetic-Hall magnetohydrodynamics (MHD) models describing the interaction of a two-fluid bulk plasma, which consists of thermal ions and electrons, with energetic, suprathermal ion populations described by Vlasov dynamics. The dynamics of the thermal components are governed by standard fluid equations in the Hall MHD limit with the electron momentum equation providing an Ohm's law with Hall and electron pressure terms involving a gyrotropic electron pressure tensor. The coupling of the bulk, low-energy plasma with the energetic particle dynamics is accomplished through the current density (current coupling scheme; CCS) and the ion pressure tensor appearing in the momentum equation (pressure coupling scheme; PCS) in the first and the second model, respectively. The CCS is a generalization of two well-known models, because in the limit of vanishing energetic and thermal ion densities, we recover the standard Hall MHD and the hybrid kinetic-ions/fluid-electron model, respectively. This provides us with the capability to study in a continuous manner, the global impact of the energetic particles in a regime extending from vanishing to dominant energetic particle densities. The noncanonical Hamiltonian structures of the CCS and PCS, which can be exploited to study equilibrium and stability properties through the energy-Casimir variational principle, are identified. As a first application here, we derive a generalized Hall MHD Grad–Shafranov–Bernoulli system for translationally symmetric equilibria with anisotropic electron pressure and kinetic effects owing to the presence of energetic particles using the PCS.

1. Introduction

The presence of energetic, suprathermal particle populations, which can modify the global plasma dynamics owing to their high kinetic energies while having small density compared with that of thermal particles, is a common feature in astrophysical and fusion plasmas. For example, such energetic particles are part of the galactic cosmic rays (Amato & Blasi 2018) originating outside the solar system and penetrate the magnetospheric plasma. Additionally, charged particles are energized and accelerated during solar flares and coronal mass ejections, owing to magnetic reconnection and shock formation (Klein & Dalla 2017). Furthermore, it is known that in the magnetospheric ring-current plasma, energetic particles coexist with particles of significantly lower energy that constitute the plasma bulk (e.g. see Daglis et al. 1999). Although, the density of the thermal particles is dominant, the high-energy content of the fast particles renders even small populations capable of significantly affecting the plasma dynamics. Analogous situations, where suprathermal particles are produced and interact with a thermal plasma bulk, occur in fusion experiments by external plasma heating mechanisms, such as neutral beam injection (NBI) and ion cyclotron resonance heating (ICRH), which accelerate hydrogen isotopes and $^{3}He$ to energies of order 1 MeV (Start et al. 1999). Also, in a burning plasma, the deuterium–tritium (D–T) fusion reactions produce particularly energetic alpha particles (3.5 MeV) that are expected to heat the plasma, and consequently their impact on the dynamics of burning plasmas cannot be neglected. Energetic particles in fusion experiments are responsible for the destabilization of the Alfvén eigenmodes (AEs) (e.g. see Chen & Zonca 2007; Wesson 2011; Todo 2019) owing to resonant wave–particle interactions occurring when there exists a particle population with velocities near the phase velocity of the Alfvén wave. These interactions lead eventually to particle losses, which prevent the energetic particles from transferring their energy to the thermal plasma and thus deteriorate the heating efficiency.

The coexistence of a thermal or cold bulk and suprathermal populations has motivated the development of several hybrid multi-scale plasma models using fluid equations to describe the bulk plasma and Vlasov, or reduced kinetic equations (drift kinetic and gyrokinetic), to describe the energetic particle dynamics. Of course, a fully kinetic description using, for example, the Maxwel–Vlasov system, contains all the micro- and macro-physics involved in such systems. However, the hybrid fluid–kinetic description significantly reduces the computational cost because it is not required to simulate the dynamics of all particle species using kinetic equations. This makes the hybrid models important tools for performing numerical simulations and studying nonlinear dynamics because resolving all the kinetic scales of large systems with complex geometries is an extremely demanding task in terms of computational resources.

The hybrid models require a set of fluid and kinetic equations that should self-consistently describe the interaction of the plasma bulk with the energetic particles and the electromagnetic fields. For the fluid equations, the most popular choice is to consider an MHD description while for the kinetic component, one may use the Vlasov or reduced kinetic theories if the magnetic field is strong and the particle magnetic moment is an adiabatic invariant. To couple the dynamics of the plasma components, there are two main strategies: one relies on coupling through the current density, called the current coupling scheme (CCS), while the other is a pressure coupling scheme (PCS), where the kinetic effects are introduced through the pressure tensors appearing in the fluid momentum equation. The kinetic-MHD model in the PCS was introduced by Cheng (1991) and the first CCS kinetic-MHD, which considered gyrokinetic particles, was introduced by Park et al. (1992) to study the nonlinear behaviour of energetic particle effects. Since then, several hybrid drift kinetic and gyrokinetic MHD models have been employed to simulate plasma dynamics containing Alfvén eigenmodes, using either the CCS or the PCS, for example Briguglio et al. (1995), Todo & Sato (1998), Park et al. (1999) and Zhu, Ma & Wang (2016).

The CCS and PCS variants of the drift kinetic and gyrokinetic MHD models were formulated using Hamiltonian variational principles and Euler–Poincaré reduction (Burby & Tronci 2017; Close, Burby & Tronci 2018), and before that, a Hamiltonian approach to the hybrid description of plasmas combining the noncanonical Poisson bracket (Morrison & Greene 1980; Morrison 2009) of ordinary MHD with the particle bracket was developed by Tronci (2010) (see also Morrison, Tassi & Tronci 2014). This approach resulted in Vlasov kinetic-MHD models in the CCS and the PCS. The Hamiltonian construction of these models opened up the possibility to employ the energy-Casimir Hamiltonian variational principle (Morrison 1998) to derive equilibrium and stability conditions for planar kinetic-MHD using CCS and PCS in the studies by Morrison et al. (2014) and Tronci, Tassi & Morrison (2015). Tronci (2010) also provides the derivation of a noncanonical Poisson bracket that correctly describes the dynamics of a standard hybrid model that treats the electrons as a fluid with zero inertia while retaining a Vlasov description for the ions (see e.g. Winske et al. 2003). This model resolves the ion kinetic scales but not those of the electron, thus saving computational resources while reproducing the structural details of the reconnection region, which may consist of thin current sheets with thickness of the order of the ion inertial length. These characteristics render the model a popular choice for studying magnetic reconnection (Hesse & Winske 1994; Le et al. 2016; Cerri & Califano 2017).

Recently, a generalized, quasineutral hybrid model, which contained an additional fluid ion component while considering an arbitrary number of kinetic species, was presented by Amano (2018). This model makes no assumptions regarding the electron mass and employs the CCS for coupling the fluid and kinetic components. The generalized model of Amano (2018) incorporates the quasineutral two-fluid (QNTF) description and the common hybrid approach of purely kinetic ions and fluid electrons with finite or vanishing electron inertia. This unifying framework enables the study of kinetic effects in a continuous manner, starting from the fluid description and proceeding to situations with several kinetic and thermal species.

In this paper, we follow the approach of Amano (2018) to derive the corresponding model in the CCS for inertia-less electrons accompanied with its Hamiltonian formulation. Starting from Hamiltonian theories, such as the Hall MHD (Holm 1987; Lingam, Morrison & Miloshevich 2015) and the Vlasov equation, which has no collision operator, one expects to find a Hamiltonian structure for the resulting hybrid model. The identification of this structure might be important for the construction of structure preserving Hamiltonian algorithms that improve the stability and the fidelity of plasma simulations (Kraus et al. 2017; Morrison 2017). Such a structure-preserving code has been recently developed for the simulation of MHD waves that interact with energetic particles in the framework of hybrid kinetic-MHD (Holderied, Possanner & Wang 2021). Moreover, we obtain a second model upon considering the method of Hamiltonian construction of Tronci (2010) for the PCS. This new model has a set of Casimir invariants, i.e. global constants of motion whose gradients are elements of the Poisson kernel (Morrison 1998), which provide an interesting coupling between the fluid and the kinetic components, as was the case for the corresponding Hamiltonian kinetic-MHD model studied by Tronci et al. (2015). This coupling leads to novel equilibrium conditions upon employing the energy-Casimir variational principle (Morrison 1998). An additional feature of our study is that we consider an anisotropic, gyrotropic electron pressure tensor. Electron pressure anisotropy is a rather ubiquitous feature in guide-field magnetic reconnection, caused by electron trapping in the parallel electric fields (Egedal et al. 2008; Le et al. 2009; Egedal, Le & Daughton 2013), and can contribute to reconnection rates even for vanishing non-diagonal elements (Cassak et al. 2015). Also, we should note that for negligible electron inertia, a gyrotropic pressure tensor can be consistently considered, as was done by Ito, Ramos & Nakajima (2007). An analogous treatment for Amano's QNTF model would require the introduction of finite-Larmor-radius and gyroviscous effects, which is not trivial in the Hamiltonian framework. The inclusion of electron inertial effects, along with non-diagonal electron pressure elements, seems to be important to reproduce fully kinetic results regarding magnetic reconnection and relevant processes within the electron diffusion region (Muñoz et al. 2018; Finelli et al. 2021). The use of a Hamiltonian approach to hybrid models with finite electron inertia would open up the possibility to construct structure-preserving algorithms for simulating a large variety of phenomena where the electron inertial effects play an important role. Although this is left for future research, we note that approaches used in the past to bridge the Hall MHD and the extended MHD Hamiltonian structures could possibly be exploited here, for example, see the works of Lingam et al. (2015) and D'Avignon, Morrison & Lingam (2016). In terms of the pressure non-gyrotropy, we note that the early approach in the context of reduced fluid modelling by Hazeltine, Hsu & Morrison (1987), which was significantly generalized by Lingam, Morrison & Wurm (2020), shows that working with the Hamiltonian structure can enable the introduction of such effects.

The rest of the paper is organized as follows: in § 2, we present the parent equations, the method for constructing the CCS and the associated Hamiltonian structure. In § 3, we derive the novel hybrid model in the PCS using Tronci's method and we compare it to the model that would have been obtained using a standard (non-Hamiltonian) approach. In § 4, we derive the translationally symmetric counterpart of the PCS Hamiltonian structure and the associated Casimir invariants. We employ also the energy-Casimir variational principle that leads to equilibrium equations and § 5 summarizes our results.

2. Current coupling scheme

2.1. Model equations

The starting point for deriving the model equations in the CCS is the same as that used by Tronci (2010) and Amano (2018), i.e. a set of multifluid equations governing the dynamics of the thermal components accompanied by the Vlasov equation for the energetic particles, which provide self-consistent closure to Maxwell's equations:

(2.1)\begin{gather} \partial_t n_s+\boldsymbol{\nabla}\boldsymbol{\cdot} (n_s{\boldsymbol V}_s)=0, \end{gather}

(2.2)\begin{gather}m_sn_s(\partial_t{\boldsymbol V}_s+{\boldsymbol V}_s\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol V}_s)= e_sn_s({\boldsymbol E}+{\boldsymbol V}_s\times{\boldsymbol B})-\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_s, \end{gather}

(2.3)\begin{gather}\partial_t f_p={-}\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f_p-\frac{e_p}{m_p} ({\boldsymbol E}+\boldsymbol{v}\times{\boldsymbol B})\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} f_p, \end{gather}

(2.4)\begin{gather}\partial_t{\boldsymbol B}={-}\boldsymbol{\nabla}\times{\boldsymbol E}, \end{gather}

(2.5)\begin{gather}\partial_t{\boldsymbol E}=\epsilon_0^{{-}1}\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-\epsilon_0^{{-}1}{\boldsymbol J}, \end{gather}

(2.6a,b)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol E}=\sigma/\epsilon_0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol B}=0, \end{gather}

(2.7a–c)\begin{gather}\sigma= \sigma_f+\sigma_k,\quad \sigma_f=\sum_s e_s n_s, \quad \sigma_k=\sum_p e_p \int \, \textrm{d}^{3}v f_p, \end{gather}

(2.8a–c)\begin{gather}{\boldsymbol J}={\boldsymbol J}_f+{\boldsymbol J}_k,\quad {\boldsymbol J}_f=\sum_s e_s n_s {\boldsymbol V}_s , \quad {\boldsymbol J}_k=\sum_p e_p \int \, \textrm{d}^{3}v f_p \boldsymbol{v} , \end{gather}

where $\boldsymbol{\mathsf{P}}_s$ is the pressure tensor of the thermal species $s, f_p=f_p({\boldsymbol x},\boldsymbol {v},t)$ is the Vlasov distribution function, i.e. the particle density in phase space $({\boldsymbol x},\boldsymbol {v})$ for the particle species $p$. In this paper, we assume that the kinetic species consist of energetic ions, e.g. populations of alpha particles or resonant ions. Note that the system (2.1)–(2.8a–c) is not fully kinetic because it has been assumed that the moment hierarchy that yields the fluid equations (2.1) and (2.2) has been truncated in view of some appropriate fluid closure that results in independent expressions for the pressure tensors $\boldsymbol{\mathsf{P}}_s$ of the fluid species. For a two-fluid plasma bulk, the subscript $s$ is $s=i,e$. Note that $\boldsymbol {v}$ is the microscopic particle velocity, and therefore,

(2.9)\begin{equation} \frac{\partial v_i}{\partial x_j}=0,\quad \forall \, i,j. \end{equation}

In the low-frequency limit, it is legitimate to impose quasineutrality in a strict manner, i.e. setting $\sigma =0$. The quasineutrality assumption renders the displacement current term negligible and the Gauss law redundant. In addition, when the electron inertial effects can be ignored but the electron pressure is comparable to the magnetic pressure, the electron momentum equation results in the following generalized Ohm's law:

(2.10)\begin{equation} {\boldsymbol E}={-}{\boldsymbol V}_e\times {\boldsymbol B}-\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e}. \end{equation}

A rigorous derivation of (2.10) requires an appropriate ordering of the various terms in the electron fluid momentum equation. We may employ an Alfvén normalization by introducing the dimensionless quantities below:

(2.11)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \tilde{\boldsymbol{\nabla}}=\ell_0 \boldsymbol{\nabla},\quad \tilde{t}=\dfrac{t}{\ell_0/{V}_A},\quad \tilde{n}_s=\dfrac{n_s}{n_0},\\ \displaystyle \tilde{B}=\dfrac{B}{B_0},\quad \tilde{{V}}_s=\dfrac{{V}_s}{{V}_A}, \quad \tilde{\boldsymbol{\mathsf{P}}}_s=\dfrac{\boldsymbol{\mathsf{P}}_s}{B_0^{2}/\mu_0},\\ \displaystyle \tilde{J}=\dfrac{J}{B_0/(\ell_0\mu_0)}, \quad \tilde{E}=\dfrac{E}{B_0{V}_A}, \end{array}\right\} \end{equation}

where ${V}_A=B_0/\sqrt {\mu _0 m_in_0}$ is the Alfvén speed, and $B_0, n_0$ and $\ell _0$ are the characteristic magnetic field, number density and length, respectively. The electron momentum equation can then be written in the following non-dimensional form:

(2.12)\begin{equation} {\boldsymbol E}={-}{\boldsymbol V}_e\times{\boldsymbol B}-d_i\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{n_e}-\frac{d_e^{2}}{d_i}\frac{\textrm{d} {\boldsymbol V}_e}{\textrm{d} t}.\end{equation}

Here, the tildes have been dropped and the parameters $d_e$ and $d_i$ are the relative electron and ion skin depths, respectively. Neglecting ${{O}} (d_e^{2})$ terms (note that for an electron-proton plasma $d_e/d_i\sim 0.023$), thus eliminating electron length scales, the electron inertial term is removed, while the electron pressure tensor survives. Then, restoring dimensions in (2.12) yields (2.10). This result stems from the assumption that the electron pressure scales as the magnetic pressure and hence it is legitimate for plasmas with high electron $\beta$. In an alternative ordering scheme, e.g. for very low $\beta$, the order of the electron pressure term would be lower than the order of the Hall term and it could thus be neglected.

As regards the form of $\boldsymbol{\mathsf{P}}_e$, because of the specific ordering described above, we may consider a gyrotropic electron pressure. To see this, let us write the electron pressure equation as it emerges from the moment hierarchy of the Vlasov equation (e.g. see Hunana et al. 2019),

(2.13)\begin{align} & m_e[\partial_t \boldsymbol{\mathsf{P}}_e+\boldsymbol{\nabla}\boldsymbol{\cdot} ({\boldsymbol V}_e\boldsymbol{\mathsf{P}}_e + \boldsymbol{\mathsf{Q}}_e)+\boldsymbol{\mathsf{P}}_e\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}_e+(\boldsymbol{\mathsf{P}}_e\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}_e)^{\top}]\nonumber\\ & \quad =e[{\boldsymbol B}\times \boldsymbol{\mathsf{P}}_e+({\boldsymbol B}\times\boldsymbol{\mathsf{P}}_e)^{\top}], \end{align}

where the superscript $\top$ denotes transpose and $\boldsymbol{\mathsf{Q}}_e$ is the electron heat flux tensor. Performing the Alfvén normalization (2.11), the electron pressure equation becomes

(2.14)\begin{align} & \frac{d_e^{2}}{d_i}[\partial_t \boldsymbol{\mathsf{P}}_e+\boldsymbol{\nabla}\boldsymbol{\cdot} ({\boldsymbol V}_e\boldsymbol{\mathsf{P}}_e + \boldsymbol{\mathsf{Q}}_e)+\boldsymbol{\mathsf{P}}_e\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}_e+(\boldsymbol{\mathsf{P}}_e\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}_e)^{\top}]\nonumber\\ & \quad =[{\boldsymbol B}\times \boldsymbol{\mathsf{P}}_e+({\boldsymbol B}\times\boldsymbol{\mathsf{P}}_e)^{\top}]. \end{align}

Neglecting electron length scales $(d_e\rightarrow 0)$, only the right-hand side of (2.13) survives, and thus $\boldsymbol{\mathsf{P}}_e$ satisfies ${\boldsymbol B}\times \boldsymbol{\mathsf{P}}_e+({\boldsymbol B}\times \boldsymbol{\mathsf{P}}_e)^{\top }=0$. A general solution of this equation is given by the gyrotropic pressure tensor

(2.15)\begin{equation} \boldsymbol{\mathsf{P}}_e=\frac{P_{e\parallel}-P_{e\perp}}{B^{2}} {\boldsymbol B}{\boldsymbol B}+P_{e\perp}\boldsymbol{\mathsf{I}},\end{equation}

where

(2.16)\begin{equation} \left. \begin{array}{c@{}} P_{e\parallel}=\boldsymbol{\mathsf{P}}_e\boldsymbol{:}{\boldsymbol b}{\boldsymbol b},\\ P_{e\perp}=\tfrac{1}{2}\boldsymbol{\mathsf{P}}_e\boldsymbol{:} (\boldsymbol{\mathsf{I}}-{\boldsymbol b}{\boldsymbol b}). \end{array}\right\} \end{equation}

Here, ${\boldsymbol b}:={\boldsymbol B}/|{\boldsymbol B}|$ and $\boldsymbol{\mathsf{A}}\boldsymbol {:}\boldsymbol{\mathsf{B}}$ indicates double contraction between the second-order tensors $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$, i.e. $\boldsymbol{\mathsf{A}}\boldsymbol {:}\boldsymbol{\mathsf{B}}=A_{ij}B_{ij}$. We consider this specific form of electron pressure throughout the rest of the paper, because it is not only legitimate in the small electron length-scale limit, but, as will be seen in the next section, facilitates also the identification of an appropriate Hamiltonian structure.

Now, using $\boldsymbol {\nabla }\times {\boldsymbol B}=\mu _0{\boldsymbol J}$ and (2.8a–c), the electron velocity can be written as

(2.17)\begin{equation} {\boldsymbol V}_e=\frac{n_i}{n_e}{\boldsymbol V}-\frac{1}{en_e}\left(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-\sum_p e_p\int \, \textrm{d}^{3}v\, \boldsymbol{v} f_p\right), \end{equation}

where ${\boldsymbol V}\equiv {\boldsymbol V}_i$. Inserting (2.17) in the Ohm's law (2.10), we obtain

(2.18)\begin{equation} {\boldsymbol E}={-}\frac{n_i}{n_e}{\boldsymbol V}\times{\boldsymbol B}+\frac{1}{en_e}(\mu_0^{{-}1}\boldsymbol{\nabla}\times {\boldsymbol B}-{\boldsymbol J}_k)\times{\boldsymbol B}-\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e}.\end{equation}

Note that owing to the energetic particle component, quasineutrality does not imply $n_i = n_e$, but $n_i=n_e-e^{-1}\sigma _k$. In view of (2.18), Faraday's law (2.4) leads to the following induction equation:

(2.19)\begin{equation} \partial_t{\boldsymbol B}=\boldsymbol{\nabla}\times\left[\frac{n_i}{n_e}{\boldsymbol V}\times{\boldsymbol B}- \frac{1}{en_e}(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-{\boldsymbol J}_k) \times{\boldsymbol B}+\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e}\right]. \end{equation}

Inserting the generalized Ohm's law (2.18) into the ion momentum equation (2.2) and the Vlasov equation (2.3), we obtain respectively

(2.20)\begin{gather} mn_i(\partial_t{\boldsymbol V}+{\boldsymbol V}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol V})= \left[\frac{n_i}{n_e}\sigma_k {\boldsymbol V} +\frac{n_i}{n_e} ({\boldsymbol J}-{\boldsymbol J}_k )\right]\times{\boldsymbol B}-\frac{n_i}{n_e} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e-\boldsymbol{\nabla} P_i, \end{gather}

(2.21)\begin{gather}\partial_t f_p+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f_p+\frac{e_p}{m_p} \left\{\left[\boldsymbol{v}-\frac{n_i}{n_e}{\boldsymbol V}+\frac{1}{en_e} ({\boldsymbol J}-{\boldsymbol J}_k)\right]\times{\boldsymbol B} -\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e} \right\}\boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{v}} f_p=0, \end{gather}

where, from now on, $m\equiv m_i$. Here, we have assumed that the thermal ions have isotropic pressure. For a magnetized plasma, a consistent consideration of thermal pressure effects for the ions within the two-fluid framework requires the inclusion of anisotropic (and in particular non-gyrotropic) pressure effects. This will be the topic of future research. Presently, it suffices for our purposes to consider an isotropic pressure $P_i$ assuming low temperature ions because thermal ions that deviate from the isotropic pressure description can be incorporated in the kinetic component described by the Vlasov equation.

The continuity equation for the ions remains unchanged, while for the electron fluid, inserting (2.17) into (2.1), we find

(2.22)\begin{equation} \partial_t n_e={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (n_i{\boldsymbol V})-\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol J}_k}{e}, \end{equation}

that is,

(2.23)\begin{equation} e\partial_t(n_i-n_e)=\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol J}_k, \end{equation}

which is an equation for electric charge conservation. The complete system of the dynamical equations consists of the ion continuity equation, the induction equation (2.19), the momentum equation (2.20), the Vlasov equation (2.21), and the electron continuity equation (2.22) or equivalently the charge conservation equation (2.23). Computing the first-order velocity moment of the Vlasov equation, one finds

(2.24)\begin{equation} \partial_t \sigma_k+\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol J}_k=0, \end{equation}

which, combined with (2.23), gives the local charge conservation. Thus, for $e(n_{i}-n_{e})+\sigma _{k}=0$ at $t=0$, the quasineutrality constraint, which is invoked in various derivations in this paper, is preserved by the dynamics.

2.2. Hamiltonian structure

Let us now follow the procedure introduced by Tronci (2010) to derive the Hamiltonian structure of the above model. The difference here is that we apply the procedure in a more general model, by considering a two-fluid bulk plasma and generalizing also for electron pressure anisotropy. The starting point for this derivation is the combination of the Hall MHD noncanonical bracket, derived by Holm (1987), with the particle Poisson bracket. This can be done by directly adding the two brackets, which are written, however, in terms of the canonical momentum density $\bar {{\boldsymbol M}}:=\rho {\boldsymbol V}+({e}/{m})\rho {\boldsymbol A}$, where $\rho :=m n_i$, and the canonical particle momenta ${\boldsymbol \pi }_p=m_p\boldsymbol {v}+e_p{\boldsymbol A}$, where ${\boldsymbol A}$ is the vector potential.

Before we proceed to this construction, let us recapitulate here some basic notions of noncanonical Hamiltonian dynamics. The Lagrange–Euler map from the Lagrangian to the Eulerian description, renders the Poisson brackets explicitly dependent on the dynamical variables of the system, say $\boldsymbol {\xi }$, i.e. it has the form

(2.25)\begin{equation} \{F,G\}=\langle F_{\xi_i}, {\mathcal{J}} _{ij}(\boldsymbol{\xi}) G_{\xi_j}\rangle, \end{equation}

where $F_\xi$ represents the functional derivative of $F$ with respect to $\xi$ and ${\mathcal {J}} (\boldsymbol {\xi })$ is the so-called Poisson operator.

This noncanonical Poisson bracket still satisfies the antisymmetry condition and the Jacobi identity:

(2.26)\begin{gather} \{F,G\}={-}\{G,F\}, \end{gather}

(2.27)\begin{gather}\{F,\{G,H\}\}+\{H,\{F,G\}\}+\{G,\{H,F\}\}=0, \end{gather}

where $F,G,H$ are functionals that are defined on the functional phase space. Owing to the explicit dependence on the dynamical variables, there exist non-trivial functionals ${\mathcal {C}}$ satisfying

(2.28)\begin{equation} \{ {\mathcal{C}} ,F\}=0,\quad \forall \, F. \end{equation}

These functionals ${\mathcal {C}}$ are called Casimirs. The equations of motion result from the following Hamilton equations:

(2.29)\begin{equation} \partial_t \xi_i=\{\xi_i,{\mathcal{H}}\}, \end{equation}

where ${\mathcal {H}}$ is the Hamiltonian of the model. Obviously, the Hamiltonian is conserved in view of the antisymmetry of the Poisson bracket and the Casimirs are constants owing to their commutative property. Such Hamiltonian structures have been identified in the context of fluid mechanics and ordinary magnetohydrodynamics (Morrison & Greene 1980; Morrison 1982,  1998), Hall MHD (Holm 1987; Lingam et al. 2015), extended-MHD (Abdelhamid, Kawazura & Yoshida 2015; Lingam et al. 2015), and Vlasov–Poisson and Vlasov–Maxwell theories (Morrison 1980) (with a correction for Vlasov–Maxwell in the papers by Weinstein & Morrison (1981), Marsden & Weinstein (1982), Morrison (1982) and a limitation to the correction pointed out by Morrison (1982), which was followed up more recently by Morrison (2013); Heninger & Morrison (2020) and then by Lainz, Sardón & Weinsten (2019)). For a more comprehensive presentation of the noncanonical Hamiltonian dynamics emerging in the Eulerian description of fluids and other continuum theories, the reader is referred to the paper by Morrison (1998). For Vlasov models, the interested reader can consult the paper by Morrison (2009).

Regarding our model, the sum of the Hall MHD and the particle brackets expressed in terms of canonical momenta is given by

(2.30)\begin{align} \{F,G\}& =\int \, \textrm{d}^{3}x\left\{\vphantom{\sum_p} \bar{\boldsymbol{M}}\boldsymbol{\cdot} (G_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{\bar{\boldsymbol{M}}}- F_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{\bar{\boldsymbol{M}}})\right.\nonumber\\ & \quad+\rho(G_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{\rho}-F_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{\rho}) -e^{{-}1}(G_{{\boldsymbol A}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_e}-F_{{\boldsymbol A}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_e})\nonumber\\ & \quad-\left.\frac{1}{en_e}(\boldsymbol{\nabla}\times{\boldsymbol A})\boldsymbol{\cdot} (F_{{\boldsymbol A}}\times G_{{\boldsymbol A}})+ \sum_p \int \, \textrm{d}^{3}{\boldsymbol \pi} f_p [\![F_{\bar{f}_p},G_{\bar{f}_p}]\!]_{{\boldsymbol \pi}_p} \right\}, \end{align}

where

(2.31)\begin{equation} [\![g,h]\!]_{{\boldsymbol \pi}_p}=\boldsymbol{\nabla} g \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\boldsymbol \pi}_p} h-\boldsymbol{\nabla} h \boldsymbol{\cdot} \boldsymbol{\nabla}_{{\boldsymbol \pi}_p} g, \end{equation}

and $\bar {f}_p$ is the distribution function expressed in terms of ${\boldsymbol \pi }_p$, i.e. $\bar {f}_p=\bar {f}_p({\boldsymbol x},{\boldsymbol \pi }_p,t)$.

The above Poisson bracket, along with the appropriate Hamiltonian, which is the sum of the Hall MHD and particle Hamiltonians expressed in terms of $\bar {\boldsymbol {M}}, {\boldsymbol \pi }$, describe correctly the dynamics of the system. However, the canonical momentum variables are not convenient because they mix up the velocity and the magnetic field potential. Ultimately, we would like to have a system of equations that treats velocity and magnetic field variables separately, i.e. like the system we derived by the standard approach in the previous subsection. To this end, following Tronci (2010) and Marsden & Weinstein (1982),we may perform the following change of variables:

(2.32a,b)\begin{equation} \bar{\boldsymbol{M}}\rightarrow {\boldsymbol V},\quad \bar{f}_p({\boldsymbol x},{\boldsymbol \pi}_p,t)=\bar{f}_p({\boldsymbol x},m_p\boldsymbol{v}+ e_p{\boldsymbol A},t)\rightarrow f_p({\boldsymbol x},\boldsymbol{v},t), \end{equation}

where

(2.33)\begin{equation} {\boldsymbol V}=\rho^{{-}1}\bar{\boldsymbol{M}}-\frac{e}{m}{\boldsymbol A}. \end{equation}

For convenience, this change of variables will be performed in two stages. First, we can write the bracket in terms of $f_p$ and then we can complete the change transforming from $\bar {\boldsymbol {M}}$ to ${\boldsymbol V}$. For the first change, we can directly use a result from Marsden & Weinstein (1982), that is, for any functional $F[{\boldsymbol A},f_p]$ and functions $g({\boldsymbol x},\boldsymbol {v}), h({\boldsymbol x},\boldsymbol {v})$, we have $F[{\boldsymbol A},f_p({\boldsymbol x},\boldsymbol {v},t)]=\bar {F}[{\boldsymbol A},\bar {f}_p({\boldsymbol x},{\boldsymbol \pi }_p,t)], g({\boldsymbol x},\boldsymbol {v})=\bar {g}({\boldsymbol x},{\boldsymbol \pi }_p), h({\boldsymbol x},\boldsymbol {v})=\bar {h}({\boldsymbol x},{\boldsymbol \pi }_p)$, and the following transformation rules must hold:

(2.34)\begin{gather} [\![\bar{g},\bar{h}]\!]_{{\boldsymbol \pi}_p}= \frac{1}{m_p}[\![g,h]\!]+\frac{e_p}{m_p^{2}}(\boldsymbol{\nabla}\times{\boldsymbol A}) \boldsymbol{\cdot} (\nabla_{\boldsymbol{v}} g\times\nabla_{\boldsymbol{v}} h), \end{gather}

(2.35)\begin{gather}\bar{F}_{{\boldsymbol A}}=F_{{\boldsymbol A}}-\sum_p \frac{e_p}{m_p}\int \, \textrm{d}^{3}v f_p\nabla_{\boldsymbol{v}} F_{f_p}, \end{gather}

where

(2.36)\begin{equation} [\![g,h]\!]:=\boldsymbol{\nabla} g \boldsymbol{\cdot} \nabla_{\boldsymbol{v}} h-\boldsymbol{\nabla} h \boldsymbol{\cdot} \nabla_{\boldsymbol{v}} g. \end{equation}

Using (2.34) and (2.35), the bracket (2.30) takes the following form:

(2.37)\begin{align} \{F,G\}& =\int \, \textrm{d}^{3}x \left\{\vphantom{\sum_p}\bar{\boldsymbol{M}}\boldsymbol{\cdot} (G_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{\bar{\boldsymbol{M}}}-F_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{\bar{\boldsymbol{M}}})\right.\nonumber\\ & \quad +\rho(G_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{\rho}-F_{\bar{\boldsymbol{M}}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{\rho}) -e^{{-}1}(G_{{\boldsymbol A}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_e}-F_{{\boldsymbol A}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_e})\nonumber\\ & \quad -\frac{1}{en_e}(\boldsymbol{\nabla}\times{\boldsymbol A})\boldsymbol{\cdot} (F_{{\boldsymbol A}}\times G_{{\boldsymbol A}})-e^{{-}1} \sum_p \frac{e_p}{m_p}f_p[\boldsymbol{\nabla} G_{n_e}\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} F_{f_p}-\boldsymbol{\nabla} F_{n_e} \boldsymbol{\cdot} \nabla_{\boldsymbol{v}} G_{f_p} ] \nonumber\\ & \quad +\frac{1}{en_e}\sum_p\frac{e_p}{m_p}f_p (\boldsymbol{\nabla}\times{\boldsymbol A})\boldsymbol{\cdot} (\nabla_{\boldsymbol{v}} F_{f_p}\times G_{{\boldsymbol A}}-\nabla_{\boldsymbol{v}} G_{f_p}\times F_{{\boldsymbol A}})\nonumber\\ & \quad -\frac{1}{e n_e}(\boldsymbol{\nabla}\times {\boldsymbol A})\boldsymbol{\cdot} \int \int \, \textrm{d}^{3}v\, \textrm{d}^{3}v' \left(\sum_p \frac{e_p}{m_p}f_p\nabla_{\boldsymbol{v}} F_{f_p}\right)\times\left(\sum_{p'} \frac{e_{p'}}{m_{p'}}f_{p'}\boldsymbol{\nabla}_{\boldsymbol{v}'} F_{f_{p'}}\right)\nonumber\\ & \quad +\left.\sum_p \int \, \textrm{d}^{3}v \frac{f_p}{m_p} \left[[\![F_{f_p},G_{f_p}]\!]+ \frac{e_p}{m_p}(\boldsymbol{\nabla}\times {\boldsymbol A})\boldsymbol{\cdot} \left(\nabla_{\boldsymbol{v}} F_{f_p}\times\nabla_{\boldsymbol{v}} G_{f_p} \right) \right]\right\}. \end{align}

Now, we can proceed by expressing this bracket in terms of ${\boldsymbol V}$ and ${\boldsymbol B}=\boldsymbol {\nabla }\times {\boldsymbol A}$, instead of $\bar {{\boldsymbol M}}$ and ${\boldsymbol A}$. Let us note first that from (2.33), an arbitrary variation of ${\boldsymbol V}$ can be written as

(2.38)\begin{equation} \delta{\boldsymbol V}=\rho^{{-}1}\delta \bar{{\boldsymbol M}}-\rho^{{-}2}\bar{{\boldsymbol M}}\delta\rho -\frac{e}{m}\delta{\boldsymbol A}. \end{equation}

Considering a functional $\bar {F}$ of $\bar {{\boldsymbol M}}, {\boldsymbol A}$ and $\rho$, and then expressing it in terms of ${\boldsymbol V}, {\boldsymbol B}$ and $\rho$, the equality $\bar {F}[\bar {{\boldsymbol M}},{\boldsymbol A},\rho ]=F[{\boldsymbol V},{\boldsymbol B},\rho ]$ must hold. Upon taking the first variation of this relation and using the chain rule for the functional derivatives, the following relations can be deduced:

(2.39)\begin{gather} \bar{F}_{\bar{{\boldsymbol M}}}=\rho^{{-}1}F_{{\boldsymbol V}}, \end{gather}

(2.40)\begin{gather}\bar{F}_\rho= F_\rho-\rho^{{-}1} \left({\boldsymbol V}+\frac{e}{m}{\boldsymbol A}\right)\boldsymbol{\cdot} F_{{\boldsymbol V}}, \end{gather}

(2.41)\begin{gather}\bar{F}_{{\boldsymbol A}}= \boldsymbol{\nabla}\times F_{{\boldsymbol B}}-\frac{e}{m}F_{{\boldsymbol V}}. \end{gather}

Substituting (2.33) and (2.39)–(2.41) to (2.37), we find the final bracket, which reads¹

(2.42)\begin{align} \{F,G\}& = \int \, \textrm{d}^{3}x\left\{\vphantom{\sum_p} (G_{\boldsymbol V}\boldsymbol{\cdot} \nabla F_\rho -F_{\boldsymbol V}\boldsymbol{\cdot} \nabla G_\rho)+ \rho^{{-}1} (\boldsymbol{\nabla}\times{\boldsymbol V})\boldsymbol{\cdot} (F_{\boldsymbol V}\times G_{\boldsymbol V})\right.\nonumber\\ & \quad +\frac{e}{m^{2}}\left(\frac{mn_e-\rho}{\rho n_e}\right){\boldsymbol B}\boldsymbol{\cdot} (F_{\boldsymbol V}\times G_{\boldsymbol V})+m^{{-}1} (G_{\boldsymbol V}\boldsymbol{\cdot} \nabla F_{n_e}-F_{\boldsymbol V}\boldsymbol{\cdot} \nabla G_{n_e})\nonumber\\ & \quad + \frac{1}{m n_e} {\boldsymbol B} \boldsymbol{\cdot} [ F_{{\boldsymbol V}}\,{\times}\,(\boldsymbol{\nabla}\,{\times}\, G_{\boldsymbol B})- G_{{\boldsymbol V}}\,{\times}\,(\boldsymbol{\nabla}\,{\times}\, F_{\boldsymbol B}) ]-\frac{1}{en_e}{\boldsymbol B}\boldsymbol{\cdot} [(\boldsymbol{\nabla}\,{\times}\, F_{\boldsymbol B})\,{\times}\,(\boldsymbol{\nabla}\,{\times}\, G_{\boldsymbol B})]\nonumber\\ & \quad +\frac{1}{en_e}\sum_p\frac{e_p}{m_p}\int \, \textrm{d}^{3}v f_p {\boldsymbol B} \boldsymbol{\cdot} \left[\vphantom{\frac{e}{m}}\nabla_{\boldsymbol{v}} F_{f_p}\times(\boldsymbol{\nabla}\times G_{\boldsymbol B})- \nabla_{\boldsymbol{v}} G_{f_p}\times(\boldsymbol{\nabla}\times F_{\boldsymbol B})\right.\nonumber\\ & \quad +\left.\frac{e}{m}(\nabla_{\boldsymbol{v}} G_{f_p}\times F_{{\boldsymbol V}}-\nabla_{\boldsymbol{v}} F_{f_p} \times G_{{\boldsymbol V}})\right]\nonumber\\ & \quad -\frac{1}{e}\sum_p \frac{e_p}{m_p}\int \, \textrm{d}^{3}v f_p (\boldsymbol{\nabla} G_{n_e}\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} F_{f_p}-\boldsymbol{\nabla} F_{n_e}\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} G_{f_p})\nonumber\\ & \quad -\frac{1}{en_e}{\boldsymbol B} \boldsymbol{\cdot} \left( \sum_{p,p'}\frac{e_p}{m_p} \frac{e_{p'}}{m_{p'}}\int\int \, \textrm{d}^{3}v \, \textrm{d}^{3}v' f_p f_{p'} \nabla_{\boldsymbol{v}} F_{f_p}\times \nabla_{\boldsymbol{v}} G_{f_{p'}}\right) \nonumber\\ & \quad +\left.\sum_p \int \, \textrm{d}^{3}v \frac{f_p}{m_p}\left[[\![F_{f_p},G_{f_p}]\!]+ \frac{e_p}{m_p}{\boldsymbol B}\boldsymbol{\cdot} (\nabla_{\boldsymbol{v}} F_{f_p}\times \nabla_{\boldsymbol{v}} G_{f_p})\right] \right\}. \end{align}

Observe, unlike the bracket (2.30), (2.42) is gauge invariant. Now, having computed the bracket of our model, we write down the Hamiltonian functional, which is the direct sum of the fluid and particle Hamiltonians, i.e.

(2.43)\begin{equation} {\mathcal{H}}=\int \, \textrm{d}^{3}x \left[\rho \frac{|{\boldsymbol V}|^{2}}{2}+\rho U(\rho)+n_e {\mathcal{U}} _e+\frac{|{\boldsymbol B}|^{2}}{2\mu_0}\right]+\sum_p\int\int \, \textrm{d}^{3}x\,\textrm{d}^{3}v m_p \frac{v^{2}}{2} f_p. \end{equation}

Here, $U(\rho )$ is the specific internal energy of the thermal ion fluid and ${\mathcal {U}} _e$ is some electron internal energy function. The dependence of this function is dictated by the nature of the electron pressure. In our case, it was first shown by Morrison (1982) (and subsequent work, e.g. Hazeltine, Mahajan & Morrison 2013) that gyrotropic pressure tensors of the form (2.15) follow from an internal energy that depends explicitly on $n_e$ and $|{\boldsymbol B}|$. If an electron internal energy has the form

(2.44)\begin{equation} {\mathcal{U}} _e= {\mathcal{U}} _e(n_e,|{\boldsymbol B}|), \end{equation}

then the following equations give the pressure tensor components $P_{e\parallel }$ and $P_{e\perp }$:

(2.45)\begin{gather} \frac{\partial {\mathcal{U}} _e}{\partial n_e}=\frac{P_{e\parallel}}{n_e^{2}}, \end{gather}

(2.46)\begin{gather}\frac{\partial {\mathcal{U}} _e}{\partial |{\boldsymbol B}|}=\frac{P_{e\perp}-P_{e\parallel}}{n_e|{\boldsymbol B}|}; \end{gather}

and hence

(2.47)\begin{gather} \frac{\delta {\mathcal{H}}}{\delta n_e}= {\mathcal{U}} _e+\frac{P_{e\parallel}}{n_e}, \end{gather}

(2.48)\begin{gather}\frac{\delta {\mathcal{H}}}{\delta {\boldsymbol B}}=\mu_0^{{-}1}(1-\gamma){\boldsymbol B}, \end{gather}

where $\gamma := ({P_{e\parallel }-P_{e\perp }})/(B^2/ \mu_0)$ is a function that measures the electron pressure anisotropy. Note that the Poisson structure (2.42) is valid for barotropic (${\mathcal {U}} _e= {\mathcal {U}} _e(n_e$)) and gyrotropic [(2.45) and (2.46)] electron pressure. The identification of an appropriate Hamiltonian structure for a generally non-gyrotropic electron pressure tensor will be pursued in a future work.

The equations of motion follow from (2.29) with the Poisson bracket (2.42), the Hamiltonian (2.43) and also using (2.47), (2.48). It can be readily seen that $\partial _t n_e=\{n_e,{\mathcal {H}}\}$ and $\partial _t\rho =\{\rho,{\mathcal {H}}\}$ are indeed (2.22) and the continuity equation for the fluid ions, respectively. The latter is

(2.49)\begin{equation} \partial_t\rho={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (\rho{\boldsymbol V}). \end{equation}

The remaining equations, i.e. the momentum, the induction and the Vlasov equation that stem from (2.29), are

(2.50)\begin{align} \partial_t{\boldsymbol V} & ={\boldsymbol V}\times\boldsymbol{\nabla}\times{\boldsymbol V} -\boldsymbol{\nabla}\left( h+\frac{{V}^{2}}{2} \right)+\frac{e}{m}\left(1-\frac{n_i}{n_e} \right){\boldsymbol V}\times{\boldsymbol B}\nonumber\\ & \quad +\frac{1}{mn_e}(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B} -{\boldsymbol J}_k)\times{\boldsymbol B}-m^{{-}1}\boldsymbol{\nabla}\left( {\mathcal{U}} _e+\frac{P_{e\parallel}}{n_e}\right)-\frac{\boldsymbol{\nabla}\times(\gamma {\boldsymbol B})}{\mu_0mn_e}\times{\boldsymbol B} , \end{align}

(2.51)\begin{align} \partial_t {\boldsymbol B} & =\boldsymbol{\nabla}\times\left[\frac{n_i}{n_e}{\boldsymbol V}\times{\boldsymbol B}- \frac{1}{en_e}(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-{\boldsymbol J}_k)\times{\boldsymbol B}\right]\nonumber\\ & \quad +\mu_0^{{-}1}\boldsymbol{\nabla}\times\left[\frac{\gamma}{en_e} ({\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol B}-\boldsymbol{\nabla}{\boldsymbol B}\boldsymbol{\cdot} {\boldsymbol B})+\frac{\boldsymbol{\nabla}\gamma\times{\boldsymbol B}}{en_e} \times{\boldsymbol B}\right], \end{align}

(2.52)\begin{align} \partial_t f_p& ={-}\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f_p-\frac{e_p}{m_p}\left[ \left(\boldsymbol{v}-\frac{n_i}{n_e}{\boldsymbol V}\right)\times{\boldsymbol B}\right.\nonumber\\ & \quad \left.+\frac{1}{en_e}\left(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-{\boldsymbol J}_k\right)\times{\boldsymbol B} -\frac{1}{e}\boldsymbol{\nabla}\left( {\mathcal{U}} _e+\frac{P_{e\parallel}}{n_e}\right)- \frac{\boldsymbol{\nabla}\times(\gamma{\boldsymbol B})}{\mu_0en_e}\times{\boldsymbol B}\right]\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} f_p. \end{align}

It is a matter of vector calculus manipulations to prove that

(2.53)\begin{align} & -\frac{1}{m}\boldsymbol{\nabla}\left( {\mathcal{U}} _e+\frac{P_{e\parallel}}{n_e}\right)- \frac{\boldsymbol{\nabla}\times(\gamma{\boldsymbol B})}{\mu_0mn_e}\times{\boldsymbol B}\nonumber\\ & \quad ={-}\frac{1}{mn_e}[\boldsymbol{\nabla} P_{e\perp}+\mu_0^{{-}1}(\gamma {\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol B} +{\boldsymbol B}{\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla}\gamma) ]={-}\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{mn_e}. \end{align}

The last equality can be easily proven by taking the divergence of the gyrotropic pressure equation (2.15). In light of (2.53), the momentum and the Vlasov equations, (2.50) and (2.52), respectively, which arise from the Hamiltonian formulation, are identical to (2.20) and (2.21). In addition, performing some manipulations on the last term of (2.51), it can be seen that

(2.54)\begin{align} & \mu_0^{{-}1}\boldsymbol{\nabla}\times\left[\frac{\gamma}{en_e} ({\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol B}-\boldsymbol{\nabla}{\boldsymbol B}\boldsymbol{\cdot} {\boldsymbol B})+ \frac{\boldsymbol{\nabla}\gamma\times{\boldsymbol B}}{en_e}\times{\boldsymbol B}\right]\nonumber\\ & \quad =\boldsymbol{\nabla}\times\left[\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e}- \frac{1}{e}\boldsymbol{\nabla} {\mathcal{U}} _e-\frac{1}{e}\boldsymbol{\nabla}\left(n_e \frac{\partial {\mathcal{U}} _e}{\partial n_e}\right)\right]= \boldsymbol{\nabla}\times\frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e}{en_e}. \end{align}

Hence, (2.51) is indeed the induction equation (2.19).

2.3. Casimir invariants

A typical procedure to identify the Casimirs of a noncanonical Poisson bracket, e.g. (2.42), is to rearrange it in the following form:

(2.55)\begin{equation} \{F,G\}=\int \, \textrm{d}^{3}x F_{u_i} \mathcal{J}_{ij} G_{u_j}, \end{equation}

where ${\mathcal {J}}$ is the Poisson operator associated with this bracket, and then seek solutions to the following system of Casimir-determining equations:

(2.56)\begin{equation} \mathcal{J}_{ij} {\mathcal{C}} _{u_j}=0,\quad i=1,\ldots,5. \end{equation}

By this procedure, we find the following Casimir invariants:

(2.57a,b)\begin{gather} {\mathcal{C}} _1= \int \, \textrm{d}^{3}x \rho,\quad {\mathcal{C}} _2 = \int \, \textrm{d}^{3}x n_e, \end{gather}

(2.58)\begin{gather} {\mathcal{C}} _3= \frac{1}{2}\int \, \textrm{d}^{3}x {\boldsymbol A}\boldsymbol{\cdot} {\boldsymbol B}, \end{gather}

(2.59)\begin{gather}{\mathcal{C}} _4= \frac{1}{2}\int \, \textrm{d}^{3}x \left({\boldsymbol V}+\frac{e}{m}{\boldsymbol A}\right)\boldsymbol{\cdot} \left(\boldsymbol{\nabla}\times{\boldsymbol V}+\frac{e}{m}{\boldsymbol B}\right), \end{gather}

(2.60)\begin{gather}{\mathcal{C}} _p= \int\int \, \textrm{d}^{3}x\,\textrm{d}^{3}v \varLambda_p(f_p). \end{gather}

This set of Casimirs, typical of magnetofluid models (e.g. Lingam, Miloshevich & Morrison 2016), has two helicity invariants, but is augmented by an additional kinetic Casimir for each particle species $p$, which involves the corresponding distribution function, while there are no cross-fluid-kinetic Casimirs. Here, $\varLambda _p(f_p)$ are arbitrary functions of their respective distribution functions $f_p$.

3. Pressure coupling scheme

3.1. Transformation of the Hamiltonian CCS

As it is well known, another method to couple the energetic and the fluid components is through the pressure tensor appearing in a center-of-mass momentum equation rather than through the current density. This can be done upon taking the first-order fluid moment of the Vlasov equations, which govern the dynamics of the energetic species, and adding it to the momentum equation of the thermal ions. The first-order fluid moment of (2.3) yields

(3.1)\begin{equation} m_p\partial_t (n_p{\boldsymbol V}_p)=e_pn_p({\boldsymbol E}+{\boldsymbol V}_p\times{\boldsymbol B})-\sum_p\boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{P}}}_p, \end{equation}

where

(3.2a–c)\begin{equation} n_p=\int \, \textrm{d}^{3}v f_p,\quad {\boldsymbol V}_p=\frac{1}{n_p}\int \, \textrm{d}^{3}v \boldsymbol{v} f_p ,\quad \tilde{\boldsymbol{\mathsf{P}}}_p=m_p\int \, \textrm{d}^{3}v \boldsymbol{v}\boldsymbol{v} f_p. \end{equation}

Adding (3.1) and (2.2) for $s=i$ and using the Ohm's law (2.18), we obtain the following momentum equation:

(3.3)\begin{equation} \partial_t{\boldsymbol M}+m\boldsymbol{\nabla}\boldsymbol{\cdot} (n_i {\boldsymbol V}{\boldsymbol V})=\mu_0^{{-}1}(\boldsymbol{\nabla}\times{\boldsymbol B})\times{\boldsymbol B} -\boldsymbol{\nabla} P_i-\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e-\sum_p\boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{P}}}_p, \end{equation}

where

(3.4)\begin{equation} {\boldsymbol M}:=mn_i{\boldsymbol V}+\sum_p m_pn_p {\boldsymbol V}_p. \end{equation}

If we replace (2.20) of the CCS model by (3.3), we obtain a totally equivalent model coupling though the fluid and the particle components through the particle species pressure tensor terms in (3.3). It is logical then to expect that this reformulation of the model would have a Hamiltonian structure resulting from some change of dynamical variables. Tronci (2010) showed that the reformulation of the Hamiltonian structure is effected by considering ${\boldsymbol M}$ instead of ${\boldsymbol V}$ as an independent dynamical variable. Performing this change of variables, one can see how the functional derivatives with respect to the old and the new sets of dynamical variables should relate. Requiring

(3.5)\begin{equation} \delta F[n_i,n_e,{\boldsymbol V},{\boldsymbol B},f_p]=\delta \tilde{F}[n_i,n_e,{\boldsymbol M},{\boldsymbol B},f_p], \end{equation}

and employing the chain rule for functional derivatives, the following equations can be deduced

(3.6)\begin{equation} \left.\begin{array}{c@{}} F_{n_e}=\tilde{F}_{n_e}, \quad F_{{\boldsymbol B}}=\tilde{F}_{{\boldsymbol B}},\quad F_{{\boldsymbol V}}=m n_i \tilde{F}_{{\boldsymbol M}},\\ F_{n_i}=\tilde{F}_{n_i}+m {\boldsymbol V}\boldsymbol{\cdot} \tilde{F}_{{\boldsymbol M}},\quad F_{f_p}=\tilde{F}_{f_p}+m_p\boldsymbol{v}\boldsymbol{\cdot} \tilde{F}_{{\boldsymbol M}}. \end{array}\right\} \end{equation}

Substituting equations (3.6) into the Poisson bracket (2.42) and using the definition (3.4) of ${\boldsymbol M}$, we find, after some algebra, the following bracket:

(3.7)\begin{align} \{F,G\}& = \int \, \textrm{d}^{3}x\,\left\{\vphantom{\frac{1}{en_e}}{\boldsymbol M}\boldsymbol{\cdot} (G_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{{\boldsymbol M}}-F_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{{\boldsymbol M}}) +n_i (G_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_i}-F_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_i})\right.\nonumber\\ & \quad + n_e(G_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_e}-F_{{\boldsymbol M}}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_e}) +{\boldsymbol B}\boldsymbol{\cdot} [F_{{\boldsymbol M}}\times(\boldsymbol{\nabla}\times{\boldsymbol B})-G_{{\boldsymbol M}}\times(\boldsymbol{\nabla}\times F_{{\boldsymbol B}})] \nonumber\\ & \quad - \sum_p\frac{e_p}{em_p}\int \, \textrm{d}^{3}v f_p(\boldsymbol{\nabla}_{\boldsymbol{v}}F_{f_p}\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_e}-\boldsymbol{\nabla}_{\boldsymbol{v}}G_{f_p}\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_e})\nonumber\\ & \quad + \frac{1}{en_e}\sum_p \frac{e_p}{m_p}\int \, \textrm{d}^{3}v f_p{\boldsymbol B}\boldsymbol{\cdot} [\boldsymbol{\nabla}_{\boldsymbol{v}}F_{f_p}\times(\boldsymbol{\nabla}\times G_{{\boldsymbol B}})-\boldsymbol{\nabla}_{\boldsymbol{v}}G_{f_p} \times(\boldsymbol{\nabla}\times F_{{\boldsymbol B}})]\nonumber\\ & \quad - \frac{1}{en_e}{\boldsymbol B}\boldsymbol{\cdot} \int \int \, \textrm{d}^{3}v\,\textrm{d}^{3}v'\, \sum_{p,p'}\frac{e_pe_{p'}}{m_pm_{p'}}f_p(\boldsymbol{v})f_{p'}(\boldsymbol{v}')\nabla_{\boldsymbol{v}} F_{f_p} \times\boldsymbol{\nabla}_{\boldsymbol{v}'}G_{f_{p'}}\nonumber\\ & \quad + \sum_p m_p^{{-}1}\int \, \textrm{d}^{3}v \, f_p \left[[\![F_{f_p},G_{f_p}]\!]+\frac{e_p}{m_p}{\boldsymbol B}\boldsymbol{\cdot} (\nabla_{\boldsymbol{v}} F_{f_p}\times \nabla_{\boldsymbol{v}} G_{f_p})\right.\nonumber\\ & \quad \left.\left.+\vphantom{\frac{e_p}{m_p}} m_p([\![F_{f_p}, \boldsymbol{v}\boldsymbol{\cdot} G_{{\boldsymbol M}}]\!]-[\![G_{f_p},\boldsymbol{v}\boldsymbol{\cdot} F_{{\boldsymbol M}}]\!])\right]-\frac{1}{en_e}{\boldsymbol B}\boldsymbol{\cdot} [(\boldsymbol{\nabla}\times F_{{\boldsymbol B}})\times(\boldsymbol{\nabla}\times G_{{\boldsymbol B}})]\right\} . \end{align}

We should also write the Hamiltonian in terms of ${\boldsymbol M}$ and then invoke Hamilton's equations to retrieve the dynamical system in this PCS formulation. The transformed Hamiltonian reads

(3.8)\begin{align} {\mathcal{H}}& = \int \, \textrm{d}^{3}x \, \left[\frac{|{\boldsymbol M}-\boldsymbol{{\mathcal{M}}}|^{2}}{2m n_i}+ n_i {\mathcal{U}} _i(n_i)+n_e {\mathcal{U}} _e(n_e,|{\boldsymbol B}|)+\frac{|{\boldsymbol B}|^{2}}{2\mu_0} \right] \nonumber\\ & \quad + \sum_p \frac{m_p}{2}\int \int \, \textrm{d}^{3}x\,\textrm{d}^{3}v f_pv^{2}, \end{align}

where

(3.9)\begin{equation} \boldsymbol{{\mathcal{M}}}:=\sum_p \int \, \textrm{d}^{3}v m_pf_p \boldsymbol{v}. \end{equation}

The functional derivative of ${\mathcal {H}}$ with respect to the new variable ${\boldsymbol M}$ is

(3.10)\begin{equation} \frac{\delta {\mathcal{H}}}{\delta {\boldsymbol M}}={\boldsymbol V}. \end{equation}

Also, one should notice that the functional derivative with respect to $f_p$ is not the same as in the previous case of the CCS, but

(3.11)\begin{equation} \frac{\delta {\mathcal{H}}}{\delta f_p}=m_p\frac{v^{2}}{2}-m_p\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol V}. \end{equation}

We can verify that Hamilton's equations are indeed the dynamical equations of the CCS model with the momentum equation being replaced by (3.3).

3.2. Conventional construction of a PCS

The Hamiltonian system described in the previous subsection is equivalent to the CCS model of § 2 because it is obtained from CCS by a mere change of variables. However, it is a commonality in the hybrid fluid–kinetic models to employ pressure coupling schemes which involve a simpler coupling of the bulk plasma with the energetic particles assuming that $n_p\ll n_i$ and ${\boldsymbol M}/\rho \sim {\boldsymbol V}$ in the dynamical equations. These assumptions result in a momentum equation governing the ion momentum density, instead of the total momentum density ${\boldsymbol M}$, and containing the divergence of a pressure tensor associated with the energetic particles (e.g. see Park et al. 1992; Kim, Sovinec & Parker 2004; Fu et al. 2006; Takahashi, Brennan & Kim 2009). A similar treatment in our case with a two-fluid, Hall-MHD bulk, results in the following set of equations:

(3.12)\begin{gather} \partial_t n +\boldsymbol{\nabla}\boldsymbol{\cdot} (n{\boldsymbol V})=0, \end{gather}

(3.13)\begin{gather}\rho (\partial_t{\boldsymbol V}+{\boldsymbol V}\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol V})={\boldsymbol J}\times{\boldsymbol B} -\sum_p\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_p-\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e-\boldsymbol{\nabla} P_i, \end{gather}

(3.14)\begin{gather}\partial_t {\boldsymbol B}=\boldsymbol{\nabla}\times\left[ {\boldsymbol V}\times {\boldsymbol B} -\frac{1}{en}{\boldsymbol J}\times {\boldsymbol B}+\frac{1}{en}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e\right], \end{gather}

(3.15)\begin{gather}\partial_t f_p+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f_p+\frac{e_p}{m_p} \left[ (\boldsymbol{v}-{\boldsymbol V})\times {\boldsymbol B} +\frac{1}{en}{\boldsymbol J}\times{\boldsymbol B}- \frac{1}{en}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e\right]\boldsymbol{\cdot} \nabla_{\boldsymbol{v}} f_p=0. \end{gather}

Here, $n=n_i=n_e$, which is the zeroth order quasineutrality condition for $n_p/n_i \ll 1$.

3.3. Hamiltonian construction of a PCS

It has been stressed in several works (Tronci 2010; Tronci et al. 2014; Burby & Tronci 2017; Close et al. 2018) that such PC models, which are widely employed in the study of energetic particle effects on plasma stability, do not conserve some energy functional as all consistent ideal plasma models do (e.g. Morrison & Greene 1980; Marsden & Weinstein 1982; Lingam et al. 2015). In addition, it has been shown that a Vlasov–MHD model in the pressure-coupling scheme exhibits spurious instabilities attributed to the lack of energy conservation (Tronci et al. 2014). In contrast, a Hamiltonian variant of the model, derived by a procedure that ensures energy conservation (Tronci 2010), does not contain these unphysical modes. Prompted by this observation, we employ the method of Tronci (2010) in our case so as to derive a simplified, yet Hamiltonian, PCS model. The idea is, instead of assuming $n_p\ll n_i$ to simplify the equations of motion, to replace the ion momentum density ${\boldsymbol M}-\boldsymbol {{\mathcal {M}}}$ by the total momentum density ${\boldsymbol M}$ in the Hamiltonian functional, implicitly assuming that ${\boldsymbol M}\approx mn_i{\boldsymbol V}$ on the Hamiltonian level. Hence, the term $|{\boldsymbol M}-\boldsymbol {{\mathcal {M}}}|^{2}/(2m n_i)$ in (3.8) becomes $|{\boldsymbol M}|^{2}/(2mn_i)$, which leads to

(3.16)\begin{align} {\mathcal{H}}& = \int \, \textrm{d}^{3}x \left[\frac{|{\boldsymbol M}|^{2}}{2mn_i}+n_i {\mathcal{U}} _i(n_i)+n_e {\mathcal{U}} _{e}(n_e,|{\boldsymbol B}|)+\frac{|{\boldsymbol B}|^{2}}{2\mu_0}\right]\nonumber\\ & \quad +\sum_p \frac{m_p}{2}\int \int \, \textrm{d}^{3}x\, \textrm{d}^{3}v f_p v^{2}\nonumber\\ & =\int \, \textrm{d}^{3}x \left[\frac{mn_i}{2}|{\boldsymbol u}|^{2}+n_i {\mathcal{U}} _i(n_i) +n_e {\mathcal{U}} _{e}(n_e,|{\boldsymbol B}|)+\frac{|{\boldsymbol B}|^{2}}{2\mu_0}\right]\nonumber\\ & \quad +\sum_p \frac{m_p}{2}\int \int \, \textrm{d}^{3}x\, \textrm{d}^{3}v f_p v^{2}, \end{align}

where ${\boldsymbol u}:={\boldsymbol M}/(m n_i)$ is a weighted sum of the ion velocities. Note that for $n_p\ll n_i$, and if ${\boldsymbol V}_p$ are comparable to ${\boldsymbol V}$, then ${\boldsymbol u}={\boldsymbol V}$ up to zeroth order in $n_p/n_i$. It is convenient for us, in terms of comparing with the results of the previous section and other studies, to write the Poisson bracket (3.7) in terms of the velocity variable ${\boldsymbol u}$. This is done in Appendix A. The Casimirs of the bracket given in Appendix A are

(3.17a,b)\begin{align} {\mathcal{C}} _1 & = \int \, \textrm{d}^{3}x n_i, \quad {\mathcal{C}} _2 = \int \, \textrm{d}^{3}x n_e, \end{align}

(3.18)\begin{align} {\mathcal{C}} _3 & = \frac{1}{2}\int \, \textrm{d}^{3}x {\boldsymbol A}\boldsymbol{\cdot} {\boldsymbol B}, \end{align}

(3.19)\begin{align} {\mathcal{C}} _4 & = \int \, \textrm{d}^{3}x\,{\boldsymbol u}\boldsymbol{\cdot} \left(\frac{1}{2}\boldsymbol{\nabla}\times{\boldsymbol u}+ \frac{e}{m}{\boldsymbol B}\right)\nonumber\\ & \quad -\sum_p\int \, \textrm{d}^{3}x\frac{m_p}{m n_i}\int \, \textrm{d}^{3}v f_p\boldsymbol{v}\boldsymbol{\cdot} \left[\boldsymbol{\nabla}\times\left( {\boldsymbol u}-\frac{1}{2m n_i} \sum_{p'}m_{p'}\int \, \textrm{d}^{3}v' f_{p'}\boldsymbol{v}'\right)+\frac{e}{m}{\boldsymbol B}\right], \end{align}

(3.20)\begin{align} {\mathcal{C}} _p & = \int \, \textrm{d}^{3}x\,\textrm{d}^{3}v \varLambda_p(f_p). \end{align}

The equations of motion that arise from Hamilton's equations with the Poisson bracket (A2) (Appendix A) and the Hamiltonian (3.16) are

(3.21)\begin{align} \partial_t n_i& ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (n_i{\boldsymbol u}), \end{align}

(3.22)\begin{align} \partial_t n_e& ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (n_e {\boldsymbol u})- \frac{1}{e}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol J}_k, \end{align}

(3.23)\begin{align} \partial_t {\boldsymbol B}& =\boldsymbol{\nabla}\times\left[{\boldsymbol u}\times{\boldsymbol B} -\frac{1}{en_e} (\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-{\boldsymbol J}_k)\times{\boldsymbol B}+ \frac{1}{en_e}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e\right], \end{align}

(3.24)\begin{align} \partial_t {\boldsymbol u}& ={\boldsymbol u}\times\boldsymbol{\nabla}\times{\boldsymbol u}-\boldsymbol{\nabla}( h_i+u^{2}/2) - \rho^{{-}1}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e\nonumber\\ & \quad -\rho^{{-}1}\sum_p \boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\mathsf{P}}}_p+\mu_0^{{-}1} \rho^{{-}1}(\boldsymbol{\nabla}\times{\boldsymbol B})\times{\boldsymbol B}, \end{align}

(3.25)\begin{align} \partial_t f_p& ={-}(\boldsymbol{v}+{\boldsymbol u})\boldsymbol{\cdot} \boldsymbol{\nabla} f_p-\frac{e_p}{m_p} \left[\boldsymbol{v}\times{\boldsymbol B}+\frac{1}{en_e}(\mu_0^{{-}1}\boldsymbol{\nabla}\times{\boldsymbol B}-{\boldsymbol J}_k)\times{\boldsymbol B}\right.\nonumber\\ & \quad -\left.\frac{1}{en_e}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}_e\right] \boldsymbol{\cdot} \boldsymbol{\nabla}_{v}f_p+\nabla_{\boldsymbol{v}} f_p\boldsymbol{\cdot} \boldsymbol{\nabla}{\boldsymbol u}\boldsymbol{\cdot} \boldsymbol{v}, \end{align}

where $\rho =m n_i$ and $h_i=(1/m) \textrm{d}[n_i \mathcal{U}_i(n_i)]/\textrm{d}n_i$ is the specific enthalpy of the thermal ion fluid. Computing the first-order velocity moment of the Vlasov equation (3.25), we find

(3.26)\begin{equation} \partial_t\sigma_k={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (\sigma_k{\boldsymbol u})-\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol J}_k, \end{equation}

which, combined with (3.21) and (3.22), yields

(3.27)\begin{equation} \partial_t \sigma={-}\boldsymbol{\nabla}\boldsymbol{\cdot} (\sigma{\boldsymbol u}), \end{equation}

which is a continuity equation implying global but not local charge conservation. However, if the plasma is initially quasineutral, i.e. $\sigma =0$ at $t=0$, then quasineutrality is ensured by (3.27) $\forall t>0$.

4. Translationally symmetric formulation and energy-Casimir equilibria

4.1. Translationally symmetric formulation

In this section, we consider simplified dynamics with all dynamical variables being invariant along a fixed straight axis in physical space. In this translationally symmetric case, the magnetic field and macroscopic velocity can be expressed in terms of five scalar Clebsch potentials. Note that the distribution functions, although translationally symmetric in space, still depend on all three microscopic velocity coordinates. Translationally symmetric models facilitate computer simulations, because the dependency on one spatial coordinate is dropped and can be considered good approximations whenever a strong guiding magnetic field directed in a fixed direction is present. Using Cartesian coordinates $(x,y,z)$ in physical space and assuming invariance along the $z$ axis, the velocity and the magnetic field can be written as

(4.1)\begin{gather} {\boldsymbol u} = u_z(x,y,t)\hat{z} +\boldsymbol{\nabla}\chi(x,y,t)\times \hat{z}+\boldsymbol{\nabla}\varUpsilon(x,y,t), \end{gather}

(4.2)\begin{gather}{\boldsymbol B} = B_z(x,y,t)\hat{z}+\boldsymbol{\nabla}\psi(x,y,t)\times\hat{z}. \end{gather}

To obtain a Hamiltonian formulation in this symmetric case, we need to translate this field decomposition in a decomposition of the vector functional derivatives in terms of functional derivatives with respect to the scalar fields $u_z,B_z,\chi,\psi,\varUpsilon$. The process of deriving these relations has been described several times before, e.g. by Andreussi, Morrison & Pegoraro (2010), Kaltsas, Throumoulopoulos & Morrison (2017), Grasso et al. (2017) and Kaltsas, Throumoulopoulos & Morrison (2018). Following the same procedure here, we find

(4.3)\begin{gather} F_{{\boldsymbol u}}=F_{u_z}\hat{z}+\boldsymbol{\nabla} F_\varOmega \times\hat{z}-\boldsymbol{\nabla} F_{w}, \end{gather}

(4.4)\begin{gather}F_{{\boldsymbol B}}=F_{B_z}\hat{z}-\boldsymbol{\nabla}(\varDelta^{{-}1}F_\psi)\times\hat{z}, \end{gather}

where $\varDelta ^{-1}$ is the inverse Laplacian operator, $\varOmega :=-\Delta \chi \hat {z}$ and $w:=\Delta \varUpsilon$. The translationally symmetric Hall MHD bracket is known from Kaltsas et al. (2017), Grasso et al. (2017), and hence we have to compute the translationally symmetric cross-kinetic-Hall MHD terms and the term accounting for electron fluid thermodynamics. This is done in Appendix B upon substituting (4.3) and (4.4) in the bracket (A2) of Appendix A. The resulting translationally symmetric bracket is

(4.5)\begin{align} \{F,G\}_{TS}& = \int \, \textrm{d}^{2}x \left\{ \vphantom{\frac{e_p}{m_p}} F_\rho \Delta G_w-G_\rho\Delta F_w\right.\nonumber\\ & \quad +\frac{n_e}{\rho}([F_\varOmega,G_{n_e}]- [G_\varOmega,F_{n_e}]+\boldsymbol{\nabla} F_w\boldsymbol{\cdot} \boldsymbol{\nabla} G_{n_e}-\boldsymbol{\nabla} G_w\boldsymbol{\cdot} \boldsymbol{\nabla} F_{n_e})\nonumber\\ & \quad +\rho^{{-}1} \varOmega ([F_\varOmega,G_\varOmega]+[F_w,G_w]+\boldsymbol{\nabla} F_w \boldsymbol{\cdot} \boldsymbol{\nabla} G_\varOmega -\boldsymbol{\nabla} F_\varOmega\boldsymbol{\cdot} \boldsymbol{\nabla} G_w)\nonumber\\ & \quad +u_z([F_\varOmega,\rho^{{-}1}G_{u_z}]-[G_\varOmega,\rho^{{-}1}F_{u_z}]+ \boldsymbol{\nabla}(\rho^{{-}1}G_{u_z})\boldsymbol{\cdot} \boldsymbol{\nabla} F_w\nonumber\\ & \quad-\boldsymbol{\nabla}(\rho^{{-}1}F_{u_z})\boldsymbol{\cdot} \boldsymbol{\nabla} G_w+\rho^{{-}1}F_\varUpsilon G_{u_z}- \rho^{{-}1}G_\varUpsilon F_{u_z}) \nonumber\\ & \quad +\psi ([F_\varOmega,\rho^{{-}1}G_\psi]-[G_\varOmega, \rho^{{-}1} F_\psi]+[F_{B_z},\rho^{{-}1}G_{u_z}]-[G_{B_z},\rho^{{-}1}F_{u_z}]\nonumber\\ & \quad +\boldsymbol{\nabla} F_w \boldsymbol{\cdot} \boldsymbol{\nabla} (\rho^{{-}1} G_\psi)-\boldsymbol{\nabla} G_w \boldsymbol{\cdot} \boldsymbol{\nabla} (\rho^{{-}1} F_\psi)+\rho^{{-}1}F_\varUpsilon G_\psi -\rho^{{-}1} G_\varUpsilon F_\psi ) \nonumber\\ & \quad +\rho^{{-}1}B_z([F_\varOmega,G_{B_z}]-[G_\varOmega,F_{B_z}]+ \boldsymbol{\nabla} F_w\boldsymbol{\cdot} \boldsymbol{\nabla} G_{B_z}-\boldsymbol{\nabla} G_w\boldsymbol{\cdot} \boldsymbol{\nabla} F_{B_z} )\nonumber\\ & \quad +\frac{1}{e}\psi([G_{B_z},n_{e}^{{-}1}F_\psi]-[F_{B_z},n_e^{{-}1}G_\psi]) -\frac{1}{en_e}B_z[F_{B_z},G_{B_z}]\nonumber\\ & \quad -\frac{1}{e}\sum_p\frac{e_p}{m_p}\int \, \textrm{d}^{3}v [n_e^{{-}1}B_z\nabla_{v_\perp} f_p \boldsymbol{\cdot} (G_{f_p} \boldsymbol{\nabla} F_{B_z}-F_{f_p}\boldsymbol{\nabla} G_{B_z})\nonumber\\ & \quad +\psi([n_e^{{-}1}F_{f_p}\partial_{v_z}f_p,G_{B_z}]-[n_e^{{-}1}G_{f_p} \partial_{v_z}f_p,F_{B_z}]\nonumber\\ & \quad +\boldsymbol{\nabla}\boldsymbol{\cdot} (n_e^{{-}1}F_\psi G_{f_p}\nabla_{v_\perp} f_p)-\boldsymbol{\nabla}\boldsymbol{\cdot} (n_e^{{-}1}G_\psi F_{f_p}\nabla_{v_\perp} f_p))]\nonumber\\ & \quad +\frac{1}{e}\sum_p\frac{e_p}{m_p}\int \, \textrm{d}^{3}v (\nabla_{v_\perp} f_p)\boldsymbol{\cdot} (F_{f_p}\boldsymbol{\nabla} G_{n_e}-G_{f_p}\boldsymbol{\nabla} F_{n_e})\nonumber\\ & \quad -\frac{1}{en_e}\sum_{p,p'} \frac{e_p}{m_p}\frac{e_{p'}}{m_{p'}}\int \int \, \textrm{d}^{3}v \, \textrm{d}^{3}v' f_p f_{p'}[B_z\langle F_{f_p},G_{f_{p'}}\rangle\nonumber\\ & \quad -\boldsymbol{\nabla}\psi\boldsymbol{\cdot} (\partial_{v_z}F_{f_p}\boldsymbol{\nabla}_{v_{{\perp}}'}G_{f_{p'}}- \partial_{v_z'}G_{f_p'}\boldsymbol{\nabla}_{v_{{\perp}}}F_{f_{p}})]\nonumber\\ & \quad +\sum_p\int \, \textrm{d}^{3}v \frac{f_p}{m_p}[[\![F_{f_p},G_{f_p}]\!]_{{\perp}}+ \frac{e_p}{m_p}(B_z\langle F_{f_p},G_{f_p} \rangle \nonumber\\ & \quad +\boldsymbol{\nabla}\psi \boldsymbol{\cdot} [\partial_{v_z}G_{f_p}\nabla_{v_\perp} F_{f_p}-\partial_{v_z}F_{f_p} \nabla_{v_\perp} G_{f_p}])\nonumber\\ & \quad +m_p([\![F_{f_p},\rho^{{-}1}(v_zG_{u_z}+\hat{z}\boldsymbol{\cdot} \boldsymbol{v}_\perp{\times}\boldsymbol{\nabla} G_{\varOmega}-\boldsymbol{v}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} G_w )]\!]_{{\perp}}\nonumber\\ & \quad-\left.\vphantom{\frac{e_p}{m_p}} [\![G_{f_p},\rho^{{-}1}(v_zF_{u_z}+\hat{z}\boldsymbol{\cdot} \boldsymbol{v}_\perp{\times}\boldsymbol{\nabla} F_{\varOmega}-\boldsymbol{v}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} F_w )]\!]_{{\perp}})]\right\}. \end{align}

Note that we have introduced a new bracket notation, namely,

(4.6)\begin{gather} {[}a,b]:= (\boldsymbol{\nabla} a\times\boldsymbol{\nabla} b)\boldsymbol{\cdot} \hat{z}, \end{gather}

(4.7)\begin{gather}\langle a,b\rangle := (\nabla_{v_\perp} a\times \nabla_{v_\perp} b)\boldsymbol{\cdot} \hat{z}, \end{gather}

(4.8)\begin{gather}{[}\![a,b]\!]_{{\perp}}:= \boldsymbol{\nabla}_\perp a\boldsymbol{\cdot} \nabla_{v_\perp} b-\boldsymbol{\nabla}_\perp b\boldsymbol{\cdot} \nabla_{v_\perp} a. \end{gather}

The translationally symmetric Hamiltonian reads as follows:

(4.9)\begin{align} {\mathcal{H}} & = \int \, \textrm{d}^{3}x \left[\frac{1}{2}\rho (u_z^{2}+|\boldsymbol{\nabla}\chi|^{2}+2[\varUpsilon,\chi]+|\boldsymbol{\nabla} \varUpsilon|^{2})\right.\nonumber\\ & \quad +\left.\frac{B_z^{2}}{2\mu_0}+\frac{|\boldsymbol{\nabla}\psi|^{2}}{2\mu_0}+ \rho U_i(\rho)+n_e {\mathcal{U}} _e(n_e,B_z,|\boldsymbol{\nabla}\psi|)\right]\nonumber\\ & \quad +\sum_p\int \, \textrm{d}^{3}x\, \textrm{d}^{3}v \frac{1}{2}m_p f_p v^{2}, \end{align}

and the functional derivatives of ${\mathcal {H}}$ with respect to the two scalars associated with the magnetic field are

(4.10)\begin{gather} \frac{\delta{\mathcal{H}}}{\delta B_z} = \mu_0^{{-}1}(1-\gamma)B_z, \end{gather}

(4.11)\begin{gather}\frac{\delta {\mathcal{H}}}{\delta\psi} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} [\mu_0^{{-}1}(1-\gamma)\boldsymbol{\nabla}\psi]. \end{gather}

Having these expressions and also computing the functional derivatives with respect to the remaining scalars, one can derive the translationally symmetric dynamical equations, and in addition, equilibrium and stability conditions. For equilibrium and stability analysis, one has to identify the families of Casimir invariants ${\mathcal {C}}$ that span the non-trivial null space of the Poisson bracket (4.5). From the Casimir determining equations, stemming from the requirement $\mathfrak {C}_{\xi _i}=0$, where $\mathfrak {C}_{\xi _i}$ appear in the following rearrangement of (4.5):

(4.12)\begin{align} \{F,G\}=\int \, \textrm{d}^{3}x \left\{F_\rho \mathfrak{C}_\rho+F_{n_e} \mathfrak{C}_{n_e}+F_{B_z} \mathfrak{C}_{B_z}+F_\psi \mathfrak{C}_\psi+F_\chi \mathfrak{C}_\chi+F_\varUpsilon \mathfrak{C}_\varUpsilon+\sum_p\int \, \textrm{d}^{3}v\,F_{f_p} \mathfrak{C}_{f_p} \right\}, \end{align}

we were able to identify the following families of Casimirs

(4.13)\begin{gather} {\mathcal{C}} _1 = \int \, \textrm{d}^{3}x n_e N(\psi) , \end{gather}

(4.14)\begin{gather}{\mathcal{C}} _2 = \int \, \textrm{d}^{3}x \rho K(\varphi), \end{gather}

(4.15)\begin{gather}{\mathcal{C}} _3 = \int \, \textrm{d}^{3}x B_z F(\psi), \end{gather}

(4.16)\begin{gather}{\mathcal{C}} _4 = \int \, \textrm{d}^{3}x \left(\varOmega+\frac{e}{m}B_z-\sum_p m_p \int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}]\right)G(\varphi), \end{gather}

(4.17)\begin{gather}{\mathcal{C}} _{p} = \int \int \, \textrm{d}^{3}x\, \textrm{d}^{3}v\,\varLambda_p(f_p), \end{gather}

where $N,K,F,G$ and $\varLambda _p$ are arbitrary functions of their respective arguments, and

(4.18)\begin{equation} \varphi:=u_z+\frac{e}{m}\psi-\rho^{{-}1}\sum_pm_p\int \, \textrm{d}^{3}v f_p v_z, \end{equation}

is a generalized ion stream function, modified owing to the presence of kinetic particle species. For $f_p\rightarrow 0$, this stream function and the entire set of Casimir invariants (4.13)–(4.17) reduce to their translationally symmetric Hall MHD counterparts (Kaltsas et al. 2017). In this fluid limit, quasineutrality implies $n_e=n_i$, because $n_p=\int \, \textrm {d}^{3}v\,f_p=0$.

4.2. Energy-Casimir equilibria

Owing to spatial symmetry, the Casimirs (4.13)–(4.17) constitute infinite families of invariants, thereby allowing for the derivation of sufficient stability conditions by the energy-Casimir method, which has been used for hybrid kinetic-MHD models by Tronci et al. (2015) and for the extended and Hall MHD models by Kaltsas, Throumoulopoulos & Morrison (2020). A preliminary step though is the definition of the stationary state that serves as the initial condition for the dynamics. The set of equilibrium equations are derived in the first step of the energy-Casimir method by setting the first-order variation of the extended Hamiltonian equal to zero. In the MHD and extended MHD models, this procedure resulted in a system of Grad–Shafranov–Bernoulli (GSB) equations. Here, we derive the corresponding GSB system for translationally symmetric hybrid kinetic-Hall MHD in the PCS. The energy-Casimir functional, or extended Hamiltonian, is given by

(4.19)\begin{equation} \mathfrak{F}={\mathcal{H}}-\sum_{i=1}^{4} {\mathcal{C}} _i-\sum_p {\mathcal{C}} _p,\end{equation}

where ${\mathcal {H}}$ is given by (4.9) and ${\mathcal {C}}$ are the Casimirs (4.13)–(4.17). The first variation of (4.19) can be written as

(4.20)\begin{align} \delta\mathfrak{F} & = \int \, \textrm{d}^{3}x \left(\vphantom{\sum_p}\mathfrak{R}_1\delta n_e+\mathfrak{R}_2 \delta\rho+\mathfrak{R}_3 \delta u_z+\mathfrak{R}_4 \delta \chi+\mathfrak{R}_5 \delta \varUpsilon\right.\nonumber\\ & \quad+\left.\mathfrak{R}_6 \delta B_z+\mathfrak{R}_7 \delta\psi +\sum_p\int \, \textrm{d}^{3}v \mathfrak{R}_p\delta f_p \right). \end{align}

Assuming independent, arbitrary variations of the scalar dynamical variables, the requirement $\delta \mathfrak {F}=0$ is equivalent to $\mathfrak {R}_i=0, i=1,\ldots,7\,$ and $\mathfrak {R}_p=0, \forall p$. These equations lead to the following equilibrium conditions:

(4.21)\begin{align} & {\mathcal{U}} _e+n_{e}^{{-}1}P_{e_{{\parallel}}}-N(\psi)=0, \end{align}

(4.22)\begin{align} & h_i(\rho)+\frac{u^{2}}{2} -K(\varphi)-\rho^{{-}1}K'(\varphi)\sum_pm_p \int \, \textrm{d}^{3}v f_p v_z \nonumber\\ & \quad -\rho^{{-}2}\left(\varOmega+\frac{e}{m}B_z -\sum_p m_p \int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}] \right)G'(\varphi) \sum_{p'}m_{p'} \int \, \textrm{d}^{3}v f_{p'}v_z \nonumber\\ & \quad -\rho^{{-}2}\sum_p m_p \int \, \textrm{d}^{3}v f_{p}[\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x},G]=0, \end{align}

(4.23)\begin{align} & -\boldsymbol{\nabla}\boldsymbol{\cdot} (\rho \boldsymbol{\nabla}\varUpsilon )+[\chi,\rho]=0, \end{align}

(4.24)\begin{align} & -\boldsymbol{\nabla}\boldsymbol{\cdot} (\rho \boldsymbol{\nabla}\chi)+[\rho,\varUpsilon]+\Delta G(\varphi)=0, \end{align}

(4.25)\begin{align} & \rho u_z-\rho K'(\varphi)-\left(\varOmega+\frac{e}{m}B_z-\sum_p m_p \int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}]\right)G'(\varphi)=0, \end{align}

(4.26)\begin{align} & \mu_0^{{-}1}(1-\gamma)B_z-F(\psi)-\frac{e}{m}G(\varphi)=0, \end{align}

(4.27)\begin{align} & \mu_0^{{-}1}\boldsymbol{\nabla}\boldsymbol{\cdot} [(1-\gamma)\boldsymbol{\nabla}\psi]+n_e N'(\psi)+B_zF'(\psi) +\frac{e}{m}\rho K'(\varphi)\nonumber\\ & \quad +\frac{e}{m} \left(\varOmega+\frac{e}{m}B_z-\sum_p m_p \int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}]\right)G'(\varphi)=0, \end{align}

(4.28)\begin{align} & \frac{1}{2}m_pv^{2}-\varLambda_p'(f_p)+m_pv_z\left[K'(\varphi)\right.\nonumber\\ & \quad +\rho^{{-}1}\left.\left(\varOmega+\frac{e}{m}B_z-\sum_p m_p\int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}]\right)G'(\varphi)\right]\nonumber\\ & \quad +\rho^{{-}1}m_p[\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x},G(\varphi)]=0. \end{align}

Equation (4.21) provides a relation between $P_{e_\parallel }$ and the variables $n_e, \psi$ and $B_z$. Upon combining (4.21) with (2.46), we find an equation for $P_{e_\perp }$ also. This one reads as follows:

(4.29)\begin{equation} P_{e_\perp}=n_e\left[|{\boldsymbol B}|\frac{\partial {\mathcal{U}} _e}{\partial |{\boldsymbol B}|} +N(\psi)- {\mathcal{U}} _(n_e,|{\boldsymbol B}|)\right]. \end{equation}

In terms of $B_z$ and $|\boldsymbol {\nabla }\psi |$, (4.29) becomes

(4.30)\begin{equation} P_{e_\perp}=n_e\left[|{\boldsymbol B}|^{2}\left(\frac{1}{B_z}\frac{\partial {\mathcal{U}} _e}{\partial B_Z}+ \frac{1}{|\boldsymbol{\nabla}\psi|}\frac{\partial {\mathcal{U}} _e}{\partial |\boldsymbol{\nabla}\psi|}\right) +N(\psi)- {\mathcal{U}} _e(n_e,B_z,|\boldsymbol{\nabla}\psi|)\right]. \end{equation}

Therefore, to fully define an equilibrium state, we need an equation of state for the electron fluid, i.e. ${\mathcal {U}} _e= {\mathcal {U}} _e(n_e,B_z,|\boldsymbol {\nabla }\psi |)$. For example, such equations, which are associated with collisionless magnetic reconnection, are provided by Le et al. (2009).

From (4.26), we readily obtain an equation that relates $B_z$ with $\psi, \varphi$ and $\gamma$,

(4.31)\begin{equation} B_z=\mu_0\frac{F(\psi)+(e/m)G(\varphi)}{1-\gamma}. \end{equation}

Equation (4.31) is in general an implicit equation for determining $B_z$ because the anisotropy function, $\gamma$, depends on $B_z$ as well, and hence $B_z$ appears nonlinearly in (4.26). We could possibly make $\gamma$ independent of $B_z$ under some special assumption, i.e. selecting the functional dependence of ${\mathcal {U}} _{e}$ on $B_z$ in a particular manner.

Now, by multiplying (4.28) with $f_p$, integrating over the velocity space and summing over the particle species, we obtain the following equation:

(4.32)\begin{align} & \sum_p\int \, \textrm{d}^{3}v \frac{m_p}{2}f_p v^{2}+\sum_p K'(\varphi) \int \, \textrm{d}^{3}v m_p f_p v_z-\sum_p\int \, \textrm{d}^{3}v f_p \varLambda_p'(f_p) \nonumber\\ & \quad +\rho^{{-}1}\left(\varOmega +\frac{e}{m} B_z -\sum_p m_p \int \, \textrm{d}^{3}v [\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}] \right)G'(\varphi)\sum_p \int \, \textrm{d}^{3}v m_p f_p v_z\nonumber\\ & \quad +\sum_pm_p \rho^{{-}1} \int \, \textrm{d}^{3}v f_p [\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x},G(\varphi)]=0. \end{align}

Combining equations (4.22) and (4.32), we find a Bernoulli equation of the form

(4.33)\begin{equation} \rho h(\rho)=\rho K(\varphi)-\frac{1}{2}\rho u^{2}-\sum_p \left[\int \, \textrm{d}^{3}v \, f_p \varLambda_p'(f_p)-\int \, \textrm{d}^{3}v \frac{m_p}{2}f_p v^{2}\right]. \end{equation}

The two Grad–Shafranov equations, one for the thermal ions and one for the electron fluid, stem from (4.25) and (4.27), respectively, using (4.23), (4.24) and (4.26), along with the definition of $\varphi$. These GS equations are

(4.34)\begin{align} & G'(\varphi)\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\frac{G'}{\rho}\boldsymbol{\nabla}\varphi\right)+ \rho\left(\varphi-\frac{e}{m}\psi\right)-\rho K'(\varphi)-\mu_0 \frac{e}{m}\frac{F(\psi)+\frac{e}{m}G(\varphi)}{1-\gamma}G'(\varphi)\nonumber\\ & \quad+\sum_p m_p \int \, \textrm{d}^{3}v (f_p v_z+G'[\rho^{{-}1}f_p,\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol x}])=0, \end{align}

(4.35)\begin{align} & \mu_0^{{-}1}\boldsymbol{\nabla}\boldsymbol{\cdot} [(1-\gamma)\boldsymbol{\nabla}\psi]+n_e N'(\psi)+\mu_0 \frac{F(\psi)+\frac{e}{m}G(\varphi)}{1-\gamma}F'(\psi)\nonumber\\ & \quad +\frac{e}{m}\rho\left(\varphi-\frac{e}{m}\psi+\rho^{{-}1} \sum_p m_p \int \, \textrm{d}^{3}v f_p v_z\right)=0. \end{align}

Note that in (4.35), the anisotropy parameter $\gamma$ appears inside the differential operator, and hence this partial differential equation might not be always elliptic (Ito et al. 2007) as is the case for an isotropic electron pressure. Finally, we need a set of equations for determining the equilibrium distribution functions $f_{p,e}$. Using (4.25), (4.28) is simplified to

(4.36)\begin{equation} \varLambda_p'(f_p)=\tfrac{1}{2}m_p v^{2}+m_p \boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol u}, \end{equation}

where we have used

(4.37)\begin{equation} {\boldsymbol u}_\perp = \rho^{{-}1}\boldsymbol{\nabla} G\times\hat{z}, \end{equation}

which can be deduced from (4.24) and from the fact that $[\boldsymbol {v}\boldsymbol {\cdot } {\boldsymbol x},G]=(\boldsymbol {\nabla } G\times \hat {z})\boldsymbol {\cdot } \boldsymbol {v}$. For invertible functions $\varLambda _p$, we take solutions to (4.36) of the form

(4.38)\begin{equation} f_{p,e}=f_{p,e}(\tfrac{1}{2} (|\boldsymbol{v}+{\boldsymbol u}|^{2}-|{\boldsymbol u}|^{2})), \end{equation}

as was the case in the paper by Morrison et al. (2014) for the planar kinetic-MHD PCS model. Here, the subscript $e$ denotes equilibrium distribution functions. With (4.36), the Bernoulli equation becomes

(4.39)\begin{equation} \rho h(\rho)=\rho K(\varphi)-\frac{1}{2}\rho u^{2}-\sum_p m_p\int \, \textrm{d}^{3}v f_p\boldsymbol{v}\boldsymbol{\cdot} {\boldsymbol u}. \end{equation}

Assuming a special functional form for $f_{p,e}$, we can, in principle, compute the velocity space integrals appearing in (4.34), (4.35) and (4.39) to obtain a Hall-MHD GSB system, modified by the presence of energetic particles. As a final remark, note that the Hall-MHD GSB equations with anisotropic electron pressure are retrieved from (4.34), (4.35), (4.39), in the limit $f_p\rightarrow 0$, as can be seen by comparison with variants of this system obtained by Kaltsas et al. (2017,  2018) for isotropic electron pressure.

5. Conclusions

We constructed two kinetic-Hall MHD models with fluid and kinetic ions and fluid electrons, neglecting electron length scales. The coupling of the kinetic and the fluid components is effected through the current density in the first case (current coupling scheme) and through the pressure tensors of the particle species in the second (pressure coupling scheme). Moreover, we consider a gyrotropic electron pressure tensor, which is legitimate in the Hall MHD limit, i.e. when neglecting electron length scales. This description, allowing for both fluid and kinetic ions, bridges ordinary Hall MHD, which can be recovered in the limit of vanishing kinetic-ion population, and the most common hybrid kinetic-ion/fluid-electron model often used in hybrid simulations. The basic structure of the Hall MHD dynamics is retained while enriched by kinetic effects which are consistently described by the Vlasov equation. The Hamiltonian structures of the models are identified using a method introduced by Tronci (2010) and, in addition, a translationally symmetric description of the PCS model is obtained along with the corresponding Hamiltonian structure and the associated Casimir invariants. These functionals are deployed in an energy-Casimir variational principle that leads to a generalized Grad–Shafranov–Bernoulli system of equilibrium equations. We note that alternative hybrid models with non-gyrotropic electron pressure could have been obtained if we had considered alternative extended MHD descriptions, e.g. that presented by Tronci (2013), where the electron mean flow inertia is neglected in the Lagrangian of the parent kinetic theory. The use of extended MHD models with finite electron inertia and non-gyrotropic electron pressure tensors for describing the bulk plasma, and also, further investigations regarding the construction of specific equilibria and the derivation of sufficient stability criteria, will be the subject of a future work.