1. Introduction
Modelling the dynamics of collisionless (or weakly collisional) plasmas at scales comparable to or smaller than the ion Larmor radius is an important issue both for laboratory and astrophysical plasmas. At the level of a kinetic description, a valuable tool is given by the gyrokinetic theory which provides a reduction of the Vlasov–Maxwell (VM) equations by focusing on phenomena with a characteristic time scale large compared with the ion gyro-period. This approach, which eliminates the dependency on the gyration angle, typically adopts, as dynamical variables, the distribution functions of the gyrocentres rather than those of the particles. Within the gyrokinetic framework, a subset of models consists of the so-called $\delta f$-gyrokinetic models, which assume the gyrocentre distribution functions of the various particle species, to be close to those of an equilibrium state. In spite of this reduction, numerical simulations of three-dimensional gyrokinetic equations in a turbulent regime (even in the
$\delta f$-gyrokinetic framework) require huge computational resources, which justifies the development of simpler (although less complete) descriptions based on gyrofluid equations governing the evolution of a finite number of moments of the gyrocentre distribution functions. The relation between gyrocentre and particle moments is well defined and can usually be computed perturbatively. As in the case of the fluid hierarchy derived from the VM equations, a gyrofluid hierarchy of equations needs to be closed, in order to obtain a gyrofluid model with a finite number of dynamical variables. An important condition to prescribe at the level of closure assumptions is the preservation, in the absence of dissipation, of the Hamiltonian character of the parent gyrokinetic equations. This guarantees that in the reduction from a gyrokinetic to a gyrofluid system, not only is no uncontrolled dissipation of the total energy introduced but also that further invariants (Casimir invariants) of the system exist, and that the dynamics takes place on hyper-surfaces in phase space where the values of these invariants are constant. Another constraint is the consistency with the linear gyrokinetic theory. In particular, this requires retaining the influence of resonant effects such as Landau damping. A closure accounting for Landau damping in the
$\delta f$-gyrokinetic approach typically introduces dissipation and thus prevents the model from being Hamiltonian. Such form of dissipation is, however, voluntarily added and the main requirement is that the model possesses a Hamiltonian structure when the dissipative terms are removed.
In this spirit, the main goal of the present paper is to construct a gyrofluid model possessing the above mentioned properties, and primarily addressed to study phenomena relevant for collisionless space plasmas. Motivated by measurements of sub-proton fluctuations in the solar wind (see Sahraoui et al. (Reference Sahraoui, Goldstein, Belmont, Canu and Rezeau2010) and Alexandrova et al. (Reference Alexandrova, Chen, Sorriso-Valvo, Horbury and Bale2013) or Bruno & Carbone (Reference Bruno and Carbone2013) for reviews), reduced fluid models have already been derived and numerically integrated to explore the dynamics of space plasmas. Kinetic Alfvén wave (KAW) turbulence was for example addressed in Boldyrev & Perez (Reference Boldyrev and Perez2012). A more general Hamiltonian reduced gyrofluid model (Passot, Sulem & Tassi Reference Passot, Sulem and Tassi2018) which, in the appropriate asymptotic limit yields the model of Boldyrev & Perez (Reference Boldyrev and Perez2012), was recently developed to simultaneously capture the three regimes of: dispersive Alfvén waves at scales larger than the sonic and/or ion Larmor radius, of KAWs at sub-ion scales and also of inertial kinetic Alfvén waves at scales comparable to the electron inertial length (Chen & Boldyrev Reference Chen and Boldyrev2017; Passot, Sulem & Tassi Reference Passot, Sulem and Tassi2017). It has been used to derive weak turbulence kinetic equations (Passot & Sulem Reference Passot and Sulem2019) and the properties of imbalanced KAW turbulence in the framework of a reduction to nonlinear diffusion equations in spectral space (Miloshevich, Passot & Sulem Reference Miloshevich, Passot and Sulem2020). Nevertheless, to the best of our knowledge, existing reduced gyrofluid models can capture parallel magnetic field fluctuations, electron inertia and ion finite Larmor radius effects, but (with the exception of the recent models of Tassi (Reference Tassi2019), which will be discussed in § 3.2.2) do not take into account a possible temperature anisotropy of the equilibrium state. It, however, turns out that collisionless space plasmas, such as the solar wind, often exhibit anisotropic distribution functions (Marsch Reference Marsch2012) that can result from various heating effects such as Landau damping (Chen, Klein & Howes Reference Chen, Klein and Howes2019) or stochastic heating (Bourouaine & Chandran Reference Bourouaine and Chandran2013; Hoppock et al. Reference Hoppock, Chandran, Klein, Mallet and Verscharen2018) or from mechanical effects such as the action of a shear flow (De Camillis et al. Reference De Camillis, Cerri, Califano and Pegoraro2016). Proton temperature anisotropies play an important role in the solar wind (Hellinger et al. Reference Hellinger, Trávníĉek, Kasper and Lazarus2006), and possibly even more so at a closer distance from the Sun (Huang et al. Reference Huang, Kasper, Vech, Klein, Stevens, Martinović, Alterman, Ďurovcová, Paulson and Maruca2020). Temperature anisotropies are usually constrained by the micro-instabilities they trigger (e.g. mirror and firehose instabilities), both for ions (Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009) and electrons (Štverák et al. Reference Štverák, Trávníček, Maksimovic, Marsch, Fazakerley and Scime2008). The range of these accessible temperature anisotropies increases as the beta parameter (corresponding to the ratio between the equilibrium thermal plasma pressure and the magnetic pressure based on the guide field) decreases. In addition, temperature anisotropies are known to affect the development of the tearing instability (Shi, Lee & Fu Reference Shi, Lee and Fu1987), and thus the stability of current sheets (Matteini et al. Reference Matteini, Landi, Velli and Matthaeus2013), which plays a major role in the turbulence evolution (Franci et al. Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017). Another important feature to account for is the coupling to ion acoustic waves, which permits the development of the parametric decay of Alfvén waves at the magnetohydrodynamic (MHD) scales at small beta (Del Zanna, Velli & Londrillo Reference Del Zanna, Velli and Londrillo2001). This effect contributes to the generation of counter-propagating waves required for the development of a turbulent dynamics at these scales. Interestingly, this parametric instability, proposed by Bowen et al. (Reference Bowen, Badman, Hellinger and Bale2018) for the generation of the solar wind compressive fluctuations, turns out to occur in a wider range of beta parameters in the presence of temperature anisotropy (Tenerani, Velli & Hellinger Reference Tenerani, Velli and Hellinger2017). A useful feature of a new gyrofluid model for space plasma studies would thus be its capability to account for equilibrium temperature anisotropies and the coupling to ion acoustic waves.
As starting point for the derivation of the model, we choose the $\delta f$-gyrokinetic equations presented in Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015) where, for the sake of simplicity, we assume an electron–proton plasma with an equilibrium state described by bi-Maxwellian distribution functions with no mean drift velocity. Such a system fully satisfies our requirements. Indeed, it is a
$\delta f$-gyrokinetic model mainly conceived for pressure-anisotropic astrophysical plasmas and, as such, it specifically accounts for equilibrium temperature anisotropy (unlike most of gyrokinetic models which consider a generic equilibrium distribution function or specialize in the case of a Maxwellian equilibrium). Also, it accounts for parallel magnetic perturbations and it has been shown to possess (at least in the limit of interest for our derivation) a Hamiltonian structure (Tassi Reference Tassi2019). As far as the number of moments to be retained in the gyrofluid model is concerned, the inclusion of Landau damping requires us to retain at least the first three moments for each particle species. This is why, under the assumption of a two-species plasma, we opt for the derivation of a six-field gyrofluid model evolving three moments, including parallel temperature fluctuations, for each species. Nevertheless, a four-field Hamiltonian reduced version will also be presented and applied. Another novelty with respect to already existing Hamiltonian gyrofluid models, is the adoption of a closure relation, referred to as quasi-static, derived from linear gyrokinetic theory in the limit of slowly evolving fields. More precisely, according to such a closure, all gyrofluid moments that are not determined by gyrofluid evolution equations, are fixed according to their expression obtained from the gyrokinetic linear theory in the limit
$|\omega /(k_z v_{{\textrm {th}}_{\parallel s}})| \ll 1$, where
$\omega$ is the frequency of a mode,
$k_z$ the component of its wave vector along the direction of the guide field and
$v_{{\textrm {th}}_{\parallel s}}$ is the thermal speed of the species
$s$, associated with the equilibrium temperature along the direction of the guide field. As such, this closure is suitable for fields that are slowly evolving (i.e. quasi-static) with respect to particles travelling at the parallel equilibrium thermal speed
$v_{{\textrm {th}}_{\parallel s}}$. The derivation of this closure relation will be presented in appendix A (see in particular (A 21) to find the expressions for the various gyrofluid moments according to the quasi-static closure). We anticipate, however, two remarkable properties that this closure possesses. The first one is that the quasi-static closure relation turns out to be compatible with a Hamiltonian structure. The second one is that it allows for exact expressions, in terms of canonical Poisson brackets, for all the nonlinear terms in the gyrofluid equations, and in particular for those involving only gyroaveraged electromagnetic fields or potentials. This is not the case, to the best of our knowledge, for the previously derived reduced gyrofluid models.
The model assumes the presence of a strong magnetic guide field and evolves, for both electrons and ions, gyrocentre density fluctuations as well as velocity and temperature fluctuations referred to the direction parallel to the guide field. A dissipative variant of the model accounting for parallel Landau damping is then formulated through a Landau-fluid modelling of the parallel heat fluxes (Hammett & Perkins Reference Hammett and Perkins1990; Hammett, Dorland & Perkins Reference Hammett, Dorland and Perkins1992). A four-field reduction assuming isothermality is presented as well. The resulting gyrofluid model retains ion and electron finite Larmor radius (FLR) corrections, electron inertia and parallel magnetic fluctuations, enables anisotropic equilibrium temperatures and does not prescribe special restrictions on the ion or electron $\beta _{e,i}$ parameters, where
$\beta _{e,i}$ indicates, for each species (
$e$ for electrons and
$i$ for ions), the ratio between the equilibrium thermal pressure and the magnetic pressure exerted by the guide field. In this respect, this gyrofluid model differs from most of the presently available gyrofluid models, which require
$\beta _{e,i} \ll 1$.
In addition to the derivation of the model, we also carry out a detailed analysis of its linearized version. Predictions of the six-field model extended with Landau damping and of its four-field Hamiltonian reduction are compared at the linear level with those of the parent gyrokinetic model, by considering the dispersion relations of KAWs and slow waves (SWs), and also by analysing the firehose or the field-swelling (Basu & Coppi Reference Basu and Coppi1982, Reference Basu and Coppi1984) instabilities. This latter instability, which requires the electron temperature transverse to the magnetic field larger than the longitudinal one, leads to a local increase in the transverse pressure, which tends to make the magnetic field ‘swell’ further locally, an effect which can be important in producing magnetic reconnection. Its properties are rather subtle, especially when considering its effect on fast modes when the disturbance propagation is nearly perpendicular to the ambient magnetic field. We thus include an appendix summarizing the results on this instability that are relevant for the discussion of the present model. As will be remarked in § 4, for the non-dissipative case, the choice of the four-field model instead of the six-field model, for the comparison with the linear gyrokinetic theory, is due to the fact that a closure at an even order as in the four-field model, provides a better agreement with the Padé approximant chosen for the electron response function in the quasi-static limit.
The investigation of the linear dispersion relations including temperature anisotropies is a first application of the model, mainly devoted to test its capability of reproducing results of the linear gyrokinetic theory in the appropriate regimes. We mention that the derivation of the model was motivated also by further physical applications of relevance to space plasmas and which will be part of subsequent works. These include the investigation of tearing instability in the presence of temperature anisotropy in a strong guide field regime, the influence of electron FLR effects and temperature anisotropy on inertial reconnection, or the effect of the coupling of Alfvén and compressible modes on the turbulence development.
The paper is organized as follows. In § 2 the gyrokinetic parent model is reviewed. In § 3 the six-field gyrofluid model is introduced and its Hamiltonian structure is presented. The two variants, corresponding to the Landau gyrofluid extension and to the Hamiltonian four-field reduction, are also described. Section 4 is devoted to the comparison of the linearized versions of the two variants of the model with other linear theories. We conclude in § 5 where we also mention the interest of the present gyrofluid model for space plasma applications. At the end of the paper three appendices are provided, presenting the derivation of the quasi-static closure relations from gyrokinetic theory, the derivation of the six-field gyrofluid model equations and a discussion of field-swelling instabilities.
2. The gyrokinetic parent model
In order to derive a gyrofluid model only based on a quasi-static closure assumption, we consider as a starting point the following set of gyrokinetics equations, which corresponds to the system provided by equations (C 58), (C 60) and (C 66)–(C 68) of Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015) when collisions and equilibrium velocities are neglected and bi-Maxwellian distribution functions are chosen as equilibrium distribution functions for all the particle species. For simplicity, we specialize to the case of a plasma consisting of two species: electrons and one species of singly ionized particles. The equations of the resulting gyrokinetic model are given by




