An introductory guide to fluid models with anisotropic temperatures. Part 1. CGL description and collisionless fluid hierarchy

We present a detailed guide to advanced collisionless fluid models that incorporate kinetic effects into the fluid framework, and that are much closer to the collisionless kinetic description than traditional magnetohydrodynamics. Such fluid models are directly applicable to modelling the turbulent evolution of a vast array of astrophysical plasmas, such as the solar corona and the solar wind, the interstellar medium, as well as accretion disks and galaxy clusters. The text can be viewed as a detailed guide to Landau fluid models and it is divided into two parts. Part 1 is dedicated to fluid models that are obtained by closing the fluid hierarchy with simple (non-Landau fluid) closures. Part 2 is dedicated to Landau fluid closures. Here in Part 1, we discuss the fluid model of Chew–Goldberger–Low (CGL) in great detail, together with fluid models that contain dispersive effects introduced by the Hall term and by the finite Larmor radius corrections to the pressure tensor. We consider dispersive effects introduced by the non-gyrotropic heat flux vectors. We investigate the parallel and oblique firehose instability, and show that the non-gyrotropic heat flux strongly influences the maximum growth rate of these instabilities. Furthermore, we discuss fluid models that contain evolution equations for the gyrotropic heat flux fluctuations and that are closed at the fourth-moment level by prescribing a specific form for the distribution function. For the bi-Maxwellian distribution, such a closure is known as the ‘normal’ closure. We also discuss a fluid closure for the bi-kappa distribution. Finally, by considering one-dimensional Maxwellian fluid closures at higher-order moments, we show that such fluid models are always unstable. The last possible non Landau fluid closure is therefore the ‘normal’ closure, and beyond the fourth-order moment, Landau fluid closures are required.