The index $s \in \{e,i \}$ adopted above indicates the particle species, so that quantities labelled with
$s$ refer to the electron species when
$s=e$ and to the ion species when
$s=i$.
In the system (2.1)–(2.4), the function $\tilde {g}_s$ is defined by

where $\tilde {f}_{s}$ is the perturbation of the gyrocentre distribution function for particles of species
$s$,
$q_{s}$ is the charge of these particles (so that
$q_e=-e$ and
$q_i=e$, with
$e$ indicating the proton charge) and
$m_s$ their mass. Furthermore,
$c$ denotes the speed of light. The bi-Maxwellian equilibrium distribution function is given by

where $n_0$ is the uniform and constant equilibrium density,
$T_{{0 }_{\parallel s}}$ and
$T_{{0 }_{\perp s}}$ are respectively the equilibrium temperatures of the particle of species s, parallel and in a plane perpendicular to an equilibrium magnetic guide field of amplitude
$B_0$, directed along the
$z$ direction of a Cartesian coordinate frame
$\{x,y,z\}$. We suppose that the spatial domain of the system corresponds to a box
$\mathcal {D}=\{ (x,y,z): -L_x \leq x \leq L_x , -L_y \leq y \leq L_y, -L_z \leq z \leq L_z \}$, with
$L_x$,
$L_y$ and
$L_z$ positive constants. All quantities of the system which depend on the spatial variables
$x$,
$y$ and
$z$ are supposed to satisfy periodic boundary conditions on the domain
$\mathcal {D}$, so that they can be expanded in Fourier series. We indicated with
$v_{\|} \in \mathbb {R}$ the velocity coordinate parallel to the guide field and with
$\mu _s=m_s v_{\perp }^2/(2 B_0) \in [0, +\infty )$ the magnetic moment of the particle of species
$s$ in the unperturbed guide field, where
$v_{\perp }$ corresponds to the velocity coordinate perpendicular to the guide field. We assume that all functions depending on
$v_{\|}$ decay to zero as
$v_{\|} \rightarrow \pm \infty$. Functions depending on
$\mu _s$ are assumed to tend to zero as
$\mu _s \rightarrow +\infty$ and to be bounded at
$\mu _s=0$. The coordinate
$t\in [0, + \infty )$ refers to time. The expressions
$\mathcal {J}_{0s}$ and
$\mathcal {J}_{1s}$ are related to the standard gyroaverage operators for the species
$s$ in Fourier space. Their definition can be introduced explicitly in the following way: adopting the notation
$\boldsymbol {x}$ to indicate a point of coordinates
$(x,y,z) \in \mathcal {D}$ and, similarly,
$\boldsymbol {k}$ to indicate a point
$(k_x, k_y, k_z)\in \mathscr {D}$, where
$\mathscr {D}$ is the lattice defined by
$\mathscr {D}=\{(2{\rm \pi} l /(2 L_x), 2 {\rm \pi}m/ (2 L_y) , 2 {\rm \pi}n / (2 L_z)) : (l,m,n) \in \mathbb {Z}^3 \}$, we can consider a function
$f : \mathcal {D} \times [0, + \infty )\rightarrow \mathbb {R}$, periodic over
$\mathcal {D}$, so that it admits the Fourier representation
$f(\boldsymbol {x} ,t)=\sum _{\boldsymbol {k} \in \mathscr {D}} f_{\boldsymbol {k}}(t) \exp (\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x})$. The action of the operators
$\mathcal {J}_{0s}$ and
$\mathcal {J}_{1s}$ on the function
$f$ is defined by


where $J_0$ and
$J_1$ indicate the zeroth- and first-order
$J$ Bessel functions, respectively,
$a_{s}=k_{\perp } \rho _{\perp _s}$ with
$k_{\perp }$ and
$\rho _{\perp _s}$ corresponding to the perpendicular wave number and the Larmor radius of the particle of species
$s$. The former is defined as
$k_{\perp }=\sqrt {k_x^2+k_y^2}$ with
$k_x=2 {\rm \pi}l/(2 L_x)$,
$k_y=2 {\rm \pi}m/(2 L_y)$ for
$l,m \in \mathbb {Z}$, while the latter is given by
$\rho _{\perp _s}=v_{\perp } / \omega _{cs}$, where
$\omega _{cs}=e B_0/ (m_s c)$ is the cyclotron frequency referred to the guide field and related to the particle of species
$s$.
The leading-order expression (up to second-order terms in the perturbations) for the magnetic field is given by

where $\widehat{\boldsymbol{z}}$ is the unit vector along the
$z$ direction,
$\tilde {A}_{\parallel }$ (referred to as magnetic flux function) corresponds to the
$z$ component of the magnetic vector potential and
$\tilde {B}_\parallel$ is the perturbation of the magnetic guide field, also referred to as parallel magnetic perturbation or parallel magnetic fluctuations. We remark that the guide field is assumed to be spatially homogeneous. This assumption is valid for a local description of space plasmas such as the solar wind, where the background magnetic field varies on scales so large that, in the local description, it can be assumed to be homogeneous. The situation would be different, for instance, in the case of tokamak plasmas. Indeed, gyrofluid models more oriented towards tokamak applications (as, for instance those of Brizard Reference Brizard1992; Snyder & Hammett Reference Snyder and Hammett2001; Scott Reference Scott2010; Waelbroeck & Tassi Reference Waelbroeck and Tassi2012; Madsen Reference Madsen2013; Keramidas Charidakos, Waelbroeck & Morrison Reference Keramidas Charidakos, Waelbroeck and Morrison2015) take into account background magnetic inhomogeneities. The set of electromagnetic quantities involved in the system is completed by the electrostatic potential
$\tilde {\phi }=\tilde {\phi }(x,y,z,t)$. In (2.2)–(2.4) we adopted the symbol
$\textrm {d} \mathcal {W}_s=(2{\rm \pi} B_0/m_s)\,\textrm {d} \mu _s\, \textrm {d}v_{\|}$ to indicate the volume element in space velocity.
The parameters $\Theta _s$ and
$\beta _{\perp _s}$ are defined by

and measure, for each species $s$, the equilibrium temperature anisotropy and the ratio between equilibrium kinetic and magnetic pressure, respectively.
Equation (2.1) is the gyrokinetic equation related to the species $s$, whereas (2.2)–(2.4) relate the gyrocentre distribution functions to electromagnetic quantities in the non-relativistic limit. In particular, (2.2) corresponds to the quasi-neutrality relation, whereas (2.3) and (2.4) derive from Ampère's law projected along directions parallel and perpendicular to the guide field, respectively.
The gyrokinetic model (2.1)–(2.4) is valid for small perturbation of the equilibrium distribution function ($\delta f$-gyrokinetic approximation) and weak variations of the fields along the direction of the guide field, the equilibrium temperature anisotropy
$\varTheta _s$ and the parameter
$\beta _{\perp _s}$ of all the species, being kept finite in this asymptotics. Further details about the derivation the regime of validity of the model and its derivation can be found in Kunz et al. (Reference Kunz, Schekochihin, Chen, Abel and Cowley2015). Its Hamiltonian structure, on the other hand, is presented in Tassi (Reference Tassi2019).
3. The gyrofluid model
We define the following gyrofluid moments:


For each particle species, the fields $\tilde {N}_s$ and
$\tilde {U}_s$ represent the fluctuations of the gyrocentre densities and parallel fluid velocities, respectively. On the other hand,
$\tilde {T}_{\parallel s}$ and
$\tilde {T}_{\perp s}$ correspond to the fluctuations of the gyrofluid temperatures defined with respect to the parallel and perpendicular gyrocentre velocities, respectively, whereas
$\tilde {Q}_{\parallel s}$ indicates the gyrocentre parallel heat flux fluctuations. In defining the parallel temperature and heat flux fluctuations, we introduced the constant
$v_{{\textrm {th}}_{\parallel s}}=\sqrt {T_{{0 }_{\parallel s}} / m_s}$, corresponding to the parallel thermal velocity associated with the species
$s$.
We intend to derive a gyrofluid model by taking moments of the gyrokinetic equation (2.1) and by imposing a closure relation derived from the quasi-static linear theory. In particular, the gyrofluid model, accounting for equilibrium temperature anisotropy, should be able, in the limit of vanishing finite Larmor radius effects, to reproduce the field-swelling instability criterion of Basu & Coppi (Reference Basu and Coppi1984). We restrict ourselves to the evolution of the first three moments referring to the parallel direction. Therefore, the resulting gyrofluid model should evolve the following six fields: $\tilde {N}_e, \tilde {N}_i, \tilde {U}_e, \tilde {U}_i, \tilde {T}_{\parallel e}$ and
$\tilde {T}_{\parallel i}$. Also, we show that the model conserves the total energy and, moreover, that it possesses a non-canonical Hamiltonian structure, as is the case for the parent gyrokinetic model (Tassi Reference Tassi2019). We refer to such model as to the six-field gyrofluid (GF6) model.
For the sake of the comparison, carried out in § 4, of the linear gyrofluid theory with the linear gyrokinetic theory we also consider an extension of GF6 accounting for a Landau-fluid closure, analogously to that discussed in Hammett & Perkins (Reference Hammett and Perkins1990), Hammett et al. (Reference Hammett, Dorland and Perkins1992), Snyder, Hammett & Dorland (Reference Snyder, Hammett and Dorland1997), Passot & Sulem (Reference Passot and Sulem2007, Reference Passot and Sulem2015) and Tassi et al. (Reference Tassi, Grasso, Borgogno, Passot and Sulem2018). This variant of the model, which we denote as GF6L, differs from GF6 for the expression of the parallel heat flux fluctuations $\tilde {Q}_{\parallel s}$. Therefore, in order to avoid some redundancy in the exposition, in the following we present the model equations leaving
$\tilde {Q}_{\parallel s}$ unspecified and we will subsequently indicate the corresponding expressions for
$\tilde {Q}_{\parallel s}$ leading to GF6 and GF6L, respectively. The closure leading to GF6L, in particular, will be given in § 3.2.
A further variant of the model, denoted as GF4, will also be considered in § 4. This model evolves only the four fields $\tilde {N}_e, \tilde {N}_i, \tilde {U}_e$ and
$\tilde {U}_i$ and represents a minimal Hamiltonian model, derived from the quasi-static closure, capable to reproduce the field-swelling instability criterion. This variant of the model will be introduced in § 3.2.
The six-field system, both in the Hamiltonian (GF6) and dissipative (GF6L) versions, can be written, in a dimensionless form, as






In (3.3)–(3.8), we adopted the following normalized quantities




where the quantities with overbars in (3.12a–d) are the dimensional spatial and time coordinates. In (3.9a,b)–(3.12a–d) the quantities (note that, according to a customary notation, in the symbols $c_{s_{\perp }}$ and
$\rho _s$, the subscript
$s$ is to indicate sonic quantities and not the particle or gyrocentre species.)

were introduced, which indicate the sound speed and the sonic Larmor radius, respectively, based on the perpendicular equilibrium temperature.
The parameter

for $s=e,i$, on the other hand, measures the ratio of the equilibrium perpendicular temperatures.
In (3.3)–(3.8) we introduced the operator $b_s= -\nabla _{\perp }^2 \rho _{\textrm {th}_{ \perp s}}^2$, with
$\nabla _{\perp }^2$ denoting the Laplacian relatively to the transverse variables and

indicating the perpendicular thermal Larmor radius associated with the species $s$.
The canonical bracket $[ \, , \, ]$, on the other hand, is defined as
$[\,f,g]=\partial _x f \partial _y g - \partial _y f \partial _x g$, for two functions
$f$ and
$g$.
The gyroaverage operators $G_{10s}$,
$G_{20s}$,
$\Gamma _{0s}$,
$\Gamma _{1s}$ in (3.3)–(3.8) can in turn be expressed in terms of the operators


which have to be intended as Fourier multipliers whose symbols are obtained by replacing $b_s$ by
$k_{\perp }^2 \rho ^2_{\textrm {th} \perp s}$. We make this statement more precise, as an example, in the case of the operator
$G_{10i}$, referred to the ion species. The expression
$G_{10i} f$ is defined by
$G_{10i} f (\boldsymbol {x} ,t)=\sum _{\boldsymbol {k} \in \mathscr {D}} G_{10} (b_i) f_{\boldsymbol {k}}(t) \exp (\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x})=\sum _{\boldsymbol {k} \in \mathscr {D}} \exp (-k_{\perp }^2 \rho _{\textrm {th}_{\perp _i}}^2/2) f_{\boldsymbol {k}}(t) \exp (\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x})$, for a function
$f$ periodic in space. Analogous expressions are valid for the other gyroaverage operators. In (3.17a,b), in particular, the symbols
$I_0$ and
$I_1$ indicate the modified Bessel functions
$I$ of order zero and one, respectively.
The expressions for the operators $G_{10s}$ and
$G_{20s}$ given in (3.16a,b) correspond to those present in Brizard (Reference Brizard1992) and follow from assuming that the perturbation of the distribution function can be written as

where $H_m$ and
$L_n$ indicate the Hermite and Laguerre polynomials, respectively, of order
$m$ and
$n$, with
$m$ and
$n$ non-negative integers. The functions
$f_{{mn}_s}$ are coefficients of the expansion and are related to the moments of
$\tilde {f}_{s}$, with respect to Hermite polynomials in
$v_{\|} / v_{{\textrm {th}}_{\parallel s}}$ and Laguerre polynomials in
$\mu _s B_0 / T_{{0 }_{\perp s}}$. Indeed, from the orthogonality properties of Hermite and Laguerre polynomials, the following relation holds:

It can be useful, in particular, to write explicitly the following relations between the lowest-order normalized moments and the quantities $f_{{mn}_s}$:



where, in (3.21a,b), we introduced the normalized gyrocentre perpendicular temperature fluctuations $T_{\perp s}=\tilde {T}_{\perp s}/T_{{0 }_{\perp s}}$.
As above anticipated, in the system (3.3)–(3.8) we temporarily left unspecified the expression for the parallel heat flux fluctuations $Q_{\|s}$. In order to obtain the model GF6, the infinite hierarchy of gyrofluid equations following from the parent gyrokinetic model (2.1)–(2.4) is closed by imposing relations obtained by computing the gyrocentre moments other than
$N_s$,
$U_s$ and
$T_{\|s}$ from a linearization of the parent gyrokinetic system about a homogeneous equilibrium, in the quasi-static limit. In the case of
$Q_{\parallel s}$, this leads (consider (3.22) and (A 21)) to