Fluid models are an extremely important tool in many areas of space physics, astrophysics and laboratory plasmas. Even though many physical systems studied in these fields are almost collisionless, where a proper kinetic description should be used, traditional fluid models with isotropic scalar pressures (temperatures), such as the usual magnetohydrodynamic description (MHD) and multi-fluid models based on this description, were extremely successful in modelling, interpreting or at least offering the first insight into many space physics phenomena. For example, it was indeed the simplified fluid approach that allowed Parker (1958b) to predict the existence of the solar wind, which was surprisingly several decades after the breakthrough discoveries in quantum mechanics and relativity. In recent years, observational studies re-sparked interest in the temperature anisotropy effects, that cannot be studied with the usual MHD fluid descriptions, and we anticipate that the interest will grow even further, once data from the Parker Solar Probe and the future Solar Orbiter missions are analysed.
The correct modelling and understanding of collisionless plasmas in a fluid framework concerns not only systems with plasma temperatures far from an isotropic state. It is sometimes forgotten that while the application of fluid theory to strongly collisional systems is intuitively obvious, the approach to collisionless (or weakly collisional) systems is not. From a linear perspective, the usual MHD description does not converge to the collisionless kinetic description even in the low-frequency long-wavelength limit (even when the kinetic distribution function is considered to be an isotropic Maxwellian), since in the absence of collisions the isotropic equation of state is never correct. Additionally, Landau damping never completely vanishes (even when electrons are hot). The situation is much more complicated from a nonlinear perspective where, on average, the effect of Landau damping might be balanced by processes called stochastic plasma echoes (Meyrand et al. 2019).
As a linear example, the ordering of the phase speeds in MHD is always slow, Alfvén, fast, i.e. v s v A v f . In contrast, in collisionless systems with a sufficiently high plasma beta, the real phase speed of the slow mode can become faster than the Alfvén mode, so that the ordering can become Alfvén, slow, fast, i.e. v A v s v f . Or in another words, the phase speeds of linear eigenmodes that are present in MHD do not hold in the collisionless (or weakly collisional) regime, and this effect exists even if the temperatures are isotropic. From a fluid perspective, the main reason for this discrepancy is that, in magnetized collisionless systems, the pressure fluctuations in the Collisionless fluid models. Part 1 5 directions parallel and perpendicular to the magnetic field lines are not equivalent, and cannot be described with a single scalar pressure equation. The pressure fluctuations have to be described with two different evolution equations for p and p ⊥ , even if the mean pressure values are equal (p (0) = p (0) ⊥ ) and no mean temperature anisotropy exists.
The simplest collisionless fluid description is the CGL fluid model -named after Chew, Goldberger & Low (1956) -and sometimes also referred to as collisionless MHD. The CGL dispersion relation is not equivalent to the MHD dispersion relation (even in the case p (0) = p (0) ⊥ ), since the pressure fluctuations still remain anisotropic and the evolution equations for p and p ⊥ remain different. In another words, in collisionless systems the distribution function is free to evolve from its initial state and to become anisotropic. By 'forcibly' prescribing only one scalar pressure, one effectively prescribes a high-collisionality regime, even if collisions are not prescribed explicitly. Formally, the MHD equations can indeed be derived from the collisionless Vlasov equation, i.e. with no explicit collisional operator. It is therefore often stated, that the MHD description is highly collisional implicitly. Moreover, as discussed for example by Kulsrud (1983), while in the presence of a magnetic field, transverse motions can in some circumstances be described by fluid-type equations, the determination of pressures as well as longitudinal motions a priori requires a kinetic description.
Here we focus on collisionless fluid models. Nevertheless, weak collisions can be incorporated easily, and calculations just yield additional terms on the right-hand sides of the parallel and perpendicular pressure and heat flux equations. For the simple Bhatnagar-Gross-Krook (BGK) collisional operator (Bhatnagar, Gross & Krook 1954), see for example Snyder, Hammett & Dorland (1997). A thorough review of anisotropic fluid models, including the collisional dynamics, was presented by Barakat & Schunk (1982). A very good discussion about collisionality of various astrophysical plasmas, such as the solar wind, interstellar medium, accretion disks and galaxy clusters, can be found for example in Schekochihin et al. (2009). Weakly collisional fluid models also seem to be applicable for modelling the upper solar photosphere and chromosphere, where a curious situation exists and the proton-proton collisional frequency is roughly equal to the proton cyclotron frequency -see for example figure 1 in Khomenko et al. (2014). We note that in this guide we use the definition of a 'collisionless fluid model' as being a fluid model that (i) is derived from the collisionless Vlasov equation with zero right-hand side; and (ii) that has two different pressure equations. Our definition therefore differs from an alternative view of for example Zank et al. (2014) and Zank (2014), where models with a scattering operator that reflects charged particle scattering by electromagnetic fluctuations on the right-hand side of the Vlasov equation are also viewed as collisionless.
Fluid modelling of collisionless plasmas is an extremely attractive subject, and an enormous amount of theoretical and numerical work was done in this field in the past. The manuscript presented here has no intention of being a proper review paper in this field. In our opinion, no satisfactory easy-to-read introductory paper exists about collisionless fluid models, and this manuscript attempts to fill such a spot. Instead of just stating the major results and discussing what was done in the past and by whom, we make a significant effort to present a (hopefully consistent) derivation of basic collisionless fluid models. On one hand, the presented calculations might be considered as too detailed in many places. On the other hand, it is exactly the relatively complicated algebra of collisionless fluid models, that makes the field difficult to enter for new researchers. The primary goal of this paper is to allow 6 P. Hunana and others researchers new to the subject to follow the algebra easily. Instead of spending years, an interested reader should be able to become comfortable with the basics of the subject in a matter of days, or perhaps a couple of weeks. The text is separated to two parts. Part 1 is dedicated to fluid models that are obtained by closing the fluid hierarchy with simple (non-Landau fluid) closures. Part 2 is dedicated to Landau fluid closures.
Here, in Part 1, § 2, we start with the detailed derivation of the pressure tensor equation. The pressure tensor p is then decomposed to its gyrotropic part p g = p bb + p ⊥ (I −bb) (also referred to as p CGL ), and non-gyrotropic part Π, the latter usually called the finite Larmor radius (FLR) corrections to the pressure tensor, or the gyroviscous stress tensor. In the highly collisional decomposition p = pI + Π, see e.g. Braginskii (1958Braginskii ( , 1965, the quantity Π is called the stress tensor. The decomposition procedure yields rigorously exact (even though still not closed) evolution equations for p and p ⊥ , see e.g. Oraevskii, Chodura & Feneberg (1968), Passot & Sulem (2004) and Goswami, Passot & Sulem (2005). Importantly, at the leading order (by neglecting the Π and also the non-gyrotropic heat flux contributions q ng ), the equations of Chew et al. (1956), hereafter referred to as CGL, are recovered. We discuss the paper by Chew et al. (1956) and point out that the paper derived the correct form of the pressure equations with the gyrotropic heat flux contributions included, however, the quantities q n , q s used in that paper are related to the usual q , q ⊥ by relations q n = q − 3q ⊥ and q s = q ⊥ . By further neglecting the heat flux contributions, the pressure equations can be written in conservative form. The resulting pressure equations became known as the CGL equations, and they can be interpreted as the conservation laws for the first and second adiabatic invariants (Kulsrud 1983;Gurnett & Bhattacharjee 2005).
Furthermore, we discuss the general equations of collisionless multi-fluid models. The pressure equations are rigorously exact, even though the system is not closed, since the FLR pressure tensor and the entire heat flux tensor are not specified; for a quick look see (2.89)-(2.92). Rewriting the system to a form where the usual CGL conservation laws are on the left-hand side, and all other contributions on the right-hand side, nicely represents the complicated plasma heating processes that can be encountered, and that are responsible for the breaking of the adiabatic invariants; see (2.100), (2.101). We also discuss the conservation of energy. For the case of periodic boundary conditions, the total conservation of energy has a very illuminating form, that can be found for example in Yang et al. (2017a,b). That formulation is obtained by considering the total 'internal' energy for each particle species E int r = 1 2 Tr p r (where the brackets represent integration over the entire spatial domain), which only reveals the total plasma heating. Here, we split the internal energy into its parallel and perpendicular components E int r = 1 2 p r , E int ⊥r = p ⊥r , and formulate the total conservation of energy with the possibility of anisotropic plasma heating, see (2.113)-(2.116).
The CGL description is analysed in great detail in § 3. We derive the CGL dispersion relation, and discuss properties of the slow, Alfvén and fast modes that are present. We verify many of the classical results of Abraham-Shrauner (1967), here written in a slightly more convenient notation by using the parallel plasma beta β , and the temperature anisotropy ratio a p = T ⊥ /T . Collisionless plasmas cannot reach arbitrarily large values of temperature anisotropy, and the linear CGL eigenmodes can indeed become unstable, with the associated instabilities referred to as the (parallel and oblique) firehose instability and the mirror instability. Similarly to MHD, the simple CGL description does not contain any dispersive effects and is technically 8 P. Hunana and others first-order tensor (here called FLR1), we recover the classical result of Yajima (1966), which is notably different from the one provided by Oraevskii et al. (1968). In the isotropic case, the FLR1 tensor is consistent with the one extracted from the stress tensor of Braginskii (1965), if the collisional terms are 'ignored'. It is noteworthy that a proper collisionless limit cannot be achieved from the stress tensor of Braginskii (1965), because expressions are proportional to τ (and also 1/τ ), where τ is time between two collisions (τ ∼ 1/ν where ν is the usual collisional frequency), so the collisionless limit τ → ∞ does not work.
We consider the Hall-CGL-FLR1 fluid model, and we provide dispersion relation for generally oblique propagation, which can be also found in . For higher-order FLR corrections, we only provide analytic dispersion relations for the parallel propagating whistler and ion-cyclotron modes, as well as the perpendicular fast mode. Nevertheless, we provide linearized, normalized and Fourier transformed equations written in the x-z plane for all the fluid models. To obtain the dispersion relation for an oblique propagation, the reader is encouraged to use analytic software such as Maple or Mathematica, or to solve the system numerically. The second-order FLR corrections (FLR2) are here defined as containing the Hall term and the time derivative ∂Π/∂t. We also consider FLR corrections with the non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ , that are here defined as FLR3. The precision of various FLR corrections can be compared nicely by considering the perpendicular fast mode in the long-wavelength limit. The comparison is especially meaningful, when written in the notation of Del Sarto, Pegoraro & Tenerani (2017), see (5.119)-(5.121). We proceed by showing that the FLR3 corrections (with the second-order non-gyrotropic heat flux vectors) indeed recover the fully kinetic dispersion relation for the fast mode in the long-wavelength (low-frequency) limit, a result reported by Mikhailovskii & Smolyakov (1985). Instead of expanding the pressure tensor equation, one can derive very precise linear FLR corrections from linear kinetic theory, which is not addressed here, and the reader is referred to papers by Passot & Sulem (2007) and Sulem & Passot (2015) and references therein.
In § 6, we investigate the parallel and oblique firehose instability. The FLR and Hall dispersive effects are crucial for the stabilization of the firehose instability at small scales, and a comprehensive discussion can be found in . That paper was essentially extracted from this guide (with many figures that we do not republish here), and an interested reader who wants to focus on the firehose instability can find further information there. Nevertheless, here we briefly investigate improvements that can be made by considering the FLR2 and FLR3 corrections, see figures 7-11. Importantly, we show that the non-gyrotropic heat flux vectors in the FLR3 model, partially reproduce the large 'bump' in the imaginary phase speed (growth rate normalized to the wavenumber), when the plasma is close to the long-wavelength limit 'hard' firehose threshold, see figure 7. The firehose instability in a fluid formalism was also investigated by Wang & Hau (2003, 2010, Schekochihin et al. (2010) and Rosin et al. (2011) and references therein.
In § 7, we derive the heat flux tensor equation. The heat flux tensor is then decomposed into its gyrotropic and non-gyrotropic parts, q = q g + q ng . The procedure yields evolution equations for the gyrotropic parallel and perpendicular heat fluxes, q and q ⊥ , that contain the tensor of the fourth-order moment r. It is emphasized that, if one wants to keep the non-gyrotropic Π contributions in the scalar heat flux equations, one needs to keep the non-gyrotropic contributions of the fourth-order moment r ng as well, since there are several possible cancellations even at the linear level. The non-gyrotropic heat flux q ng can be further decomposed to the non-gyrotropic heat Collisionless fluid models. Part 1 9 flux vectors S ⊥ , S ⊥ ⊥ and tensor σ . The detailed algebra of the non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ , i.e. how to express them through lower-order moments, is presented in appendix D. The first-order expressions are obtained at the nonlinear level and the second-order expressions at the linear level. We do not address how to decompose the tensor σ through lower-order moments. Such calculations require a complicated 'inversion' procedure for a third-rank tensor (b × σ ) S , and an advanced reader is referred to Ramos (2005).
In § 8, we consider the fourth-order moment r which is a tensor of fourth rank, r ijkl . The moment is again decomposed to its gyrotropic and non-gyrotropic parts, r = r g + r ng and the gyrotropic part has three scalar components, r , r ⊥ and r ⊥⊥ . We show that for a bi-Maxwellian distribution function, the gyrotropic components can indeed be evaluated as r = 3p 2 /ρ, r ⊥ = p p ⊥ /ρ and r ⊥⊥ = 2p 2 ⊥ /ρ. This constitutes a 'normal' closure, a name suggested by Chust & Belmont (2006). By using a similar procedure to the one provided by Grad (1949) for dilute gases, it is possible to express the nongyrotropic bi-Maxwellian r ng through a combination of p g and Π, see e.g. Oraevskii et al. (1968).
In § 9, a dispersion relation of a fluid model closed with the bi-Maxwellian 'normal' closure is provided for generally oblique propagation, and we call this model second-order CGL (CGL2). We specifically focus on the mirror instability, since this simple fluid model (without any Landau damping) corrects the erroneous 1/6 factor in the 'hard' mirror threshold found in the basic CGL description, a result also reported by Dzhalilov, Kuznetsov & Staude (2011). The mirror instability is not addressed to a higher level of sophistication in this guide. However, we note that capturing the mirror growth rate (when the threshold is crossed) sufficiently well requires Landau fluid models (Snyder et al. 1997) and the stabilization at small scales requires FLR corrections (Passot & Sulem 2007;Sulem & Passot 2015). We also provide the dispersion relation of the Hall-CGL2 fluid model. Finally, the CGL2 model can be simplified by considering slow-dynamics regime and constructing generalized isothermal closure that is called the 'static' closure, yielding the simplest fluid model that captures the correct mirror threshold (Constantinescu 2002;Chust & Belmont 2006;Passot et al. 2006).
In § 10, we consider the bi-kappa distribution function. We show that the closure at the fourth-order moment is constructed by r = α κ 3p 2 /ρ, r ⊥ = α κ p p ⊥ /ρ and r ⊥⊥ = α κ 2p 2 ⊥ /ρ, where the coefficient α κ = (κ − 3/2)/(κ − 5/2), and the closure is valid for κ > 5/2. We call this closure and the associated fluid model 'BiKappa', and we discuss its dispersion relations. Even though the linear modes are generally different in this fluid model than in the CGL2 fluid model, we show that the 'hard' thresholds for the parallel and oblique firehose instability, and for the highly oblique mirror instability, are not affected by and are independent of the κ value. The Hall-BiKappa fluid dispersion relation is also provided. We also provide the first-order non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ . We do not calculate the r ng for the bi-kappa distribution, and therefore we do not provide the second-order non-gyrotropic heat flux vectors.
In § 11, we discuss the core differences between the usual fluid models and kinetic theory. Namely, we discuss why the usual fluid models do not contain collisionless damping mechanisms, such as Landau damping, regardless of the order to which the fluid hierarchy is developed. The effect of Landau damping is present in the collisionless Vlasov equation, and the crucial difference between the usual fluid hierarchy and kinetic calculations just lies in the technique of how the Vlasov equation is integrated over the velocity space. We introduce preliminary ideas as to how the 10 P. Hunana

and others
Landau fluid closures will be constructed. For example, for closures performed at the fourth-order moment instead of the 'normal' closure, one needs to consider perturbations around this state, and prescribe r = 3p 2 /ρ + r , r ⊥ = p p ⊥ /ρ + r ⊥ and r ⊥⊥ = 2p 2 ⊥ /ρ + r ⊥⊥ . The deviations r , r ⊥ , r ⊥⊥ will be calculated from linear kinetic theory in Part 2, by performing Landau fluid closures.
In § 12, we derive the evolution equation for a general nth-order fluid moment (a tensor with 3 n components). We then consider fluid models in the simplified one-dimensional (1-D) geometry that can be viewed as an electrostatic case (or as a propagation along the mean magnetic field), and that are closed at a general nth-order level by a Maxwellian fluid closure. A dispersion relation is obtained, which for n > 4 always yields some solutions that are unstable. It is therefore concluded that the last non-Landau fluid closure is the 'normal' closure and that for n > 4, Landau fluid closures are required. This surprising result, first reported in Hunana et al. (2018), serves as motivation for Part 2, which is a detailed guide to Landau fluid closures.

Pressure tensor equation
Collisionless plasmas are described by the Vlasov equation, which in CGS units reads ∂f r ∂t and which describes how a distribution function f r (x, v, t) evolves in time. The r is the index of species and r = p for protons, r = e for electrons, etc. The q r is the particle charge, m r the particle mass, c the speed of light, E the electric field vector and B the magnetic field vector. The species index r can sometimes be confusing in lengthy tensor calculations with multiple indices and for clarity we will often drop it and reintroduce it when required. To derive the fluid equations, we need to integrate (perform an averaging) at each spatial point over the 'kinetic' velocity v. It is important to realize that the distribution function just describes the probability of finding a particle with velocity v at the position x, t and that the 'kinetic' velocity v entering the distribution function f (x, v, t) is a completely independent variable from x, t, i.e. ∂v i ∂t = 0; ∂v i ∂x j = 0.
(2.2) Also, the magnetic and electric fields B(x, t), E(x, t) are macroscopic quantities that do not depend on v and can be moved outside of velocity integrals over d 3 v, or in another words ∂B i /∂v j = 0 and ∂E i /∂v j = 0. The definitions of the fluid moments are Collisionless fluid models. Part 1 11 where we have omitted the tensor product notation that is sometimes written down explicitly as uu = u ⊗ u and in the index notation (uu) ij = u i u j . If two vectors or tensors are next to each other without an operator between them, a tensor product is always assumed. The number density n is a scalar, the fluid velocity u is a vector, the pressure tensor p is a tensor of second rank (3 × 3 matrix), the heat flux tensor q is a tensor of third rank (with 3 × 3 × 3 components), the tensor r is of fourth rank (with 3 × 3 × 3 × 3 components), etc. Directly from the definitions, it is obvious that all fluid tensors must be symmetric in all of their indices, i.e. p ij = p ji , q ijk = q ikj = · · · . The fluctuating velocity is defined as and should not be confused with the speed of light c.
The second important concept that is used in deriving the fluid hierarchy, is the use of the usual Gauss-Ostrogradsky (divergence) theorem. The divergence theorem is used in velocity space and, written in a form that is typically encountered when calculating the fluid hierarchy, it reads (2.9) The A is a general nth-order tensor, f is a distribution function, the left-hand side is a 3-D integral calculated over the entire velocity volume V and the right-hand side is a surface integral calculated over a boundary of that volume dS =n dS, wheren is a unit normal vector to the local surface area pointing outwards. When such an integral is encountered in the fluid hierarchy, the result is always assumed to be zero. The volume integrals are from v = −∞ to v = ∞ in each velocity component and the integration on the right-hand side of (2.9) is therefore performed over the velocity surface area at infinity. The necessary (but technically not sufficient) condition for the integral to be negligible is lim v→∞ f (v) = 0. (2.10) Since the area dS ∼ v 2 dv, the sufficient condition for the integral to vanish, can be estimated more precisely as When calculating the fluid hierarchy, the encountered expressions are always A(v) ∼ v n , where n is a positive integer. Then, for example, for a Maxwellian distribution f (v) ∼ e −v 2 the limit (2.11) is always zero for all n. Even for a slower converging distribution f (v) ∼ e −|v| , the limit is always zero. For distribution functions proportional to a power law, for example f (v) ∼ (v 2 ) −(κ+1) as for a kappa distribution, the situation is more restrictive, with some minimum required values of κ. When we will calculate the nth-order fluid moment ( § 12), we will see that two surface integrals are encountered, one with A(v) ∼ v n and one with A(v) ∼ v n+1 . Therefore, for a kappa distribution, the strict condition for lim v→∞ (v 2 ) −(κ+1) v n+1 v 2 = 0 yields the requirement κ > (n + 1)/2. For example, many fluid models discussed in this guide are closed at the fourth-order moment r, so n = 4, which implies the requirement κ > 5/2. At first sight, the required limit (2.11) can be considered only a technical 12 P. Hunana and others mathematical detail since, physically, one can argue that particles with enormously large energies will not be measured/observed, and some physical mechanism that is responsible for a cutoff of the distribution function at finite energies can be usually assumed. Or in another words, even observational studies that fit the data with a kappa distribution do not assume that the fit is valid all the way up to v → ∞, and usually a cutoff is implicitly assumed. However, once we calculate the fourth-order moments for a kappa distribution, we will see that, for example r = α κ 3p 2 /ρ, where α κ = (κ − 3/2)/(κ − 5/2), so the restriction κ > 5/2 returns, and has to be applied regardless. When the limit (2.11) is satisfied, the neglect of the surface integrals (2.9) is therefore based on solid theoretical principles, and it is not an approximate or an ad hoc choice, the surface integrals are really zero. Nevertheless, for complete clarity in the upcoming fluid hierarchy calculations, we will differentiate between expressions that are zero exactly, and between the surface integrals (2.9) that are zero asymptotically, by using = 0 and → 0.
We start directly with the derivation of the pressure tensor equation, since a detailed derivation of the density and momentum equations can be found in many books. To derive the pressure tensor equation, it is possible to multiply the Vlasov equation by mc i c j or mv i v j (another possibility is mc i v j ). Here, we will use the first choice, which is slightly easier to present in detail, because the second choice requires the use of density and momentum equations to cancel several terms. It is useful to derive the following identities (2.14) We write down the derivatives with respect to time ∂/∂t and velocity ∂/∂v i explicitly, but we abbreviate the derivative with respect to spatial coordinates as ∂/∂x i ≡ ∂ i . We will need Collisionless fluid models. Part 1 (2.20) where in the last identity we used the result klm δ lk = 0 since the Levi-Civita tensor ijk is antisymmetric and the Kronecker δ ij is a symmetric tensor. Integrating the first term of the Vlasov equation yields The second term of the Vlasov equation yields To go back and forth between the index notation and the vector notation in a fully consistent matter, it is important to establish some conventions that were not required in simple fluid models. One important convention that is typically used is that the divergence of a tensor is meant to be through its first component, i.e. for a general tensor X ijk...n , the divergence ∇ · X means ∂ i X ijk...n . This convention is the reason why in the above expression we used that q ijk = q kij and also the ordering of u and p inside ∇ · (up) reflects that in the index notation we have ∂ k (u k p ij ).
The third term of the Vlasov equation calculates

P. Hunana and others
Finally, the fourth term of the Vlasov equation calculates In the expression above, one encounters a vector product of a vector with a matrix, B × p. This operator might appear unusual at first, but it is just a generalization of a vector product of two vectors, it is defined as (B × p) ij = ikl B k p lj and the result is a matrix. Similarly, a vector product between a vector and a tensor of third rank is defined as (B × q) ijk = ilm B l q mjk and for a tensor of nth rank (B × X ) i...jk = ilm B l X m...jk (i.e. the second index in is the vector B index and the third index in is the first index of tensor X ). However, a definition of a vector product where a matrix is applied on a vector reads (p × B) ij = jkl p ik B l , and in general p × B = −(B × p T ) T , but since the p is here symmetric, p × B = −(B × p) T . The vector product (cross-product) is addressed further in appendix B. The term (2.24) is sometimes rewritten as Combining all the results together 1 + 2 + 3 + 4 = 0, yields the entire pressure tensor equation as One notices that, for example, the pressure tensor equation of Kulsrud (1983), equation (57), contains a minus sign typo in the last term, which should be written as (2.25). That it is indeed a typo and not a definition of p × B is evident from his subsequent equation (59). It is useful to introduce a symmetric operator that acts on a matrix A according to which yields a more compact form of the pressure tensor equation Collisionless fluid models. Part 1 15 2.1. Pressure tensor decomposition By introducing a unit vector in the direction of the local magnetic fieldb = B/|B| and the gyrofrequency Ω = qB 0 /(mc), the last term in the pressure tensor equation (2.26) can be rewritten as (2.29) At very low frequencies ω Ω (and at very long spatial scales) this term will dominate and must be equal to zero. The simplest possibility that makes this term to be equal to zero is p iso = pI, where p is the scalar pressure and I is the unit matrix. In index notation p iso ij = pδ ij , and we can verify that There is however a much more general solution that makes the term (2.29) equal to zero, which is where the 'g' stands for gyrotropic. For quick analytic calculations it is sometimes easier to use p g = (p − p ⊥ )bb + p ⊥ I since we have to deal withbb only once instead of twice. We need to verify that the term (2.29) indeed disappears, (2.33) It is easy to see that the parallel and perpendicular pressures can be extracted from the gyrotropic pressure tensor matrix by performing double contractions according to p = p g :bb, and p ⊥ = p g : The double contraction (the double dot product, represented by the colon) is a very useful operator that is frequently encountered in higher-order fluid models. For two matrices it is defined as A : B = A ij B ij and yields a scalar. Sometimes, the double contraction is defined as A : B = A ij B ji , which for symmetric matrices is equivalent to the previous definition. As a reminder, the usual matrix product of two matrices is (A · B) ij = A ik B kj and yields a matrix. It is useful to write down the following identities that involve the double contraction,bb :bb = 1; (2.36) I :bb = 1; (2.37) (I −bb) :bb = 0; (2.38) I : I = 3; (2.39) (I −bb) : (I −bb)/2 = 1.

P. Hunana and others
can now be verified easily. Even more revealing is to apply the decomposition (2.35) directly to the definition of the entire pressure tensor (2.5), which yields where we have used the fact that the magnitude of the 3-D velocity vector can be decomposed as |v − u| 2 = |v ⊥ − u ⊥ | 2 + (v − u ) 2 , or equivalently |c| 2 = |c ⊥ | 2 + c 2 . A distribution function is called isotropic when the direction of the velocity vector v does not matter, and the distribution function depends only on magnitude |v|, as for example the Maxwellian distribution e −|v| 2 . A distribution function is called gyrotropic when the direction of the perpendicular velocity vector v ⊥ does not matter, and the function depends only on |v ⊥ |, as for example the bi-Maxwellian distribution e −α v 2 e −α ⊥ |v ⊥ | 2 . In another words, the gyrotropic distribution function is isotropic only in its transverse velocity components. The same vocabulary is used for the fluid moments, and a fluid moment is called gyrotropic when it involves integrals over |c ⊥ | 2i , i = 0, 1, 2, 3. . . . For the pressure tensor, there are only two possibilities, represented by the parallel and perpendicular pressures (2.41), (2.42).
The part of the pressure tensor extracted by the decomposition (2.35) is therefore called the gyrotropic pressure p g , sometimes also called the CGL pressure p CGL . The gyrotropic approximation is sufficient at very long spatial scales, where the Larmor radius (= the gyroradius) of particles gyrating around the magnetic field is small, and the non-gyrotropic contributions become negligible. However, at sufficiently small spatial scales comparable to the gyroradius, the non-gyrotropic pressure contributions become significant and we represent these through a tensor Π, that is called the finite Larmor radius corrections to the (gyrotropic) pressure tensor, or sometimes non-gyrotropic or gyroviscous stress tensor. The entire pressure tensor is thus decomposed according to (2.43) and for clarity we write down the decomposition explicitly in matrix notation In the momentum equation, the pressure tensor enters through its divergence, which is useful to break down into its components. Since ∇ · (bb) =b(∇ ·b) + (b · ∇)b, and ∇ · (p ⊥ I) = ∇p ⊥ , the divergence of the pressure tensor calculates as By further using the first requirement, the second requirement can be rewritten as By using the pressure decomposition (2.43) in the pressure tensor equation (2.26), and using the fact that the gyrotropic part of pressure satisfiesb × p g + (b × p g ) T = 0, the pressure tensor equation can be rewritten as We can now derive the time dependent equations for parallel and perpendicular pressures by applying the usual double contractions to this pressure tensor equation.

Parallel pressure equation
We calculate the double contraction withbb term by term. We will need the following identities ∂ ∂t (bb) :bb = 0; ∂ k (bb) :bb = 0. (2.51) The first term calculates as

P. Hunana and others
The second term calculates as (2.53) The third term calculates as The fourth term calculates similarly as (2.56) and the final fifth term becomes The entire equation for the parallel pressure therefore reads Using the symmetric operator (2.27) and the definition for the convective derivative d/dt ≡ ∂/∂t + u · ∇, the equation is written in a simple form as which corresponds to (26) of Passot & Sulem (2004), equation (9) of Goswami et al. (2005) (see also Oraevskii et al. (1968), Passot & Sulem (2007) and Passot, Sulem & Hunana (2012)).

Perpendicular pressure equation
It is possible to either apply the double contraction : (I −bb)/2 to the pressure tensor equation to directly obtain the equation for ∂p ⊥ /∂t, or to apply the trace (the double contraction : I) to obtain the equation for ∂p /∂t + 2∂p ⊥ /∂t, and then subtract the previously obtained equation for ∂p /∂t. Both approaches are equivalent, since p ⊥ = (Tr(p) − p )/2. In another words, directly calculating the trace of the pressure tensor yields where we have used that Tr Π = Π ii = 0. Applying the trace operator to the pressure tensor equation (2.50) term by term yields for the first term for the second term the fourth term is equivalent to the third term even though we kept the last expression in transpose form so that we can use the symmetric operator after we combine all the terms. Finally the fifth term Combining all the terms, using the parallel pressure equation (2.59) and dividing by 2 yields the time-dependent equation for the perpendicular pressure that reads The result corresponds to, for example, equation (25) of Passot & Sulem (2004), equation (8) of Goswami et al. (2005). It is important to note that the parallel 20 P. Hunana and others and perpendicular pressure equations are exact equations. These equations are valid regardless of what kind of higher-order closure will be adopted later. It is useful to rewrite the term Π : (d/dt)(bb) slightly so that we can directly use the induction equation and eliminate the time derivative. The term can be rewritten as (2.69) We will need the following identity for the time derivative of the unit vector (see later in the text), which, when used in the above expression together withb · Π ·b = 0, allows us to write which is an expression of general validity since no specific form of the induction equation was assumed yet. In the next section we will calculate the decomposition of the heat flux tensor q and we will show that, if only the gyrotropic heat flux components q , q ⊥ are considered, the heat flux terms entering the pressure equations readb (2.73) Therefore, by splitting the heat flux into gyrotropic and non-gyrotropic parts, q = q g + q ng , the pressure equations read Terms with the non-gyrotropic heat flux q ng can be further split into non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ and the heat flux tensor σ , which is addressed in appendix D, see (D 124), (D 126). Terms containing the FLR pressure tensor Π are usually called the FLR stress forces, and these forces can generate complicated plasma heating processes. The FLR stress forces can generate both parallel and perpendicular plasma heating that can be determined only from fully nonlinear numerical simulations, since all the terms disappear at the linear level. This process is also sometimes called 'stochastic heating'. The importance of stochastic heating was shown by , who performed 1.5-dimensional numerical simulations with the sophisticated FLR Landau fluid model of Passot & Sulem (2007), Passot et al. (2012).
2.4. On the paper by Chew et al. (1956) The very influential paper by Chew et al. (1956) that is typically cited for the CGL equations (with zero heat fluxes) actually also derived the pressure equations when the gyrotropic heat fluxes are considered. It is important to emphasize that, even though the scalar pressures p n , p s in the notation of that paper are equal to the usual pressures p , p ⊥ , this is not the case for the heat fluxes. The gyrotropic heat flux tensor in that paper is decomposed into components q n , q s according to (Chew et al. 1956, equation (34)) q g ijk = q nbibjbk + q s (δ ijbk + δ jkbi + δ ikbj ). (2.76) In contrast, the usual gyrotropic heat flux decomposition to q , q ⊥ reads (see later in the text) q i.e. the last term is not present in their decomposition. However, this does not mean that their pressure equations (Chew et al. 1956, equations (31), (32)) are incorrect. It only implies that their equations are written in terms of q n = q − 3q ⊥ ; q s = q ⊥ . (2.78) By neglecting the FLR corrections Π and q ng , the parallel and perpendicular pressure equations (2.74), (2.75) simplify to dp dt (2.80) By using (2.78), the heat flux contributions in the parallel pressure equation above can be rewritten as which agrees with the parallel pressure equation obtained by Chew et al. (1956), equation (31). The pressure equations of Chew et al. (1956) are therefore correct, but it has to be remembered that q n is not the usual parallel heat flux q = m (v − u ) 3 f d 3 v. 1 Actually, the authors do not even call the quantities q n , q s heat flux components, but components of a 'pressure-transport' tensor.
To conclude this subsection, we will go slightly ahead and use (3.51), (3.52) (without the zero right-hand sides) that are the usual CGL pressure equations and 22 P. Hunana and others valid only at long spatial scales when the Hall term in the induction equation is neglected. Under this assumption, the pressure equations (2.79), (2.80) are rewritten as d dt It is often forgotten that the paper by Chew et al. (1956)  /ρ 3 = const. and p ⊥ /(ρ|B|) = const., see e.g. Mjølhus (2009). Additionally, similarly to the definition of (appropriately normalized) entropy in MHD, s = ln(p/ρ γ ), it is useful to define parallel and perpendicular entropy in the CGL model, according to (2.86) The total entropy, s + s ⊥ = ln(p 1/3 p 2/3 ⊥ ρ −5/3 ), is equivalent to the MHD entropy when p = p ⊥ , see for example (9)-(11) in Abraham-Shrauner (1967).

Physical meaning of CGL equations
The exact derivation of the CGL equation (2.84), (2.85) will be completed in the next section. Here, we briefly discuss the physical meaning of these equations. The CGL equations (2.85), (2.84) can be interpreted as the conservation of the first and second adiabatic invariants, as is nicely discussed for example by Kulsrud (1983) and by Gurnett & Bhattacharjee (2005). Considering many particles with velocity components v and v ⊥ with respect to the mean magnetic field, the parallel and perpendicular pressures can be estimated as where the brackets represent an average over all particles. To complete the estimate, we need expressions for v 2 and v 2 ⊥ , which come from the conservation of the adiabatic invariants. The first adiabatic invariant is the conservation of magnetic moment µ of a particle that is gyrating periodically around a mean magnetic field. The second adiabatic invariant, sometimes also called the longitudinal invariant, is associated with a particle bouncing periodically between two magnetic mirror points.

23
The conserved quantity is the integral over the parallel momentum of a particle J ≡ b a mv dl, where the magnetic mirrors are located at 'a' and 'b' and the integral is performed along the magnetic field line. The distance between the two magnetic mirrors at 'a' and 'b' can be labelled as L. The conservation of the first and second adiabatic invariants is therefore Since we are considering only non-relativistic particles, the particle mass 'm' is also a constant. The conservation of magnetic moment µ therefore implies v 2 ⊥ ∼ B, and using this in (2.87) yields p ⊥ ∼ Bρ, recovering (2.85). The CGL equation (2.85) therefore corresponds to conservation of magnetic moment, i.e. the first adiabatic invariant. The conservation of the second adiabatic invariant J is more tricky, since we need to somehow estimate the non-intuitive length L between the two magnetic mirrors, and no magnetic mirrors are explicitly assumed, since we are just dealing with somewhat random magnetic field lines. The length L can be nevertheless estimated from two fundamental physical principles. Consider a magnetic flux tube (a deformed cylinder) with cross-sectional area A and length L. The conservation of the total magnetic flux through the area A implies AB = const., which yields an estimate for the area A ∼ 1/B. The second principle is the conservation of the total mass of particles that are completely trapped in that flux tube and that cannot escape, mnAL = const., yielding L ∼ 1/(ρA). The use of A ∼ 1/B yields the final estimate for the non-intuitive length L ∼ B/ρ. The conservation of the second adiabatic invariant J therefore implies v ∼ 1/L ∼ ρ/B. Using this result in (2.87) implies p ∼ ρ 3 /B 2 , recovering (2.84). The CGL equation (2.84) corresponds to conservation of the second adiabatic invariant J, and the equation is valid for particles that are completely trapped and therefore cannot carry a heat flux. As we discuss later, the consideration of the Hall term, the heat flux and the FLR stress forces (the stochastic heating), leads to the breaking of these two adiabatic invariants. Nevertheless, exact conservation laws can still be derived, and the very simple CGL equations (2.84), (2.85) have to be modified with non-zero right-hand sides.
In Part 2 of our guide (see § 4.1), we show that the conservation of magnetic moment is very useful for understanding the form of the perturbed distribution function f (1) in the gyrotropic limit. By performing calculations in the laboratory reference frame, we show that to obtain the correct f (1) in the gyrotropic limit actually requires the usual complicated kinetic integration around the unperturbed orbit (see the Appendix of Part 2) and only then the gyrotropic limit can be imposed. In contrast, performing calculations in the guiding-centre reference frame allows one to prescribe the conservation of magnetic moment from the beginning, and obtain the correct f (1) in the gyrotropic limit perhaps more intuitively. Conservation of adiabatic invariants is important in many areas of space physics and astrophysics. For example, conservation of adiabatic invariants is used to construct models that describe particles trapped inside of magnetic islands and that are being accelerated during magnetic island contraction and merging (Drake, Swisdak & Fermo 2013;Zank et al. 2014).
2.6. Exact equations of anisotropic multi-fluid models Before we discuss solutions of specific fluid models, it is beneficial to summarize the most general equations of multi-fluid models that were derived with no simplifications 24 P. Hunana and others at all. By reintroducing the species index r, the equations that were directly obtained by integrating the collisionless Vlasov equation read where the convective derivative d r /dt = ∂/∂t + u r · ∇. In the above equations, the pressure tensor for each species was just decomposed into gyrotropic contributions p r , p ⊥r and the rest of the pressure tensor (the FLR pressure tensor Π r ), according to which is a rigorous decomposition not introducing any simplifications. The heat flux q r contributions are here left at the most general level without any simplifications either. The equations are accompanied by Maxwell's equations, that for now we write down here with no simplifications to emphasize that no Maxwell equations were used in deriving the above system, (2.95) Equations (2.89)-(2.95) represent an exact (even though not closed) kinetic Vlasov-Maxwell system formulated in fluid variables, and the description is valid for any general distribution function and for all considered spatial and temporal scales. All advanced collisionless multi-fluid models have to be based in one way or another on these equations. One can concentrate on the proton dynamics or one can concentrate on the electron dynamics and simplify these equations accordingly for these spatial scales. Naturally, one can consider complicated multi-fluid descriptions composed of many particle species. For a new reader who just jumped straight to this section, the symmetric operator 'S' acts on a matrix A ij according to A S ij = A ij + A ji . If one does not like the symmetric operator in the expressions above, the symmetric operator is actually not necessary, since where we also got rid of the double contractions. Note that Π ik = Π ki and also p ik = p ki . Equations (2.89)-(2.92) contain two very important catches -the entire heat flux tensor q r and the non-gyrotropic pressure tensor Π r are not specified yet. In this moment, the entire class of collisionless fluid models separates to many possible classes and sub-classes, depending on how precisely one wants to evaluate the heat flux tensor q r and the FLR pressure tensor Π r . There are several complications how to correctly model these two quantities, but before we discuss this, we want to point out that by using the simple general identities (3.46), (3.47) (that are nothing more than calculating derivatives and using the density equation for each species separately), the above pressure equations can be rewritten to the following form d r dt Again, no simplifications were introduced and no Maxwell equations were used yet.
The equations are exact. The equations represent the complicated plasma heating processes that are responsible for anisotropic plasma heating, and that are in general very difficult to classify. The left-hand sides of (2.100)-(2.101) are the familiar CGL equations that represent conservations of the first and second adiabatic invariants, and all the expressions on the right-hand sides break these adiabatic invariants. On the right-hand side, we have the ∂B/∂t that couples various species together through Maxwell's equations, and that also introduces the Hall term responsible for the simplest dispersive effects. The heat flux tensor q r (a tensor of third rank), is decomposed into its gyrotropic and non-gyrotropic (FLR) parts q r = q g r + q ng r . At large scales, the gyrotropic heat flux can be viewed as a gateway for the simplest form of collisionless damping mechanism, known as Landau damping, that can be further separated to a 'pure' electrostatic Landau damping, and its magnetic analogue, the transit-time damping -see Part 2 of the text, subsection 'Coulomb force and mirror force (Landau damping and transit time damping)'. Nevertheless, fluid models that contain the gyrotropic heat flux fluctuations, but that do not contain any Landau damping, can be constructed as well (see CGL2 model later in the text). The FLR pressure tensor Π r introduces further dispersive effects and therefore also modifies the collisionless damping rates, as well as the heat flux q r modifies the FLR pressure corrections Π r . The quantities q r and Π r are therefore generally coupled. Importantly, the FLR pressure tensor Π r introduces complicated turbulent plasma heating processes, referred to as stochastic heating. Since these FLR stress forces completely disappear at the linear level in the above pressure equations, they can be explored only with fully 26 P. Hunana and others nonlinear numerical simulations (Π r still enters at the linear level in the momentum equation and modifies the dispersion relations).
From a fluid perspective, there are several major difficulties how to correctly model the quantities Π r and q r . The heat flux tensor q r is described by an infinite hierarchy of higher-order fluid moments, and one needs to find an appropriate way, how to close the system. The FLR pressure tensor Π r is described by equations that are implicit, and the FLR corrections have to be typically expanded on temporal and spatial scales in order to obtain tractable expressions that can be used for numerical simulations. This expansion will naturally restrict the area of validity of such fluid models to those scales considered. For example, for the proton-electron plasma where the spatial and temporal scales are largely separated because of the mass ratio m p /m e = 1836, simple first-order FLR corrections expanded around the proton scales kρ i , ω/Ω p will long lose their validity at the electron scales. The most complicated Landau fluid models that use linear kinetic theory to evaluate Π and that contain the Bessel functions (Passot & Sulem 2007;Passot et al. 2012;Sulem & Passot 2015), have technically no restriction for wavenumbers k ⊥ ρ i , but the first-order frequency restriction is there nevertheless.

Conservation of energy
The equations (2.100)-(2.101), that correctly split the entire plasma heating into parallel and perpendicular heating, are indeed quite complicated. To gain further insight, let us briefly consider a case in which we are not interested in parallel and perpendicular heating (whose net effect can possibly be zero), and we are only interested in the heating of the entire system. Such formulations are sometimes used to interpret plasma heating in fully kinetic simulations (Yang et al. 2017a,b). By summing together the pressure equations (2.91), (2.92), we are interested in the evolution equation ∂(p + 2p ⊥ )/∂t, or in another words, we are interested in ∂Tr p/∂t. Summing the equations together (actually the trace of the entire pressure tensor equation was calculated a few pages back, that is how the p ⊥ equation was derived) yields ∂ ∂t (p r + 2p ⊥r ) + ∇ · ((p r + 2p ⊥r )u r ) + ∇ · (Tr q r ) + 2((p g r + Π r =p r ) · ∇) · u r = 0. (2.102) Note that Tr ∇ · q r = ∇ · (Tr q r ), and also Tr Π r = 0. By integrating over the entire spatial volume of the plasma, one can define the internal energy (or thermal energy) for each particle species as By considering a special case of periodic boundary conditions, i.e. typical numerical simulations of turbulence in a periodic box, one can use the Gauss-Ostrogradsky theorem (2.9), now in a spatial domain, and for periodic boundary conditions all divergence operators in (2.102) vanish, yielding (2.104) Equations (2.102), (2.104) are equivalent to (4) and (7) of Yang et al. (2017b). Here the equations are just written in such a way that one can explicitly see the pressure components p , p ⊥ , Π, but the equations are equivalent. Following Yang et al. (2017b), one can also easily derive the evolution equations for the 'kinetic' energy and the electromagnetic energy For the kinetic energy, by first applying u r · to the momentum equation (2.90) u r · ρ r ∂u r ∂t + ρ r u r · ∇u r = −u r · (∇ · p r ) + q r n r u r =j r · E, (2.106) which further implies and for the electromagnetic energy 1 8π (2.108) Therefore, by integrating over the entire volume and assuming periodic boundary conditions, the total conservation of energy can be expressed as (Yang et al. 2017b) which beautifully clarifies how the energy can be transferred. Nevertheless, we point out that it is exactly the splitting into the parallel and perpendicular heating that is so complicated, since even with periodic boundary conditions, only a few terms in the pressure equations (2.91), (2.92) vanish. Considering anisotropic heating, and splitting the internal (thermal) energy into parallel and perpendicular parts where ( p r · ∇) · u r = p rb · ∇u r ·b − p ⊥rb · ∇u r ·b + p ⊥r ∇ · u r + (Π r · ∇) · u r . (2.117) The anisotropic plasma heating is obviously a very complicated process.

CGL description
The CGL model is the simplest possible fluid model that incorporates anisotropic temperatures. In contrast to MHD, CGL contains two pressure equations. For very good MHD reviews, from the perspective of turbulence in the solar wind, and the use of MHD in modelling heliospheric and astrophysical plasmas, we recommend the reviews by Goldstein, Roberts & Matthaeus (1995), Tu & Marsch (1995), Zank (1999), Zhou, Matthaeus & Dmitruk (2004) and Bruno & Carbone (2013). The CGL model is obtained by using the parallel and perpendicular pressure equations (2.74), (2.75) and prescribing zero heat flux q = q ⊥ = 0, q ng = 0, and by neglecting the FLR pressure tensor Π = 0. The pressure equations greatly simplify to the form (3.2) By using the notation for the convective derivative d/dt = ∂/∂t + u · ∇, the pressure equations can be rewritten as dp dt Note that the CGL pressure equations are fully nonlinear and that the nonlinearities are actually of fourth order, i.e. CGL pressure equations are more nonlinear than the usual MHD pressure equation dp/dt + γ p∇ · u = 0, where the nonlinearity is only Collisionless fluid models. Part 1 29 of second order. At first sight, expressionsb · ∇u ·b might appear to be a little unusual when written without any brackets. From MHD, one is familiar with nonlinear expressions of the type u · ∇u, which are also usually written without any brackets and can be interpreted in two equivalent forms, either as (u · ∇)u (a scalar u · ∇ applied on a vector u; or u · (∇u) (a vector u multiplied by a matrix ∇u). In a similar fashion, the expressionb · ∇u ·b is meant to be interpreted as (b · ∇u) ·b orb · (∇u) ·b; but importantly, it is not meant to be interpreted asb · ∇(u ·b), because the derivative is meant to act only on u in this case. For direct numerical simulations, the fourth-order nonlinearity explicitly readŝ +b y (b x ∂ y u x +b y ∂ y u y +b z ∂ y u z ) +b z (b x ∂ z u x +b y ∂ z u y +b z ∂ z u z ), (3.5) and is of course a scalar. So far, we have not made any assumptions about particle species, and have only integrated the Vlasov equation and derived the density, momentum and scalar pressure equations, which can be done for each species separately. By re-introducing the species index r, where r = p for protons, r = e for electrons etc., the n-fluid CGL-type equations become ∂p r ∂t + ∇ · (u r p r ) + 2p rb · ∇u r ·b = 0; (3.8) ∂p ⊥r ∂t + ∇ · (u r p ⊥r ) + p ⊥r ∇ · u r − p ⊥rb · ∇u r ·b = 0. (3.9) The above CGL pressure equations are valid regardless of the form of the induction equation ∂B/∂t and are therefore very general and valid for a wide range of CGL-type n-fluid models. It is advisable to keep them in this form when considering direct numerical simulations or solving linear dispersion relations. The above equations are accompanied by Maxwell's equations. Specifically, (i) Gauss's law ∇ · E = 4πρ c where the total charge density ρ c = r q r n r , (ii) ∇ · B = 0, (iii) Faraday's induction equation where the total current j = r q r n r u r . In many areas of space physics and astrophysics it is often very useful to eliminate very high frequency effects occurring at frequencies higher than the plasma frequency. This is achieved by prescribing local charge neutrality (so called quasi-neutrality since the plasma is still ionized) by ρ c = 0 and by neglecting the second term in Ampère's law (1/c)(∂E/∂t) known as Maxwell's displacement current (or Maxwell's correction to Ampère's law). This eliminates high-frequency effects such as for example plasma oscillations, known as Langmuir waves. This assumption is also equivalent to assuming that the Alfvén speed is much less than the speed of light V A /c 1. It is important to emphasize that Gauss's law ∇ · E = 4πρ c is replaced by the quasi-neutrality condition ρ c = 0 and no requirement is imposed on ∇ · E (other than that it is small). For further discussion concerning the charge neutrality and ∇ · E, see for example Braginskii (1965) page 264, Webb et al. (2007) equation (7) and Passot & Sulem (2007) (3.12) and this form of Maxwell's equations is considered in most of the fluid models that we discuss, regardless of future complexities arising from FLR corrections, heat flux, Landau damping and so on. The exception is the section in Part 2 'Electron Landau damping of the Langmuir mode', where the displacement current has to be retained.
Here, we further focus only on a plasma consisting of protons and electrons, with charges q p = e and q e = −e, but models with more species can of course be considered. Charge neutrality implies that the proton and electron number densities are equal, i.e. n p = n e = n. By using the definition of current j = enu p − enu e , the electron velocities are related to the proton velocities by which after using (3.12) yields (3.14) By using the momentum equation (3.7) for electrons, the electric field can be expressed as The last term represents the effects of electron inertia and is for example responsible for the effect that at small spatial scales the frequency ω of the parallel whistler mode converges to the electron cyclotron frequency Ω e instead of increasing to infinity. When considering only relatively low frequencies, it is very useful to neglect the electron inertia term. 2 This yields somewhat simpler analytic expressions, but most importantly, it provides a great benefit for direct numerical simulations since the time step is not restricted by the requirement to fully resolve the electron motion. This of course restricts the validity of a model to frequencies that are sufficiently smaller than the electron cyclotron frequency. By neglecting the electron inertia and using (3.14), the electric field can be expressed as Using this expression in the proton momentum equation and in the induction equation eliminates the electric field from the entire system. The system of equations therefore reads (3.20) and is accompanied by the parallel and perpendicular pressure equations (3.8), (3.9) for the proton and electron species. In the induction equation (3.20), the first term is the familiar term from the MHD induction equation and considering only this first term yields the classical CGL model (where in addition the electrons are prescribed to be cold to eliminate p e from the momentum equation). The second term in the induction equation represents the Hall term, and fluid models containing this term are generally referred to as Hall-CGL fluid models. The third term represents the electron pressure contributions and in CGL plasmas the term is in general non-zero, since In contrast, by prescribing the electrons to be isotropic and isothermal (a typical assumption used in hybrid simulations) with pressure p e = p e I, p e = nT (0) e , so that ∇ · p e = ∇p e = T (0) e ∇n, eliminates the pressure term since ∇ × [(1/n)∇n] = ∇ × [∇ ln(n)] = 0. By prescribing the electrons to be isotropic but not isothermal, with pressure p e = p e I, the term initially does not disappear and it is equal to (3.22) The term only disappears after the equation of state is prescribed, for example for p e ∼ n γ . The term (3.22) is sometimes called the 'battery' term (Biermann 1950;Kulsrud et al. 1997;Khomenko et al. 2017), since it is argued that for any deviations from an ideal equation of state (caused for example by shocks), the term can produce magnetic field fluctuations even if there is no magnetic field initially.
3.1. Normalized equations and definitions As with MHD and Hall-MHD, it is often useful to work with normalized equations. The density, velocity, magnetic field and pressure are normalized with respect to ρ 0 , u 0 , B 0 , p (0) . The length is normalized with respect to (for now) arbitrary x 0 and the time with respect to t 0 = x 0 /u 0 . The normalized quantities are then

P. Hunana and others
Obviously n = ρ. After dropping the tilde, the normalization yields the same density and pressure equations, whereas the normalized momentum and induction equation become (3.26) where the Alfvén speed The momentum equation implies that it is beneficial to choose the (so far unspecified) normalizing velocity u 0 to be the Alfvén speed u 0 = V A . If the second and third terms are neglected in the induction equation (i.e. when the usual MHD induction equation is used), the CGL system is independent of length scale x 0 , similarly to the MHD system. If the Hall term is considered, it is beneficial to choose the (so far unspecified) normalizing length scale x 0 to be the ion inertial length (3.27) Therefore by choosing x 0 = d i , the normalizations that we use in all fluid models in this paper read The normalization also has the advantage that, when transferred to Fourier space, one obtains the dispersion relations for normalized wavenumber k and frequency ω as It is useful to define the parallel and perpendicular thermal speeds (3.30) and the definition of parallel and perpendicular temperatures are T = p /n and T ⊥ = p ⊥ /n. We use the usual convention with the Boltzmann constant k B = 1. 3 The parallel and perpendicular plasma beta are defined according to and later we will often use the abbreviation a p for the proton temperature anisotropy ratio Collisionless fluid models. Part 1

33
The normalized momentum and induction equations therefore read (tilde are dropped) (3.34) Another useful quantity that we will use later in the text is the proton Larmor radius (the gyroradius) ρ i , defined according to We note that normalizations do not have to be done with respect to mean values, and normalizations with respect to fluctuating (turbulent) quantities are used for example in the nearly incompressible models of Zank et al. (2017) and Zank & Matthaeus (1993).

Classical CGL model with cold electrons
Here we want to consider the simplest possible CGL model for the proton species. We assume the electrons to be massless and cold (p e = 0), and we also neglect the Hall term in the induction equation. Since only proton species are present, for simplicity we can drop the proton index 'p'. Let us work in physical units for a moment. The classical CGL model with cold electrons therefore reads (written in physical units)  (3.42) and the result is of course a scalar. Similarly for the spatial and convective derivative (3.43) Later we will need (3.44) Now we can directly calculate d dt (3.45) and, by using the density equation dρ/dt = −ρ∇ · u and the identity (3.43), we obtain d dt (3.46) The above equation is completely general since we did not use the induction equation so far and we will use this equation in the next section when we consider the Hall-CGL model. In a similar way one derives a completely general identity d dt (3.47) Now, we use the induction equation Equations (3.51)-(3.52) verify that the CGL pressure equations (3.41) written in the 'conservative' form and the directly derived CGL equations (3.38)-(3.39) are equivalent. We will see that the right-hand side of the usual CGL equations (3.41) is non-zero if the Hall term is considered, and that it is further modified by the heat flux contributions and by the FLR stress forces.
3.3. Linearized CGL equations To obtain the dispersion relations, the equations need to be linearized and transformed to Fourier space. Equations (3.36)-(3.40) are linearized with respect to mean values where the magnetic field B 0 is assumed to be in the z-direction. We assume that the mean value of the velocity is zero (the zero mean velocity value u (0) = 0 should not be confused with the normalizing velocity as in (3.23), which is later chosen to be V A ). The termŝ b · ∇u ·b are linearized asb · ∇u ·b lin = ∂ z u z and the linearized system reads ∂ρ ∂t The divergence of the pressure tensor is 58) and the first step of linearization yields, by components, which can be conveniently written in matrix notation as (3.62) Notice that, since ∂ i |B| =b · (∂ i B), the derivatives of unit vectors are calculated according to which by components reads where, importantly, ∂ zbz lin = 0. The divergence of the pressure tensor in the liner approximation is therefore expressed as where the second term contributes if the temperature anisotropy is present, i.e. when The entire set of linearized CGL equations reads Additionally, the 3 components of the induction equation are not independent, but satisfy the ∇ · B = 0 constraint. In order to find the linear dispersion relations, without any loss of generality, we consider wave propagation in the x-z plane and assume that there is no variation with respect to the y-coordinate (i.e. one assumes that all expressions containing ∂ y are zero; alternatively, one can consider propagation in the y-z plane, where all expressions containing ∂ x would be zero). By exploring the above system, it is noteworthy that the density equation does not play any role in determining the dispersion relation, since no other equation contains the variable ρ.
3.4. Alfvén mode On considering (3.68)-(3.74) written in the x-z plane with ∂ y = 0, it is apparent that the equations for u y and B y decouple from the entire system and applying ∂/∂t to the equation for u y yields which is a wave equation describing the propagation of (generally oblique) Alfvén waves in the CGL description. The same wave equation can be obtained for the component B y . For isotropic mean pressures p (0) = p (0) ⊥ the equation is equivalent to the MHD description with the usual Alfvén speed where a p is the temperature ratio (3.32), the Alfvén wave equation can be rewritten as Considering a wave propagating obliquely in the x-z plane with the wavevector k = (k ⊥ , 0, k ) at an angle θ with respect to the z-direction (the B 0 direction), the parallel and perpendicular wavenumbers are defined as k = k cos θ; k ⊥ = k sin θ , (3.78) and the transformation to Fourier space is performed according to (see appendix A) The wave equation (3.77) is then which yields the dispersion relation for the Alfvén waves It is important to emphasize that, similarly to MHD, the Alfvén mode propagates in all the oblique directions, with a frequency ω that is proportional to the projection of the wavenumber k in the direction of B 0 , i.e. k = k cos θ. If the expression under the square root in (3.81) becomes negative, i.e. if 1 + β (a p − 1)/2 < 0, the frequency ω is imaginary and the Alfvén mode becomes unstable. This instability, that affects the obliquely propagating Alfvén mode, is known as the oblique firehose instability. The necessary (but not sufficient) condition for the instability to develop 38 P. Hunana and others is a p = T (0) ⊥ /T (0) < 1. The firehose instability criterion can be rewritten in many forms, the most common being The first expression directly implies that, if β 2, no firehose instability can exist. The firehose instability therefore exists only for a relatively high plasma beta β > 2 and only if the parallel temperature is higher than the perpendicular temperature. Importantly, the firehose instability criterion in the CGL model is equivalent to the one found from linear kinetic theory in the long-wavelength limit (for example, Rosenbluth (1956), Chandrasekhar, Kaufman & Watson (1958, Parker (1958a) and Vedenov & Sagdeev (1958)). Alfvén wave propagation only affects components u y and B y which are decoupled from the rest of the system, and yields the other (eigenvector) components u x , u z , B x , B z as zero. The u y component does not enter the density equation and Alfvén mode fluctuations therefore do not produce any density fluctuations or pressure fluctuations. The Alfvén mode is therefore incompressible. It is also useful to define the magnetic compressibility χ b (k , k ⊥ ) in Fourier space, that measures the relative ratio of magnetic energy in the parallel component versus the total magnetic energy (3.83) Since the Alfvén mode produces only B y fluctuations, with B x = B z = 0, the magnetic compressibility for this mode is for all the propagation directions. Sometimes the magnetic compressibility is defined as a ratio of parallel and perpendicular energies |B z | 2 /(|B x | 2 + |B y | 2 ), this quantity however has a slight disadvantage in that the values can become very large (and actually equal to infinity), which makes plotting of this quantity problematic. In contrast, the magnetic compressibility as defined in (3.83), has a nicely bounded range of values between 0 and 1. A similar quantity is defined for the ratio of the energy in velocity fluctuations (3.85) and the Alfvén mode has χ u = 0 for all propagation directions, since only the perpendicular velocity component u y is non-zero.

Slow and fast modes
By applying ∂/∂t to the equations for ∂u x /∂t and ∂u z /∂t and using the evolution equations for B x , B z , p , p ⊥ one obtains The easiest way to solve this coupled system is to transform to Fourier space, which yields (3.88) The system has non-trivial solutions only if the determinant of the matrix is zero, yielding the CGL dispersion relation at fourth order in frequency ω that contains the slow and fast modes One can proceed in several ways with the above dispersion relation.
In the specific cases of strictly parallel (k ⊥ = 0) and strictly perpendicular (k = 0) propagation directions it is actually easier to work directly with the system (3.88), because the off-diagonal components in the matrix are zero and the dynamics in the u x and u z components decouples. For strictly parallel propagation, one solution is in the u x component, with the dispersion relation identical to the Alfvén wave (3.81). Since the magnetic and velocity fluctuations are only in the u x , b x components, the compressibility χ b = 0 and χ u = 0 for this mode. The second solution for parallel propagation is in the u z component, and it is called the sound mode, or the ion-acoustic mode, with dispersion relation ω = ±k 3p (0) /ρ 0 = ±k V A 3β /2. The parallel ion-acoustic mode has χ u = 1 and χ b = 0. Considering strictly perpendicular propagation, the waves in the u z component vanish with the solution ω = 0. The waves in the u x component are fast magnetosonic waves with dispersion relation , and χ u = 0 and χ b = 1. It is also possible to define the parallel and perpendicular sound speeds C = 3p (0) /ρ 0 ; ⊥ /ρ 0 , which yields the parallel acoustic mode dispersion ω = ±k C and the perpendicular fast mode dispersion ω = ±k ⊥ V 2 A + C 2 ⊥ . 4 4 Sometimes the parallel and perpendicular sounds speeds are defined as C = p (0) /ρ 0 ; C ⊥ = p ⊥ /ρ 0 , which yields the parallel acoustic mode dispersion ω = ±k √ 3C and the perpendicular fast mode dispersion

P. Hunana and others
In the general case of oblique propagation, we choose to use the normalized frequencies and wavenumbers (3.29), and rewrite the dispersion relation as It is sometimes useful to introduce the parallel Alfvén phase speed v A , and its since v 2 A can be used to write the CGL dispersion relation in a shorter form, as done later in (3.193). The solution of the polynomial (3.90) is simply ω 2 = (A 2 ± A 2 2 − 4A 0 )/2, and it is obvious that the coefficient A 2 is always positive, since both the plasma beta and the temperature anisotropy ratio must always be positive. It is also possible to show that the discriminant under the square root A 2 2 − 4A 0 0 (proof shown a few lines below). The solution with the minus sign is called the slow mode and the solution with the plus sign is called the fast mode. For the slow mode, ω 2 becomes negative (and ω complex) when A 0 < 0, implying that the slow mode becomes unstable when where we suppress the tildes on the wavenumbers since the above condition is valid for both k and k forms of wavenumbers. The above condition implies that there are two competing instabilities. The highly parallel propagation with k k ⊥ yields an instability threshold equivalent to (3.82) and the associated instability is called the parallel firehose instability. The highly oblique propagation with k ⊥ k yields the mirror instability and the threshold in the CGL description reads 5 1 + a p β − 1 6 a 2 p β < 0; 1 6 It is well documented that the mirror instability threshold in the CGL description is not equivalent to the one found in linear kinetic theory. In particular, the CGL mirror threshold contains a factor of 6 error (Abraham-Shrauner 1967;Tajiri 1967;Kulsrud 1983;Ferrière & André 2002), since the correct threshold obtained from kinetic theory reads 1 + a p β − a 2 p β < 0; The threshold (3.93) implies that the necessary condition for the mirror instability to exist in the CGL description is a p = T (0) ⊥ /T (0) > 6, whereas in kinetic theory the necessary condition is a p > 1. A very good discussion about the physical mechanism of the mirror instability and the role of resonant particles can be found in Southwood & Kivelson (1993) and Kivelson & Southwood (1996). The correct mirror instability threshold is obtained when using the normal closure discussed later on in the text or, more simply, its quasi-static limit (Constantinescu 2002;Chust & Belmont 2006;Passot et al. 2006). In the opposite case where ions are cold, a similar electron temperature anisotropy instability is obtained which is called the field-swelling instability, see e.g. Basu & Coppi (1984).
Returning to the oblique dispersion relation, we show that indeed A 2 2 − 4A 0 0. By using (3.78) the wavenumber k can be pulled out of the A 0 , A 2 coefficients and the dispersion relation can be rewritten in terms of a phase speed as .97) and (after staring for a sufficiently long time) it can be shown that which is obviously always non-negative. The dispersion relation for the slow and fast waves in the CGL description can be written in the compact form which is equivalent to the dispersion relation of Abraham-Shrauner (1967), equation (30), here written in the convenient units of β and a p = T (0) ⊥ /T (0) (for a direct comparison with the units used in that paper, the relations are β = 2S and a p β = 2S ⊥ ). The magnetic compressibility of the slow and fast modes can be calculated easily directly from the induction equations (3.72). The slow and fast modes generate only components B x and B z with the component B y = 0. The transformation of (3.72) to Fourier space yields −ωB x = B 0 k u x , and −ωB z = −B 0 k ⊥ u x , and by applying the magnitude operator yields the magnetic compressibility For quasi-parallel propagation angles the magnetic compressibility χ b is therefore not a good indicator for mode recognition, since all three modes have χ b close to zero. The magnetic compressibility is, however, an excellent tool for oblique propagation angles, because the Alfvén mode has χ b = 0 regardless of the propagation direction, whereas for the slow and fast modes χ b increases with θ towards χ b = 1. We will leave discussion about the CGL velocity eigenvector (the ratio χ u ) to a later subsection since the discussion is a bit technical, and we instead focus now on an interesting effect that can be directly obtained by considering the CGL dispersion relations. 3.6. Slow mode can become faster than Alfvén mode There is a very interesting phenomenon present in the CGL description that is also observed in the kinetic description but that is absent in the MHD description. For a sufficiently high plasma beta, the slow mode phase speed can be larger than the Alfvén mode phase speed. To obtain the critical plasma beta when this happens, one can write the dispersion relation for the Alfvén wave as ( ω/ k) 2 A = (1 − β /2 + a p β /2) cos 2 θ ≡ A A and using this expression with the dispersion relation (3.95), the critical condition when the slow mode speed matches the Alfvén speed is It takes quite a few algebraic operations to simplify the result, the easiest way is to collect terms with different powers of β and then simplify. A very valuable approach available now is to use mathematical software such as Maple or Mathematica for guidance, and for example the Maple command simplify(expr,trig) immediately simplifies the result, which can then be of course verified with pencil and paper. The result is The expression is naturally satisfied for cases of strictly parallel and strictly perpendicular propagation. For θ = π/2, both the slow mode and the Alfvén mode Collisionless fluid models. Part 1 43 have zero frequency. For θ = 0, the slow mode speed can be only smaller or equal to the speed of the Alfvén mode (never higher), depending on another condition that will be easier to clarify later when specifically considering the parallel propagation limit. Now, on which side of the threshold (3.101) does the speed of the slow mode become faster than the speed of the Alfvén mode? Working with inequalities is a bit trickier. By assuming that both modes have real frequencies (and are not firehose or mirror unstable), the condition is derived to be − sin 2 θ cos 2 θ 1 + β 3 2 a p − 2 + β 2 a p 3 4 a p − 2 0, (3.102) which agrees with the result obtained by Abraham-Shrauner (1967), equation (42). The expression (3.102) reveals that there is a condition that can be satisfied for all directions of propagation and that reads This is a quadratic equation in β which can be solved easily for a given a p (or solved for a p for a given β ). To visualize this equation, the threshold is plotted in figure 1 as a black dashed line. A particularly useful result is obtained for isotropic temperatures (a p = 1). In this case the inequality (3.103) becomes 1 − β /2 − β 2 5/4 0 and the critical plasma beta at which the slow mode speed becomes larger than the Alfvén mode speed is Thus, this effect exists in the CGL description even when the temperatures are isotropic.
3.7. CGL parallel propagation When numerically solving the CGL dispersion relations for finite (non-zero) propagation angles, the slow and fast mode expressions always get the correct ordering with ω s ω f . It is sometimes confusing about how the dispersion relations hold in the limit of parallel propagation θ → 0 which we will address here. In this limit, the dispersion relations (3.99) for the slow and fast modes can be directly evaluated as (3.105) Now, because √ a 2 = |a|, for further calculation one needs to determine the sign of the expression under the square root. For 1 + 1 2 a p β − 2β > 0, i.e. when β < 2/(4 − a p ), the above equations yield so the slow mode is the sound/acoustic mode and the fast mode coincides with the Alfvén mode. However, for 1 + 1 2 a p β − 2β < 0, i.e. when β > 2/(4 − a p ), the dispersion relations yield

P. Hunana and others
so the slow mode coincides with the Alfvén mode and the fast mode is the sound/acoustic mode. In this region, the slow mode phase speed will always stay equal to the Alfvén mode phase speed, regardless of further increases of β . In the special case when 1 + 1 2 a p β − 2β = 0, which is equivalent to all three modes propagate with the same phase speed ( ω/ k) 2 Asf = 3 2 β . The criterion (3.108) can be satisfied only for β 1/2 and in the limit β → ∞ the asymptote is a p = 4, so the possible range of a p is [0, 4]. This curve is plotted in figure 1 as a green dotted line. The curve separates two regions. In the first region β < 2/(4 − a p ) and v s < v A . In the second region β > 2/(4 − a p ) and v s = v A . The interesting result to remember is that, for isotropic temperatures (a p = 1) and parallel propagation, the critical plasma beta at which the slow mode speed matches the Alfvén mode speed is β crit (0) 2 3 .
( 3.109) 3.8. CGL with T ⊥ = T is not equivalent to MHD It is important to emphasize that, even for isotropic temperatures a p = 1, the CGL dispersion (3.99) is not equivalent to the MHD dispersion for the slow and fast waves. The MHD dispersion is usually written in the form (see (C 20) which nicely shows that the expression under the square root is always greater or equal to zero. The MHD sound speed is defined as C 2 s = γ p (0) /ρ 0 , where γ = 5/3. Rewritten with the parameters we use here C 2 s /V 2 A = γβ /2 = 5β /6, the MHD dispersion for the slow and fast waves reads In the CGL description with isotropic temperatures (a p = 1) the dispersion relation (3.99) for slow and fast waves reads (3.115) We wrote down all major possibilities on purpose, since we wanted to clearly show that, regardless of the choice of β and γ , the MHD dispersion relation cannot match the CGL dispersion relation simultaneously for all propagation directions θ . One could even play with an abstract idea to make the MHD γ dependent on the angle θ , and possibly β . Even then, there is no obvious way to modify the MHD dispersion relation so that it behaves like the CGL model. Additionally, the CGL description leads to the development of pressure anisotropy even when the initial pressure is isotropic.
3.9. MHD velocity eigenvector Before we proceed with the discussion of the CGL velocity eigenvector (which can appear to be a bit tricky), for the sake of clarity it is useful to consider the simpler MHD case first. By following the steps outlined above, it is easy to show that again the Alfvén mode dynamics in the u y decouples from the dynamics in the other directions and the velocity eigenvector matrix for the slow and fast waves in MHD is This yields the MHD dispersion relations for the slow and fast modes in the form of (3.110) or alternatively (3.111). It possible to either express the ratio u x /u z or u z /u x . We want to calculate the velocity ratio χ u and we choose the first possibility that allows us to express that for the slow and fast modes (since u y = 0) (3.118) The ratio u x /u z can be found for example from the second row of the matrix as

P. Hunana and others
or equivalently from the first row of the matrix as . (3.120) The possible range of the ratio u x /u z is from −∞ to +∞ and the possible range for χ u is from 0 to 1. Solutions can be easily divided into two categories -depending on the behaviour for parallel propagation -and are plotted in figure 2. For the sake of clarity, we only consider propagation with positive wavenumbers k and k ⊥ , so the possible range of angle θ is from 0 to π/2 and sin θ 0, cos θ 0. Consider first the case C s < V A . For parallel propagation, the dispersion relation (3.111) directly yields that for the fast mode (ω/k) 2 f | θ=0 = V 2 A and for the slow mode (ω/k) 2 s | θ=0 = C 2 s . Using the expression (3.119) then yields that for the fast mode (u x /u z ) f | θ=0 = +∞ and for the slow mode (u regardless of the actual value of C s /V A < 1. The second case C s > V A is straightforward, since for parallel propagation everything is just turned around; the fast mode dispersion relation is (ω/k) 2 f | θ=0 = C 2 s and the slow mode dispersion relation is regardless of the actual value of C s /V A > 1. The perpendicular propagation limit is even easier since the split into two categories is not required and for both cases the dispersion (3.110) yields (ω/k) 2 f | θ=π/2 = V 2 A + C 2 s and (ω/k) 2 s | θ=π/2 = 0, implying that for the fast mode (u x /u z ) f | θ=π/2 = +∞ and that for the slow mode (u x /u z ) s | θ=π/2 = 0. All the χ u (θ ) curves therefore 'end' at whether the ratio C s /V A < 1 or C s /V A > 1. The behaviour of the χ u (θ ) curves is plotted in figure 2. In MHD, the ratio C 2 s /V 2 A is sometimes used as the definition of plasma beta, which we will call here β MHD ≡ C 2 s /V 2 A . As can be seen in figure 2, for β MHD approaching 1, the curves display a 'singular' behaviour, and it is difficult to guess what will happen when β MHD reaches the value of 1 exactly. This special case is addressed in the next subsection. There are also two special cases, one for C s V A or β MHD 1 and one for C s V A or β MHD 1. Both limits C s V A and C s V A are easily obtained by rewriting the dispersion relation (3.110) to a form (3.124) In both limits the second term under the square root is small and the expression can be expanded as √ 1 − = 1 − /2, yielding the slow and fast mode dispersion relations which further yield dispersion relations (3.126), (3.128). Note that the limit C s V A can be easily obtained directly from the velocity matrix (3.117) by neglecting the offdiagonal terms that decouple the two modes. In contrast, the limit C s V A is not that obvious and one needs to calculate the determinant. In the limit β MHD 1, the MHD dispersion relation simplifies to The expression (3.119) yields for the fast mode (u x /u z ) f (θ ) = ∞ and the expression (3.120) yields for the slow mode (u (3.127) In the opposite limit β MHD 1, the MHD dispersions simplify to (3.128) In this limit, the expression (3.119) yields (u that describes the behaviour of the fast mode for large beta values and the complete and these two curves are plotted in figure 2 as dashed lines. 3.9.1. MHD special case V A = C s Consider the very peculiar case of parallel propagation θ = 0 with V A = C s . The velocity eigenvector matrix (3.117), when evaluated exactly at these values reads It is not possible to determine the ratio u x /u z and the eigenvector must be found by taking the limit θ → 0. For V A = C s the dispersion relation for the slow and fast modes As before, for the sake of simplicity we consider only positive wavenumbers, so the propagation angle θ is in the range [0, π/2] and sin θ 0, cos θ 0. Evaluating the velocity matrix (3.117) at V A = C s only, it is possible to pull the sin θ out of the system so that which nicely demonstrates that for θ → 0 the system approaches (3.130). Nevertheless, we can now calculate the velocity ratio (for example from the second row) as FIGURE 2. MHD velocity ratio χ u (θ) = |u | 2 /|u| 2 for the slow mode (red lines) and the fast mode (blue lines) for different ratios of C 2 s /V 2 A = β MHD . Panel (a) is for β MHD < 1, and the thickness of the lines increases with β MHD when approaching the critical value of β MHD = 1. The plotted lines have β MHD = 0.3; 0.6; 0.8; 0.9; 0.99. Panel (b) is for β MHD > 1 and the thickness of the lines decreases with β MHD when going away from the critical value of β MHD = 1. The plotted lines have β MHD = 1.01; 1.1; 1.2; 2.0; 10; 10 3 . In both figures, the dotted lines are the critical case β MHD = 1 with solutions (3.134). In (b) the dashed lines represent the limit β MHD 1 with solutions (3.129) and the dashed lines nicely overlap with a thin solid line solutions obtained for β MHD = 10 3 . The Alfvén mode has χ u (θ) = 0 regardless of the plasma beta and is not plotted. ratio χ sf u (0) = 1/2 for both the fast and slow modes. The expression (3.132) can of course be directly obtained by using the result (ω/k) 2 = V 2 A ± V 2 A sin θ in the expression (3.119) and by cancelling the sin θ in the nominator and denominator. The expression (3.132) is valid for any angle θ . However, for the slow mode it becomes singular at θ = π/2, and has to be rearranged by multiplying and dividing with 1 + sin θ (or equivalently, use the first row of the matrix). To have useful expression for the entire range of θ, we therefore split the expressions for the slow and fast modes as which yields, in this peculiar case of V A = C s , that the ratio χ u for the slow and fast modes can be expressed as These solutions are plotted in figure 2 as dotted lines.
3.10. CGL velocity eigenvector To find the ratio for velocity components χ u in the CGL description, it is useful for a moment to work with abbreviated symbols, and we write the matrix (3.88) for the u x and u z components as where the coefficients are (3.136) The determinant must be zero, which implies ( ω/ k) 4 − ( ω/ k) 2 (a + c) + ac − b 2 = 0, and has a solution (3.137) The last step was actually used in obtaining (3.98) which is otherwise quite difficult to see. Using this result in the above matrix, the eigenvector matrix becomes The ratio of u x /u z can be expressed (for example from the second row) as (3.140) or alternatively (from the first row) as (3.141) The behaviour of these solutions is plotted in figure 3, where we plot solutions for isotropic temperatures a p = 1. Similarly to MHD, two classes of solutions can be characterized based upon the behaviour for parallel propagation θ = 0, depending if the quantity a − c (3.139) evaluated at θ = 0 is positive or negative. The quantity (a − c)| θ=0 = 1 + 1 50 P. Hunana and others this quantity represents the green dotted curve in figure 1 that separates the system into two regions. For (a − c)| θ=0 > 0, i.e. in the region β < 2/(4 − a p ), for the slow mode (use (3.141)) the ratio (u x /u z ) s (0) = 0 and for the fast mode (use (3.140)) the ratio (u x /u z ) f (0) = +∞. In this region, all the χ u (θ ) curves start at in the region β > 2/(4 − a p ), for the slow mode (use (3.140)) the ratio (u x /u z ) s (0) = +∞ and for the fast mode (use (3.141)) the ratio (u In this region, all the χ u (θ) curves start at For perpendicular propagation θ = π/2 the quantity (a − c)| θ=π/2 = 1 + a p β is always positive and the separation into two regions is not necessary. For the slow mode (u x /u z ) s (π/2) = 0 and for the fast mode (u x /u z ) f (π/2) = +∞. All the χ u (θ ) curves, regardless of the value of β therefore end at This summarizes the general behaviour of the solutions. However, as can be seen from figure 3, there is a special case when the parameters approach β = 2/(4 − a p ), i.e. when (a − c)| θ=0 = 0, that will be discussed in the next subsection. One can again consider the two limiting cases when β is small or large. The β 1 case yields dispersion relations 6 Using expression (3.141) for the slow mode yields (u x /u z ) s (θ ) = 0 and expression (3.140) for the fast mode yields (u x /u z ) f (θ ) = ∞, implying that, in low beta limit, the χ u (θ ) solutions are (3.148) and the constant solutions approaching this limit are visible in the plot of figure 3. The opposite β 1 limit in the CGL description is not particularly revealing since there is no obvious way to simplify the expansion, even for isotropic temperatures. For a p = 1, the β 1 expansion yields One can directly work with the dispersion relation (3.137) and for β 1 and a p bounded and less than an extremely large values, the parameter (a − c) is positive and (a − c) ∼ 1, on the other hand the parameters b 1 and c 1, so an expansion of the square root is possible as (3.145) The fast mode ( ω/ k) 2 f a + (b 2 /a) a, and the slow mode ( ω/ k) 2 s c − (b 2 /a) c, which yields low β dispersion relations. With 'less confidence' the same result can be directly obtained from the velocity matrix (3.135) by neglecting the off-diagonal terms.
where we neglected a term proportional to 1/β . The expression results in the correct parallel limit ( ω/ k) 2 f (0) = 3 2 β , ( ω/ k) 2 s (0) = 1 and the correct perpendicular limit ( ω/ k) 2 f (π/2) = 1 + β , ( ω/ k) 2 s (π/2) = 0. If the perpendicular fast mode limit is sufficient to be ( ω/ k) 2 f (π/2) = β (which is acceptable since β 1), only the terms proportional to β can be retained and the fast mode dispersion relation can be simplified to The expression (3.140) yields the velocity ratio of the fast mode as that is independent of β . For the slow mode wave the dispersion relation cannot be further simplified because if we collect only terms proportional to β , the resulting expression will be non-zero everywhere except exactly at θ = 0 (and also θ = π/2). For θ = 0, the entire contribution of ( ω/ k) 2 s (0) = 1 comes from terms in (3.149) that are not proportional to β . Nevertheless, the analytic velocity ratio (3.152) for the fast mode is sufficiently simple and to obtain analytic χ s u (θ ) for the slow mode, it is actually possible to cheat a little by realizing that the slow and fast mode curves are always symmetric and that χ s u (θ) + χ f u (θ) = 1. Our final large beta results therefore are Solutions are plotted in figure 3 as dashed lines and fit the solutions of the full dispersion relations with β = 10 3 very nicely.
3.10.1. CGL special case β = 2/(4 − T ⊥ /T ) In this peculiar case when all 3 modes propagate in the parallel direction with the same phase speed, one again faces complications when calculating the velocity eigenvector, as with the MHD peculiar case when V A = C s . This CGL special case can be encountered only for β 1/2 because the temperature anisotropy must be a p 0. We start with the velocity eigenvector matrix (3.138) and for the special case considered 1 + 1 2 a p β − 2β = 0, the expression (3.139) simplifies to a − c = (4β − 1) sin 2 θ . One can already see that for strictly parallel propagation the entire matrix will be zero and it will not be possible to evaluate the eigenvector. Similarly to the MHD case, the sin θ can be pulled out of the matrix or alternatively one can use (3.140), (3.141) and calculate the limit. For the special case considered here, the parameter b = 1 2 (4β − 2) sin θ cos θ. Considering again for simplicity only positive wavenumbers, the expression (3.140) yields in this particular case FIGURE 3. CGL velocity ratio χ u (θ) = |u | 2 /|u| 2 for the slow mode (red lines) and the fast mode (blue lines) for different values of β , with isotropic temperature T ⊥ = T (or the temperature anisotropy ratio a p = 1). Panel (a) is for β < 2/3, and the thickness of the lines increases with β when approaching the critical value of β = 2/3. The plotted lines have β = 0.5; 0.6; 0.64; 0.66; 0.666. Panel (b) is for β > 2/3 and the thickness of the lines decreases with β when going away from the critical value of β = 2/3. The plotted lines have β = 0.67; 0.7; 0.8; 1.0; 10, 10 3 . In both figures, the dotted lines are the critical case β = 2/3 with solutions (3.155). In (b) the dashed lines represent the limit β 1 with solutions (3.153), (3.154) and the dashed lines nicely overlap with a thin solid line solutions obtained for β = 10 3 . The Alfvén mode has χ u (θ) = 0 regardless of the plasma beta and is not plotted.
The expression can be evaluated at θ = 0, which yields for the slow mode (u x /u z ) s (0) = −1 and for the fast mode (u x /u z ) f (0) = +1, implying that both modes start at For the perpendicular propagation the discussion is unchanged from the previous general case and yields χ s u (π/2) = 1 and χ f u (π/2) = 0.
3.11. Empirical models with polytropic indices γ , γ ⊥ The differences between MHD and CGL can be further clarified with polytropic indices. The linearized MHD pressure equation reads One can introduce the concept of polytropic indices γ and γ ⊥ to the CGL pressure equations (3.73), (3.74) and write where for the CGL model γ = 3 and γ ⊥ = 2. 7 The reasoning for such a construction is to enable a smooth transition between two extreme cases -the CGL model on one side, and the isothermal model γ = 1 and γ ⊥ = 1 on the other side. The important observation is that in the parallel pressure equation the γ acts only on the parallel velocity component u z , and in the perpendicular pressure equation the γ ⊥ acts only on the perpendicular velocity components u x , u y . In MHD, the γ acts on all three components u x , u y , u z . The values of polytropic indices are directly related to the number of degrees of freedom i, through the well-known relation γ = (i + 2)/i. In the MHD description the number of degrees of freedom is i = 3 and so γ = 5/3. In contrast, the CGL description can be interpreted as being composed of a (strongly coupled) mixture of one-dimensional and two-dimensional dynamics, where for the parallel direction i = 1 and γ = 3 and for the perpendicular direction i = 2 and γ ⊥ = 2. The apparently subtle, but nevertheless fundamental difference between (3.157) and (3.158), (3.159) is the core of many differences between the MHD and CGL descriptions, such as for example the effect of the slow mode becoming faster than the Alfvén mode. The two systems remain different even if the mean CGL pressures (or temperatures) are prescribed to be isotropic p (0) = p (0) ⊥ and equal to the mean MHD pressure p (0) . This means that the MHD does not match the CGL description even if the distribution function is isotropic, as for example Maxwellian.
We briefly consider models with polytropic indices γ , γ ⊥ that are based on the linearized pressure equations (3.158), (3.159). If these pressure equations are 'un-linearized', i.e. if one performs a 'reverse engineering' procedure to obtain nonlinear equations (a procedure that is possible to do because of the known guidance by the CGL model), the nonlinear equations read dp dt which can be rewritten as d dt (3.162) Empirical fluid models of this type were studied for example by  and  with the motivation that spacecraft observations in the Earth's magnetosheath often show behaviour that cannot be explained by the CGL values of γ = 3, γ ⊥ = 2 or by the isothermal values of γ = 1, γ ⊥ = 1, but can be fitted by the values that lie somewhere in between (and actually quite close to the isothermal case). For a direct comparison with the CGL model, the equations (3.162) can be rewritten to a form d dt (3.163) 7 One can then define the parallel and perpendicular sound speeds as C = γ p (0) /ρ 0 and C ⊥ = ⊥ /ρ 0 , since in this model the parallel acoustic mode will have a dispersion ω = ±C k and the perpendicular fast mode ω = ±k ⊥ V 2 A + C 2 ⊥ , as can be easily verified.

54
P. Hunana and others d dt (3.164) or to an interesting alternative form Obviously, for CGL values of γ = 3, γ ⊥ = 2 the right-hand sides disappear and one obtains the classical CGL expressions. The right-hand sides of the above expressions can be compared with equations of Chew et al. (1956) equations (2.82), (2.83) which contain the gyrotropic heat flux contributions. Therefore, the deviations from the CGL values of γ = 3, γ ⊥ = 2 can be considered to offer a very simple empirical model (perhaps too simple) that shows how the heat flux contributions modify the CGL dynamics. The dispersion relations for this fluid model can be found in  and also in Abraham-Shrauner (1973), who studied an even more general empirical model with four polytropic indices (essentially, the γ , γ ⊥ acting on |B| and ρ in (3.162) can be further split to two independent parameters). Interestingly, complex polytropic indices have been suggested to model Landau damping (Belmont & Mazelle 1992). Considering the linearized CGL equations (3.68)-(3.72) and linearized polytropic pressure equations (3.158)-(3.159) written in the x-z plane yields the dispersion relations for this model. The Alfvén wave is still decoupled from the rest of the system and the Alfvén wave dispersion relation is unaffected by the different pressure equations, therefore also yielding the same oblique firehose instability threshold. What is modified is the dispersion relation of the slow and fast waves. The velocity equations read (3.168) and the transformation to Fourier space yields the velocity matrix Calculating the determinant yields a dispersion relation for the slow and fast modes in the form Solutions for parallel propagation (k ⊥ = 0) are the usual CGL Alfvén mode and the ion-acoustic mode ω = ± (γ /2)β k , and for perpendicular propagation (k = 0) the fast mode ω = ± 1 + (γ ⊥ /2)a p β k ⊥ . Similarly to the CGL description it can be shown that A 2 2 − 4A 0 0 because of the trick (3.137), implying that the slow mode becomes unstable when A 0 < 0, i.e. when The two competing instabilities are again the parallel firehose instability and the highly oblique mirror instability. The parallel firehose instability has the same threshold as obtained previously. However, the mirror instability threshold reads which is consistent with the result of , their equation (14), written in the form 8 (3.176) For γ = 3, γ ⊥ = 2, the condition (3.175) is naturally equivalent to the CGL mirror threshold. However, the condition offers a temptingly simple result to keep the perpendicular polytropic index unchanged γ ⊥ = 2 and modifies the parallel polytropic index to γ = 1/2. In this way, the mirror threshold matches the correct kinetic threshold (3.94). Nevertheless, this is just an interesting empirical concept, it fixes the highly oblique mirror threshold but at the same time moves the slow and fast mode dispersion relations in other directions, for example the parallel sound speed in such a model is C 2 = 1 2 p (0) /ρ 0 , which appears to be unrealistically low. Such a 'tuning' of models with free parameters -only in one particular direction and only at the level of the dispersion relations -can be sometimes useful for interpretation of observational data, however, such models will usually lead to unphysical situations. We will not consider models with polytropic indices further in this text. Instead, to obtain a better match with kinetic description, it is advisable to use the correct CGL values of γ = 3, γ ⊥ = 2 and focus on fluid models that derive the correct heat flux contributions q and q ⊥ . 56 P. Hunana and others 3.12. CGL and radially expanding solar wind Similarly to MHD, the CGL fluid model can be very useful for understanding the evolution of temperature in the radially expanding solar wind. A good discussion can be found for example in Matteini et al. (2007), Matteini et al. (2013) and references therein. Let us first consider the MHD case. Combining the pressure equation dp/dt + γ ∇ · u = 0 with the density equation dρ/dt + ρ∇ · u = 0 yields which in a steady state (or co-moving with the convective derivative) implies T ∼ ρ γ −1 . By considering simple radially expanding solar wind, one prescribes that the density ρ evolves with the heliocentric distance r as ρ ∼ 1/r 2 , and by using γ = 5/3, the ideal MHD model yields (in the absence of any dissipation and associated turbulent heating) In contrast, the CGL equations (3.41) are rewritten as which in a steady state implies For the magnetic field it is possible to consider profiles of roughly |B| ∼ 1/r 2 , or |B| ∼ 1/r (one can be more precise and consider the Parker spiral profile). The first magnetic field profile yields |B| ∼ 1/r 2 ; ⇒ T ∼ const.; T ⊥ ∼ 1/r 2 ; The second magnetic field profile yields There is also a curious case when the magnetic field profile is prescribed to be |B| ∼ r −4/3 , i.e. a profile between two cases of 1/r and 1/r 2 , that yields (3.183) and the temperature anisotropy stays constant. The estimations are of course valid only until a firehose threshold or a mirror threshold is reached. For example, evolution according to (3.181) will lead the system to a firehose threshold, and evolution according to (3.182) will lead the system to a mirror threshold. In general, if we fix the density profile to ρ ∼ 1/r 2 , any magnetic field profile steeper than r −4/3 will evolve the system towards a firehose threshold, and a magnetic field profile shallower than r −4/3 will evolve the system towards a mirror threshold. The estimations do not reveal the presence of the firehose instability or the mirror instability in the CGL model, since the solutions are written in a steady state, and concern only evolution of mean temperatures (i.e. the solutions should be written with T (0) , T (0) ⊥ ). To find an instability, one of course needs to consider evolution equations for fluctuating variables and to analyse the associated dispersion relations, since it is the fluctuations that become unstable. Further discussion can be found in .
From a linear perspective, the CGL model has to be accompanied by the Hall term and FLR corrections, so that the instabilities are stabilized at small spatial scales. However, considering nonlinear evolution with finite amplitudes of fluctuations, even the non-dispersive CGL model is stabilized. For example, considering parallel propagation, Tenerani, Velli & Hellinger (2017) has shown that the firehose instability criterion reads (see (9) in that paper) i.e. for infinitely small amplitudes of the magnetic field, the usual firehose criterion is obtained. However, when the CGL system reaches firehose-unstable regime and the amplitude of the fluctuations starts to grow, the firehose criterion is modified according to (3.184). Therefore, the amplitude of the fluctuations cannot grow without bounds, and beyond some critical amplitude, the system becomes again stabilized. The condition (3.184) nicely demonstrates how the CGL system is stabilized during nonlinear evolution, once the firehose instability is reached. For further discussion, see Tenerani et al. (2017), .
To finish this subsection, the heuristically generalized CGL model with polytropic indices γ , γ ⊥ yields that, in the steady state, the mean temperatures evolve according to and such heuristic models can yield a wide array of possible temperature profiles.
3.13. CGL model -summary The usual 'combining' of evolution equations (by performing ∂/∂t and substituting one equation into another) is very useful to gain insight into the linear eigenmodes of a considered system. However, for more advanced fluid models, such a procedure can become very analytically involved, and it is much easier to calculate the determinant of the entire system with analytic software, such as Maple or Mathematica. Therefore, for more advanced fluid models we often do not write down the final dispersion relation for arbitrary propagation angle (which might be uneconomically large to write down), and we only provide normalized, linearized in the x-z plane and Fourier transformed equations. For easy comparison with more advanced fluid models, it is beneficial to write the final CGL equations in the same form. It is useful to work in normalized units, and we drop writing the tilde. The normalized, linearized in the x-z plane, and Fourier transformed CGL equations read 58 P. Hunana and others where we have introduced the normalized parallel Alfvén speed v 2 A ≡ 1 + (β /2)(a p − 1) (the tilde on v A is dropped too). The determinant can be easily calculated and factorized in Maple, yielding the CGL dispersion relation in normalized units (tilde are dropped) (3.195) which of course agrees with the previously obtained CGL solutions for the Alfvén mode (3.81), and the slow and fast modes (3.90).

Hall-CGL model
Here we consider a slightly generalized CGL model when the Hall term is included in the induction equation. Concerning the isotropic Hall-MHD model, a good comparison with solutions of kinetic theory can be found for example in Howes (2009). Before we proceed with the discussion of dispersion relations, we want to point out how the Hall term and the electron pressure contributions and the electron inertia contributions modify the usual 'CGL conservation equations' for protons. The general electric field was already written down in (3.15). For a moment it is beneficial to consider this most general form including the electron inertia and the electron pressure contributions and separate the electric field into two parts (4.2) The notation E H is perhaps slightly misleading since only the first term in (4.2) is usually called the Hall electric field. Note that the first term in E H still has to be used in the proton momentum equation even for the simplest CGL (or MHD) models. The induction equation then reads which is then rewritten with the use of the convective derivative as The last expression is then used in the completely general equations (3.46), (3.47), which yields d dt (4.7) By cancelling pressure equations (3.3), (3.4), that are also the pressure equations of the Hall-CGL model, one obtains that for a general Hall-CGL model the electric field E H modifies the conservation laws according to d dt implying the Hall term breaks the first and second adiabatic invariants.

Hall-CGL model with cold electrons
Here, we want to study the simplest Hall-CGL model for the proton species. We neglect the electron inertia and we make the electrons cold with p e = 0. The full model is described by the usual CGL equations (3.36)-(3.39), but the induction equation now reads The induction equation is linearized as By using (3.67), direct calculation of the second term yields (4.12) and by using c/(4πen 0 ) = V 2 A /(Ω p B 0 ), the Hall-term contributions are (4.13) 60

P. Hunana and others
The entire linearized induction equation reads (4.14) which when written in the x-z plane (all ∂ y = 0) yields These equations replace the ∂B/∂t equations in the linearized CGL system written in the x-z plane. And the trouble is immediately apparent. In contrast to the CGL system, we are no longer able to separate the Alfvén mode which was before exclusively in the B y , u y components. Now, all the components are coupled and all three modes, the Alfvén mode, the slow mode and the fast mode are coupled. This is not surprising, since the same situation is in the Hall-MHD. In this case, we have no other choice and we have to calculate the dispersion relation for all three modes. The dispersion relation can be solved analytically only for special cases and, in general, the dispersion relation has to be solved numerically. Let us quickly consider the two special cases of parallel and perpendicular propagation that can be solved analytically.

Parallel propagation
For propagation parallel to the mean magnetic field, all ∂ x = 0 and the familiar sound mode/ion-acoustic mode decouples in the u z , p components as with the familiar dispersion relation ω 2 = C 2 k 2 . The dispersion relation of this mode does not depend on the Hall term. However, the other two modes are now coupled through If the last terms in the B-field equations are neglected, i.e. if the Hall term is neglected, the two modes decouple and one naturally recovers the CGL parallel propagating Alfvén modes, both propagating with the phase speed Applying ∂ t to the B-field equations yields which when transformed to Fourier space read (4.26) and the zero determinant requirement implies The solutions of the dispersion relation (quadratic polynomial in ω 2 ) is simply where, importantly, k 2 = |k | was used. The solution (4.28) is the whistler mode, and the solution (4.29) is the ion-cyclotron mode. For isotropic temperatures (a p = 1), the result is equal to the Hall-MHD dispersion relation. At long wavelengths (k → 0) the solutions become non-dispersive and connect to the usual CGL Alfvén modes. For parallel propagation, the Hall term is therefore responsible for splitting the two modes that otherwise would propagate at the same speed, while it does not influence the ion-acoustic mode. Or in another words, using the vocabulary of slow, Alfvén, fast, depending on plasma β and temperature anisotropy a p , the Hall term either splits the parallel propagating Alfvén and fast mode, or the slow and Alfvén mode. At very short spatial scales (k 1), the limits are 9 k → ±∞ : (4.30) 62

P. Hunana and others
For isotropic temperatures a p = 1 the frequency of the ion-cyclotron mode (that should be actually called the proton-cyclotron mode) converges towards the proton-cyclotron frequency ±Ω p . For anisotropic temperatures, the frequency of the mode starts to be β dependent and will converge to frequencies that can be higher or lower than Ω p . The frequency of the whistler mode is not bounded and technically goes to infinity, which is a consequence of neglecting the electron inertia. Importantly, both modes are stable for high wavenumbers. The solutions (4.28), (4.29) are analytically correct, but similarly to the Hall-MHD, the solutions can be rewritten to a more useful form, by realizing that the dispersion equation (4.27) can be decomposed as 10 The first bracket yields two solutions (4.33) and the second bracket yields two solutions It is straightforward to check that (ω − ω W,1 )(ω − ω W,2 ) = 0 yields the whistler mode solution (4.28), and (ω − ω IC,1 )(ω − ω IC,2 ) = 0 yields the ion-cyclotron mode solution (4.29). Often, solutions (4.32)-(4.35) are written in an alternative form where the absolute value on k is removed (which obviously can be done), and the solutions are just slightly re-arranged. Additionally, for more complicated fluid models, we want to write solution in one line instead of four lines. In this example, we just write a shortcut that two solutions are with another two solutions obtained by substituting ω with −ω.
10 One can double check the result by working with quantities

Firehose instability
Now we can discuss the firehose instability. If the quantity 1 + (β /2)(a p − 1) > 0, then all the solutions (4.32)-(4.35) have frequencies that are purely real for the entire range of wavenumbers. For isotropic temperatures a p = 1, the dispersion relations do not depend on the value of β , which is a consequence of neglecting the FLR pressure corrections, as will be discussed later, and all solutions are equivalent to the Hall-MHD dispersion relations. Solutions (4.32)-(4.35) with a p = 1 are plotted in figure 4(a), where we plotted both positive and negative wavenumbers and positive and negative frequencies to clearly understand how the modes are connected to each other. The magenta and cyan lines are plotted as dashed lines, so that the growth rates in the firehose-unstable regime (plotted in figure 5) are clearly visible. For anisotropic temperatures a p = 1, the solutions are naturally β dependent. To show solutions in a firehose unstable regime, one has to choose β > 2 since the instability does not exist for lower values of β . In figure 4(b), we have chosen a value of β = 4 and varied the temperature anisotropy from a p = 2 down to a critical value of a p = 0.5. In this regime, all the modes have purely real frequencies.
Importantly, if the temperature anisotropy parameter a p is further decreased, the frequencies of all the modes will become complex numbers (with a real and imaginary part) in the region where which can be viewed as a modified Hall-CGL firehose instability threshold that is now length-scale dependent. At long spatial scales (k → 0), the above criterion is naturally equivalent to the CGL firehose criterion. The term proportional to k 2 in (4.37) is always positive and as k increases, after some critical wavenumber the frequency 64 P. Hunana and others  (a), has an imaginary part of the frequency which is very close to zero for the entire range of k . For this reason this mode is not plotted in (b), where instead a mode with a p = 0.4 is plotted. A mode that has Im ω > 0 is unstable and growing in time. A mode that has Im ω < 0 is stable and damped. According to the figure, whistler modes with ω r > 0 are unstable, and ion-cyclotron modes with ω r < 0 are unstable. It is assumed that √ −1 = +i. The solutions for ω r < 0 are here non-causal, and (b) should be re-plotted with causal solutions (4.43)-(4.46), so that whistler modes are always unstable, and ion-cyclotron modes stable. of modes will become purely real. The Hall term is therefore responsible for the stabilization of the firehose instability at sufficiently high wavenumbers (small spatial scales). The solutions in this firehose-unstable regime are plotted in figure 5, where again β = 4 and the temperature anisotropy is varied from a p = 0.49 to the lowest possible value of a p = 0.0. Figure 5(a) shows real frequencies and the (b) shows imaginary frequencies.
It is of interest to determine which mode becomes unstable. According to our Fourier transform (see appendix A), a mode with positive imaginary frequency is unstable and growing, and a mode with negative imaginary frequency is stable and damped. If the firehose threshold (4.37) is satisfied, for some range of k the second terms in solutions (4.32)-(4.35) become purely imaginary and the first terms stay purely real. Consider one of the four possible quadrants in figure 5(a), for example the quadrant with k > 0; ω r > 0. Up to some critical k where the Hall-CGL firehose criterion is satisfied, both the blue mode and the cyan mode propagate with the same real frequency ω r = +k 2 V 2 A /(2Ω p ). The differences are in the imaginary frequency, see figure 5(b), where in this quadrant the blue mode is unstable and the cyan mode is stable (and damped). After a critical k (for example for a p = 0 it is k V A /Ω p = 2), both modes become purely real and evolve with different real frequencies, here the blue mode is obviously the whistler mode and the cyan mode is the ion-cyclotron mode. The same situation is in the quadrant k < 0; ω r > 0. We can therefore conclude that for positive real frequencies, ω r > 0, the whistler mode is unstable and the ion-cyclotron mode is stable and damped. Now we could say that considering both positive and negative wavenumbers and frequencies is redundant, and that we are finished. However, from a perspective of this guide, we consider useful to show that one might easily come into contradictions, if negative frequencies are not handled with sufficient care. Importantly, according to figure 5, for negative real frequencies ω r < 0, the whistler mode is stable and the ion-cyclotron mode is unstable. This is indeed contradictory. Keeping the k > 0, a change from positive to negative ω r represents a change of direction of propagation along B 0 . One can argue, that in a homogeneous system (e.g. excluding propagation in a stratified fluid), a change of direction of propagation should not yield a change of stability. Additionally, kinetic theory and simulations show, that it is the whistler mode that is firehose unstable for both ω r > 0 and ω r < 0, and the ion-cyclotron mode is stable. It is important to emphasize that, from a perspective of solutions (4.32)-(4.35), at the range of wavenumbers where the firehose instability is present, the whistler and ion-cyclotron modes propagate (in a given quadrant) with the same real frequency. Moreover, both modes have the same polarization of electric and magnetic field (see § 4.2.3). Both modes are completely degenerate, and not distinguishable. Thus, how the modes are continued for higher wavenumbers once the firehose instability disappears, is not obvious (since √ −1 can be possibly both +i or −i). One should consider all the possible solutions and verify, which solutions are physically plausible. Therefore, after the firehose threshold is crossed, the analytic solutions (4.32)-(4.35) are not correctly separated. Obviously, more 'physical information' is needed to distinguish between the ion-cyclotron and whistler modes in the firehose-unstable regime, and to determine which solutions are causal and which are non-causal.

Generalization to causal (correct) solutions
The mind boggling puzzle as why the Hall-CGL solutions (4.32)-(4.35) are in contradiction to kinetic theory and to common sense (a discrepancy found by one of the referees of this text), was beautifully solved by Paul Cally. The easiest way to obtain correct solutions is to introduce a very small 'causal' dissipation into the system, and remove it later. Alternatively, it should be possible to obtain the same result by using Laplace-Fourier transforms on temporal-spatial scales and considering an initial value problem.
We introduce causality, by adding a tiny amount of dissipation on the right-hand sides of momentum equations (4.19)-(4.20), in the following form where is a dimensionless real small positive number, for example = 10 −8 . The dissipation will be later completely removed by the limit lim →0 + . The V 2 A /Ω p is there so that the parameter is indeed dimensionless (since here we usually normalize with respect to V A and Ω p , but other normalizations can be of course considered). The factor of 2 is there for pure convenience, so that the final dispersion relations (4.43)-(4.46) are of the utmost beauty, and that we do not have to later redefine . The momentum equations are accompanied by the induction equations (4.21)-(4.22). The reader is encouraged to verify the calculations by switching to normalized variables − 1), we continue in physical units.

P. Hunana and others
By adding the dissipation, the dispersion relation (4.31) is modified to the following form It is useful to define the Hall-CGL firehose threshold and similar generalizations can be done for models with FLR corrections. Then, since is small, it is easy to show that the non-causal solutions (4.32)-(4.35) are 'corrected' by the tiny dissipation to the following causal form The reader is encouraged to plot these solutions, for example with = 10 −8 , and verify that figure 5(b), is indeed corrected, i.e. that the whistler is unstable for both ω r > 0 and ω r < 0. We want to further address the limits. For clarity, let us first consider the following simplest case lim where one can further suppress the small real /2 on the right-hand side. The limit (4.47) provides an answer, when √ −1 is equal to +i or −i. The result (4.47) can be easily generalized, and limits in expressions (4.43)-(4.46) are calculated according to (4.48) Obviously, the whistler solutions (4.44), (4.45) are unstable, and the ion-cyclotron solutions (4.43), (4.46) are damped. Only now we can conclude with confidence, that the parallel firehose instability as an instability of the whistler mode. We note that if instead of the momentum equations, we added the dissipation (diffusivity) to the induction equations, the solutions (4.43)-(4.46) are recovered with → − . In other words, by adding dissipation to the induction equations, the whistler mode would be always damped and the ion-cyclotron mode unstable. The result makes sense with respect to kinetic theory, since the ion-cyclotron resonances enter through the velocity moment (by integration over the perturbation of the distribution function f (1) ), or more precisely through the current j = r q r n r u r , see Part 2.

Polarization
The polarization of the magnetic field can be easily checked, actually without plugging any dispersion relations into the matrix (4.26), since the matrix can be represented simply as It is useful to define the polarization as the angle between the complex B x and B y components, which is calculated as Arg(B y /B x ). If the result is smaller than zero, the wave is left polarized. If the result is larger than zero, the wave is right polarized. The determinant of the matrix (4.49) is a 2 − b 2 = 0, so a = ±b. The case a = −b (which is directly equal to the first bracket in (4.31) with solutions (4.32), (4.33)) yields that B y = −iB x , so that (4.50) The second case a = +b (which is equal to the second bracket in (4.31) with solutions (4.34), (4.35)) yields that B y = +iB x , so that Similarly, one can calculate the polarization of electric field E, and (when not specified that the B field is considered), the words left/right polarized are usually meant for the E field. The ratio E y /E x can be easily calculated, for example from equations (4.74), (4.75) (written there in the normalized form), which when multiplied by ω and use of momentum equations, yields for the special case of parallel propagation considered here (4.52) When plotting many solutions, it is useful to additionally divide the angle by π, and define the B and E field polarizations as (4.53) Electric field polarization defined this way was used for example by Hunana et al. (2013) to investigate properties of (highly oblique) kinetic Alfvén waves in various fluid models and the kinetic description; and also by Camporeale & Burgess (2017), who compared various linear modes in hybrid, gyrokinetic and fully kinetic descriptions. 11 When P E < 0 the wave is called left polarized and when P E > 0, the wave is called right polarized. For the first group (4.50) this yields E y /E x = −i and 68 P. Hunana and others P E = −1/2 < 0 (left polarized wave), and for the second group (4.51) E y /E x = +i and P E = +1/2 > 0 (right polarized wave). Finally, it is of interest to note another possible definition of polarization, that is defined with respect to positive real frequency, i.e. the electric field polarization (4.53) is redefined to (4.54) A similar definition of polarization with sign(ω r ) is used for example in the book by Gary (1993), page 91. With such a definition, the whistler mode is always right polarized, both in the firehose-stable and firehose-unstable regimes, since for the whistler mode ((4.34), blue line) the real frequency remains ω r > 0, and for the whistler mode ((4.33), magenta line) ω r < 0. The situation is more confusing for the ion-cyclotron mode, where the firehose-stable and firehose-unstable regimes have to be considered separately. (i) In the firehose-stable regime, both ion-cyclotron modes are left polarized, since for the mode ((4.32), green line) ω r > 0, and for the mode ((4.35), cyan line) ω r < 0. (ii) In contrast, in the firehose-unstable regime, both ion-cyclotron modes are right polarized, since for the mode ((4.32), green line) ω r < 0, and for the mode ((4.35), cyan line) ω r > 0.

Simplest ion-cyclotron resonances
When talking about 'resonances', one often expects to see expressions with a denominator that becomes zero. Generally speaking, the simplest 'singular' expressions representing resonances just arise from a technique used in deriving the dispersion relations. If one eliminates the electric field E from the beginning and expresses the fluid model through a matrix multiplied by a vector (ρ, u x , u y , u z , b x , b y , b z , p , p ⊥ ), leads to dispersion relations that are not singular explicitly. In contrast, expressing the model through a matrix multiplied by a vector (E x , E y , E z ), yields equivalent dispersion relations where the resonances are shown explicitly. The dispersion equation (4.31) can be easily rewritten as , (4.55) which, in the cold plasma limit (β = 0) or in this case also for isotropic temperatures (a p = 1), simplifies to making it more evident that we have the simplest ion-cyclotron resonances here too. Indeed, for k → ∞ the frequency ω → Ω p (first case) and ω → −Ω p (second case). For a plasma with finite β and non-isotropic temperatures, the k → ∞ limit yields more general ion-cyclotron resonances (4.30).
To better understand what effects are present and what effects are absent in the Hall-CGL model, we plot fully kinetic solutions obtained with the WHAMP code in figure 6. The proton temperature is prescribed to be isotropic T ⊥ /T = 1 and the electrons are prescribed to be cold with T e /T p = 0 (in the WHAMP code, the value is chosen to be 10 −8 ). In the WHAMP code it is necessary to prescribe the ratio of (4.60) 70 P. Hunana and others ∂p ⊥p ∂t + ∇ · (p ⊥p u p ) + p ⊥p ∇ · u p − p ⊥pb · ∇u p ·b = 0. (4.61) Linearizing in normalized units is even easier, because (temporarily introducing back the tilde for clarity) p (0) p = 1; p (0) ⊥p = a p ; ρ 0 = 1; B 0 = 1. Linearizing these equations, writing them in the x-z plane and Fourier transforming yields (dropping also the proton index 'p') −ωρ + k ⊥ u x + k u z = 0; (4.62) where the normalized v 2 A = 1 + (β /2)(a p − 1). Alternatively, it is sometimes illuminating to work with the general induction equation where the electric field is used (4.73) and in this case with cold electrons the electric field components are (4.74) E y = u x + ik B y ; (4.75) E z = 0. (4.76) The induction equation (4.71)-(4.73) with the above electric field components is of course equivalent to the induction equation (4.66)-(4.68). All equations are normalized, the tilde are dropped, and the electric field (4.78) (1 + a p β ) + 3 2 β 1 + a p β − 1 6 a 2 p β + k 2 k 2 β 3 2 k 2 + a p k 2 ⊥ ; (4.79) where the wavenumber k 2 = k 2 + k 2 ⊥ . Perhaps the nicest form is to keep the CGL contributions on the left-hand side (so the terms can be factorized to form the CGL dispersion relation) and put the Hall term contributions on the right-hand side, which yields the Hall-CGL dispersion relation Solutions for the parallel propagation (k ⊥ = 0) are of course the whistler and the ioncyclotron mode with other two solutions obtained by exchanging ω with −ω. By bringing back the tilde to the normalized solution (4.84), the result (4.36) in physical units is easily recovered. The other parallel solution is the ion-acoustic mode ω = ± 3β /2k . For perpendicular propagation (k = 0), the fast mode solution is ω = ±k ⊥ 1 + a p β , which is equivalent to CGL since the Hall term disappears. We investigated solutions of the Hall-CGL model in much greater detail in our paper 'On the Parallel and Oblique Firehose Instability in Fluid Models' ), a paper which was essentially pulled out of this guide. We do not want to repeat the discussion with many associated figures here, and a reader who is further interested in the Hall-CGL firehose instability can find much more information in that paper.

FLR corrections to the pressure tensor
The finite Larmor radius corrections to the pressure tensor can be derived at several levels of approximation. The pressure tensor equation (2.50) can be rewritten aŝ This equation is exact and since p = p g + Π, the FLR tensor Π is described implicitly.

Fully nonlinear FLR corrections
If preservation of all nonlinearities is desired, this implicit equation can be further rewritten by the following procedure that can be found for example in Hsu et al. (1986), Passot & Sulem (2004) and Ramos (2005) as a brief note. By applyingb× to the left-hand side, the first term calculates as where we have used that ikl lrs = δ ir δ ks − δ is δ kr . The second term on the left-hand side of (5.1) calculates quite differently and we will need an identity 72

P. Hunana and others
Let us also remind ourselves that in the index notation Π ss = Tr Π = 0 andb kbl Π kl = Π :bb = 0. The second term calculates as The entire left-hand side therefore reads Performing a usual matrix (single) dot product of this result with matrixbb yields This further yields that and adding this result together with its transpose implies where we have used that the Π ij = Π ji . One therefore derives that the FLR tensor can be extracted from the left-hand side of (5.1) by performing an 'inversion' procedure We can introduce matrix κ that will represent the right-hand side of (5.1) aŝ By performing operations (5.9) in this equation yields the FLR tensor in the following form The expression (5.12) is completely general since no simplifications were introduced, the expression is exact. However, the equation is still implicit, since on the righthand side the tensor κ contains the full pressure tensor p = p g + Π. Nevertheless, the equation is extremely useful, once the expansion of the right-hand side is performed, as we will do below. Note that the brackets are not unique, since where the last choice leaves the brackets unspecified (since both options are allowed), similarly tob · ∇u. Alternatively, one can get rid of the operator ·I (which acting on a matrix yields the same matrix), and split the result into two parts The definition is of course motivated by the pressure decomposition (2.46). By applying this operator to (5.10), the left-hand side remains the same, yielding The solution (5.12) for the FLR tensor therefore rewrites The subtraction of scalar pressure equations is motivated by the observation that when working in the linear approximation directly with (5.10), i.e. without performing the inversion procedure, the scalar pressure equations have to be subtracted at the end. However, by performing the inversion procedure, the scalar pressure equations are subtracted 'automatically' during calculations, and it is actually easier to work directly with (5.12) instead of (5.16). Therefore, let us work with (5.12). The leading-order nonlinear FLR corrections are obtained by neglecting the non-gyrotropic contributions Π and q ng on the right-hand side of (5.12), i.e. by prescribing Therefore, the only non-gyrotropic quantities that were retained on the right-hand side of (5.12), are the perpendicular components of the velocity u. The slightly unappealing factors B 0 /|B| in all the expressions can be removed, by redefining the cyclotron frequency with |B| instead of B 0 , so that Ω |B| = e|B|/(mc). Let us calculate (5.17) step by step. For the brevity of the calculations, it is useful to introduce convective derivative d/dt = ∂/∂t + u · ∇ and work with

P. Hunana and others
We want to present complete derivation with the least amount of algebra involved, and we first want to obtain the form of Schekochihin et al. (2010). Let us first write each term in κ (1) in the index notation. By using in the following form and grouping similar terms in κ (1) together yields The first line in the above expression is left unchanged, since it is the final result, defined below as a matrix W ij = p ⊥ (∂ i u j + ∂ j u i ) + ∂ j (q ⊥bi ) + ∂ i (q ⊥bj ), that can be written as W = p ⊥ ∇u + ∇(q ⊥b ) + (p ⊥ ∇u + ∇(q ⊥b )) T = (p ⊥ ∇u + ∇(q ⊥b )) S . However, calculating the rest of the expression according to (5.17) simplifies several terms, and for example all terms in the third line of (5.27) vanish. The terms in the third line containingb ibj vanish after applyingb×, since (b × (bb)) ij = iklbkblbj = 0. The terms in the third line containing δ ij vanish after applying the symmetric operator since (5.28) We therefore need to concentrate only on the second line of (5.27), and slightly rewrite that one. For example, the first term in the second line of (5.27), the term d/dt(b ibj ), calculates as The second term in the second line of (5.27), the termb ib · ∇u j , vanishes after applyingb×, since The third term in the second line of (5.27), the termb jb · ∇u i however does not vanish. Note that b jb · ∇u i = ((b · ∇u)b) ij and the third term calculates as The fourth term in the second line of (5.27), calculates in the same fashion than the second and third terms sinceb Collecting all the results together yields the FLR tensor where W is a matrix and w is a vector. The result is equivalent to (5)-(8) of Schekochihin et al. (2010), derived in the Appendix of that paper in a somewhat simpler way. In Schekochihin et al. (2010), instead of our notation for matrix W and vector w, the symbols S and σ are used, where S = W and σ = −w. We wanted to avoid the use of S and σ , since these symbols are typically used when decomposing 76 P. Hunana and others the heat flux tensor. For a reader who just jumped to this result, the symmetric operator is defined as A S ij = A ij + A ji . For clarity, for any vector w and any matrix W , the following identities hold The last identity follows from an identity for general matrix A, that reads The result (5.33)-(5.35) can be slightly simplified. The matrix W in the index notation reads and by applyingb× eliminates only the last term implying that the matrix W could possibly be redefined to However, now the matrix is not symmetric. Note that the brackets in the expression (5.33) are not unique, and one can for example pull theb× out, so that the expression reads Let us multiply the last term in (5.40) by ·(I + 3bb), it calculates as Obviously, it is beneficial to move this term to the vector w. Therefore, redefining W and w, the solution reads Now the matrix W is again symmetric, and w is just a vector. This solution with matrix (5.44) is slightly nicer than the solution with matrix (5.34). Alternatively, one can write the solution in its full form The nonlinear solution for the FLR tensor Π can be further re-arranged. For example, to obtain the form reported by Ramos (2005), one needs to separate a part that is obtained by performing Π ·b. The solution (5.43) in the index notation reads and direct calculation yields where it was useful to define the new vector h ≡ W ·b + w, by adding the quantity W ·b = p ⊥ ((∇u) ·b +b · ∇u) + q ⊥b · ∇b (5.50) to the vector w. By separating the Π ·b part, the solution (5.43) for Π is therefore re-arranged to the following form (5.51) (5.52) (5.53) The solution (5.51) has a nice advantage in that it is decomposed with respect to the magnetic field lines. In other words, projecting Π along the field lines by performing Π ·b cancels the first term, and directly yields vectorb × h/Ω |B| , that could be defined as vector Π z .

Employing a non-dispersive induction equation
Considering the long-wavelength low-frequency limit, the solution (5.45) for vector w and the solution (5.53) for vector h can be further simplified, by using the usual 78 P. Hunana and others non-dispersive CGL/MHD induction equation, that can be written in the following form (see later in the text) (5.54) By further applyingb× to the induction equation, the last term on the right-hand side disappearsb and vectors (5.45), (5.53) simplify to The result (5.57) can be further rewritten with vorticity ω = ∇ × u, by using identity Results (5.51), (5.52), (5.58) are equivalent to equations (49)-(51) of Ramos (2005) (see also (11)-(13) of Macmahon (1965)). Further evaluation of nonlinear FLR corrections to higher precision is out of the scope of this simple guide, and the interested reader is referred to Ramos (2005). We will continue in a much simpler linear approximation.

Simplest FLR corrections (FLR1)
We could continue by evaluating the obtained nonlinear results in the linear approximation. However, the nonlinear FLR calculations were quite complicated and it is very useful to clarify how we derive the simplest FLR corrections directly in the linear approximation. Let us start again, and work with the general equation (5.1). To derive the simplest finite Larmor radius corrections to the pressure tensor, one has to cancel the Π contributions on the right-hand side of equation (5.1) and also neglect the heat flux contributions, which yieldŝ (5.59) Considering complete linearization (so that we can easily go to Fourier space), this is equivalent to performing expansion to the first order in frequency ω/Ω and to the first order in wavenumber r L k. We do not want to perform complete linearization just yet. However, to move further, we introduce simplification by not evaluating the FLR corrections along the magnetic field lines, but along the ambient magnetic field B 0 that is prescribed to be in the z-direction, and thereforeb = (0, 0, 1) ≡b 0 . This simplification, typically used in numerical simulations, is appropriate if the magnetic field lines are not too distorted with respect to the mean magnetic field. Here, we will distinguish between complete linearization and evaluation alongb 0 . However, it is important to note that the procedure of evaluating some terms withb 0 , while Collisionless fluid models. Part 1 79 other nonlinearities are kept, is of course not a fully consistent procedure. With the approach presented here, the only fully consistent procedure that simplifies the general equation (5.1), is complete linearization. Additionally, as discussed by Passot et al. (2014), neglecting magnetic field line distortion may in some instances lead to spurious instabilities. Let us evaluate (5.59) alongb 0 . On the left-hand side we calculate each component in the matrix separately. Since (b × Π) ij = iklbk Π lj , it is obvious that non-zero contributions are obtained only ifb k =b z = 1, so the index k = z, and (b × Π) ij = izl Π lj . Now, for example, the xx component calculates as (b × Π) xx = xzl Π lx and the only non-zero contribution is obtained for index l = y, which yields (b × Π) xx = xzy Π yx = − xyz Π yx = −Π yx = −Π xy , where in the last step we used that the FLR tensor is symmetric. The entire left-hand side The gyrotropic pressure tensor is of course The first term on the right-hand side of (5.59) calculates as 62) and using the matrix notation The second term in (5.63) contributes because of the induction equation ∂B/∂t = −c∇ × E. For now, we will not use the induction equation and keep our FLR calculations in the above form. Let us continue with the second term on the right-hand side of (5.59). When completely linearized this term does not contribute, but for now we will keep it and combine it with the first term into a convective derivative ∂p g /∂t + u · ∇p g = dp g /dt, so that ∂/∂t in (5.63) is replaced by d/dt. The remaining terms on the right-hand side of (5.59) are straightforward to evaluate since no derivatives ofb are encountered, and the third term p g ∇ · u|ˆb =b 0 is equal to matrix (5.61) multiplied by ∇ · u. The last two terms of (5.59) are together evaluated as 2p ⊥ ∂ y u y ; p ⊥ ∂ y u z + p ∂ z u y p ∂ z u x + p ⊥ ∂ x u z ; p ∂ z u y + p ⊥ ∂ y u z ; 2p ∂ z u z   .

P. Hunana and others
Now we need to subtract the scalar CGL pressure equations obtained previously. Sincê b · ∇u ·b|ˆb =b 0 = ∂ z u z , the scalar pressure equations read dp dt By cancelling terms coming from these two equations (affecting only the diagonal components), the entire equation (5.59) rewrites The identity Π :bb = 0 evaluated forb =b 0 implies that Π zz = 0 and the identity Tr Π = 0 further implies that Π xx = −Π yy . The components of the pressure tensor therefore read Finally, the time derivative of the unit vector The general induction equation reads Collisionless fluid models. Part 1 81 which evaluates in theb =b 0 approximation as These expressions are used in the FLR tensor components Π xz and Π yz , and the FLR tensor reads However, since the FLR tensor was derived with a precision to only first order in wavenumber k, to make everything consistent, the Hall contributions in the induction equation are usually (but not always) neglected. We therefore have db x /dt = ∂ z u x and db y /dt = ∂ z u y which yields the leading-order FLR pressure tensor evaluated along the magnetic fieldb = (0, 0, 1) in the form which is equivalent to equation (5) of Yajima (1966). Perhaps surprisingly, the result is also equivalent to the highly collisional equations (2.20)-(2.24) of Braginskii (1965), after one prescribes p = p ⊥ and after terms containing the collisional time τ are 'ignored' (τ ∼ 1/ν where ν is the usual collisional frequency). A proper collisionless limit τ → ∞ cannot be obtained from the expressions of Braginskii (1965), since the results are proportional to τ (and also 1/τ ). Furthermore, evaluation alongb = (0, 0, 1) basically means that the magnetic field terms of the system were linearized and other parts of the system were not. The expansion procedure preserved some nonlinearities, but other nonlinearities that are of the same order were implicitly neglected. Therefore, the only fully consistent procedure is to evaluate the FLR tensor in the linear approximation, yielding Notably, the tensor is very different from (3.10)-(3.13) of Oraevskii et al. (1968). The FLR corrections to the pressure tensor (these ones or more complicated ones) were used in a number of direct numerical simulations, see for example Borgogno et al.  −ωρ + k ⊥ u x + k u z = 0; (5.75) where v 2 A = 1 + (β /2)(a p − 1), and the first-order FLR tensor Π = Π (1) is The dispersion relation for the Hall-CGL-FLR1 fluid model can be written as where the polynomial on the right-hand side (5.91) The solution for strictly parallel propagation is For strictly perpendicular propagation the solution is and the result is consistent with equation (17) of Yajima (1966), after neglecting the electron inertia in that paper (by ω 0 → ∞). The strictly perpendicular propagation is an excellent example to point out the deficiencies of the Hall-CGL-FLR1 fluid model, or actually the CGL-FLR1 model, since the Hall contributions are zero for k = 0. The result (5.93) is also equal to equation (34) of Del Sarto et al. (2017), who discuss in detail that, compared to solutions of a truncated kinetic Vlasov system, and also of a fluid model with more precise FLR corrections, the solution (5.93) actually has the opposite sign in front of the FLR term ∼k 2 ⊥ . For further refinement with FLR2, see (5.117) and the subsequent discussion.

CGL-FLR1 fluid model (no Hall term)
If the Hall term is neglected, the general dispersion relation of the CGL-FLR1 fluid model reads where the CGL parameters A 2 , A 0 are (5.87) and the FLR1 polynomial on the righthand side has the form

P. Hunana and others
where v 2 A = 1 + (β /2)(a p − 1). For strictly parallel propagation, the solution for the whistler and ion-cyclotron waves reads (5.98) For strictly perpendicular propagation the solution is equivalent to (5.93) since the Hall contributions vanish.
5.5. Second-order FLR corrections (FLR2) If more precise FLR corrections are desired, it is relatively easy to increase the precision to the second order in frequency ω/Ω while keeping the precision in wavenumber kr L at first order. The pressure tensor equation (5.1) is simplified tô ( 5.99) For clarity, let us first consider the case with the heat flux q = 0. Compared to the firstorder FLR corrections in the previous section, we just have one more matrix present, (∂/∂t)Π. As before, the FLR tensor must be symmetric and also satisfy conditions Π zz = 0 and Π yy = −Π xx . Following the previous derivation, the system (5.66) now includes also the time derivative of the FLR tensor on the right-hand side. Let us write down the system directly in the linear approximation. The system reads There are only four FLR components that have to be considered and which can be rewritten separately as Collisionless fluid models. Part 1

85
The FLR components could be considered as independent fluid quantities described by the time-dependent equations as written above. However, from a linear perspective this introduces additional linear modes, and from a direct numerical simulation perspective it is preferred not to introduce an additional four time-dependent equations (nevertheless, numerical simulations with the full pressure tensor equation are sometimes performed, see e.g. Wang et al. (2015), Del Sarto, Pegoraro & Califano (2016) and Ng et al. (2017)). The above equations are therefore typically simplified by expanding the FLR tensor to the first and second orders as Π = Π (1) + Π (2) and the time derivative of the second-order tensor is neglected, i.e. (∂/∂t)Π (2) = 0 and the system reads where we also used the linearized induction equation in the last two equations. Now, there are two possibilities to handle the Hall-term contributions coming from the induction equation. It is possible to either define Π (1) to be equal to set (5.74) and move the Hall term to Π (2) , or it is possible to keep the Hall term in Π (1) that is equal to linearized set (5.72). The first choice, i.e. if Π (1) is set, (5.74) yields the second-order components The second choice, i.e. when Π (1) is equal to linearized set (5.72), yields The difference between the two approaches therefore is that the first approach neglects the time derivative of the Hall term. It is noted that the second approach is not necessarily better or more accurate than the first approach. For example, for parallel propagation, the frequency of the whistler mode seems to be better described by the first approach and the solutions are closer to kinetic theory.

Hall-CGL-FLR2 fluid model
The second-order FLR tensor (5.109) is normalized and written in Fourier space according to where for cold (and massless) electrons in the x-z plane We purposely wrote the equations with E H so that the generalization to a more elaborate electric field will be straightforward. Note that the Hall contributions (5.112) in the FLR tensor completely vanish for isotropic temperatures a p = 1, regardless of the form of E H , i.e. regardless if the electrons are cold or not.
To write down the solution for parallel propagation (k ⊥ = 0) it is useful for brevity to introduce the following quantity 113) and the parallel propagating whistler and ion-cyclotron modes factorize as (5.114) and the two solutions are with two further solutions corresponding to −ω on the left-hand side. It is perhaps useful to rewrite expression (5.114) to the 'ion-cyclotron' resonance form  of Del Sarto et al. (2017), obtained by solving the pressure tensor equation without the expansion to Π (1) and Π (2) , and by expanding for small wavenumber afterwards. The difference between the FLR1 solution (5.93) and the FLR2 solution (5.118) is perhaps even more clear, when written in physical units for ω and k, and by using the gyroradius ρ i = v th⊥ /Ω p , in the following form . As discussed in that paper, the FLR1 model not only fails to capture the correction to the Alfvén velocity, it also introduces a correction to the thermal velocity with the incorrect sign. The kinetic result can be obtained for example from (2.10) of Mikhailovskii & Smolyakov (1985), when written for cold electrons (β e = 0), and rewritten to our notation (thermal speed in that paper does not have a factor of 2, so (ρ 2 i ) (MS) = ρ 2 i /2). Coming back to the generally oblique propagation, the second approach to writing the FLR2 tensor, i.e. when the time derivative of the Hall term is not neglected and kept in the first-order tensor Π (1) with components yields the second-order tensor The resulting dispersion relation for strictly parallel propagation is slightly different

P. Hunana and others
with explicit solutions For isotropic temperatures (a p = 1) the result is equivalent to (5.115) since the Hallterm contributions for the FLR tensor disappear for all propagation directions. Also, for perpendicular propagation, the dispersion relation is naturally equal to (5.117).

FLR corrections with gyrotropic heat flux
The equation (5.99) contains the divergence of the gyrotropic heat flux, ∇ · q g , which in the index notation calculates as (5.126) Evaluation along magnetic fieldb =b 0 = (0, 0, 1) yieldsb k ∂ k → ∂ z and the matrix is evaluated as (∇ · q g ) xx = ∂ z q ⊥ + q ⊥ (∇ ·b + 2∂ xbx ); (5.127) (∇ · q g ) xy = q ⊥ (∂ ybx + ∂ xby ) = (∇ · q g ) yx ; (5.128) (5.129) (∇ · q g ) yz = (q − 2q ⊥ )∂ zby + ∂ y q ⊥ + q ⊥ ∂ ybz = (∇ · q g ) zy ; (5.130) (∇ · q g ) yy = ∂ z q ⊥ + q ⊥ (∇ ·b + 2∂ yby ); (5.131) where in the last 'zz' component we also used ∂ zbz = 0, but for brevity we left the divergence ∇ ·b = ∂ xbx + ∂ yby intact in all expressions. Note that (∇ · q g ) yy = −(∇ · q g ) xx , which is not a problem as will be described below. The pressure tensor equation is here approximated as and evaluated with respect tob =b 0 = (0, 0, 1). All components were already evaluated and the calculation is straightforward, however, we need to subtract from diagonal components the scalar pressure equations that contain the gyrotropic heat flux contributions and that evaluated alongb 0 read dp dt The requirement again is that Π zz = 0 and Π yy = −Π xx . Direct calculation yields Note that after subtraction of the perpendicular pressure equation, the xx and yy components are again anti-symmetric, i.e. (∇ · q g ) yy − ∂ z q ⊥ − 2q ⊥ ∇ ·b = −q ⊥ (∂ xbx − ∂ yby ), and the system indeed satisfies that Π yy = −Π xx . Also, considering the zz component, after subtraction of the parallel pressure equation all the terms cancel, i.e. (∇ · q g ) zz − ∂ z q − (q − 2q ⊥ )∇ ·b = 0 and the system indeed satisfies Π zz = 0. Finally, expressing the above system in the linear approximation, with the assumption that mean heat flux values are zero, i.e. q (0) = 0, q (0) ⊥ = 0, the only gyrotropic heat flux contributions that remain are ∂ x q ⊥ and ∂ y q ⊥ , in equations (5.138), (5.139). The FLR tensor is again separated to Π (1) + Π (2) and similarly to the Hall-term contributions, it is again a matter of choice if the heat flux contributions are pushed to Π (2) or kept in Π (1) . We prefer the first choice, i.e. when Π (1) is equivalent to (5.74), which yields the second-order tensor Note that we could have derived the contributions ∂ x q ⊥ and ∂ y q ⊥ in a much quicker way, if the expression (5.126) was linearized from the beginning. However, we wanted to demonstrate that if nonlinear FLR corrections are considered, the heat flux nonlinearities can significantly complicate the dynamics, even if simplified evaluation alongb 0 is performed. Right now, we cannot use these FLR corrections to obtain a dispersion relation, since we have no means to determine a closure for the gyrotropic heat flux q ⊥ . We would have to consider CGL2 fluid model or more complicated Landau fluid models, that contains evolution equations for q and q ⊥ . Additionally, the perpendicular propagating fast mode, as well as the parallel propagating ion-cyclotron and whistler modes, are not influenced by the gyrotropic heat flux. Instead, as the last step, we will consider contributions of the non-gyrotropic heat flux.

FLR corrections with non-gyrotropic heat flux vectors (FLR3)
It is possible to further increase precision of the FLR corrections by considering nongyrotropic heat flux contributions. Detailed calculations with the non-gyrotropic heat flux vectors S ⊥ and S ⊥ ⊥ are presented in appendix D. At the linear level, these nongyrotropic heat flux vectors contribute to the linearized pressure equations (5.142) and also to the linearized equations for the FLR tensor where the gradients are meant to be further linearized. Nevertheless, the above form is useful to point out that the non-gyrotropic heat fluxes are proportional to the gradients of the temperature. For clarity, partially linearized expressions are also written down by components in appendix D, see (D 24), (D 25) and (D 61), (D 62). It is important to emphasize that the expressions (5.148), (5.149) were derived for a bi-Maxwellian distribution function. As discussed later in the text, this is achieved by prescribing closures for the gyrotropic fourth-order moment in the form r = 3p 2 /ρ + r , r ⊥ = p p ⊥ /ρ + r ⊥ , and r ⊥⊥ = 2p 2 ⊥ /ρ + r ⊥⊥ . In (5.148), (5.149), we additionally neglected the perturbations r (since right now we do not want to consider models where these perturbations are evaluated from linear kinetic theory). Therefore, the FLR corrections with the non-gyrotropic heat flux, that we call here FLR3, are distribution function specific. Nevertheless, as discussed later in the text, very similar closures can be also obtained for the bi-kappa distribution function, just with additional coefficients α κ = (κ − 3/2)/(κ − 5/2), and expressions similar to (5.148), (5.149) can be derived.
The second-order vectors are also derived in appendix D, this time directly in the linear approximation, equations (D 88), (D 89) and the vectors read As a reminder, Π ij = Π ji . For clarity, we wrote the expressions above explicitly in the index notation. The non-gyrotropic heat flux vectors are perpendicular tob 0 = (0, 0, 1), and the non-zero components are for index k = x, y. The expressions can be written in a more elegant vector form, by defining vector Π z ≡ (Π xz , Π yz , Π zz = 0), vector ∇ ⊥ = (∂ x , ∂ y , 0) and matrix (5.154) which are equations (53) and (54) of Passot et al. (2012). Normalizing the equations (dropping the tilde), Fourier transforming, and writing them in the x-z plane (with ∂ y = 0) yields where the terms with Hall electric field components are specified in (5.112). At this moment we do not consider higher-order fluid models with gyrotropic heat fluxes such as CGL2 or Landau fluids (see later in the text) and to have a closed model, here we prescribe q ⊥ = 0, q = 0. It is useful to briefly explore what fluid models are obtained, if we decide to keep only the first-order vectors S (1) ⊥ , S ⊥(1) ⊥ , and ignore the secondorder corrections S (2) ⊥ , S ⊥(2) ⊥ . By considering only the S (1) ⊥ , S ⊥(1)

P. Hunana and others
either put them into to Π (1) or Π (2) . However, it can be shown that both options yield dispersion relation for the perpendicular fast mode (long-wavelength, in physical units) Rather amusingly, in comparison to the FLR2 result (5.120), the wrong sign in the correction of the thermal speed is back! 5.9. Hall-CGL-FLR3 fluid model Obviously, if matching with kinetic theory for the fast perpendicular mode is of upmost importance (for small wavenumbers), we have no other choice and we have to keep the second-order non-gyrotropic heat flux contributions S (2) ⊥ , S ⊥(2) ⊥ . We call this FLR pressure tensor as FLR3, and as stated previously, the FLR3 pressure tensor is distribution function specific. This model is very important to us and especially for the following discussion of the firehose instability. In order to be completely clear on which equations we are solving and also so that our results can be easily reproduced, let us state the entire model in all of its beauty. The linearized, normalized, Fourier transformed equations written in the x-z plane read where each component of the (FLR) pressure tensor Π, and each component of the (FLR) heat flux vectors S ⊥ , S ⊥ ⊥ , is separated to The components of the Π tensor are given by where the Hall electric field contributions for cold massless electrons read (5.163) and the components of the non-gyrotropic (FLR) heat fluxes are given by Importantly, the Π (2) expressions contain both the first-and second-order contributions from the heat flux vectors S ⊥ , S ⊥ ⊥ . To explicitly see where the closure for the gyrotropic heat fluxes q = 0, q ⊥ = 0 was performed, we scratched the terms involving these quantities. Additionally, we also scratched contributions from (S ⊥ ⊥ ) (1) x , which are here zero at the linear level.
Let us check the solution of the Hall-CGL-FLR3 fluid model for the perpendicular propagation (k = 0), which in this case reads (5.165) As a 'sanity check', since the solution stays always positive for all the wavenumbers and we do not have any instability, the model appears to be good. For small wavenumbers, the expansion yields ω 2 = k 2 ⊥ 1 + a p β − k 2 ⊥ a p β 1 8 + 5 16 a p β + · · · , (5.166) and in physical units Therefore, the Hall-CGL-FLR3 fluid model finally matches the analytic result from kinetic theory! We note that if the Π contributions are neglected in the second-order heat flux expressions (5.164), instead of the correct −5/16 FLR correction to the thermal speed, one obtains −7/16. The necessity to keep the second-order heat flux contributions to recover the kinetic result for the perpendicular fast mode was discovered by Mikhailovskii & Smolyakov (1985). Checking the solution for the parallel propagating (k ⊥ = 0) ion-cyclotron and whistler modes yields the following dispersion relation

P. Hunana and others
where the quantities v b = β (1 − a p /2) and v 2 A = 1 + (β /2)(a p − 1). We will use this dispersion relation to study the firehose instability. For a general oblique propagation direction, the dispersion relation is obviously way too large to write down, and we recommend using analytic software such as Maple or Mathematica.

Moving the Hall term to Π (1)
For completeness, just in case we want to investigate later the influence of the Hall term, moving it to Π (1) yields the following dispersion relation for the parallel propagation (5.169) 6. Parallel and oblique firehose instability 6.1. Parallel propagation For clarity, it is useful to summarize all 3 major models that describe the parallel propagating ion-cyclotron and whistler modes. By prescribing k ⊥ = 0 in the equations of the Hall-CGL-FLR3 fluid model, the model greatly simplifies. The parallel magnetic field B z = 0, the ion-acoustic mode decouples, and both the first-and second-order contributions to Π xx = 0, Π xy = 0, S ⊥ ⊥ = 0. The normalized, Fourier transformed equations written in the x-z plane read where the components of the non-gyrotropic (FLR) pressure tensor Π are given by x , (6.2) and the components of the non-gyrotropic (FLR) heat flux S ⊥ read When the entire Π is neglected, yields the Hall-CGL solution Neglecting the Π (2) contributions yields the Hall-CGL-FLR1 solution (6.5) and neglecting the heat flux S ⊥ yields the Hall-CGL-FLR2 solution Finally, the full model yields Hall-CGL-FLR3 solution (6.7) The quantities v b and v 2 A used in the solutions are defined as Importantly, the solutions are written here for ω (and not ω 2 ), and all the models of course yield 4 solutions, the other two solutions are obtained by substituting ω with −ω. For a reader who just jumped to this section, and is confused on how the solutions were split, see § 4.2 where the Hall-CGL model is discussed with final equation (4.36). The solutions can be written in various forms, and one possibility is to keep |k | in front of the square roots. The (parallel and oblique) firehose instability was studied in detail by , who focused on the Hall-CGL and Hall-CGL-FLR1 models. One of the conclusions reached in that paper was that the main reason for the relatively large discrepancy between the usual fluid models and the kinetic description is the appearance of a huge 'bump' in the imaginary phase speed, when close to the firehose threshold. The situation is demonstrated in figure 7, where kinetic solutions obtained by the WHAMP code (blue solid lines) are compared to solutions of the Hall-CGL-FLR2 model (a, blue dashed lines) and the Hall-CGL-FLR3 fluid model (b, black dashed lines). The plasma β = 4, so the long-wavelength firehose threshold, that we call the 'hard threshold', is at a p = 0.5. The temperature anisotropy is varied from a p = 0 to a p = 0.501. For solutions with the Hall-CGL and Hall-CGL-FLR1 models see figures 3 and 4 of . In the WHAMP code the value of a p = 0 cannot be used, and a p = 10 −4 was chosen instead. Also, since here we concentrate on proton dynamics, the influence of electrons in the WHAMP code was eliminated by making the electrons cold (with T e /T p = 10 −8 ). FIGURE 8. Same parameters as in figure 7, but the growth rate is plotted, and a linear scale is used for the x-axis. Notice the excellent precision of the FLR3 model for small wavenumbers up to kd i = 0.1-0.2. In comparison to kinetic theory, the fluid solutions are stabilized much more 'rapidly'. Nevertheless, the value of the maximum growth rate, and the wavenumber where the maximum growth rate is achieved, is surprisingly close to kinetic theory. This is an excellent result for a fluid model, which does not contain collisionless ion-cyclotron damping.

P. Hunana and others
Very surprisingly, figure 7 shows that the large 'bump' is reproduced by the FLR3 model, and the precision is quite good! The situation is further analysed in figure 8, where instead of the imaginary phase speed, the growth rate is plotted. The same conclusion is obtained. For the temperature anisotropy a p = 0.501, all simpler fluid models (Hall-CGL, FLR1, FLR2) are fully stable for all ranges of wavenumbers. In contrast, the FLR3 model still develops a strong firehose instability. Therefore, it is the non-gyrotropic heat flux S ⊥ of the FLR3 model that is responsible for the 'bump' (for the case of strictly parallel propagation). The appearance of the 'bump' can be understood analytically, by evaluating the fluid dispersion relations exactly at the 'hard' firehose threshold a p = 1 − 2/β , where the quantities v 2 A = 0 and v b = 1 + β /2. At the hard firehose threshold, the expression under the square root of the FLR2 solution (6.6) can be factorized as (6.9) implying that the FLR2 model is always stable for all values of k and β . Such a factorization cannot be achieved for the FLR3 model, and the solution (6.7) still becomes unstable at some range of wavenumbers, where which is the (strictly parallel) firehose instability criterion of the FLR3 model. We note that the model can become unstable also for the temperature anisotropy a p > 1. This can be perhaps considered as some remnant of the ion-cyclotron anisotropy instability, but we did not study the situation further since the instability should not be reproduced correctly. The parallel firehose instability for high plasma beta value β = 100 is shown in figure 9. The left figure is from , and it shows solutions of the Hall-CGL model (blue dashed lines) and of the Hall-CGL-FLR1 model (green dashed lines). The kinetic solutions are blue solid lines. It is shown that the FLR corrections are crucial for the correct stabilization of the firehose instability. Figure 9(b) shows refinement with FLR2 (blue dashed lines) and FLR3 tensors (black dashed lines). Obviously, for high plasma beta values, the maximum growth rate of the (strictly) parallel firehose instability is captured very precisely by the Hall-CGL-FLR3 model. The maximum growth rate (for parallel propagation) can be easily found analytically only for the Hall-CGL and Hall-CGL-FLR1 models. For models with the FLR2 and FLR3 tensors it is easier to find the maximum growth rate numerically. For example, let us consider the Hall-CGL model. Assuming the firehose-unstable regime, the frequency of the Hall-CGL solution (6.4) can be split into ω = ω r + iω i , where the imaginary part ω i = k −v 2 A − k 2 /4. Then, calculating ∂ω i /∂k = 0, yields the wavenumber k max and the maximum growth rate γ max in the following form Similarly, the Hall-CGL-FLR1 model yields . 6.2. Oblique propagation Here, we briefly investigate the parallel and oblique firehose instability for oblique propagation directions. Contour plots of the growth rate in the k − θ plane are shown in figure 10. The plasma β = 4 and the temperature anisotropy a p = 0.49, so all the fluid models are in the firehose-unstable regime. We compare four different fluid models, Hall-CGL (a), and enhancement with FLR1 (b), FLR2 (c) and FLR3 (d). The FLR2 solution shows an additional instability at higher wavenumbers that we did not study further. Note the large differences between the four solutions. Importantly, the FLR3 solution shows large enhancement of the growth rate. We use the WHAMP code for kinetic calculations, and we do not provide a contour plot for the kinetic theory. Instead, we plot the growth rate for several propagation angles, so that the comparison with the contour plots can be made easily. The figure (e) is the FLR3 model, and ( f ) is the kinetic result. It is noted that, similarly to other fluid models with higher-order FLR corrections, the FLR3 model can develop secondary instabilities at scales below the proton gyroscale.

Hellinger's contours for the Hall-CGL-FLR3 model
Here, we prescribe a fixed value for the maximum growth rate, γ max = 10 −3 ; 10 −2 ; 10 −1 , and plot solutions in the β − a p plane, which is shown in figure 11. (a,c) Show the parallel firehose instability, and (b,d) show the oblique firehose instability. Panels (a,b) are plotted with the usual logarithmic scales, and (c,d) with linear scales. Additionally, only solutions with β < 6 are shown at the bottom panels. Solutions for the Hall-CGL-FLR3 fluid model are black dashed lines, and kinetic solutions are solid blue lines (a,c) and solid red lines (b,d). The kinetic solutions were provided to us by P. Hellinger (private communication) and are identical to solutions shown in figure 1 of Hellinger et al. (2006). It is shown that the γ max = 10 −3 ; 10 −2 contours clearly lie below the hard firehose threshold, i.e. the firehose instability in the FLR3 fluid model develops at some range of wavenumbers, even if the model is stable in the long-wavelength limit. Importantly, the kinetic contours of Hellinger et al. (2006) were not calculated for cold electrons, but for isotropic electrons with β e = 1. This does not matter for the parallel propagation, since the isotropic electron pressure does not influence the dispersion relations for the parallel propagating whistler and ion-cyclotron modes. Only the effect of electron inertia will enter, however, the effect should be negligible at the scales considered here. Unfortunately, the solution for the parallel firehose instability cannot be further improved by currently developed fluid models. The oblique case is different, since even isotropic electron pressure will influence the dynamics. The contours for the oblique firehose instability should be therefore recalculated by using a proper proton-electron two fluid model. The oblique case could be further improved, by considering FLR contributions from the 100 P. Hunana and others FIGURE 11. (a,b) Solutions for the parallel (a,c) and oblique (b,d) firehose instability in the β − a p plane for a prescribed maximum growth rate γ max = 10 −3 ; 10 −2 ; 10 −1 . Solid blue and red lines are kinetic solutions from Hellinger et al. (2006). Black dashed lines are solutions of the Hall-CGL-FLR3 model. The magenta line is the 'hard' (long-wavelength limit) firehose threshold. (c,d) Same as top panels, but with linear scales for both axes, and only showing results for β < 6. non-gyrotropic heat flux tensor σ that were neglected here (contributions for parallel propagation are zero). For the oblique case, one can also use higher-order fluid models with the gyrotropic heat flux fluctuations q , q ⊥ (such as the CGL2 model discussed later), or even Landau fluid models. The main reason for the relatively large discrepancies found in the contours for the parallel firehose instability is demonstrated in figure 12.

Heat flux tensor equation
Multiply the Vlasov equation by mc i c j c k and integrate over the velocity space. Naturally, the other possibilities are to multiply by mc i c j v k , mc i v j v k or mv i v j v k , although none of them is more revealing and we will use the first choice. It is convenient to define the symmetric operator S that acts on a tensor of third rank according to  when there is no instability in the long-wavelength limit. Nevertheless, by increasing the value of a p beyond 0.5, the growth rate in the FLR3 model quickly falls off, and the instability disappears for a p > 0.5317. In contrast, kinetic theory (blue solid lines) develops a strong firehose instability even at a p = 0.6.
i.e. the operator represents all possible cyclic permutations. For the first term we will need the following identities (7.3) and the entire first term of the integrated Vlasov equation calculates as For the second term we will need identities m cccvf d 3 v = r + qu;  (7.8) and the second term calculates as Note again that the divergence of a tensor operates through its first component Similarly, if the operator acts on a tensor from the right-hand side, the most natural way is to define that it operates through the last component, i.e. [q · ∇] ij = q ijl ∂ l . In the expressions above [q · ∇u] ijk = q ijl ∂ l u k and [pu · ∇u] ijk = (pu) ijl ∂ l u k = p ij u l ∂ l u k . The last term can also be rewritten as p ij (u · ∇u) k , where the expression (u · ∇u) k = u l ∂ l u k is familiar from MHD. For the third term we will need identities ∂c i ∂v l = δ il ; (7.10) ∂ ∂v l (c i c j c k ) = δ il c j c k + c i δ jl c k + c i c j δ kl , (7.11) and the entire third term calculates as For the fourth term we will need identities (7.16) and the entire fourth term calculates as Combining all the results together 1 + 2 + 3 + 4 = 0, the entire heat flux tensor equation obtained by direct integration of the Vlasov equation reads Now we will need to use 3 different momentum equations that will cancel various terms. We will also need identity We need to multiply the momentum equation for ∂u i /∂t by p jk , the equation for ∂u j /∂t by p ki and the equation for ∂u k /∂t by p ij . All the three momentum equations that we want to subtract from the heat flux equation can be written together as (7.20) and because of identity (7.19), this is equivalent to Note that, because of the symmetric operator, it does not matter if the tensor p is applied to the momentum equation from the left or right, since all the expressions are symmetric in this regard, i.e. for example [(∇ · p)p] S = [p(∇ · p)] S . By subtracting (7.21) from (7.18), the final heat flux tensor equation reads

P. Hunana and others
By defining the cyclotron frequency vector Ω = qB/mc, this equation identifies with equation (A 5) in Chust & Belmont (2006). By using the scalar cyclotron frequency defined with respect to |B| as Ω = q|B|/mc, the heat flux tensor equation is written explicitly in the index notation ∂ ∂t q ijk + ∂ l (r lijk + u l q ijk ) + q ijl ∂ l u k + q jkl ∂ l u i + q kil ∂ l u j + Ωb l ( ilm q mjk + jlm q mki + klm q mij ) − 1 mn (p ij ∂ l p lk + p jk ∂ l p li + p ki ∂ l p lj ) = 0, (7.23) which is equivalent to (4) in Goswami et al. (2005) with the term (b × q) S ijk = ilmbl q mjk + jlmbl q mki + klmbl q mij = −b l ( iml q jkm + jml q ikm + kml q ijm ).
(7.24) 7.1. Heat flux tensor decomposition Considering only the gyrotropic part of the heat flux tensor (3 × 3 × 3) cube, the heat flux is decomposed into the scalar parallel and perpendicular heat flux components q , q ⊥ according to (7.25) and which in the index notation reads This part of the pressure tensor represents only the gyrotropic part, and the full heat flux decomposition can be written as q = q g + q ng . In the heat flux tensor equation, the term B × q is proportional to the cyclotron frequency Ω = qB 0 /mc and the equation rewrites as The situation is now similar to the previously studied pressure tensor. At long spatial scales (low frequencies ω), this term will dominate and the gyrotropic contribution must be equal to zero At first look, it is not that obvious that the decomposition (7.25) satisfies this equation. It is however possible to verify that indeed = q ⊥br ( irjbk + irkbj ); (7.29) Putting all terms together (7.32) and we see that all three parts of the symmetric operator are required to make this term equal to zero. Therefore, the exact heat flux tensor equation reads Since this should be an introductory text, we do not want to be bothered right now with the complicated algebra of the non-gyrotropic heat flux q ng . For the clarity of the presented material, here we separate the non-gyrotropic heat flux into a separate term Q ng , and write the heat flux tensor equation in the following form We will address the non-gyrotropic heat flux contributions in appendix D.
In a similar fashion to the pressure decomposition (2.35), we are going to frequently apply double contractions withbb and (I −bb)/2. Applying these operators to the heat flux tensor q ijk , must yield quantities that are vectors. It is therefore logical to define parallel and perpendicular heat flux vectors S ≡ q :bb; S ⊥ ≡ q : (I −bb)/2. (7.36) The scalar parallel and perpendicular heat flux components q , q ⊥ (which are the only parts that are gyrotropic) are obtained by further projecting these heat flux vectors along the magnetic field lines, i.e. by performing ·b, so that q = (q :bb) ·b; q ⊥ = (q : (I −bb)/2) ·b. (7.37) Briefly considering only the gyrotropic heat flux (7.26) on the right-hand side, it is easy to verify that the parallel decomposition indeed works For the perpendicular decomposition, it is useful to specifically calculate components Tr q g = q g : I = q b + 2q ⊥b ; (7.42) q g : (I −bb)/2 = (q g : I − q g :bb)/2 = (q b + 2q ⊥b − q b )/2 = q ⊥b ; (7.43) (7.44) which verifies that the perpendicular decomposition (7.37) is satisfied for q g . Now we apply the decomposition (7.37) directly to the definition of the entire heat flux tensor (2.6), which yields (7.46) and the perpendicular decomposition The decomposition (7.37) obviously works for the entire heat flux q, as well as for the gyrotropic part q g . Similarly to the pressure decomposition, this further yields the required properties that the non-gyrotropic heat flux q ng must satisfy. By using (7.37), the entire heat flux decomposition reads Obviously, it would be useful to introduce a triple-contraction operator and instead of (q :bb) ·b to write something like q . . .bbb, but we do not want to introduce new notations. Also, an alternative and perhaps prettier expression is to move the ·b to the left-hand side, as (q :bb) ·b =b · q :bb, but we will keep the first choice. Applying :bb and ·b to the above equation yields (q :bb) ·b = (q :bb) ·b + (q ng :bb) ·b, (7.52) implying the first requirement for the non-gyrotropic heat flux (q ng :bb) ·b = 0. (7.53) Similarly, the second requirement is obtained by applying : (I −bb) and ·b, yielding (q ng : (I −bb)) ·b = 0. (7.54) By using the first requirement, the second requirement simplifies to (q ng : I) ·b = Tr q ng ·b = 0, (7.55) where obviously the trace and ·b operators commute. The two requirements in the index notation read q ng ijkbibjbk = 0; q ng iikbk = 0. (7.56) Instead of decomposing q = q g + q ng , an alternative and very useful decomposition of the entire heat flux tensor reads q = S + σ , (7.57) with the requirement σ :bb = 0 and σ : (I −bb) = 0 (or equivalently σ : I = 0). The heat flux vectors (7.36) therefore satisfy S = S :bb and S ⊥ = S : (I −bb)/2. The heat flux vectors contain both gyrotropic and non-gyrotropic contributions. Since the gyrotropic contributions are obtained by projecting these vectors along the magnetic field lines q = S ·b, q ⊥ = S ⊥ ·b, it is useful to introduce the following decomposition The vectors S ⊥ , S ⊥ ⊥ are referred to as the non-gyrotropic heat flux vectors, and their algebra is addressed in appendix D. Here, we only state that the q ng can be decomposed into vectors S ⊥ , S ⊥ ⊥ and tensor σ according to The entire heat flux tensor σ is of course non-gyrotropic. Now we need to verify the heat flux contributions (2.72), (2.73) that we used in the pressure equations (2.74), (2.75). To calculate the heat flux contributions to the pressure equations, we will need (7.60) (7.62) and the contributions arê

P. Hunana and others
The left-hand side of the above equation can be also written naturally as (∇ · q g ) : (I − bb)/2, which is consistent with an alternative derivation of the perpendicular pressure equation, where instead of doing trace and subtracting the parallel pressure equation, one can directly perform : (I −bb)/2.

Parallel heat flux equation
Now we apply the decomposition (7.37) to the heat flux tensor equation (7.34) step by step. We consider only gyrotropic heat flux components. Starting with the equation for the parallel heat flux, the first term calculates as the second term calculates as Note that in all of the expressions here, the operator ·b can be naturally moved to the left asb·. The third term calculates as For the final fourth term we will need (7.70) and the fourth term calculates as Combining all the terms yields the equation for the scalar parallel heat flux where Q ng ≡ (Q ng :bb) ·b. The heat flux equations contain the fourth-order moment r ijkl , which we will consider in the next section and at this stage it is unspecified.

Perpendicular heat flux equation
We will apply the trace operator to (7.34) and also perform ·b step by step. In the index notation, we are basically multiplying the entire equation by δ ijbk . The first term calculates The second term ∇ · (r + uq g ) calculates as The components r liik = r iikl are nothing else but the trace of the moment r [Tr r] kl = δ ij r ijkl = r iikl . (7.81) The third term in (7.34) calculates as where we have used that q The fourth term in (7.34) calculates as [p(∇ · p)] S ijk = p ij ∂ l p lk + p jk ∂ l p li + p ki ∂ l p lj ; (7.85) Tr[p(∇ · p)] S ijk = p ii ∂ l p lk + p ik ∂ l p li + p ki =p ik ∂ l p li = p ii ∂ l p lk + 2p ik ∂ l p li ; (7.86) Collisionless fluid models. Part 1 111 and we have to calculate each term separately. Since p ii = p + 2p ⊥ and by using the already calculated equation (7.70) forb k ∂ l p lk and identities p ikbk = p b i + Π ikbk and Π ikbibk = 0 one obtains (7.89) and the final result of (7.87) is Collecting all the terms, one obtains (7.91) and subtracting the parallel heat flux equation (7.74) and dividing by two yields the perpendicular heat flux equation where Q ng ⊥ ≡ (Q ng : (I −bb)/2) ·b. The term containing r in the above equation can be rewritten to many possible forms, for example :bb] ·b =b · (∇ · r) : (I −bb)/2.

P. Hunana and others
7.4. Scalar heat flux equations continued The parallel and perpendicular heat flux equations (7.74), (7.92) contain expressions for the fourth-order moment r, and even though we will consider this moment in detail in the next section, here we want to finish the derivation of the scalar heat flux equations, and we write down the required expressions. Similarly to the pressure tensor and the heat flux tensor, the fourth-order moment can be decomposed into its gyrotropic and non-gyrotropic parts, r = r g + r ng . The gyrotropic part r g is decomposed according to (8.10), and it can be shown (see later in the text), that direct calculation yieldsb which yields the parallel heat flux equation (7.96) and the perpendicular heat flux equation Exact nonlinear expressions for Q ng and Q ng ⊥ that represent contributions from the non-gyrotropic heat flux q ng are calculated in appendix D, see (D 110), (D 120). The scalar heat flux equations (7.96) and (7.97) are completely general at this stage, since no distribution function has been prescribed yet, and no simplification has been introduced. These equations are exact.
Nevertheless, equations (7.96), (7.97) are very complicated, and at this stage, it is beneficial to simplify them. One possibility is to cancel all the non-gyrotropic contributions Π, q ng and r ng , and we will study such fluid models later. Another possibility is to keep only those non-gyrotropic terms that have some non-zero contribution at the linear level. This eliminates the last term in the second line of (7.97) that is proportional to [. . .] · Π ·b. Also, by following the derivations in appendix D, it is easy to show that the terms Q ng , Q ng ⊥ (that represent the non-gyrotropic heat flux q ng ), do not contribute at the linear level. Therefore, the heat flux equations simplify

P. Hunana and others
terms entering the heat flux equations directly calculate as (7.106) The use of these expressions in (7.102), (7.103) cancels various terms, and the scalar heat flux equations read (7.108)

Fourth-order fluid moment
The algebra of the fourth-order fluid moment r = m ccccf d 3 v can be quite intimidating at first, since the moment is a 4-D cube (3x3x3x3). We are not going to derive the time evolution equation for this moment step by step, nevertheless, later in the text we derive the evolution equation of the nth-order fluid moment X (n) , see (12.16), and therefore, for n = 4, the equation reads The symmetric operator 'S' is here defined as Here, we are not that interested in the evolution equation for r, we just want to clarify its decomposition, that we need in the heat flux equations. The definition of gyrotropy means that the integral has to be evaluated only over combinations of (v − u ⊥ ) ≡ c and |v ⊥ − u ⊥ | 2 ≡ c 2 ⊥ . For the fourth moment r, there are obviously only 3 possibilities: c 4 , c 2 c 2 ⊥ and (c 2 ⊥ ) 2 , and these gyrotropic components will be called r , r ⊥ and r ⊥⊥ . We have already seen that the double contractions withbb and (I −bb)/2, were very useful operators to extract the gyrotropic components for the lower-order moments, p and q. To obtain any scalar quantity from the fourth moment, we obviously need to apply two double contractions. There are 3 possibilities: we can applybb twice, we can applybb and (I −bb)/2 or we can apply (I −bb)/2 twice. Not surprisingly, these double contractions with r indeed extract the 3 possible gyrotropic parts, as it is easy to verify that We are now ready to guess how to write the decomposition of the gyrotropic fourth-order moment. Motivated by the previous decompositions, it obviously has to be something in the form of where we still did not determine how the symmetric operator acts here. Importantly, the symmetric operator 'Sym' is not equivalent to the symmetric operator 'S' that cycles all the indices around, equation (8.2). The fluid hierarchy obviously needs two symmetric operators, one unique 'S' that is used to derive the evolution equation of a given fluid moment X (n) , and one non-unique 'Sym' that is used for the decomposition of that fluid moment. The determination of how 'Sym' acts here is not that obvious. Nevertheless, one can consider in how many ways one can extract the gyrotropic components from r ijkl . For r , one does two double contractions with (bb). The possible choices are (bb) ij (bb) kl ; (bb) ik (bb) jl and (bb) il (bb) jk , however, all of these choices are equivalent. To obtain the r ⊥ , we perform double contractions with (bb) and (I −bb)/2, where the last operator contains a function δ ij . How many different delta functions we can obtain from 4 indices i, j, k, l? There are 4 2 = 4!/2!(4 − 2)! = 6 different possibilities: and all of them are non-equivalent. The symmetric operator acting in the second term therefore has 6 components Lastly, to obtain the r ⊥⊥ component, one applies two double contractions with (I −bb)/2 that results in combinations of δ ij δ kl . How many possibilities do we have? Obviously, there are 6 possibilities for the first delta function, which, by pairing with the other delta function in a way that indices are not repeated, yields together 6 possibilities δ ij δ kl ; δ ik δ jl ; δ il δ jk ; δ jk δ il ; δ jl δ ik ; δ kl δ ij . However, 3 possibilities have an equivalent pair, δ ij δ kl = δ kl δ ij ; δ ik δ jl = δ jl δ ik and δ il δ jk = δ jk δ il and there are only 3 116

P. Hunana and others
non-equivalent combinations. The symmetric operator acting on the last term in (8.6) can be determined to be

[(I −bb)(I −bb)]
Sym ijkl The factor 1/2 is actually not that obvious, and one needs to verify that the decomposition indeed satisfies (8.5). The entire fourth-order gyrotropic moment r g is decomposed in the index notation according to It might be tempting to rearrange this decomposition to a more compact form nevertheless, in actual calculations we find the form (8.10) to be more useful. It is important to verity that the decomposition (8.10) really works. A straightforward calculation yields which is consistent with (8.3). The r ⊥ component calculates as (8.17) which is consistent with (8.4). The r ⊥⊥ component calculates as  (8.22) which is consistent with (8.5). The decomposition (8.10) indeed works. Now we need to verify expressions (7.94), (7.95) that were used in the scalar heat flux equations. The first expression calculates as (8.26) and the second expression calculates similarly as which verifies (7.94), (7.95). Furthermore, by exploring the ∂r/∂t equation (8.1), at frequencies that are much smaller than the gyrofrequency, the gyrotropic part of r should satisfy (8.32) in the same way that the gyrotropic parts of p and q satisfied this requirement. By using the gyrotropic (8.10), it is easy to calculate for example that and by adding together all the 4 representations of the 'S' operator, one can indeed verify that all the terms cancel, yielding (8.32). By decomposing the entire fourthorder moment into its gyrotropic and non-gyrotropic parts r = r g + r ng , (8.34) the evolution equation (8.1) can therefore be rewritten as

P. Hunana and others
8.1. Non-gyrotropic r ng Basic properties of the non-gyrotropic tensor r ng can be easily determined with a similar procedure as those we used for the non-gyrotropic pressure Π and the nongyrotropic heat flux q ng . The decomposition of the entire fourth-order moment is (8.36) and the meaning of the 'Sym' operators was specified by (8.8), (8.9). By using definitions (8.3), (8.4), (8.5) for r , r ⊥ , r ⊥⊥ in the expression above, the full decomposition reads The last two properties can be also written as Tr r ng :bb = 0, Tr Tr r ng = 0. Finally, writing all 3 properties in the index notation One would assume that the non-gyrotropic r ng can be evaluated by rewriting (8.35) as (8.42) and then expanding the right-hand side similarly to what we did for the non-gyrotropic Π (and also the non-gyrotropic heat flux vectors in appendix D). However, the equation does not seem to yield anything useful. Instead, one needs to consider a specific example of a bi-Maxwellian distribution function.

Bi-Maxwellian distribution
Here we consider a special case of a bi-Maxwellian distribution function It is useful to use the fluctuating velocity c = v − u, where c 2 = c 2 + c 2 ⊥ and c 2 ⊥ = c 2 x + c 2 y . For brevity of the calculations we introduce where here the thermal speeds are meant to be spatially dependent, i.e. they are written with T and not with T (0) . The bi-Maxwellian distribution therefore reads It is useful to remind ourselves of the following one-dimensional integrals (8.46) and that the general integral with x n , where n is an integer, reads The double factorial (n − 1)!! = 1 · 3 · 5 · · · (n − 1), with (−1)!! = 1. The identity (8.47) is easily obtained by performing differentiation ∂/∂α of the first result in (8.46), and this technique is useful in the calculation of the fluid moments.
To become more familiar with the bi-Maxwellian distribution, it is useful to verify whether integrals over f 0 indeed yield the expected fluid moments. Since the macroscopic velocity u is independent of v, changing from variable v to c just yields d 3 v = d 3 c (similarly to the substitution y = x + const. that yields dy = dx). We use the notations d 3 c = dc d 2 c ⊥ and d 2 c ⊥ = dc x dc y . Also, because the integrals are evaluated from −∞ to ∞, substitution from v to c does not change these bounds. Integration in velocity space therefore yields For simplicity we have integrated in Cartesian coordinates (instead of perhaps the more elegant cylindrical coordinates), which has an additional benefit in that we can stop writing the integral bounds, since the integrals are always from −∞ to ∞. By using the bi-Maxwellian f 0 it is very easy to verify that

P. Hunana and others
Continuing with the heat flux components yields It is important to emphasize that the zero heat flux values were obtained by assuming some prescribed equilibrium f 0 , here bi-Maxwellian. A crucial principle used in constructing fluid models with higher-order moments is that the specific f 0 is assumed only for the last retained moment. Indeed, if one closes the hierarchy by prescribing q = 0, q ⊥ = 0, the CGL fluid model is obtained. However, to go higher in the fluid hierarchy, evolution equations for q and q ⊥ cannot be eliminated, even if the hierarchy will be eventually closed, for example at the fourth-order moment level, by assuming a bi-Maxwellian f 0 . Or in another words, the zero heat flux values were obtained for a distribution function that is strictly f 0 , however, fluctuations/perturbations around f 0 are still allowed. Therefore, to obtain the next available model after the CGL, the heat flux equations must be retained, and a closure is performed on the fourth-order moment. Only this last moment is calculated by assuming a specific distribution function, that is strictly f 0 . For a bi-Maxwellian, the gyrotropic fourth-order moments are calculated according to Therefore, for a bi-Maxwellian, the following closure can be constructed 56) and this closure is known as the 'normal' closure, a name suggested by Chust & Belmont (2006).

8.3.
Bi-Maxwellian non-gyrotropic r ng contributions Let us consider a specific case of a bi-Maxwellian distribution function. Then, an expansion procedure analogous to the one developed by Grad (1949) for rarefied gases can be used, where the distribution function is expanded in a series of Hermite polynomials. By keeping only the first few terms, one can express the fourth-order moment through the second-order (pressure) moments. The procedure is written down in Oraevskii et al. (1968), and the entire fourth-order moment can be decomposed in the following way: see (2.31) in Oraevskii et al. (1968) (where a small typo on the left-hand side is present, where instead of q αβγ , there should be q αβγ δ ), equation (24) in Goswami et al. (2005), equation (9) in Passot & Sulem (2007) where the non-gyrotropic contributions read (8.58) and where terms Π ij Π lk + Π ik Π jl + Π il Π jk were neglected. By using p g ij = p b ibj + p ⊥ (δ ij −b ibj ), it can be indeed shown by direct straightforward calculation that the gyrotropic part (8.59) This agrees with the decomposition (8.10) valid for a general distribution function, after one specifies that for a bi-Maxwellian distribution r = 3p 2 /ρ, r ⊥ = p p ⊥ /ρ and r ⊥⊥ = 2p 2 ⊥ /ρ (last term in (8.10) contains r ⊥⊥ /2). Now it is possible to calculate the non-gyrotropic contributions r ng , by using (8.58). It is useful to pre-calculate the trace of (8.58) that reads (8.60) The two terms that enter the parallel and perpendicular heat flux equations (7.96),  (8.64) so the term entering the q ⊥ equation (7.97) reads Note that (b · ∇Π) :bb = −Π : (b(b · ∇)b) S , and the term does not contribute at the linear level. Of course, many nonlinear terms were neglected, and technically, the only fully consistent procedure is to explicitly write those terms at the linear level, in the formb The results show that, if the Π contributions are kept in the heat flux equations, the r ng contributions cannot just be straightforwardly neglected, since the Π and r ng contributions partially cancel out.

Bi-Maxwellian fluid model, second-order CGL (CGL2)
In the previous section we have seen that for a bi-Maxwellian distribution function, one can close the fluid hierarchy by prescribing the 'normal' closure It is useful to summarize the fully nonlinear model that we want to consider here. We consider only proton species, and we make the electrons massless and cold. We simplify the pressure equations, by neglecting the FLR stress forces (which enter only nonlinearly), and we keep only those non-gyrotropic contributions in the pressure and heat flux equations that have a non-zero linear contributions. The nonlinear model reads (9.4) ∂p ∂t + u · ∇p + p ∇ · u + 2p b · ∇u ·b + ∇ · (q b ) − 2q ⊥ ∇ ·b + ∇ · S ⊥ = 0; (9.5) ∂p ⊥ ∂t + u · ∇p ⊥ + 2p ⊥ ∇ · u − p ⊥b · ∇u ·b + ∇ · (q ⊥b ) + q ⊥ ∇ ·b + ∇ · S ⊥ ⊥ = 0; (9.6) ∂q ∂t To develop a vocabulary, let us first consider the above fluid model with all the non-gyrotropic fluctuations neglected, i.e. when all the terms with Π, S ⊥ , S ⊥ ⊥ are neglected. Additionally, let us neglect the Hall term in the induction equation. Such a fluid model is non-dispersive, and represents a generalization of the non-dispersive CGL model. Even though such a fluid model was not explicitly considered by Chew, Goldberger and Low, the pressure equations that include the gyrotropic heat flux contributions q , q ⊥ , were indeed given by Chew et al. (1956). As discussed previously, the authors just used a different notation with q n = q − 3q ⊥ and q s = q ⊥ , but their pressure equations contain the gyrotropic heat flux. The authors did not consider evolution equations for the gyrotropic heat fluxes, equations (9.7), (9.8). Nevertheless, the name CGL is so well recognized by the community, i.e. the name CGL is essentially recognized as the collisionless MHD, that we suggest calling the non-dispersive fluid model that used the above gyrotropic heat flux equations, 'the second-order CGL', abbreviated as 'CGL2'. The name CGL2 has a nice advantage in that the abbreviations for many models considered previously can be easily and naturally generalized. The CGL2 fluid model does not contain any dispersive effects and it is length-scale invariant, similarly to the CGL and MHD models. If the Hall term in the induction equation is considered, this yields the 'Hall-CGL2' model. Considering also the first-order FLR corrections to the pressure tensor (FLR1) or the second-order corrections (FLR2), yields the 'Hall-CGL2-FLR1' and the 'Hall-CGL2-FLR2' fluid models. Finally, considering the non-gyrotropic heat flux fluctuations S ⊥ , S ⊥ ⊥ , yields the Hall-CGL2-FLR3 model. Fluid models with the Hall term neglected can be easily abbreviated as 'CGL2-FLR1', 'CGL2-FLR2' etc. Even though, if the dispersive effects of the FLR corrections are considered, there is really no reason to neglect the simple Hall term, and the use of such fluid models is discouraged. Obviously, the name 'CGL2' is extremely beneficial and the very natural for the classification of fluid models, and we will use this name henceforth.
It is important to correctly normalize the heat flux, and the normalization is according to q ,⊥ = q ,⊥ /(p (0) V A ), so both the parallel and perpendicular heat fluxes are normalized in the same way (similarly, both p and p ⊥ are normalized with respect to p (0) ). By dropping the tilde, the scalar pressure equations (9.5), (9.6) remain unchanged, and the normalized nonlinear heat flux equations read Now we need to linearize the system, where one prescribes q (0) = 0 and q (0) ⊥ = 0. Once normalized, linearized, transformed to Fourier space and written in the x-z plane, the model reads (dropping the tilde everywhere) −ωρ + k ⊥ u x + k u z = 0; (9.11) 124 P. Hunana and others As previously, v 2 A = 1 + (β /2)(a p − 1).

CGL2 dispersion relation
Neglecting all the dispersive effects, prescribing Π = 0, S ⊥ = 0, S ⊥ ⊥ = 0, and eliminating the Hall term in the induction equation (i.e. terms proportional to ik 2 in the above induction equation), yields the CGL2 model. By exploring the system, it is obvious that the Alfvén mode separates from the system in the u y , B y components. The dispersion relation of the Alfvén mode in the (non-dispersive) CGL2 model is therefore the same as in the CGL model, and reads (the normalization tildes are dropped) The 'hard' threshold of the oblique firehose instability is therefore the same as in the CGL model, and correctly reproduced. The entire dispersion relation for the CGL2 fluid model can be written in the following form

25)
A 2 = k 4 β 2 k 2 9 4 v 2 A + 3 8 β + k 2 ⊥ 9 4 1 + a p β − 5 9 a 2 p β ; (9.26) The CGL2 fluid model therefore contains 5 forward and 5 backward propagating modes. The oblique Alfvén mode is separated as stated above, and is incompressible. The remaining 4 forward and 4 backward waves are generally compressible and coupled together through the eighth-order polynomial in ω (fourth-order polynomial in ω 2 ). Let us ignore for a moment the distinction between the forward and backward modes, since these will always be symmetric in ω. The CGL2 model contains 5 waves/modes. The oblique Alfvén mode is separated, and the remaining 4 compressible modes are in general strongly coupled. In contrast to the simpler CGL (and MHD) model, that contains only 2 compressible modes (slow and fast), we therefore have 2 'new' modes, that do not have an analogy in the usual CGL and MHD descriptions. There is no standardized vocabulary on how the modes should be named, perhaps the expression 'higher-order modes' is the most appropriate. For the highly oblique propagation, considering the mirror instability, sometimes the expression 'mirror modes' is used. Nevertheless, since none of the 4 compressible modes match the slow and fast CGL dispersion relations, the distinction between modes becomes blurry. The reader has to get used to the fact that, by considering higher-order fluid moments, we perhaps came closer to the kinetic theory, which admits an infinite number of modes that are difficult to classify, unless a specific situation is considered. The difference between the CGL2 model and the kinetic description is that, unless a firehose or mirror instability threshold is reached, the fluid modes are un-damped. For strictly parallel propagation (k ⊥ = 0), in addition to the parallel Alfvén mode (9.22) ω = ±k v A that was already separated, the solution of the eighth-order polynomial contains another Alfvén mode ω = ±k v A , which was of course expected and is equivalent to the CGL model, since the components u x , u y , B x , B y de-couple for parallel propagation. The remaining solutions are (9.30) Solutions (9.28) and (9.29) can be obtained by considering a pure 1-D geometry (where p ⊥ and q ⊥ disappear since the kinetic velocity v ⊥ disappears) and considering only fluctuations in the quantities ρ, u z , p , q . The q heat flux fluctuations in the CGL2 fluid model are therefore responsible for 'splitting' of the CGL ion-acoustic mode ω = ±k 3β /2 into two modes (9.28) and (9.29).
Considering solutions for strictly perpendicular propagation (k = 0), the only nonzero mode is the fast mode (9.31) and the dispersion relation is unchanged from the CGL model (for perpendicular propagation the q , q ⊥ contributions to the equations for p , p ⊥ naturally vanish and the CGL2 model is equivalent to the CGL model).
9.2. Mirror instability A very interesting direction of propagation for this fluid model is the highly oblique limit, k ⊥ k , since the CGL2 fluid model contains the correct mirror instability threshold. The correct mirror threshold can be already seen in the quantity A 0 , where, in contrast to the usual CGL model, the factor 1/6 disappears in the last term ∼ a 2 p β . In the highly oblique limit, the fast mode (9.31) can be separated from the coupled eighth degree polynomial (actually fourth degree polynomial in ω 2 ) in the following 126 P. Hunana and others way. Importantly, one cannot prescribe a completely perpendicular propagation k = 0, and a highly oblique limit k ⊥ k , or k → 0, has to be considered. The polynomial coefficients can be approximated as A 6 = k 2 ⊥ (1 + a p β); (9.32) A 4 = k 2 β k 2 ⊥ 7 2 1 + a p β − 1 7 a 2 p β ; (9.33) A 2 = k 4 β 2 k 2 ⊥ 9 4 1 + a p β − 5 9 a 2 p β ; (9.34) An alternative approximation is to use k ⊥ → k in the above expressions. The coefficient A 6 does not contain any k and is very large compared to A 4 , A 2 , A 0 that all contain k . The solutions will inevitably be one fast mode, and 3 very slow modes. In this specific case, the fast mode can be quickly separated from the eighth-order polynomial in (9.23) with a neat trick (9.36) The trick can be verified by multiplying both brackets term by term, dividing by −A 6 and observing that (A 4 , A 2 , A 0 )/A 6 is negligible. The fast mode is always obtained correctly in this way (when the A 6 is much larger than all the other coefficients), however, the dispersion relation for the 3 slow modes is correct only in this specific case when the system is non-dispersive, and all 3 slow modes have constant phase speeds. A proper limit and expansion should be always checked. Nevertheless, in the highly oblique limit, the fast mode ω = ±k ⊥ 1 + a p β is separated, and the rest of the dispersion is of sixth order in ω and contains 3 'slow' modes. By cancelling out the k 2 ⊥ that is common in all of the remaining coefficients (or k 2 ), an obvious substitution is offered that readsω = ω/(k β ), and that yields a dispersion equation of third order inω 2 in the form −(1 + a p β )(ω 2 ) 3 + 7 2 1 + a p β − 1 7 a 2 p β (ω 2 ) 2 − 9 4 1 + a p β − 5 9 a 2 p β (ω 2 ) Now, at exactly the mirror threshold a 2 p β = 1 + a p β , the last term disappears and we have one solutionω 2 = 0, and two solutions (9.38) that are both positive. Slightly below or beyond the mirror threshold, we can prescribe a 2 p β = (1 + a p β )(1 + ), where is a small parameter, meaning that for < 0, the system is slightly below the expected mirror threshold, i.e. stable. For > 0, the system is slightly beyond the expected mirror threshold, i.e. unstable. The polynomial transforms to It is possible to show that this polynomial has a discriminant 12 > 0, implying the existence of 3 distinct real roots. Now, when is small, the two solutions (9.38) will remain almost the same and change only slightly, importantly, the solutions will remain positive. Solutions of a cubic polynomial (x − x 1 )(x − x 2 )(x − x 3 ) = 0 form the last coefficient in that polynomial, −x 1 x 2 x 3 , here equal to 3 8 , and here x 1 and x 2 are positive. Therefore, < 0 (slightly below the threshold) implies x 3 > 0; and > 0 (slightly beyond the threshold) implies x 3 < 0. The negative solution (in (ω 2 )) represents the mirror instability, which finishes the analytic proof that the (highly oblique) mirror threshold in CGL2 model is indeed a 2 p β = 1 + a p β . The CGL2 fluid model therefore contains the same mirror instability threshold, as is found in kinetic theory.
The conclusion can be double checked by numerically exploring solutions of (9.39) for several values of , and it is indeed possible to conclude that, for < 0 (slightly below the threshold), the signs of the solutions areω 2 = +, +, +; and that, for > 0 (slightly beyond the mirror threshold), the signs of the solutions areω 2 = +, +, −. The solution with the minus sign represents the mirror instability. It is noted that the Landau fluid description allows for the construction of a fluid model with quasi-static heat fluxes, that also correctly reproduces the 'hard' mirror instability threshold . Considering Landau fluid models, the mirror instability was numerically investigated in detail for example by Passot & Sulem (2007), Passot et al. (2012) and Sulem & Passot (2015). As a suggestion for possible future work, it would be beneficial to numerically calculate 'Hellinger's' contours (see § 6.3) with the prescribed maximum growth rate for the mirror instability. If one uses Landau fluid models with sufficiently precise FLR corrections, the results should match the kinetic contours very precisely.

Hall-CGL2 dispersion relation
Considering the Hall term, the dispersion relation of the Hall-CGL2 fluid model reads = k 2 k 2 ω 2 ω 6 − ω 4 β 7 2 k 2 + a p k 2 ⊥ + ω 2 β 2 k 2 9 4 k 2 + 3a p k 2 ⊥ − k 4 β 3 3 8 k 2 + a p k 2 ⊥ , (9.40) where the left-hand side represents the CGL2 dispersion relation (9.23) and the righthand side is the Hall-term contribution. For strictly parallel propagation, the solutions are the usual Hall-CGL whistler and ion-cyclotron waves ω = ±k 2 /2 + k v 2 A + k 2 /4, ω = ±k 2 /2 − k v 2 A + k 2 /4, that are accompanied by the CGL2 solutions (9.28), (9.30) and (9.29). For strictly perpendicular propagation, the solution is (9.31) since the Hall term vanishes. The dispersion relation of the Hall-CGL-FLR1 model was already quite large to write down, and we do not provide the final dispersion relation for the Hall-CGL2-FLR1 model, instead we recommend working with the full system of linearized equations and solving it numerically or using analytic software such as Maple.
12 Discriminant of a cubic polynomial ax 3 https://www.cambridge.org/core/terms. https://doi.org/10.1017/S0022377819000801 128 P. Hunana and others 9.4. 'Static' closure -generalized isothermal closure There exists even simpler fluid closure that recovers the correct mirror threshold, called the 'static' or 'quasi-static' closure (Constantinescu 2002;Chust & Belmont 2006;Passot et al. 2006). By using heat flux equations (9.7), (9.8) and considering the quasi-static regime with ∂q /∂t = 0, ∂q ⊥ /∂t = 0, (9.41) where the FLR pressure tensor Π was neglected in the second equation. It is useful to define ∂ ≡b · ∇, which represents the gradient along the magnetic field lines, and to use the identity (9.42) further yielding (9.43) Equation (9.43) describes the evolution of the parallel and perpendicular temperature fluctuations along the magnetic field lines, and it can be verified that the solution is where a p = T (0) ⊥ /T (0) . Solution (9.44) is referred to as the 'static' closure and it is, for example, equivalent to (19), (20) of Passot et al. (2006), and (2) of Constantinescu (2002). For isotropic mean temperatures a p = 1, equation (9.44) yields T ⊥ = T (0) ⊥ = T (0) and both parallel and perpendicular temperatures are isothermal. A similar result is obtained for general a p with |B| = B 0 , which yields T ⊥ = T (0) ⊥ . The closure (9.44) can be therefore viewed as a generalization of the isothermal closure in the presence of temperature anisotropy and variations of magnetic field strength. The closure is used directly in the momentum equation (9.3), and the time-dependent pressure equations (9.5), (9.6) are disregarded. Therefore, even though the closure was derived by considering fourth-order moments, the resulting fluid model is actually simpler than the CGL2 and CGL models. As discussed for example by Passot et al. (2006), the closure prescribes a limitation for the minimal value of |B|/B 0 > 1 − 1/a p that is required to prevent temperature singularity. Or in other words, by separating |B| = B 0 + B (1) the requirement reads B (1) > −1/a p , meaning that the fluctuating B (1) cannot be too negative. Mirror instability is often associated with nonlinear structures in the form of magnetic holes and humps, and the requirement implies that the magnetic holes cannot be too deep. To easily analyse the dispersion relations, it is useful to write the closure at the linear level, in the following form The result can be obtained by linearizing (9.44) or (9.43), and at the linear level ∂ i |B| lin a p (1 − a p )B z . Considering the simplest non-dispersive model (by neglecting the Hall term and Π), the dispersion relation of this fluid model reads (9.46) In comparison with the dispersion relation (3.193) of the adiabatic CGL model, the 'static' closure eliminates the erroneous 1/6 factor in the A 0 coefficient and yields the correct mirror threshold. The erroneous 1/6 factor in the CGL model can be therefore interpreted as a result of inadequacy of adiabatic closures in the very slow-dynamics context, such as the mirror instability. For completeness, solutions for parallel propagation (k ⊥ = 0) are ω = ±v A k and ω = ± (β /2)k .
For perpendicular propagation (k = 0) the solution is ω = ± 1 + a p β − a 2 p (β /2)k ⊥ . Obviously, somewhere beyond the mirror threshold, the perpendicular fast mode experiences unphysical instability, which can be interpreted as a result of inadequacy of 'static' (isothermal) closures in the fast-dynamics context.
To clearly demonstrate that variations of magnetic field in the 'static' closures (9.45) are indeed crucial in recovering the mirror instability, let us neglect them and quickly consider strictly isothermal closure p = nT (0) ; p ⊥ = nT (0) ⊥ , which yields the following dispersion relation As can be seen from the A 0 coefficient, this model does not recover the mirror instability. For parallel propagation the solutions are the same as for the 'static' closure, nevertheless, the perpendicular fast mode is now always stable with ω = ± 1 + a p (β /2)k ⊥ . As a double check, after prescribing a p = 1, the dispersion relation (9.47) is equal to (9.46), since the B z contributions in the 'static' closure disappear. Also, we considered models with polytropic indices γ , γ ⊥ in § 3.11, and the dispersion relation (9.47) is consistent with (3.171), after prescribing γ = 1, γ ⊥ = 1.

Bi-kappa fluid model (BiKappa)
10.1. Bi-kappa distribution function Often in the solar wind and in space plasma physics, distribution functions with elongated high-energy 'suprathermal' tails are observed, that cannot be modelled with the bi-Maxwellian distribution. A slightly more general distribution function is the bi-kappa distribution, that contains one free parameter κ, see for example 130 P. Hunana and others and references therein. The κ parameter can be different in the parallel and transverse directions (see e.g. Basu (2009), dos Santos, Ziebell & Gaelzer (2015). Here, we consider only one free κ parameter. The distribution has the following form where the generalized thermal speeds and Γ is the usual gamma function. In the limit κ → ∞, the generalized thermal speeds are the usual thermal speeds and the distribution function is bi-Maxwellian, since lim Why the generalized thermal speeds are defined this way will become clear from later calculations of the moments, and the choice of the power law index −(κ + 1), instead of perhaps the more logical −κ, is attributed to purely historical reasons. Calculating second-order (pressure) integrals with this distribution requires for convergence κ > 3/2, and this value is also the minimum value for the generalized thermal speeds (10.2) to have real values. Sometimes, in the definition of the isotropic kappa distribution, a value of Γ (3/2) = √ π/2 is used and π 3/2 = 2π √ π/2 = 2πΓ (3/2). Again, for brevity of the calculations, it is useful to introduce α = 1/(κθ 2 ), α ⊥ = 1/(κθ 2 ⊥ ), and rewrite the bi-kappa distribution (10.1) in the following form , a > 1 2 . (10.5) For the critical value of a = 1/2, the integral 1/ √ 1 + x 2 dx = arcsinh(x), and since lim x→±∞ arcsinh(x) = ±∞, the integral diverges. For an isotropic kappa distribution f (c 2 ), the three-dimensional integrals can be conveniently evaluated by transformation to spherical coordinates (10.6) In this case, the important integral is and by a direct substitution For an isotropic kappa distribution therefore , (10.9) which explains the normalization factors. For an anisotropic bi-kappa distributions f (c 2 , c 2 ⊥ ), the integration is slightly more complicated because, in contrast to a bi-Maxwellian distribution, the integration cannot be fully separated into two independent integrations over parallel and perpendicular velocity components, and has to be performed successively. It is possible to use a cylindrical coordinate system and integrate However, here we keep the old fashioned Cartesian coordinates, and at the expense of a slightly bit more algebra (i.e. one more integration), we will integrate over d 3 c = dc d 2 c ⊥ = dc dc x dc y . A small added benefit is that we do not have to follow whether the integrals are ∞ −∞ or ∞ 0 , and we can stop writing the boundaries. The most important integral is a slight generalization of (10.5) in the form where for simplicity b, α are assumed to be positive constants. So one can calculate (1 + α y y 2 + α z z 2 ) −(a−1/2) , (10.12) where, importantly, the 'power law' changed from −a to −(a − 1/2). Performing a successive three-dimensional integral therefore yields (10.13) and integrating the bi-kappa distribution over the velocity space ; κ > 1 2 , (10.14) 132

P. Hunana and others
which verifies that the normalization constants in the bi-kappa definition (10.4) are indeed correct, and that f 0 d 3 c = n. For later calculations, it is useful to write down two partial integrals when the integration is done over dc and d 2 c ⊥ that read and (10.20) Continuing with the calculation of velocity moments, it is easy to show that cf 0 d 3 c = 0, which verifies nu = vf 0 d 3 v. Continuing with the higher-order moments, it is possible to calculate those from table integrals (a great help for double checking is an analytic software like Maple or Mathematica) or similarly for a bi-Maxwellian distribution to use a trick with the differentiation with respect to α as where we have also used a property of the gamma function xΓ (x) = Γ (x + 1), and the last integral requires κ > 1 2 . A slightly more general integral is which can be used to evaluate the following integrals (10.25) Since we like to double check everything, the above integral can be also calculated as For the calculation of the parallel pressure, we first integrate over d 2 c ⊥ so we can use result (10.20) and then over dc where we have used (κ − 3 2 )Γ (κ − 3 2 ) = Γ (κ − 1 2 ). Calculating the perpendicular pressure, we here first integrate over dc so we can use result (10.17) and then over

P. Hunana and others
The last two calculations clarify the reasoning behind the definition of the bi-kappa generalized thermal speeds θ , θ ⊥ (here rewritten with the notation α and α ⊥ ), where one starts with a bi-kappa distribution with unspecified α and α ⊥ , and calculating the parallel and perpendicular pressure integrals, yields the required forms α −1 = (κ − 3 2 )(2T /m); α −1 ⊥ = (κ − 3 2 )(2T ⊥ /m). Then, the split to α −1 = κθ 2 ; α −1 ⊥ = κθ 2 ⊥ is dictated by the requirement that for κ → ∞, the generalized thermal speeds converge to the usual thermal speeds, and which also yields that in this limit the distribution is bi-Maxwellian.
The bi-kappa distribution was specified with perhaps a somewhat obscure power law −(κ + 1) instead of the more logical −κ, which is attributed to purely historical reasons. If one wants to define an anisotropic power-law distribution, a first guess would be which then yields normalization constants Then, calculating the parallel and perpendicular pressures yields the requirement α −1 ,⊥ = (κ − 5 2 )(2T ,⊥ /m). Further writing α −1 ,⊥ = κθ 2 ,⊥ yields a distribution function with generalized thermal speeds (10.32) The convergence of the pressure integrals requires κ > 5 2 , which is also the limiting case for the thermal speeds to be real. In the limit κ → ∞ the generalized thermal speeds converge to the usual thermal speeds and the distribution converges to the bi-Maxwellian. To double check that the distribution (10.31) is really correct, one can of course substitute κ → κ + 1, which fully recovers the original bi-kappa distribution (10.1), since κθ 2 → (κ + 1)(1 − 5/2(κ + 1))(2T/m) = κ(1 − 3/2κ)(2T/m). From now on, let us continue calculations with the original bi-kappa distribution (10.1), that uses the −(κ + 1) power-law index.
10.2. Bi-kappa fluid closure Now that we are familiar with the bi-kappa distribution, we are ready to calculate higher-order moments. The heat fluxes are zero and easy to calculate, since both q , q ⊥ are anti-symmetric in c , and by a direct calculation As discussed previously for the bi-Maxwellian distribution, the heat fluxes must be kept in the general form, since we want to consider closure performed at the fourthorder moment level. To calculate the fourth-order fluid moments, we will need the following integrals The r moment then calculates as (10.39) where in the last step we used that Γ (x)(x + 1) 2 /Γ (x + 2) = (x + 1)/x with x = κ − 5/2. The r ⊥ moment calculates as (10.40) and finally the r ⊥⊥ moment calculates as

P. Hunana and others
To summarize, the following closure can be constructed where for the bi-kappa distribution (10.1), the α κ coefficient reads It is noteworthy that all three fourth-order moments are just multiplied by the same constant α κ . Obviously, in the limit κ → ∞, α κ → 1, and the 'normal' bi-Maxwellian closure is obtained. As suggested for example by Chust & Belmont (2006), a very large class of distribution functions can be accounted for by considering closures where the 3 proportionality constants α , α ⊥ , α ⊥⊥ have to be determined, once a specific distribution function is prescribed. Such an unspecified closure can account for a very wide class of distribution functions. Also, a much more complicated distribution functions can be modelled by considering multi-fluid models. For example, one could consider a plasma consisting of 3 different electron species, that would describe the core, halo and the tail/strahl parts of the electron distribution function. Further information about the core, halo and strahl components of the electron distribution function can be found for example in Vocks et al. Here, we consider closure (10.42), (10.43). Nevertheless, in the following calculations one can we keep the value of α κ unspecified. Therefore, the results are actually valid for a much larger class of fluid models, and not just for a bi-kappa distribution. Let us completely neglect the non-gyrotropic Π and r ng in the heat flux equations (7.98) and (7.99). We need to calculate (10.46) and the nonlinear heat flux equations for the bi-kappa fluid model read Obviously, for α κ = 1 (i.e. for κ → ∞), the bi-Maxwellian heat flux equations are obtained. The model is now closed, and it is accompanied by (9.2)-(9.6). The model is nonlinear and numerical simulations for a chosen (and fixed) value of κ can be considered. Even though the κ coefficient is fixed, similarly to bi-Maxwellian nonlinear simulations, the temperature anisotropy is free to develop, and complicated plasma heating processes can be studied (we recommend adding the FLR pressure stress forces to the scalar pressure equations). We abbreviate the non-dispersive fluid model that uses the above heat flux equations as 'BiKappa'. If the Hall term is considered in the induction equation, the 'Hall-BiKappa' fluid model is obtained. Similarly, considering FLR1 and FLR2 corrections to the pressure tensor, yields the 'Hall-BiKappa-FLR1' and 'Hall-BiKappa-FLR2' fluid models. The heat flux equations are normalized, linearized, written in the x-z plane and Fourier transformed according to where we have used the same plasma beta definition as always β /2 = p (0) /(V 2 A ρ 0 ). The heat flux equations are accompanied by the linearized system (9.11)-(9.19).
To consider FLR3 corrections, one needs to derive non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ for the bi-kappa distribution. For the first-order vectors the derivation is straightforward, and by following the calculations in appendix D, it is easy to show that the nonlinear Here, we neglected perturbations r, since we will not consider Landau fluid models for bi-kappa distribution. For α κ = 1, the results are of course equivalent to the bi-Maxwellian expressions (D 37), (D 60). The second-order heat flux vectors are not addressed here, since we would have to derive r ng for the bi-kappa distribution.
The equation (11.3) is actually not that much different from (11.1), the Vlasov equation is just fully linearized and written in Fourier space. However, in kinetic calculations, the (11.3) is not directly integrated in this form. The crucial technique that reveals the presence of Landau damping in the Vlasov equation is that the equation (11.3) is first divided by (ω − k z v z ), and an expression for f (1) is obtained. Only then is the equation integrated over the velocity space, yielding The integral contains a 'singularity' in the denominator, and it is not obvious how to calculate this integral. Such singularities are called wave-particle resonances, and can have a more general form ω − k z v z ± lΩ, where Ω is the cyclotron frequency of the considered species, and l is an integer. The resonance for l = 0 is called the Landau resonance, and the other resonances for l = 0 are called cyclotron resonances. It is the presence of these resonances in integrals such as (11.4), that yields collisionless damping mechanisms, the Landau damping and the cyclotron damping. In Part 2 of our guide, we will see that the integral actually cannot be 'calculated', i.e. the integral cannot be expressed through elementary functions, even if f 0 is a simple Maxwellian. The integral is just expressed through the famous plasma dispersion function Z(ζ ), where ζ = ω/(|k z |v th ), since the plasma dispersion function was specifically developed to describe this integral. Skipping many technical (but very important) details, such as the definition of Z(ζ ), for a Maxwellian f 0 the density integral is simply expressed as (11.5) One can define function R(ζ ) = 1 + ζ Z(ζ ), that is called the plasma response function. Now, for example, considering a proton-electron plasma with T (0) e = T (0) p at wavelengths that are much longer than the Debye length, yields the following dispersion relation and numerical solution R(ζ p ) + R(ζ e ) = 0; ⇒ ζ p = ω |k z |v thp = ±1.457 − 0.627i. (11.6) The negative imaginary part represents the Landau damping, and the effect comes from the Landau resonance in the integral (11.4). Obviously, integrating the Vlasov equation directly in the form (11.1) or (2.1) neglects the wave-particle resonances in velocity space, which clarifies why traditional fluid models do not contain collisionless damping mechanisms. It is f (1) -the fluctuations/perturbations around f 0 -where the wave-particle resonances are written down explicitly, that yields the collisionless damping mechanisms in kinetic theory. To introduce collisionless damping into fluid models, we have no other choice, and we have to calculate the hierarchy of moments again, this time by integrating over the 'kinetic' f (1) . We find it useful to introduce the following vocabulary: (i) Fluid moments = integrals over a general (unspecified) distribution function f , or over a specific f 0 . Evolution equations for fluid moments are obtained by direct integration of the Vlasov equation written in a form that does not explicitly contain wave-particle resonances. The resulting expressions can be nonlinear or linear and contain only fluid variables, i.e. expressions do not contain Z(ζ ).

P. Hunana and others
(ii) Kinetic moments = integrals over f (1) (perturbations f (1) = f − f 0 , where a specific equilibrium f 0 is assumed), where the wave-particle resonances in velocity space are written down explicitly. The resulting expressions are linear and contain Z(ζ ), which is not a fluid variable. Importantly, we want to keep the entire nonlinear fluid hierarchy discussed so far, including the nonlinear heat flux equations. Therefore, considering bi-Maxwellian distribution f 0 , instead of the 'normal' closure, the fourth-order moment is decomposed in the following form (11.7) The first terms represent the 'normal' closure, i.e. when an exact bi-Maxwellian f 0 is prescribed. The tilde components represent fluctuations/perturbations around f 0 , and a good notation would also be δr. The decomposition (11.7) yields the following form of the nonlinear heat flux equations Similarly to decomposition (11.7), 'kinetic' corrections with wave-particle resonances can also be considered for the non-gyrotropic r ng , which is not addressed in this guide. Neglecting those contributions, the heat flux equations are equivalent for example to equations (14), (15) of Passot & Sulem (2007), and equations (9), (10) of Passot et al. (2012). In Part 2 of this guide, we will calculate the 'kinetic moments' for r , r ⊥ , r ⊥⊥ , by integrating over f (1) . In the 3-D electromagnetic geometry, we will use a more sophisticated f (1) in the gyrotropic limit, which describes both Landau damping and its magnetic analogue, the transit-time damping. The expressions for the entire 'kinetic' hierarchy will contain Z(ζ ), or actually R(ζ ), which is not a fluid variable. Nevertheless, we will see that by introducing the concept of the Padé approximation to R(ζ ), that sometimes a rare possibility arises, when the 'kinetic' moment r can be expressed through lower-order 'kinetic' moments. Such a closure, where the last retained kinetic moment (in this case r) is expressed through lower-order moments, in a way such that the kinetic Z(ζ ) function is completely eliminated and that the closure is valid for all ζ values, will be called a Landau fluid closure. For the impatient reader, we will reveal that one of the possibilities is going to be (written in Fourier space) (11.10) (11.11) r ⊥⊥ = 0. (11.12) 144

P. Hunana and others
Again, because of the symmetric operator, it does not matter if X (n) is applied to the momentum equation from the left or right. The subtraction also conveniently eliminates the electric field, however it introduces ∇ · p. We also need an identity (12.15) where in the last step we emphasized that it does not matter if the B× applies only to u or to the entire tensor uX (n−1) , since the vector product operates only through the first component of the tensor anyway. The final evolution equation reads ∂ ∂t (12.16) and it is valid for n 2. Evaluation at m = 3 recovers the heat flux tensor equation (7.22) and evaluation at m = 2 recovers the pressure tensor equation (2.28). Note that X (4) = r, X (3) = q, X (2) = p, however, and very importantly, according to our definition (12.1), X (1) = 0 and it is not equal to the velocity u, which is defined according to (2.4). Additionally, X (0) = ρ. The momentum equation can still be recovered by evaluating (12.13) at n = 1, and by dividing by ρ (where naturally the symmetric operator does not have any influence on a vector, and u S i = u i ). Similar equation to our (12.13) was also obtained for example by Waelbroeck (2010), equation (22), even though we did not verify whether it is consistent with ours since X (n) in that work is defined with v instead of our c.
We want to always decompose the tensor X (n) into its gyrotropic and non-gyrotropic parts, and to have space to write the 'g' and 'ng', let us move the index (n) down, so that X (n) = X g (n) + X ng (n) . (12.17) By generalizing the results we have seen for the gyrotropic p g , q g , r g , the same result is obtained for the nth moment, and at low frequencies the gyrotropic part must satisfy (12.18) which can be indeed viewed as a definition of gyrotropy. The general (12.16) therefore rewrites as (12.19) yielding an implicit equation for the non-gyrotropic part in the form (12.20) Ideally, an 'inversion procedure' should be found for the left-hand side, so that an equation for X ng (n) can be found, and then by expanding the right-hand side (for 150 P. Hunana and others (12.48) Continuing to higher orders, the hierarchy of deviations of Maxwellian moments therefore evaluates as The equations are written in physical units, and these results will be useful in Part 2, where we will calculate Landau fluid closures for deviations of various moments.
12.2. Impossibility of going beyond CGL2 without Landau fluid closures One can analyse the dispersion relations easily in physical units, but since in this Part 1 we used normalized units with k = kV A /Ω p , ω = ω/Ω p and β = v 2 th /V 2 A in almost all dispersion relations, let us rewrite the fluid hierarchy to normalized units, so that everything feels more familiar. The X (5) is normalized with p 0 V 3 A , the X (6) with p 0 V 4 A and X (n) with p 0 V n−2 A . By dropping the normalization tilde as usual (but obviously not for even moments such as r, that should be perhaps called δr), the equations in Fourier space read −ωρ + ku z = 0; (12.55) −ω X (6) + kX (7) − 45 4 β 2 kq = 0. (12.61) Now, using just first 3 equations (12.55)-(12.57), and closing the system with a closure q = 0, yields the CGL model, with solution CGL : ω = ±k β 3 2 , (12.62) the familiar CGL ion-acoustic mode. Going higher in the fluid hierarchy and using the first 4 equations (12.55)-(12.58), and closing the system with r = 0 (which is equivalent to prescribing r = 3p 2 /ρ), yields the CGL2 model. The dispersion relation is ω 4 − 3βk 2 ω 2 + 3 4 β 2 k 4 = 0, with solutions CGL2 : ω = ±k β 3 2 + 3 2 ; (12.63) the already obtained (9.28), (9.29). The heat flux fluctuations in the CGL2 model therefore did 'split' the CGL ion-acoustic mode to two modes. One might assume that this will always be the case when going higher and higher in the fluid hierarchy.
With the next example we demonstrate that this does not happen. Let us consider fluctuations in r, and use the first 5 equations (12.55)-(12.59) with the closure X (5) = 0. We might call this model CGL3. The dispersion relation reads ω 4 = 15 4 β 2 k 4 , and the solutions are Importantly, the last two modes are imaginary, and one is unstable. Therefore, such a fluid model does not make any physical sense. Why did this happen? For example, expressing everything through pressure fluctuations, and putting the r fluctuations on the right-hand side yields ω 4 − 3βk 2 ω 2 + 3 4 β 2 k 4 p = k 2 ω 2 r; (12.66) = k 3 ωX (5) + −3βk 2 ω 2 + 9 2 β 2 k 4 p.
(12.67) Therefore, if r = 0, so when (12.66) is used, the CGL2 model is obtained. However, when fluctuations in r are considered, so when (12.67) is used, the term −3βk 2 ω 2 cancels on both sides, with the resulting dispersion relation ω 4 − 15 4 β 2 k 4 p = k 3 ωX (5) , (12.68) and by prescribing closure X (5) = 0, yields the dispersion relation of the CGL3 model. Now it is easy to go higher in the fluid hierarchy. By considering the first 6 equations (12.55)-(12.60), the dispersion relation reads and prescribing closure X (6) = 0 (which is equivalent to prescribing X (6) = 15p 3 /ρ 2 ), yields a model that we can call CGL4, with numerical solutions  All the models beyond CGL2 contain modes that are unstable.

Double check
Let us double check our results, to verify that we calculated everything correctly. As stated previously, a closure should be performed only at the last retained fluid moment, and this is especially true when nonlinear numerical simulations are performed. 13 The CGL and CGL2 systems were already discussed at great length. Let us check the CGL3 dispersion relation, and let us work directly in physical units. The system of equations that we have in mind reads (12.80) and the closure is performed by setting X (5) = 0. Importantly, the fourth-order moment r is kept undetermined, and it is not separated to its 'core' and r, nor any specific form of a distribution function is prescribed yet (other than the mean values of odd moments X (l) 0 = 0). Calculating dispersion relation of this system (a matrix multiplied by vector (ρ, u z , p, q, r)) yields result ω 4 − 5k 4 r 0 /ρ 0 = 0, which for r 0 = 3p 2 0 /ρ 0 recovers the CGL3 dispersion relation. Our calculations were therefore done right. Similarly, double checking the CGL4 dispersion, the system (12.80) is supplemented with (12.33), and a closure is performed with X (6) = 15p 3 /ρ. Again, the r is left untouched, and dispersion relation of this system recovers the CGL4 result. Checking the CGL5 dispersion relation, the system is supplemented with (12.34), and a closure X (7) = 0. Again, the X (6) is kept untouched and generally undetermined. Calculating the dispersion relation yields ω 6 − 7k 6 X (6) 0 /ρ 0 = 0, and prescribing Maxwellian X (6) 0 = 15p 3 0 /ρ 2 0 recovers the CGL5 result. Going higher in the hierarchy, and using closure X (8) = 105p 4 /ρ 3 recovers the CGL6 results. And finally, keeping all the variables up to X (8) undetermined and using closure X (9) = 0 yields dispersion relation ω 8 − 9k 8 X (8) 0 /ρ 0 = 0, which for Maxwellian X (8) 0 = 105p 4 0 /ρ 3 0 recovers the CGL7 result. Or, in general, without yet performing a closure, and going higher and higher in the hierarchy step by step, a straightforward calculation yields CGL : −ω 2 + 3 where the validity can be proven by induction and using (12.37), (12.38). CGLl model (or CGL of lth order) is therefore defined as performing a closure on X (l+2) moment. In the equations (12.81)-(12.89) we used the word CGL, even though no specific distribution function (nor any closure) was prescribed yet.
density and momentum equations, it is possible to derive an evolution equation for a general nth-order fluid moment (n 2), see (12.16). Then, evolution equations for the pressure tensor and heat flux tensor can be obtained by simply evaluating (12.16) at n = 2 and n = 3.
• Even though the fluid hierarchy is infinite, the reformulation of the kinetic description to a fluid formalism is not necessarily complete, since the usual calculations neglect wave-particle resonances in velocity space. These wave-particle resonances, that are implicitly present in the Vlasov equation, are responsible for collisionless damping mechanisms, such as Landau damping, transit-time damping and cyclotron damping.
• The presence (or absence) of collisionless damping mechanisms in a fluid hierarchy is determined by the type of closure selected to truncate that fluid hierarchy. We differentiate between two classes of closures: a) 'classical' (non-Landau fluid) closures, that neglect wave-particle resonances, and b) Landau fluid closures, that account for Landau wave-particle resonances and associated Landau damping and transit-time damping, addressed in Part 2.
• There are currently no known closures that account for cyclotron wave-particle resonances and associated cyclotron damping. Nevertheless, there is a priori no reason why such fluid closures can not be found in the future, at least for the simplified case of electromagnetic propagation along the magnetic field (slab geometry).
• Considering 'classical' closures, the hierarchy of fluid moments therefore does not contain collisionless damping mechanisms, regardless of the order to which the hierarchy is developed. Or in another words, collisionless damping is beyond all orders in the classical hierarchy of fluid moments.
• In fact, it is impossible to go beyond the fourth-order moment with classical bi-Maxwellian fluid closures X (n) = X (n) 0 + X (n) , where X (n) 0 represents bi-Maxwellian value of X (n) and its perturbation X (n) = 0, since all the fluid models contain higherorder modes that are unstable. This is perhaps the most surprising result discussed in Part 1. For a detailed proof, see § 12.2, where a 1-D geometry is considered and a bi-Maxwellian closure is used at the nth-order fluid moment, which yields dispersion relation (12.96). The same dispersion relation will be valid in a 3-D geometry for the propagation along the magnetic field, where additional de-coupled modes will be present as well. For n > 4, all the fluid models contain unphysical instabilities. The same result is expected for other distribution functions.
• From a collisionless (or weakly collisional) perspective, the long-wavelength low-frequency limit of a distribution function is not necessarily an isotropic Maxwellian, but a general distribution function that is gyrotropic, i.e. isotropic only in its perpendicular velocity components. Fluid moments are therefore typically decomposed into their gyrotropic and non-gyrotropic parts. For example 156 P. Hunana and others the pressure tensor p = p g + Π, the heat flux tensor q = q g + q ng , the fourth-order moment r = r g + r ng . The non-gyrotropic parts are also referred to as the FLR corrections, since they represent deviations from gyrotropy at small spatial scales, when the Larmor radius (gyroradius) is not infinitely small.
• A general nth-order fluid moment contains 1 + int[n/2] scalar gyrotropic moments, where the function 'int' means integer part.
• The pressure tensor is second-order fluid moment and contains two gyrotropic moments, p and p ⊥ . Therefore, any fluid description of collisionless (or weakly collisional) plasmas, requires two separate evolution equations for p and p ⊥ . Importantly, from a linear perspective the evolution equations remain different even when the distribution function is isotropic, i.e. even when the mean pressure values are equal (p (0) = p (0) ⊥ ), since pressure fluctuations in the directions parallel and perpendicular to the local magnetic field remain anisotropic. We often use proton pressure (temperature) anisotropy coefficient a p = T (0) ⊥ /T (0) = p (0) ⊥ /p (0) .
• The simplest fluid model that describes a collisionless plasma in an adiabatic regime is the CGL description (model) -named after Chew-Goldberger-Law (Chew et al. 1956). The model is non-dispersive and obtained by a closure with zero heat flux. The CGL description should be viewed as collisionless magnetohydrodynamics (collisionless MHD), since the usual MHD description with isotropic scalar pressure is highly collisional implicitly.
• The CGL pressure tensor is decomposed with respect to the direction of the local magnetic field lines according to p g = p bb + p ⊥ (I −bb). It is possible to switch to the reference frame of the magnetic field lines, where the pressure tensor is a diagonal matrix. However, it should not be forgotten that in a such a reference frame, imposing an external (mean) magnetic field is not straightforward. For nonlinear numerical simulations of turbulence and plasma heating, the laboratory reference frame is strongly recommended.
• The two pressure equations of the CGL model can be interpreted as conservation laws for the first and second adiabatic invariants, see § 2.5. The first adiabatic invariant is a conservation of the magnetic moment of a particle that is periodically gyrating around a mean magnetic field. The second adiabatic invariant is a conservation of the average parallel momentum of a particle that is completely trapped and periodically bouncing inside of a magnetic bottle.
• Without performing (yet) any kind of closure and leaving the heat flux tensor q and FLR pressure tensor Π unspecified, it is possible to derive rigorously exact evolution equations for p and p ⊥ , see (2.91)-(2.92).
• Rewriting the pressure equations to an alternative form, equations (2.100)-(2.101), shows that the adiabatic invariants in the CGL model are broken by the inclusion of the Hall term in the induction equation, gyrotropic heat flux, non-gyrotropic FLR corrections to the pressure and heat flux, as well as by coupling of various species together. All of these contributions therefore yield very complicated anisotropic plasma heating processes, that can be studied only by nonlinear numerical simulations.
• The anisotropic plasma heating does not simplify much even in the case of periodic boundary conditions (i.e. when the system can be viewed as completely isolated) and when averaging over the entire spatial volume is performed and expressed as a conservation of total energy, see (2.113)-(2.116). Only when the anisotropic heating is neglected, i.e. when only total plasma heating studied, the conservation of energy significantly simplify, see (2.109)-(2.111).
• By using the concept of polytropic indices, the pressure equations in the ideal CGL model can be interpreted as having γ = 3 and γ ⊥ = 2. Consequently, the CGL dispersion relation is in general not equal to the MHD dispersion relation, even when a p = 1. The exception is the Alfvén mode, which for a p = 1 propagates with the same phase speed ω/k = V A cos θ in both CGL and MHD models, for all the propagation directions θ .
• Polytropic indices can be further related to the number of degrees of freedom i through γ = (i + 2)/i, which for the CGL model yields i = 1 and i ⊥ = 2. The CGL pressure equations can be therefore viewed as being composed of strongly coupled 1-D and 2-D dynamics, whereas the dynamics in the MHD model with γ = 5/3 and i = 3 can be viewed as isotropically three-dimensional.
• In the CGL model, the Alfvén mode propagates with the phase speed ω/k = V A cos θ 1 + (β /2)(a p − 1), which matches the (collisionless) kinetic theory in the long-wavelength limit. Consequently, for a p = 1, the Alfvén mode in the MHD model matches the kinetic theory as well. The CGL result is very useful when numerically solving kinetic dispersion relations (for example with the WHAMP solver), since the Alfvén mode can be easily identified at long wavelengths.
• In contrast to MHD, where the ordering of phase speeds is always v s v A v f (slow, Alfvén, fast), in collisionless plasmas the oblique slow mode can become faster than the oblique Alfvén mode. In general, oblique slow and fast modes in the CGL model do not necessarily match the kinetic theory in the long-wavelength limit. Nevertheless, the effect when v s > v A is present in the CGL model and exists even for a p = 1, see § 3.6. The effect is very important for the parallel firehose instability.
• Considering strictly perpendicular propagation (θ = 90 • or k = 0), the fast mode in the CGL model propagates with the phase speed ω/k = V A 1 + a p β , which matches the kinetic theory in the long-wavelength limit. The expression can be rewritten as (ω/k) 2 = V 2 A + v 2 th⊥ = V 2 A + 2(p (0) ⊥ /ρ 0 ), and can be contrasted with the MHD result (ω/k) 2 = V 2 A + 5 3 (p (0) /ρ 0 ). The factor of 2 in the CGL result comes from γ ⊥ = 2, which shows that the adiabatic assumption of the CGL model is appropriate for the perpendicular fast mode.
• Contributions of the Hall term completely disappear for strictly perpendicular propagation.
• It is more difficult to analytically show where the adiabatic CGL value γ = 3 is appropriate, since for example the parallel propagating ion-acoustic (sound) mode is strongly Landau damped in kinetic theory (unless the electrons are hot), which also modifies its real frequency. The ion-acoustic mode in the CGL model therefore does not match kinetic theory. One needs to consider Landau fluid models with proton and electron species, which is addressed in Part 2.

P. Hunana and others
• Nevertheless, when electrons are hot (T e T i ) and the ion-acoustic mode is Landau damped only weakly, the real phase speed obtained from kinetic theory in the long-wavelength limit reads ω/k = √ (T e + 3T i )/m i . Thus, in this example the kinetic theory is matched by a fluid model where protons are described adiabatically (with the CGL value γ = 3), while electrons are isothermal (with γ = 1). Also, a good example for γ = 3 is the real frequency of the Langmuir mode, ω 2 = ω 2 pe + 3(T (0) e /m e )k 2 , which is addressed in Part 2.
• For anisotropic temperatures (a p = 1), some modes can become unstable. A threshold of an instability that is obtained in the long-wavelength (non-dispersive) limit is called a 'hard' threshold. The CGL model contains 3 instabilities, the oblique firehose instability, the parallel firehose instability and the mirror instability. The mirror instability is not described correctly.
• The oblique firehose instability is an instability of the oblique Alfvén mode. Since the Alfvén mode is described correctly by the CGL model, its instability threshold 1 + (β /2)(a p − 1) < 0 matches the kinetic theory. The instability requires a p < 1 and (in the long-wavelength limit) β > 2. Since the CGL model is non-dispersive, here the instability growth rate only has a simple cos θ dependence. Once dispersive effects are considered, the instability reaches maximum growth rate at some wavenumber and angle θ that is oblique, see figure 10.
• The parallel firehose instability is an instability of the whistler mode, which can be shown by considering the Hall-CGL model. However, considering usual (non-causal) analytic solutions (4.32)-(4.35), one arrives at a contradiction that for ω r > 0 the whistler mode is unstable, and that for ω r < 0 the whistler mode is stable. The contradiction arises, because at the range of wavenumbers where the firehose instability exists, the ion-cyclotron and whistler modes are completely degenerate and distinguishing between them loses sense. The problem is resolved by introducing a small (∼ ) causal dissipation into the momentum equations, which is later removed by the limit lim →0 + . The procedure yields causal dispersion relations (4.43)-(4.46), where the whistler mode is firehose unstable for both positive and negative ω r , and the ion-cyclotron mode is stable.
• The parallel firehose instability occurs for quasi-parallel (small θ ) propagation directions, with the maximum growth rate at θ = 0 • , see figure 10. In the non-dispersive CGL model it corresponds to the instability of the slow mode. Its threshold obtained for θ = 0 • is equivalent to the oblique firehose instability (since v s = v A ), and matches the kinetic theory. At first it might sound surprising that the slow mode connects to the whistler mode once dispersive effects are considered. However, considering θ = 0 • , β > 2, a p < 1 yields that the fast mode is the ion-acoustic (sound) mode and v s = v A , see figure 1. After dispersive effects are introduced, which split the Alfvén and slow mode to the ion-cyclotron and whistler modes, it is clarified that the instability is associated with the whistler mode. For quasi-parallel propagation directions the situation is more obvious, since v s > v A (see figure 1), and the slow mode clearly becomes the whistler mode.
• In the non-dispersive CGL model, all 3 instabilities are non-propagating, i.e. purely growing with zero real frequency. Once dispersive effects are introduced, the parallel firehose instability becomes propagating.
• The mirror instability requires a p > 1 and usually develops for highly oblique propagation angles. In the CGL model it corresponds to the instability of the slow mode. Its threshold obtained in the highly oblique (but not completely perpendicular) limit reads 1 + a p β − 1 6 a 2 p β < 0, and with respect to kinetic theory the 1/6 factor is erroneous. The erroneous 1/6 factor shows that the adiabatic CGL closure is not appropriate for the very slow dynamics of the mirror instability.
• The simplest fluid closure that recovers the correct threshold of the mirror instability is the 'static' closure, see § 9.4. This closure is derived from the 'normal' closure in a simplified quasi-static approximation, and represents a generalization of the isothermal closure (in the presence of temperature anisotropy and variations of magnetic field strength). The 'static' closure is very useful since it clarifies that the highly oblique very slow dynamics of the slow mode (and of the associated mirror instability) is better described by a generalized isothermal closure than the adiabatic CGL closure.
• However, the 'static' closure modifies the dynamics of the perpendicular fast mode, that now does not match the kinetic theory. Additionally, the fast mode now experiences an unphysical instability, even though the instability is beyond the mirror threshold, and should perhaps not play a role in numerical simulations. The erroneous result shows that the generalized isothermal closure is not appropriate for the relatively fast dynamics of the perpendicular fast mode, where the adiabatic CGL closure is appropriate.
• Heuristically, the correct mirror threshold can be also obtained by keeping the CGL value γ ⊥ = 2 and modifying the parallel polytropic index to γ = 1/2. Such a closure does not alter the dynamics of the perpendicular fast mode, since the phase speed depends only on γ ⊥ (and not on γ ). Thresholds of parallel and oblique firehose instability are not altered either. Nevertheless, the γ = 1/2 is even below the isothermal value γ = 1, the closure is heuristic, and it is mentioned only as a curiosity.
• Obviously, to correctly capture both the slow dynamics of the highly oblique mirror instability, as well as the fast dynamics of the perpendicular fast mode, neither adiabatic nor (generalized) isothermal closures are sufficient. Considering classical (rigorously derived) closures, one has no other choice but to keep the correct adiabatic CGL values γ = 3, γ ⊥ = 2 unaltered, and instead, break the (long-wavelength) adiabaticity by going higher in the fluid hierarchy. One needs to consider fluid models with evolution equations for the heat flux tensor q, that has two gyrotropic components, q and q ⊥ . The system is closed at the fourth-order moment r, that has 3 gyrotropic components, r ; r ⊥ ; r ⊥⊥ . To perform a closure, a specific distribution function (or a class of distribution functions) has to be assumed.
• Considering a bi-Maxwellian distribution function, the closure r = 0 is known as the 'normal' closure (see above). We call this non-dispersive bi-Maxwellian fluid model the 'second-order CGL', abbreviated as CGL2. The abbreviation CGL2 seems beneficial, since the name CGL is associated with collisionless magnetohydrodynamics, now just taken to one higher order. Additionally, by showing that all classical bi-Maxwellian models beyond CGL2 do not make physical sense, fluid models can be easily classified as based on CGL or CGL2 descriptions. For direct comparison of the two systems, see dispersion relations 160 P. Hunana and others of the CGL model (3.193) and the CGL2 model (9.23). The Hall term introduces the simplest dispersive effects; see the dispersion relations of the Hall-CGL model (4.81), and the Hall-CGL2 model (9.40).
• By comparing the dispersion relations, it is obvious that the gyrotropic heat flux fluctuations in the CGL2 model modify the general dynamics (the mean heat flux values are zero q (0) = q (0) ⊥ = 0). One of the exceptions is the perpendicular fast mode, which is not altered by the gyrotropic heat fluxes, and therefore described correctly. Other exceptions are the oblique Alfvén mode (the oblique firehose instability), and once dispersive effects are introduced also the parallel propagating ion-cyclotron and whistler modes (the parallel firehose instability).
• Two additional evolution equations typically create two additional modes. For example, for the case of highly oblique CGL2 propagation, one obtains the Alfvén mode, the fast mode, and 3 'slow' modes. One of these slow modes is responsible for the mirror instability, and is sometimes called the mirror mode. The CGL2 mirror threshold matches the kinetic theory. The CGL2 model is therefore the simplest classical model that correctly captures (in the long-wavelength limit) both the mirror instability threshold, as well as the perpendicular fast mode. Nevertheless, only the 'hard' mirror threshold is correctly recovered, and capturing the (long-wavelength) mirror instability growth rate requires Landau damping.
• If a closure for the bi-kappa distribution is performed at a higher nth-order moment (which is presumably not possible with classical closures since it is expected that for n > 4 all models contain unphysical instabilities), the closure will be restricted to κ > (n + 1)/2, which increases with n. The requirement for the minimum κ value is there, since the power law in the bi-kappa distribution is assumed all the way till infinite velocities, and the requirement is necessary for the convergence of the velocity integrals. Therefore, by going higher and higher in the bi-kappa fluid hierarchy, the minimum value of κ increases, and in the limit κ → ∞ the required distribution function converges towards bi-Maxwellian.
• The dispersion relation for the BiKappa model is given by (10.53). A general α κ coefficient is used, so the dispersion relation is valid for a much larger class of distribution functions that can be closed with an analogous closure with just one α κ coefficient. The dispersion relation shows that the general dynamics is of course altered by the value of the α κ coefficient. Nevertheless, few exceptions are: the oblique Alfvén mode (oblique firehose threshold), the perpendicular fast mode, the parallel firehose threshold and the mirror threshold. When dispersive effects are considered, and figures similar to figure 10 created, the growth rates for oblique propagation will of course be affected by the α κ value. The Hall-BiKappa dispersion relation is given by (10.69).
• The non-gyrotropic (FLR) corrections Π to the pressure tensor were studied in § 5. The tensor Π is described by the pressure tensor equation implicitly. To prevent introducing several new independent fluid variables for components of the tensor Π, each with its own evolution equation, expansion on temporal and spatial scales is required. Correct evaluation of the FLR tensor with respect to magnetic field lines is quite cumbersome, and at the leading order (in temporal and spatial scales) one obtains (5.43)-(5.45), or equivalently (5.51)-(5.53). The FLR tensor Π is therefore often evaluated in the linear approximation, which is appropriate when the magnetic field lines are not too distorted.
• Even though the FLR pressure tensor Π has a zero trace, its contributions should not be called 'off diagonal', since the diagonal components are non-zero. In the collisionless case considered here with the mean field in the z-direction, only Π zz = 0, and Π xx = −Π yy = 0. In the collisional case (which was not considered) also Π zz = 0.
• For the purpose of comparing various contributions, we differentiate between FLR1, FLR2 and FLR3 corrections, even though the classification can be a bit blurry and there are various possibilities. Perhaps the best definition that we used later is as follows. The classical FLR1 tensor Π contains only velocity gradients, see (5.74). The FLR2 tensor also contains ∂Π/∂t and the Hall term from the induction equation, see (5.109). It can also contain gyrotropic heat flux contributions, see (5.140). The FLR3 contains non-gyrotropic heat flux vectors S ⊥ , S ⊥ ⊥ , described in § 5.8, with the detailed algebra presented in appendix D. It can also contain the nongyrotropic heat flux tensor σ , which we neglected. For complete clarity in analysing the dispersion relations, the entire Hall-CGL-FLR3 model (linearized in the x-z plane, normalized and Fourier transformed) is written down in § 5.9, see (5.160)-(5.164).
• Considering the strictly perpendicular fast mode, kinetic theory yields the leadingorder FLR corrections (in the long-wavelength limit) to the phase speed in the following form (ω/k) 2 = V 2 A (1 − 1 8 k 2 ρ 2 i ) + v 2 th⊥ (1 − 5 16 k 2 ρ 2 i ). The kinetic result is reproduced by the CGL-FLR3 model (the Hall contributions are zero). Importantly, both the first-and second-order non-gyrotropic heat flux vectors have to be retained.
• If the second-order non-gyrotropic heat flux contributions are neglected in the FLR3 model, yields solution for the perpendicular fast mode (5.159). In such a model, the correction to the Alfvén speed is captured correctly, however, the correction to the thermal speed is + 1 16 k 2 ρ 2 i instead of − 5 16 k 2 ρ 2 i , i.e. the correction has the wrong sign. This is surprising, since the solution of the FLR2 model, equation (5.120), has a correction to the thermal speed − 1 16 k 2 ρ 2 i , i.e. at least the sign is correct. Finally, the FLR1 model yields solution (5.119), and does not capture any correction to the Alfvén speed, additionally, the correction to the thermal speed has the wrong sign as well.
• We used the Hall-CGL-FLR3 model to investigate the parallel and oblique firehose instability, and compare it with results of the FLR2 and FLR1 models, see figures 7-11. It is shown that the growth rates of the parallel and oblique firehose instability are strongly enhanced by the non-gyrotropic heat flux vectors of the FLR3 model, see figure 10 (see the scales of the colour bars).

P. Hunana and others
• In general, when the maximum growth rates are sufficiently large (let us say ω i /Ω ∼ 0.1) and not tiny (such as 0.001), the FLR3 model reproduces the kinetic results with unexpectedly good accuracy. See for example bottom panels of figure 10. Not only is the value of the maximum growth rate for the oblique firehose instability (red lines) approximately the same in the FLR3 model and kinetic theory, the maximum growth rate is also reached approximately for the same angle of propagation and the same wavenumber.
• The oblique firehose instability can be (in general) reproduced with better accuracy than the parallel firehose instability. Nevertheless, for very high β values, the parallel firehose instability is reproduced very accurately, see figure 9.
• Importantly, we show that the non-gyrotropic heat flux vectors in the FLR3 model, partially reproduce the large 'bump' in the imaginary phase speed (growth rate normalized to the wavenumber), when the plasma is close to the long-wavelength limit 'hard' firehose threshold, see figure 7. The result clearly shows that, similarly to kinetic theory, fluid models can develop the firehose instability at small spatial scales, even when they are stable in the long-wavelength limit. Or in other words, fluid models can become unstable even before the 'hard' firehose threshold is reached. This unexpected result is perhaps the second most surprising result discussed in Part 1.
integral can be very confusing. We therefore recommend using the first decomposition, where waves with ω i < 0 are damped. This choice yields that Fourier transformations in the x-z plane (with B 0 in the z-direction) are performed according to a shortcut The transformation (A 1) is technically the inverse/backward Fourier transform F −1f (k, ω). The forward Fourier transform reads The location of the normalization constants 1/(2π) is an ad hoc choice, one just needs to be consistent in using them, especially when calculating convolutions.

Appendix C. MHD dispersion relation
We assume that many people reading this text will have some previous experience with MHD. It is therefore beneficial to obtain the MHD dispersion relation, so that it will be clear how more advanced fluid models are treated in this text. The MHD fluid model reads ∂ρ ∂t + ∇ · (ρu) = 0; (C 1) where γ = 5/3. As in advanced fluid models, it is beneficial to normalize the speed with respect to the Alfvén speed V A , and normalize the length with respect to the ion inertial length d i = V A /Ω p , i.e. by using normalizations (3.23)-(3.28), which in Fourier space yields k = kd i , ω = ω/Ω p . Normalizing the MHD equations and dropping the tilde yields unchanged density, induction and pressure equations, and the momentum equation reads One can decide how to rewrite p (0) /(ρ 0 V 2 A ), either by introducing the usual MHD sound speed C 2 s = γ p (0) /ρ 0 , or by introducing the plasma beta where we use the Boltzmann constant k B = 1, see the footnote 3 after the definition of thermal speeds (3.30). By specifying the mean magnetic field to be in the z-direction, the normalized equations are linearized according to ∂ρ ∂t + ∇ · u = 0; (C 7)

P. Hunana and others
The 'inversion' procedure for the left-hand side of this is equation is very complicated, and will not be addressed here. Nevertheless, the inversion procedure exists, and an interested reader can check equation (43) of Ramos (2005). The leading-order q ng (first order in frequency and wavenumber) can be obtained by making the quantities q, r, p on the right-hand side gyrotropic, and we want to solve The heat flux tensor decomposition is q = q g + q ng , or alternatively q = S + σ , where σ :bb = 0, σ : I = 0. Since Terms on the right-hand side of (D 2) calculate as (D 10) (∇ · r g ) ijkbibj =b kb · ∇r + (r − 3r ⊥ )(b · ∇b k +b k ∇ ·b) −b kb · ∇r ⊥ + ∂ k r ⊥ ; (D 11) (q ng · ∇u) S ijk δ ijbk = (q ng ijl ∂ l u k + q ng jkl ∂ l u i + q ng kil ∂ l u j )δ ijbk = (q ng : I) l (∂ l u k )b k + 2q ng ilkbk ∂ l u i = (S ⊥ + 2S ⊥ ⊥ ) · (∇u) ·b + 2(q ng ·b) : ∇u.