The model GF6 is thus obtained by inserting the relation (3.23) into the system (3.3)–(3.8). Details on the derivation of the closure relations originated from the quasi-static assumption can be found in appendix A. A remarkable property of this closure is that it leads to the annihilation of all the contribution of the higher-order moments in the gyrofluid equations.
We find it useful to provide also a reformulation of the evolution equations (3.3)–(3.5), which should help putting in evidence the physical nature of the terms contributing to the evolution of the various fields. Equations (3.3)–(3.5) can indeed be rewritten as



where the parallel gradient operator $\nabla _{\parallel s}$ is defined, for each species
$s$, by

for a function $f$. From the formulation (3.24)–(3.26) it emerges that the gyrocentre density, parallel momentum and parallel temperature fluctuations are all advected, in the perpendicular plane, by the incompressible velocity field

Such velocity field includes a first contribution, associated with $G_{10s} \phi$, which corresponds to the usual
$\boldsymbol {E}\times \boldsymbol {B}$ drift (with E indicating the electric field based on the gyroaveraged electrostatic potential), ubiquitous in low-
$\beta$ gyrofluid models such as those discussed in Waelbroeck & Tassi (Reference Waelbroeck and Tassi2012), Snyder & Hammett (Reference Snyder and Hammett2001), Keramidas Charidakos et al. (Reference Keramidas Charidakos, Waelbroeck and Morrison2015) and Waelbroeck, Hazeltine & Morrison (Reference Waelbroeck, Hazeltine and Morrison2009). When higher
$\beta$ values are allowed, however, the perpendicular advection acquires a further contribution due to the parallel magnetic perturbations, as it transpires from (3.28). We remark that, similarly to the
$\boldsymbol {E}\times \boldsymbol {B}$ contribution, also the latter contribution does not vanish in the limit
$b_s \rightarrow 0$ of negligible FLR corrections. This is a consequence of the fact that, in the presence of parallel magnetic perturbations, gyrocentre density fluctuations differ from particle density fluctuations, even in the absence of FLR corrections (Brizard Reference Brizard1992).
From the continuity equation (3.24), we see that the evolution of gyrocentre density fluctuations has also a source in the last term on the left-hand side, which is due to the compressibility of the gyrocentre mean velocity along the direction of the magnetic field.
From (3.25) we see that parallel momentum, in addition to be advected by $\boldsymbol {u}_{\perp _s}$, evolves due to the term
$-\textrm {sgn} (q_s)E_{\parallel s}$, where
$E_{\parallel s}$ is defined by

Such a term represents the force exerted by the gyroaveraged electric field, along the direction of the magnetic field. A further source for the parallel momentum is due to the term $\nabla _{\parallel s} \mathcal {P}_s$, where

This term is associated with the parallel component of the divergence of an anisotropic pressure tensor.
The parallel temperature equation (3.26) has, on its left-hand side, the same structure of the continuity equation (3.24). We just remark the presence of the coefficient $2$ multiplying
$\nabla _{\parallel s} U_s$. This coefficient, of course, follows directly from taking the second order moment, in Hermite polynomials for the normalized parallel velocity, of the gyrokinetic equation (2.1), as discussed in appendix B. However, as pointed out by Keramidas Charidakos et al. (Reference Keramidas Charidakos, Waelbroeck and Morrison2015), in the presence of background magnetic curvature such coefficient, in general, has to be adjusted in order to obtain a Hamiltonian structure. Finally, we consider the term on the right-hand side of (3.26), associated with the parallel heat flux. If the expression for
$Q_{\parallel s}$ is chosen according to the quasi-static closure, i.e. imposing the relation (3.23), this term vanishes and the system, as will be shown in § 3.1, is Hamiltonian. On the other hand, if the Landau fluid closure (3.35) is chosen, this term acts as a sink and the system is not energy conserving.
3.1. Hamiltonian structure of GF6
The quasi-static closure (3.23) leading to GF6, allows the resulting model to be cast in Hamiltonian form. In particular, it can be verified by direct computation that, when the electromagnetic fluctuations $\phi$,
$A_{\parallel }$ and
$B_\parallel$ can be expressed in terms of the dynamical variables
$N_e$,
$N_i$,
$M_e$ and
$M_i$ (where we introduced the short-hand notation
$M_s=(m_s/m_i)U_s+\textrm {sgn}(q_{s})G_{10s} A_{\parallel }$ to indicate the parallel canonical momenta), by making use of the relations (3.6)–(3.8), the evolution equations (3.3)–(3.5) complemented by (3.23), can be written in the Hamiltonian form

with $s=e,i$. In (3.31a–c)
$H$ is the Hamiltonian functional

and $\{ \, , \, \}$ is a non-canonical Poisson bracket given by

for two functionals $F$ and
$G$. For details about the non-canonical Hamiltonian formulation of fluid models one can refer, for instance, to Morrison (Reference Morrison1998). In (3.33) the subscripts on the functionals indicate functional derivatives. In order to verify the formulation (3.31a–c) it is convenient to remark that, from (3.32), one obtains

In order to derive the relations (3.34a–c) we made use of the formal symmetry of the operators $G_{10s}$ and
$G_{20s}$, i.e.
$\int \textrm {d}^3x \, f G_{10s} g = \int \textrm {d}^3 x \, g G_{10s} f$ and
$\int \textrm {d}^3x \, f G_{20s} g = \int \textrm {d}^3 x \, g G_{20s} f$ for two functions
$f$ and
$g$, as well as of the formal symmetry of the linear operators in terms of which one can express
$\phi , B_\parallel$ and
$A_{\parallel }$ in terms of
$N_s$ and
$M_s$ through (3.6)–(3.8). We did not provide the explicit expression for such operators which can, however, be obtained considering the representation in Fourier series of the fields involved, following the procedure discussed in Tassi (Reference Tassi2019).
The Hamiltonian functional (3.32) is a conserved quantity for the dynamics and corresponds to the total energy. In (3.33), the sum of all the terms with $s=e$ is a Poisson bracket in its own right. Similarly, all the terms with
$s=i$ form a Poisson bracket. The sum of these two contributions is a direct sum of Poisson brackets which is in turn a Poisson bracket verifying in particular the Jacobi identity. The Poisson brackets referring to the electron and ion quantities correspond to those already discussed in other Hamiltonian reduced fluid models. The reader can in particular refer to Tassi (Reference Tassi2015) and Keramidas Charidakos et al. (Reference Keramidas Charidakos, Waelbroeck and Morrison2015) for the verification of the Jacobi identity for brackets of such form and for a discussion of the corresponding Casimir invariants. The model GF6, in the two-dimensional limit when the dependence on the
$z$ coordinate is suppressed, can also be cast in the form of a system of advection of equations for Lagrangian invariants, as is the case for several other reduced fluid and gyrofluid models (see, e.g. Schep, Pegoraro & Kuvshinov Reference Schep, Pegoraro and Kuvshinov1994; Waelbroeck et al. Reference Waelbroeck, Hazeltine and Morrison2009; Waelbroeck & Tassi Reference Waelbroeck and Tassi2012; Grasso & Tassi Reference Grasso and Tassi2015; Keramidas Charidakos, Waelbroeck & Morrison Reference Keramidas Charidakos, Waelbroeck and Morrison2015; Tassi Reference Tassi2015, Reference Tassi2019).
Remark The gyroaverage operators $G_{10}$ with a form different from (3.16a,b) have been proposed in the literature (see, e.g. Dorland & Hammett Reference Dorland and Hammett1993; Snyder & Hammett Reference Snyder and Hammett2001; Scott Reference Scott2010), when the expansion (3.18) is not assumed, and are frequently adopted. In particular,
$G_{10} (b_s)=\varGamma _0^{1/2} (b_s)$ was shown to provide better agreement with the linear theory at large
$b_s$ (Dorland & Hammett Reference Dorland and Hammett1993). We point out, however, that many important features of the model (3.3)–(3.8), such as the total energy conservation and the Hamiltonian structure, are guaranteed whatever the form of the operators
$G_{10}$ and
$G_{20}$ of (3.16a,b) is, provided these operators are linear and formally symmetric, in the sense defined above. This is in particular the case with
$G_{10} (b_s)=\varGamma _0^{1/2} (b_s)$. These issues are also discussed by Mandell, Dorland & Landreman (Reference Mandell, Dorland and Landreman2018).
3.2. Variants of the model
3.2.1. Six-field model with Landau closure (GF6L)
The variant GF6L of the six-field gyrofluid model, accounting for Landau damping, corresponds to the system (3.3)–(3.8) with $Q_{\parallel s}$ given by

In (3.35) we introduced the constant $\alpha _s=(2/{\rm \pi} )^{1/2}(m_i/m_s)^{1/2}$. The operator
$\mathcal {L}$ holds for the Landau damping operator. Its modelling in the nonlinear regime is discussed in Tassi et al. (Reference Tassi, Grasso, Borgogno, Passot and Sulem2018). In the linear approximation, it reduces to the negative of the Hilbert transform in the direction of the ambient magnetic field (here taken in the
$z$ direction). The presence of this Landau operator in reduced fluid models breaks the Hamiltonian structure by violating energy and Casimir conservation (Tassi et al. Reference Tassi, Grasso, Borgogno, Passot and Sulem2018; Grasso et al. Reference Grasso, Borgogno, Tassi and Perona2020). Its purpose, on the other hand, is to introduce terms that allow the linear dispersion relation of the gyrofluid model to reproduce that of the parent gyrokinetic model.
3.2.2. Four-field model with quasi-static closure (GF4)
The second variant GF4 is obtained by retaining the evolution equations for $N_s$ and
$M_s$ and imposing the quasi-static closure on the parallel temperature fluctuations
$T_{\parallel s}$. Considering (3.21a,b) and (A 21), this amounts to setting

The resulting model reads





and corresponds to taking (3.3), (3.4), (3.6), (3.7) and (3.8) of GF6 with $T_{\parallel s}=0$. Its Hamiltonian structure is given by the Hamiltonian

and by the Poisson bracket

If one neglects electron FLR effects (i.e. $b_e \rightarrow 0$), parallel magnetic perturbations (i.e.
$B_\parallel =0$), equilibrium temperature anisotropies (i.e.
$\varTheta _e=\varTheta _i=1$) and sets
$G_{10} (b_i)=\varGamma _{0}^{1/2} (b_i)$ (i.e. takes the alternative form of the ion gyroaverage operator mentioned in the above remark in § 3.1), GF4 reduces to the Hamiltonian gyrofluid model of Waelbroeck & Tassi (Reference Waelbroeck and Tassi2012), the latter taken in the limit of vanishing magnetic curvature and equilibrium density gradients. We note, however, that, once that the quasi-static closure relations are determined, as in (A 21), all the terms in GF4 (and likewise for GF6), are determined exactly. In particular no approximations of the gyroaverage operators (unlike, for instance, in Scott Reference Scott2010; Waelbroeck & Tassi Reference Waelbroeck and Tassi2012) are carried out. Terms involving gyroaverage operators are determined exactly also in Brizard (Reference Brizard1992), but without making use of the quasi-static closure. As a result, in Brizard (Reference Brizard1992), nonlinear terms involving more than one gyroaverage operator do not result in having the single canonical bracket structure (as is the case in GF4, for instance, with the term
$\textrm {sgn} (q_{s}) [ G_{10s} \phi ,G_{10s} A_{\parallel }]$ appearing in the second line of (3.38)) and which is crucial for determining the Lie–Poisson Hamiltonian structure. This property, which follows from the quasi-static closure, differentiates the models presented in our paper also from its closest predecessors, i.e. the gyrofluid models constructed with the technique recently presented in Tassi (Reference Tassi2019). The latter gyrofluid models, in fact, are also Hamiltonian and account for equilibrium temperature anisotropy, but adopt a different closure. Namely, all gyrofluid moments involving finite powers of the magnetic moment
$\mu _s$ (e.g. the perpendicular temperature fluctuations) are set equal to zero. This allows for a Hamiltonian structure but the terms involving gyroaverage operators are not all determined exactly from taking the moments of the gyrokinetic equations. For this reason, we think that, for situations where the quasi-static assumption
$|\omega /(k_z v_{{\textrm {th}}_{\parallel s}})| \ll 1$ is satisfied, the models presented in this paper might be preferable to the models described in Tassi (Reference Tassi2019).
4. Comparison with the linearized gyrokinetic parent model
4.1. KAW dispersion relation
We first discuss the KAW dispersion relation as predicted by the gyrofluid models, with and without Landau damping (GF6L and GF4 respectively), in comparison with the predictions of the linearized parent model (Kunz et al. Reference Kunz, Abel, Klein and Schekochihin2018) which identifies with the low-frequency kinetic theory described in Kuznetsov, Passot & Sulem (Reference Kuznetsov, Passot and Sulem2012). This latter dispersion relation involves the plasma response function $R(\zeta _s)$ of the particles of species
$s$, which is related to the corresponding plasma dispersion function
$Z(\zeta _s)$ by
$R(\zeta _s) = 1 + \zeta _s Z(\zeta _s)$ with
$\zeta _s = \omega /(2^{1/2}|k_z|v_{\textrm {th}_{\| s}})$. Different Padé approximants
$R_{ij}(\zeta _e)$ (for which we follow the notations of Hunana et al. Reference Hunana, Tenerani, Zank, Goldstein, Webb, Velli and Adhikari2019) are used to estimate the electron response function.
In order to test the validity of the quasi-static assumption that affects the form of the FLR terms, independently from the effects resulting from Landau damping, we are led to compare the prediction of the four-field gyrofluid model GF4 with the gyrokinetic dispersion relation obtained by choosing for the electron plasma response function, the function $R_{21}(\zeta _e)=1/(1-2\zeta _e^2)$ (model denoted GKNL). This choice directly results from the assumption
${T_{{\|e}}}=0$ (see (B7) of Passot & Sulem Reference Passot and Sulem2007) and ignores electron Landau damping.Footnote 1 Predictions of the GF6L model will, on the other hand, be compared with the full gyrokinetic dispersion relation, referred to as GK. In the latter description, we use for the electron Padé approximant the function
$R_{86}(\zeta _e)$, which shows an excellent agreement with the exact response function. For the ions, we use in all the cases the Padé approximant
$R_{20}(\zeta _i)=1/(1-2\zeta _i^2-\textrm {i}{\rm \pi} ^{1/2}\zeta _i)$ (or
$R_{21}(\zeta _i)$ in the absence of Landau damping). Higher-order Padé approximants give almost identical results.
Because we are making use of several acronyms to refer to the different adopted models, we summarized them in table 1. In all the examples presented in this section, $\varTheta _i=1$. Introducing the ion to electron parallel temperature ratio at equilibrium
$\tau _{\|} = {T_{0_{\|i}}}/{T_{0_{\|e}}}=\tau _{\perp _ i}\varTheta _e/{\varTheta _i}$ and denoting by
$\alpha$ the angle between the wave vector and the guide field (propagation angle), we compare, in figure 1(a), the (real) normalized KAW frequency
$\omega /(k_z c_A)$ for
$\beta _{\perp e}=1$,
$\tau _{\|}=1$,
$\varTheta _e=1$,
$\alpha =89^{\circ }$ as a function of
$k d_i$ (where
$k=\sqrt {k_{\perp }^2+k_z^2}$ and
$d_i$ is the ion inertial length defined by
$d_i = c_A/\omega _{ci}$, with
$c_A =\sqrt {B_0/(4{\rm \pi} m_i n_0)}$ the Alfvén velocity), calculated using GF4 (black diamond symbols), with the GKNL prediction (green solid line). An excellent agreement is found, the slight deviation appearing as
$kd_i>40$ probably resulting from the failure of the quasi-static approximation used in the calculation of all the moments starting from the temperature. For the same plasma parameters, the case with Landau damping is displayed as a red solid line for GK and blue diamond symbols for GF6L, with a good agreement up to scales
$kd_i \lessapprox 20$. The damping rate is displayed with the same symbols in figure 1(b). Interestingly, the prediction of GF6 (not shown) departs from GKNL significantly at all the scales. Indeed, in GF6, the parallel temperature fluctuations that obey a dynamical equations with zero heat flux do not approach a quasi-static dynamics. Such a non-dissipative odd-order closure rather fits an adiabatic regime (see footnote in this section). Landau damping is requested to ensure convergence to a quasi-static regime. In that sense, GF4 is preferable to GF6 for addressing a non-dissipative problem, the interest of the latter model being mainly to provide a framework where Landau damping can be supplemented to an otherwise Hamiltonian description.

Figure 1. (a) Normalized real part of the KAW frequency $\omega /(k_z c_A)$ for
$\beta _{\perp e}=1$,
$\tau _{\|}=1$,
$\varTheta _e=1$,
$\alpha =89^{\circ }$, as a function of
$k d_i$, calculated using GF4 (black diamond symbols) and GF6L (blue diamond symbols) models (respectively without and with Landau damping) and with the gyrokinetic dispersion relations GKNL (green solid line) and GK (red solid line). (b) Damping rate from the GK theory (red solid line) and from the GF6L model (blue diamond symbols). (c) Normalized growth rate of the firehose instability as a function of
$k d_i$, for
$\beta _{\perp e}=0.5$,
$\tau _{\|}=10^{-3}$,
$\varTheta _e=0.16$,
$\alpha =89^{\circ }$, as predicted by GF4 and GF6L models, together with GKNL and GK dispersion relations, using the same graphic conventions.
Table 1. Acronyms and corresponding models.

4.2. Firehose instability
Figure 1(c) displays for $\beta _{\perp e}=0.5$,
$\tau _{\|}=10^{-3}$,
$\alpha =89^{\circ }$,
$\varTheta _e=0.16$, the growth rate of the firehose instability as a function of
$kd_i$ predicted by GF4 and GF6L (black and blue diamond symbols respectively), together with the GKNL and GK dispersion relations (green and red solid lines). In all the cases, the agreement is excellent, the quenching of the instability being in particular well reproduced.
The agreement found between the KAW dispersion relations predicted by the gyrofluid models and the gyrokinetic theory, even when limited to scales such that $kd_i \lessapprox 20$, and despite a large value of
$\zeta _i = \sqrt {2\varTheta _e/(\beta _{\perp _e}\tau _{\|})}\;\omega /(|k_z|c_A)$, can be attributed to the fact that, at least within the linear theory, ion acoustic and kinetic Alfvén waves remain essentially decoupled. The influence of the ion closure relation on the KAW properties thus remains limited. Deviations from the gyrokinetic theory at small scales are mostly due to the fact that when reaching these scales
$\zeta _e=\sqrt {m_e/m_i}\sqrt {2\varTheta _e/\beta _{\perp _e}}\;\omega /(|k_z|c_A)$ becomes non-negligible.
4.3. SW dispersion relation and field-swelling instability
Figure 2 concerns a similar comparison in the case of the field swelling instability discussed in appendix C, for $\beta _{\perp e}=1$,
$\tau _{\|}=1$,
$\alpha =89^{\circ }$,
$\varTheta _e=2.2$ (panel a). For both the cases with and without Landau damping, the relatively good agreement found at large scale between the gyrofluid theories and the gyrokinetic ones deteriorates at smaller scales. The stabilization scale is, however, correctly captured. This discrepancy is associated with a value of
$\zeta _i$, which is not small enough (in this case,
$\zeta _e$ remains reasonably small). We show in figure 2(b) with
$\varTheta _e=2.41$, again with
$\beta _{\perp e}=1$ and
$\alpha =89^{\circ }$, that a much better agreement can be found when ions are hotter (
$\tau _{\|}=50$), which prescribes a small
$\zeta _i$. The case with cold ions (
$\tau _{\|}=10^{-5}$), for which the ion dynamics is decoupled, is displayed in figure 3, showing an even better agreement. Figure 3(a) displays the real part of the slow wave for
$\varTheta _e=1$, obtained with the GKNL model (solid green line) or the GF4 model (diamond symbols). Figure 3(b) shows the growth rate of the field-swelling instability for
$\varTheta _e=2.01$ (keeping unchanged the other parameters).

Figure 2. Normalized growth rate of the field-swelling instability versus $kd_i$ for
$\alpha =89^{\circ }$ and
$\beta _{\perp e}=1$. (a) Predictions of GF4 and GF6L, compared with those of GKNL and GK respectively, for
$\tau _{\|}=1$ and
$\varTheta _e=2.2$. (b) Prediction of GF4 compared with that of GKNL, for
$\tau _{\|}=50$ and
$\varTheta _e=2.41$. Same graphic conventions as in figure 1 are used.

Figure 3. (a) Normalized slow-wave frequency versus $kd_i$ in regimes where the ions are cold (
$\tau _{\|}=10^{-5}$), with
$\beta _{\perp e}=1$ and
$\alpha =89^{\circ }$, as predicted by the GF4 model and the GKNL dispersion relation. (b) Normalized growth rate of the swelling instability for
$\tau _{\|}=10^{-5}$ and
$\varTheta _e=2.01$. Same conventions as in figure 1 are used.
5. Conclusion
We derived a six-field Hamiltonian gyrofluid model, referred to as GF6, retaining the gyrocentre density, the parallel velocity and temperature fluctuations for each species, under the sole assumption that all the other gyrocentre moments are calculated from the quasi-static linear kinetic theory. Such an assumption on the closure turns out to yield exact expressions for all the terms of the model, without, in particular, requiring approximated expressions for the terms involving gyroaverage operators. Nonlinear terms involving more than one gyroaverage operator, in particular, appear in the form of a single canonical bracket, which naturally lets the model fit in the class of Hamiltonian models with a Lie–Poisson structure. The model accounts for equilibrium temperature anisotropy and also retains both ion and electron FLR corrections, electron inertia and parallel magnetic fluctuations. In a variant of the model (GF6L) parallel Landau damping is retained through a Landau-fluid modellization of the gyrocentre parallel heat fluxes. A second variant of the model (GF4) is obtained by prescribing parallel isothermality, which still falls in the frame of the quasi-static closure and allows for a Hamiltonian formulation. The comparison of the dispersion relations of KAWs and SWs predicted from GF6L or GF4, with those derived from the parent gyrokinetic theory where the plasma response function is replaced by a Padé approximant, provides an estimate of the maximal transverse wavenumber beyond which the phase velocity of the corresponding wave is too large compared with the electron (in the former case) or ion (in the latter case) parallel thermal velocity for consistency with a closure condition based on a quasi-static assumption. It turns out that the agreement extends to transverse scales significantly smaller than the ion Larmor radius in the case of KAWs, mostly because, at least at the linear level, SWs and KAWs are essentially decoupled, making the influence of the ion closure relation on the KAW properties relatively limited. This situation contrasts with the case of the SWs for which the dispersion relation is accurately reproduced only at scales larger than a significant fraction of the ion Larmor radius. Under these conditions, the model reproduces the instabilities induced by temperature anisotropy, such as firehose or field-swelling instabilities. It should nevertheless be noted that, as it assumes small perturbations of an equilibrium state, the model does not permit evolution of the mean temperatures, an effect usually considered as contributing efficiently to the saturation of these instabilities. The subcritical nonlinear regime is, however, expected to be accurately described. The model will in particular be most useful for studying the coupling of KAWs with SWs which can generate large-scale parametric decay instabilities at small $\beta _e$, a regime especially relevant in the regions of the solar wind relatively close to the Sun explored by space missions such as Parker Solar Probe or Solar Orbiter.
In general, to the best of our knowledge, our gyrofluid model is the only one, at the present moment, possessing the following features, which could make it a valuable tool for local investigations of basic plasma phenomena of interest for space plasmas: it accounts for equilibrium temperature anisotropies as well as parallel magnetic perturbations; it reproduces, in a rather wide range of values of parameters, compatible with the quasi-static closure, quantitative features of known kinetic linear dispersion relations; the model equations, and in particular the terms involving FLR corrections, are calculated exactly, unlike other gyrofluid models which adopt truncations or approximations of such terms; it possesses a Hamiltonian structure.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of closure relations from the gyrokinetic linear theory
We consider the linearization of the gyrokinetic system (2.1)–(2.4) about the equilibrium state $\tilde {g}_s=0$ (or, equivalently,
$\tilde {f}_{s}=0$), with
$\tilde {\phi }=\tilde {A}_{\parallel }=\tilde {B}_\parallel =0$. The resulting linear system can be written in the form




where we adopted the same notation with the tilde symbol that we used in (2.1)–(2.4) to indicate the dynamical variables $\tilde {f}_{s}$ of the linearized system and the field perturbations
$\tilde {\phi }, \tilde {A}_{\parallel }$ and
$\tilde {B}_\parallel$.
We introduce the following Fourier series representation:


with $\omega \in \mathbb {C}$ indicating the complex frequency.
For any given $\boldsymbol {k} \in \mathscr {D}$, from (A 1) we obtain the relation

where $\tilde {\zeta _s}=\omega /(k_z v_{{\textrm {th}}_{\parallel s}})$.
We consider now the quasi-static limit $\vert \tilde {\zeta _s} \vert \ll 1$. In this limit, the relation (A 7) reduces to

In the following, we make use of the relation (A 8), derived from the linear theory, in order to determine the closure relations to insert in the hierarchy of nonlinear gyrofluid equations. For this purpose, we adopt the Hermite–Laguerre expansion of (3.18) for the perturbation of the distribution function in the linearized system and we write

with

in Fourier representation. From (A 9), (A 10), using the orthogonality relations for Hermite and Laguerre polynomials, one obtains

Inserting the relation (A 8) into (A 11) and using the orthogonality relations for Hermite polynomials, one has

where the operators $G_{1n}$ and
$G_{2n}$ are defined by


with

Explicit expressions for the operators $G_{1n}$ and
$G_{2n}$ can be found by computing the integrals in (A 13) and (A 14), which yields


In order to obtain the expressions (A 16)–(A 17a,b) use was made of the orthogonality of Laguerre polynomials as well as of the relations (Szegö Reference Szegö1975)


where $L_n^{(1)}$ is a generalized Laguerre polynomial. Making use of the Fourier representations (A 10), (A 5a,b) and (A 6a,b) for
$f_{{mn}_s}$,
$\tilde {\phi }$ and
$\tilde {B}_\parallel$, respectively, one can deduce from (A 12) the relation

or, equivalently,

where we also made use of the normalization (3.11a–c) for $\phi$ and
$B_\parallel$. The operators
$G_{1ns}$ and
$G_{2ns}$ are defined, consistently with the definition of
$G_{10s}$ and
$G_{20s}$ given in § 3, by
$G_{1ns} f (\boldsymbol {x} ,t)=\sum _{\boldsymbol {k} \in \mathscr {D}} G_{1n} (b_s) f_{\boldsymbol {k}}(t) \exp (\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x})$ and
$G_{2ns} f (\boldsymbol {x} ,t)=\sum _{\boldsymbol {k} \in \mathscr {D}} G_{2n} (b_s) f_{\boldsymbol {k}}(t) \exp (\textrm {i} \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x})$, for a function
$f$ and
$n\geq 0$ (in the specific case of the linear dispersion relation, the dependence on
$t$ is provided by the factor
$\textrm {e}^{-\textrm {i} \omega t}$, but when the relations (A 21) are used as closures for the nonlinear models, the dependence on
$t$ is of course left arbitrary).
The relations (A 21) descending from the quasi-static assumption, are adopted as closure relations in GF6 and GF6L for all the moments involved in the model, except for $N_s$,
$M_s$,
$T_{\parallel s}$, which are derived by solving the evolution equations (3.3)–(3.5). The parallel heat flux fluctuations
$Q_{\parallel s}$, on the other hand, are determined, as already mentioned, again by a quasi-static closure for GF6 ((3.23) which follows from (A 21) for
$m=3, n=0$) or by the Landau closure (3.35) for GF6L. The closure (3.36) adopted for GF4, is again a quasi-static closure following from (A 21) when
$m=2$ and
$n=0$.
Appendix B. Derivation of the model equations
The gyrofluid system (3.3)–(3.8) is derived from the parent gyrokinetic system (2.1)–(2.4) by making use of the expansion (3.8), applied to the perturbations of the distribution functions. In order to obtain a closed system with a finite number of equations, such expansion is constrained in the following way. The moments $f_{00_s}, f_{10_s}$ and
$f_{20_s}$ (or, equivalently, by virtue of (3.20a,b)–(3.21a,b), the gyrofluid densities, parallel velocities and temperatures
$N_s$,
$U_s$ and
$T_{\parallel s}$), for each species
$s$, are determined by evolution equations obtained by making the product of all the terms of the gyrokinetic (2.1) with the zero-, first- and second-order Hermite polynomial in the variable
$v_{\|}/ v_{{\textrm {th}}_{\parallel s}}$ and integrating over the velocity volume element
$\textrm {d} \mathcal {W}_s$. For GF6L, the parallel heat flux
$Q_{\parallel s}$ is determined by the relation (3.35). All the other moments, on the other hand, are assumed to be given by the relations (A 20), or, in normalized form, by (A 21), obtained from the linear theory in the quasi-static limit. With this prescription, the expansion (3.18) becomes

where we made use of the fact that $H_0 (v_{\|} / v_{{\textrm {th}}_{\parallel s}})=1$,
$H_1(v_{\|} / v_{{\textrm {th}}_{\parallel s}})=v_{\|} / v_{{\textrm {th}}_{\parallel s}}$,
$H_2 (v_{\|} / v_{{\textrm {th}}_{\parallel s}})=v_{\|}^2/v_{{\textrm {th}}_{\parallel s}}^2-1$ and
$H_3(v_{\|} / v_{{\textrm {th}}_{\parallel s}})=v_{\|}^3 / v_{{\textrm {th}}_{\parallel s}}^3 - 3 v_{\|} / v_{{\textrm {th}}_{\parallel s}}$.
Inserting the expansion (B 1) into (2.5) and (2.1), and integrating over $\textrm {d} \mathcal {W}_s$ one obtains

where we made use of the definitions (A 13) and (A 14) and of the orthogonality of Hermite polynomials. We remark, at this point, that the sum of the last term in the first line of (B 2) with the first term on the second line of (B 2) yields zero, because of the antisymmetry of the canonical bracket $[ \, , \, ]$. We thus conclude that the quasi-static closure (A 20) has the remarkable property of annihilating, in the continuity equation, all the contributions associated with the moments
$f_{0n_s}$, for
$n \geq 1$. By virtue of this cancellation, from (B 2) we obtain

Applying to (B 3) the normalization (3.9a,b)–(3.12a–d), one obtains (3.3).
In order to derive (3.4) we point out first that, with the help of the identities




we obtain the relation (see also Brizard Reference Brizard1992)



and, by an analogous procedure, the relation

Upon multiplying (2.1) by $(1/n_0) v_{\|} / v_{{\textrm {th}}_{\parallel s}}$ and integrating over
$\textrm {d} \mathcal {W}_s$ one obtains, adopting the expansion (B 1) as well as the relations (A 13), (A 14), (B 8) and (B 9), the following equation:

Also in this case, the quasi-static closure leads to a remarkable cancellation. Indeed, among the nonlinear terms involving only electromagnetic quantities (i.e. $\tilde {\phi }$,
$\tilde {A}_{\parallel }$ and
$\tilde {B}_\parallel$) all those containing gyroaverage operators
$G_{1ns}$ and
$G_{2ns}$, with
$n \geq 1$, vanish. As a result, (B 10) reduces to

Equation (3.4) is then obtained from (B 11) after applying the normalization (3.9a,b)–(3.12a–d).
Equation (3.5) is obtained upon multiplying (2.1) by $(1/n_0) (v_{\|}^2/ v_{{\textrm {th}}_{\parallel s}}^2 -1)$ and integrating over
$\textrm {d} \mathcal {W}_s$. Making use of the expansion (B 1) and of the orthogonality of Hermite polynomials one obtains

Adopting the normalization (3.9a,b)–(3.12a–d), (3.5) follows from (B 12).
Equations (3.6), (3.7) and (3.8) follow from (2.2), (2.3) and (2.4), respectively, upon inserting the expansion (B 1), evaluating the integrals and applying the normalization (3.9a,b)–(3.12a–d). With regard to the evaluation of the integrals and the derivation of the equations in the form (3.6), (3.7) and (3.8), we remark that, in addition to (A 13)–(A 14), the following relations are of use:



Appendix C. Field-swelling instabilities
C.1. Dispersion relation at the MHD scales
In this appendix, we first provide a simple derivation of the dispersion relation for fast and slow modes at MHD scales in the presence of temperature anisotropy, starting from the kinetic-MHD equations (see equations (37), (38b), (44a–c) (46)–(48) from Kulsrud (Reference Kulsrud, Galeev and Sudan1983), or equivalently equations (1)–(8) from Snyder et al. (Reference Snyder, Hammett and Dorland1997)).
The transverse velocity can be decomposed into compressible and solenoidal parts by writing

One immediately gets (see e.g. equations (48), (52) and (56) of Passot & Sulem (Reference Passot and Sulem2006) where FLR corrections are neglected)


where $p_{\perp }$ and
$p_{\|}$ denote the total (ion plus electron) perpendicular and parallel pressure fluctuations, with the subscript
$0$ referring to the equilibrium values. Furthermore,
$\rho _0$ denotes the equilibrium plasma density.
In (C 2), the perpendicular pressure fluctuations are given by the drift-kinetic theory, as expressed, e.g., by equation (27) of Snyder et al. (Reference Snyder, Hammett and Dorland1997), in the form

together with

where ${q}_r=e$ for ions and
$-e$ for electrons, respectively. The potential
$\psi$ is defined in terms of the parallel electric field by
$E_z = -\partial _z \psi$. These expressions for the perpendicular pressure and density perturbations identify with the large-scale limit of formulas given in appendix B of Passot & Sulem (Reference Passot and Sulem2007). Quasineutrality requires the equality of the electron (
$n_e$) and ion (
$n_i$) number-density fluctuations, which prescribes

Plugging the expressions for the pressures in (C 2)–(C 3) thus leads to

We note that ${\zeta ^2_i={\omega ^2}/{2\tau _{\|} \omega ^2_s}}$ with
$\omega _s= k_z c_{s_{\|}}$ where
${c_{s_{\|}}=\sqrt {{T_{{0 }_{\parallel e}}}/{m_i}}}$ is the sound speed based on the parallel electron temperature.
C.2. Link with the fluid theory
In order to make a link with fluid theory as performed in Basu & Coppi (Reference Basu and Coppi1984), we note that keeping all the terms in the ion density and parallel velocity equations (in particular the time derivatives), is equivalent to expanding $R(\zeta _i)$ for
$\zeta _i$ large, i.e. in the adiabatic limit where
$R(\zeta _i)\approx -1/(2\zeta _i^2)$. It thus follows that the ions cannot be hot, at least outside an angular boundary layer near the transverse direction.
In addition to adiabatic ions, let us also assume a quasi-static limit for the electrons, which leads to assume $\zeta _e \to 0$ and thus
$R(\zeta _e)\approx 1$. Equation (C 7) is then rewritten

For cold ions ($\tau _{\|}=\tau _{\perp _i}=0$), we get

which is also rewritten in the form of equation (44) of Basu & Coppi (Reference Basu and Coppi1984)

The slow mode is obtained when $\omega \sim \omega _s$, while the fast mode corresponds to
$\omega \gg \omega _s$.
C.3. The field-swelling instabilities
C.3.1. Slow mode
As discussed in Basu & Coppi (Reference Basu and Coppi1984), it follows from (C 10) that the slow mode becomes unstable when

As $\varTheta _e$ is increased from 1, the phase velocity decreases and becomes zero at
$1+1/\beta _{\perp e}$. As
$\varTheta _e$ is further increased, a pair of purely imaginary complex conjugate roots appears, leading to the so-called slow-mode swelling instability.
C.3.2. Fast mode
According to the theory of Basu & Coppi (Reference Basu and Coppi1984) that assumes $\zeta _e$ very small, the fast mode becomes unstable when

The validity conditions require in particular that $(k_z/k_{\perp })^2>(m_e/m_i)(1/\beta _{\| e})$ in order to ensure that
$\zeta _e$ is small enough.
C.4. Instability growth rate
Let us now consider the full dispersion relation, given by (C 7), with a general electron response function. The instabilities are illustrated in figure 4 which displays, as a function of the perpendicular to parallel electron temperature anisotropy $\varTheta _e$, the imaginary part of
$r=\omega /kc_A$ for the unstable mode, solution of (C 7) with
$\beta _{\perp e}=1$ and cold ion temperatures, for propagation angle
$\alpha = 89^{\circ }$ (panel a) and
$\alpha = 80^{\circ }$ (panel b), when using four different approximations for the plasma response functions. The first one, referred to as BC84 (black solid line), uses
$R(\zeta _i)= -1/(2\zeta _i^2)$ and
$R(\zeta _e)=1$ and corresponds to the asymptotic solution of Basu & Coppi (Reference Basu and Coppi1984), the second one, denoted DK (red solid line), uses for the electrons the Padé approximant
$R_{20}(\zeta _e)=1/(1-2\zeta _e^2-i{\rm \pi} ^{1/2}\zeta _e)$ (the use of
$R_{42}$ leads to almost identical results) and for the ions the same large-
$\zeta _i$ limit as in the BC84 model, since the ions are cold. The third model, called DKNL (green solid line), differs from the previous one by the fact that electron Landau damping is suppressed (using
$R_{21}(\zeta _e)=1/(1-2\zeta _e^2$)). Superimposed diamond symbols refer to the prediction of GF4, taken in the large-scale limit. It is legitimate to consider this model, derived in the limit
$\zeta _i\to 0$, in the regime of cold ions, because in this case the dynamics is insensitive to the closure assumption. Figure 4(c) corresponds to the case with
$\tau _{\perp _i}=1$ for
$\alpha =89^{\circ }$. It shows that the GF4 model is still approximately valid even in a situation where
$\zeta _i$ is not small.

Figure 4. Imaginary part of $r=\omega /kc_A$ for the unstable mode, solution of (C 7), for propagation angles
$\alpha =89^{\circ }$ (a) and
$80\,^{\circ }$ (b) as a function of the perpendicular to parallel electron temperature anisotropy
$\varTheta _e$ in the case
$\beta _{\perp e}=1$ and cold ion temperatures (leading to take
$R(\zeta _i) = -1/(2\zeta _i^2)$) for various approximations of
$R(\zeta _e)$: BC84 (black solid line) uses
$R(\zeta _e) = 1$, DK (red solid line) uses
$R_{20}(\zeta _e)$ and DKNL (green solid line) uses
${ R_{20}}(\zeta _e)$ (no electron Landau damping). Superimposed diamond symbols refer to predictions of GF4 model taken in the large-scale limit. Panel (c) corresponds to DKNL and GF4 models with
$\alpha =89^{\circ }$ in the case
$\tau _{\perp i}=1$.
Figure 5 displays, as a function of $\varTheta _e$, the positive real part of the three roots of (C 7) (referred to as min, int and max in increasing order of magnitude, displayed in turquoise, magenta and brown colours respectively), in the case of DKNL model (solid lines) or for the DK model (dash-dotted lines), with the predictions of the GF4 model superimposed as diamond symbols, again for
$89^{\circ }$ (panel a) and
$80^{\circ }$ (panel b) propagation angles and cold ions. For the
$89^{\circ }$ angle (which falls outside the range of admissible angles for the BC84 model), the destabilization of the slow mode for
$\varTheta _e>2$ is almost similar for the four models. Within the BC84 model, the fast mode becomes unstable for
$\varTheta _e>4$: the real part of the associated root vanishes (not shown) and the instability growth rate (black solid line in figure 4a) increases rapidly. The slow mode reappears when
$\varTheta _e>4$. The DKNL model displays a very different behaviour, whereby the fast mode (brown solid line in figure 5a) remains almost unchanged for the whole range of values of
$\varTheta _e$. The slow mode (turquoise solid line in the same panel) disappears for
$\varTheta _e>2$ and reappears on the intermediate branch (magenta solid line in figure 5a) for
$\varTheta _e>4$. The intermediate root for
$\varTheta _e <4$ is one of the many plasma modes that coexist with the usual slow and fast modes of fluid theory; it is here the only extra mode for the present choice of the Padé approximant (
$R_{21}(\zeta _e)$, leading to what we called DKNL)). The instability that continues to exist for
$\varTheta _e >4$ (its growth rate corresponds to the green solid line in figure 4a) is associated with the destabilization of this extra plasma mode (while both the fast and slow modes continue to exist and remain stable). The only small difference observed in the presence of Landau damping is that the intermediate mode becomes purely imaginary in a range of values of
$\varTheta _e$ between 2 and 4 (see figure 5a, dash-dotted line). We also note that the use of the
$R_{42}$ Padé does not change the roots associated with the slow and fast modes, but only those associated with the extra damped plasma modes (not shown). We conclude that in an angular boundary layer close to
$90^{\circ }$, the fast mode is always stable, and the instability that continues to exist for
$\varTheta _e>4$ is of the same nature as the slow-mode swelling instability, i.e. its growth rate tends to zero as
$\alpha$ approaches
$90^{\circ }$.

Figure 5. Positive real part of the three roots (denoted min (turquoise), int (magenta) and max (brown) in increasing order of magnitude for DKNL (solid lines) and DK (dash-dotted lines), together with predictions of GF4 (diamond symbols), versus $\varTheta _e$, for
$\alpha = 89^{\circ }$ (a) and
$\alpha = 80^{\circ }$(b), in the case of cold ions.
For $\alpha =80^{\circ }$, the value of
$\zeta _e$ is sufficiently small for the approximation
$R(\zeta _e)= 1$ to be valid. In this case, the behaviour of the slow-mode swelling instability is similar to the previous case, the main difference affecting the fast mode. Its phase velocity for
$\varTheta _e<4$ now corresponds to the intermediate branch, the one with the largest real part corresponding to the extra plasma mode which, in the presence of Landau damping, is heavily damped (damping rate not shown). For
$\varTheta _e>4$, the real frequency of the fast mode vanishes and this mode becomes unstable, as predicted in Basu & Coppi (Reference Basu and Coppi1984). The slow mode reappears on this intermediate branch with a very good match between the GF4 and DK as well as DKNL predictions.
Complementary information is presented in figure 6 which displays the growth rate of the unstable mode as a function of the angle $\alpha$ for the three models described above for
$\varTheta _e=3$ (panel a) and
$\varTheta _e=4.5$ (panel b). For the case corresponding to the slow-mode swelling instability of Basu & Coppi (Reference Basu and Coppi1984) (
$\varTheta _e=3$), the three models are very similar (the case without Landau damping is actually almost identical to the BC84 model so that both curves are superimposed). The GF4 model gives very similar growth rates for angles between
$80^{\circ }$ and
$90^{\circ }$ but its predictions deviate for smaller angles. For
$\varTheta _e=4.5$ the fast-mode instability as predicted by BC84 displays a growth rate proportional to
$k$ up to angles
$\alpha =90^{\circ }$. This fast-mode instability is only recovered with DKNL (and DK) at oblique angles, the deviation with BC84 starting to be significant for
$\alpha > 75^{\circ }$. For this value of
$\varTheta _e$, the slow mode is always stable in the GF4 model.
(i) Influence of warm ions: if one considers the fast mode at an angle close (but not equal) to
$90^{\circ }$, one can assume warm ions and at the same time
$R(\zeta _i)= 0$. Taking also
$R(\zeta _e)=1$, and using (C 7), one gets the dispersion relation (20) from Pokhotelov & Onishchenko (Reference Pokhotelov and Onishchenko2014), where the authors show that ions are stabilizing.
(ii) The case with ion temperature anisotropy: with finite (but isotropic) ion temperatures, other modes are present, but the one which becomes first unstable when electron temperature anisotropy is increased is still the slow mode. The case where the instability is driven by ion temperature anisotropy (so-called classical mirror instability with isotropic electrons) is in contrast different since the instability originates from the extra mode associated with the finite ion temperature fluctuations and not from the slow mode which continues to exist and to be stable above the mirror threshold. An interesting point is that the mirror mode, usually thought of being non-propagating, originates from one of the damped propagating ‘ion temperature modes’. These modes only become non-propagating for a large enough ion temperature anisotropy.

Figure 6. Growth rate of the unstable mode as a function of the propagation angle $\alpha$ for the BC84 (black), DK (red), DKNL (green) and GF4 (diamond symbols) models in the case
$\varTheta _e=3$ (a) and
$\varTheta _e=4.5$ (b). For
$\varTheta _e=3$, BC84 and DKNL curves are superimposed.
Editor Hong Qin thanks the referees for their advice in evaluating this article.